1 /* $OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:18 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* expl.c 20 * 21 * Exponential function, 128-bit long double precision 22 * 23 * 24 * 25 * SYNOPSIS: 26 * 27 * long double x, y, expl(); 28 * 29 * y = expl( x ); 30 * 31 * 32 * 33 * DESCRIPTION: 34 * 35 * Returns e (2.71828...) raised to the x power. 36 * 37 * Range reduction is accomplished by separating the argument 38 * into an integer k and fraction f such that 39 * 40 * x k f 41 * e = 2 e. 42 * 43 * A Pade' form of degree 2/3 is used to approximate exp(f) - 1 44 * in the basic range [-0.5 ln 2, 0.5 ln 2]. 45 * 46 * 47 * ACCURACY: 48 * 49 * Relative error: 50 * arithmetic domain # trials peak rms 51 * IEEE +-MAXLOG 100,000 2.6e-34 8.6e-35 52 * 53 * 54 * Error amplification in the exponential function can be 55 * a serious matter. The error propagation involves 56 * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ), 57 * which shows that a 1 lsb error in representing X produces 58 * a relative error of X times 1 lsb in the function. 59 * While the routine gives an accurate result for arguments 60 * that are exactly represented by a long double precision 61 * computer number, the result contains amplified roundoff 62 * error for large arguments not exactly represented. 63 * 64 * 65 * ERROR MESSAGES: 66 * 67 * message condition value returned 68 * exp underflow x < MINLOG 0.0 69 * exp overflow x > MAXLOG MAXNUM 70 * 71 */ 72 73 /* Exponential function */ 74 75 #include <float.h> 76 #include <openlibm.h> 77 78 #include "math_private.h" 79 80 /* Pade' coefficients for exp(x) - 1 81 Theoretical peak relative error = 2.2e-37, 82 relative peak error spread = 9.2e-38 83 */ 84 static long double P[5] = { 85 3.279723985560247033712687707263393506266E-10L, 86 6.141506007208645008909088812338454698548E-7L, 87 2.708775201978218837374512615596512792224E-4L, 88 3.508710990737834361215404761139478627390E-2L, 89 9.999999999999999999999999999999999998502E-1L 90 }; 91 static long double Q[6] = { 92 2.980756652081995192255342779918052538681E-12L, 93 1.771372078166251484503904874657985291164E-8L, 94 1.504792651814944826817779302637284053660E-5L, 95 3.611828913847589925056132680618007270344E-3L, 96 2.368408864814233538909747618894558968880E-1L, 97 2.000000000000000000000000000000000000150E0L 98 }; 99 /* C1 + C2 = ln 2 */ 100 static const long double C1 = -6.93145751953125E-1L; 101 static const long double C2 = -1.428606820309417232121458176568075500134E-6L; 102 103 static const long double LOG2EL = 1.442695040888963407359924681001892137426646L; 104 static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L; 105 static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L; 106 static const long double huge = 0x1p10000L; 107 #if 0 /* XXX Prevent gcc from erroneously constant folding this. */ 108 static const long double twom10000 = 0x1p-10000L; 109 #else 110 static volatile long double twom10000 = 0x1p-10000L; 111 #endif 112 113 long double 114 expl(long double x) 115 { 116 long double px, xx; 117 int n; 118 119 if( x > MAXLOGL) 120 return (huge*huge); /* overflow */ 121 122 if( x < MINLOGL ) 123 return (twom10000*twom10000); /* underflow */ 124 125 /* Express e**x = e**g 2**n 126 * = e**g e**( n loge(2) ) 127 * = e**( g + n loge(2) ) 128 */ 129 px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */ 130 n = px; 131 x += px * C1; 132 x += px * C2; 133 /* rational approximation for exponential 134 * of the fractional part: 135 * e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) 136 */ 137 xx = x * x; 138 px = x * __polevll( xx, P, 4 ); 139 xx = __polevll( xx, Q, 5 ); 140 x = px/( xx - px ); 141 x = 1.0L + x + x; 142 143 x = ldexpl( x, n ); 144 return(x); 145 } 146