1141296Sdas 2141296Sdas/* @(#)e_exp.c 1.6 04/04/22 */ 32116Sjkh/* 42116Sjkh * ==================================================== 5141296Sdas * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 72116Sjkh * Permission to use, copy, modify, and distribute this 8141296Sdas * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 13176451Sdas#include <sys/cdefs.h> 14176451Sdas__FBSDID("$FreeBSD$"); 152116Sjkh 162116Sjkh/* __ieee754_exp(x) 172116Sjkh * Returns the exponential of x. 182116Sjkh * 192116Sjkh * Method 202116Sjkh * 1. Argument reduction: 212116Sjkh * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 222116Sjkh * Given x, find r and integer k such that 232116Sjkh * 24141296Sdas * x = k*ln2 + r, |r| <= 0.5*ln2. 252116Sjkh * 26141296Sdas * Here r will be represented as r = hi-lo for better 272116Sjkh * accuracy. 282116Sjkh * 292116Sjkh * 2. Approximation of exp(r) by a special rational function on 302116Sjkh * the interval [0,0.34658]: 312116Sjkh * Write 322116Sjkh * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 33141296Sdas * We use a special Remes algorithm on [0,0.34658] to generate 34141296Sdas * a polynomial of degree 5 to approximate R. The maximum error 352116Sjkh * of this polynomial approximation is bounded by 2**-59. In 362116Sjkh * other words, 372116Sjkh * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 382116Sjkh * (where z=r*r, and the values of P1 to P5 are listed below) 392116Sjkh * and 402116Sjkh * | 5 | -59 41141296Sdas * | 2.0+P1*z+...+P5*z - R(z) | <= 2 422116Sjkh * | | 432116Sjkh * The computation of exp(r) thus becomes 442116Sjkh * 2*r 452116Sjkh * exp(r) = 1 + ------- 462116Sjkh * R - r 47141296Sdas * r*R1(r) 482116Sjkh * = 1 + r + ----------- (for better accuracy) 492116Sjkh * 2 - R1(r) 502116Sjkh * where 512116Sjkh * 2 4 10 522116Sjkh * R1(r) = r - (P1*r + P2*r + ... + P5*r ). 53141296Sdas * 542116Sjkh * 3. Scale back to obtain exp(x): 552116Sjkh * From step 1, we have 562116Sjkh * exp(x) = 2^k * exp(r) 572116Sjkh * 582116Sjkh * Special cases: 592116Sjkh * exp(INF) is INF, exp(NaN) is NaN; 602116Sjkh * exp(-INF) is 0, and 612116Sjkh * for finite argument, only exp(0)=1 is exact. 622116Sjkh * 632116Sjkh * Accuracy: 642116Sjkh * according to an error analysis, the error is always less than 652116Sjkh * 1 ulp (unit in the last place). 662116Sjkh * 672116Sjkh * Misc. info. 68141296Sdas * For IEEE double 692116Sjkh * if x > 7.09782712893383973096e+02 then exp(x) overflow 702116Sjkh * if x < -7.45133219101941108420e+02 then exp(x) underflow 712116Sjkh * 722116Sjkh * Constants: 73141296Sdas * The hexadecimal values are the intended ones for the following 74141296Sdas * constants. The decimal values may be used, provided that the 752116Sjkh * compiler will convert from decimal to binary accurately enough 762116Sjkh * to produce the hexadecimal values shown. 772116Sjkh */ 782116Sjkh 79226596Sdas#include <float.h> 80226596Sdas 812116Sjkh#include "math.h" 822116Sjkh#include "math_private.h" 832116Sjkh 842116Sjkhstatic const double 852116Sjkhone = 1.0, 862116SjkhhalF[2] = {0.5,-0.5,}, 872116Sjkho_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ 882116Sjkhu_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ 892116Sjkhln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 902116Sjkh -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ 912116Sjkhln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 922116Sjkh -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ 932116Sjkhinvln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ 942116SjkhP1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 952116SjkhP2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 962116SjkhP3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 972116SjkhP4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 982116SjkhP5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 992116Sjkh 100176373Sdasstatic volatile double 101251024Sdashuge = 1.0e+300, 102176373Sdastwom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/ 1032116Sjkh 10497407Salfreddouble 105117912Speter__ieee754_exp(double x) /* default IEEE double exp */ 1062116Sjkh{ 107176074Sbde double y,hi=0.0,lo=0.0,c,t,twopk; 10817141Sjkh int32_t k=0,xsb; 1092116Sjkh u_int32_t hx; 1102116Sjkh 1112116Sjkh GET_HIGH_WORD(hx,x); 1122116Sjkh xsb = (hx>>31)&1; /* sign bit of x */ 1132116Sjkh hx &= 0x7fffffff; /* high word of |x| */ 1142116Sjkh 1152116Sjkh /* filter out non-finite argument */ 1162116Sjkh if(hx >= 0x40862E42) { /* if |x|>=709.78... */ 1172116Sjkh if(hx>=0x7ff00000) { 1182116Sjkh u_int32_t lx; 1192116Sjkh GET_LOW_WORD(lx,x); 1208870Srgrimes if(((hx&0xfffff)|lx)!=0) 1212116Sjkh return x+x; /* NaN */ 1222116Sjkh else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 1232116Sjkh } 1242116Sjkh if(x > o_threshold) return huge*huge; /* overflow */ 1252116Sjkh if(x < u_threshold) return twom1000*twom1000; /* underflow */ 1262116Sjkh } 1272116Sjkh 1282116Sjkh /* argument reduction */ 129141296Sdas if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 1302116Sjkh if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 1312116Sjkh hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 1322116Sjkh } else { 133141296Sdas k = (int)(invln2*x+halF[xsb]); 1342116Sjkh t = k; 1352116Sjkh hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 1362116Sjkh lo = t*ln2LO[0]; 1372116Sjkh } 138226596Sdas STRICT_ASSIGN(double, x, hi - lo); 139141296Sdas } 1402116Sjkh else if(hx < 0x3e300000) { /* when |x|<2**-28 */ 1412116Sjkh if(huge+x>one) return one+x;/* trigger inexact */ 1422116Sjkh } 1432116Sjkh else k = 0; 1442116Sjkh 1452116Sjkh /* x is now in primary range */ 1462116Sjkh t = x*x; 147176074Sbde if(k >= -1021) 148176074Sbde INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0); 149176074Sbde else 150176074Sbde INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0); 1512116Sjkh c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 152141296Sdas if(k==0) return one-((x*c)/(c-2.0)-x); 1532116Sjkh else y = one-((lo-(x*c)/(2.0-c))-hi); 1542116Sjkh if(k >= -1021) { 155176074Sbde if (k==1024) return y*2.0*0x1p1023; 156176074Sbde return y*twopk; 1572116Sjkh } else { 158176074Sbde return y*twopk*twom1000; 1592116Sjkh } 1602116Sjkh} 161238722Skargl 162238722Skargl#if (LDBL_MANT_DIG == 53) 163238722Skargl__weak_reference(exp, expl); 164238722Skargl#endif 165