1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */ 2 /* 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 4 * 5 * Permission to use, copy, modify, and distribute this software for any 6 * purpose with or without fee is hereby granted, provided that the above 7 * copyright notice and this permission notice appear in all copies. 8 * 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 11 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 12 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 13 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 14 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 15 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 16 */ 17 /* 18 * Exponential function, long double precision 19 * 20 * 21 * SYNOPSIS: 22 * 23 * long double x, y, expl(); 24 * 25 * y = expl( x ); 26 * 27 * 28 * DESCRIPTION: 29 * 30 * Returns e (2.71828...) raised to the x power. 31 * 32 * Range reduction is accomplished by separating the argument 33 * into an integer k and fraction f such that 34 * 35 * x k f 36 * e = 2 e. 37 * 38 * A Pade' form of degree 5/6 is used to approximate exp(f) - 1 39 * in the basic range [-0.5 ln 2, 0.5 ln 2]. 40 * 41 * 42 * ACCURACY: 43 * 44 * Relative error: 45 * arithmetic domain # trials peak rms 46 * IEEE +-10000 50000 1.12e-19 2.81e-20 47 * 48 * 49 * Error amplification in the exponential function can be 50 * a serious matter. The error propagation involves 51 * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ), 52 * which shows that a 1 lsb error in representing X produces 53 * a relative error of X times 1 lsb in the function. 54 * While the routine gives an accurate result for arguments 55 * that are exactly represented by a long double precision 56 * computer number, the result contains amplified roundoff 57 * error for large arguments not exactly represented. 58 * 59 * 60 * ERROR MESSAGES: 61 * 62 * message condition value returned 63 * exp underflow x < MINLOG 0.0 64 * exp overflow x > MAXLOG MAXNUM 65 * 66 */ 67 68 #include "libm.h" 69 70 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 71 long double expl(long double x) 72 { 73 return exp(x); 74 } 75 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 76 77 static const long double P[3] = { 78 1.2617719307481059087798E-4L, 79 3.0299440770744196129956E-2L, 80 9.9999999999999999991025E-1L, 81 }; 82 static const long double Q[4] = { 83 3.0019850513866445504159E-6L, 84 2.5244834034968410419224E-3L, 85 2.2726554820815502876593E-1L, 86 2.0000000000000000000897E0L, 87 }; 88 static const long double 89 LN2HI = 6.9314575195312500000000E-1L, 90 LN2LO = 1.4286068203094172321215E-6L, 91 LOG2E = 1.4426950408889634073599E0L; 92 93 long double expl(long double x) 94 { 95 long double px, xx; 96 int k; 97 98 if (isnan(x)) 99 return x; 100 if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */ 101 return x * 0x1p16383L; 102 if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */ 103 return -0x1p-16445L/x; 104 105 /* Express e**x = e**f 2**k 106 * = e**(f + k ln(2)) 107 */ 108 px = floorl(LOG2E * x + 0.5); 109 k = px; 110 x -= px * LN2HI; 111 x -= px * LN2LO; 112 113 /* rational approximation of the fractional part: 114 * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2)) 115 */ 116 xx = x * x; 117 px = x * __polevll(xx, P, 2); 118 x = px/(__polevll(xx, Q, 3) - px); 119 x = 1.0 + 2.0 * x; 120 return scalbnl(x, k); 121 } 122 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 123 // TODO: broken implementation to make things compile 124 long double expl(long double x) 125 { 126 return exp(x); 127 } 128 #endif 129