11dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 21dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* @(#)e_exp.c 1.6 04/04/22 */ 31dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 41dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 51dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 61dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 71dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Permission to use, copy, modify, and distribute this 81dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * software is freely granted, provided that this notice 91dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * is preserved. 101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#ifndef lint 141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic char rcsid[] = "$FreeBSD: src/lib/msun/src/e_exp.c,v 1.10 2005/02/04 18:26:05 das Exp $"; 151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#endif 161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* __ieee754_exp(x) 181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Returns the exponential of x. 191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Method 211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 1. Argument reduction: 221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Given x, find r and integer k such that 241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x = k*ln2 + r, |r| <= 0.5*ln2. 261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Here r will be represented as r = hi-lo for better 281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * accuracy. 291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2. Approximation of exp(r) by a special rational function on 311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * the interval [0,0.34658]: 321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Write 331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * We use a special Remes algorithm on [0,0.34658] to generate 351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * a polynomial of degree 5 to approximate R. The maximum error 361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * of this polynomial approximation is bounded by 2**-59. In 371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * other words, 381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * (where z=r*r, and the values of P1 to P5 are listed below) 401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * and 411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * | 5 | -59 421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * | 2.0+P1*z+...+P5*z - R(z) | <= 2 431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * | | 441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The computation of exp(r) thus becomes 451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2*r 461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * exp(r) = 1 + ------- 471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * R - r 481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * r*R1(r) 491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * = 1 + r + ----------- (for better accuracy) 501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2 - R1(r) 511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * where 521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2 4 10 531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * R1(r) = r - (P1*r + P2*r + ... + P5*r ). 541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 3. Scale back to obtain exp(x): 561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * From step 1, we have 571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * exp(x) = 2^k * exp(r) 581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Special cases: 601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * exp(INF) is INF, exp(NaN) is NaN; 611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * exp(-INF) is 0, and 621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * for finite argument, only exp(0)=1 is exact. 631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Accuracy: 651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * according to an error analysis, the error is always less than 661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 1 ulp (unit in the last place). 671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Misc. info. 691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * For IEEE double 701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * if x > 7.09782712893383973096e+02 then exp(x) overflow 711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * if x < -7.45133219101941108420e+02 then exp(x) underflow 721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Constants: 741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The hexadecimal values are the intended ones for the following 751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * constants. The decimal values may be used, provided that the 761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * compiler will convert from decimal to binary accurately enough 771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * to produce the hexadecimal values shown. 781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math.h" 811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math_private.h" 821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic const double 841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectone = 1.0, 851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjecthalF[2] = {0.5,-0.5,}, 861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projecthuge = 1.0e+300, 871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projecttwom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ 881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projecto_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ 891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectu_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ 901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ 921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ 941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectinvln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ 951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectP1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectP2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectP3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectP4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source ProjectP5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 1001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectdouble 1031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project__ieee754_exp(double x) /* default IEEE double exp */ 1041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 1051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double y,hi=0.0,lo=0.0,c,t; 1061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project int32_t k=0,xsb; 1071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t hx; 1081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project GET_HIGH_WORD(hx,x); 1101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project xsb = (hx>>31)&1; /* sign bit of x */ 1111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx &= 0x7fffffff; /* high word of |x| */ 1121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* filter out non-finite argument */ 1141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx >= 0x40862E42) { /* if |x|>=709.78... */ 1151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx>=0x7ff00000) { 1161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t lx; 1171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project GET_LOW_WORD(lx,x); 1181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(((hx&0xfffff)|lx)!=0) 1191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return x+x; /* NaN */ 1201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 1211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x > o_threshold) return huge*huge; /* overflow */ 1231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x < u_threshold) return twom1000*twom1000; /* underflow */ 1241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* argument reduction */ 1271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 1281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 1291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 1301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 1311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project k = (int)(invln2*x+halF[xsb]); 1321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = k; 1331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 1341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lo = t*ln2LO[0]; 1351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x = hi - lo; 1371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else if(hx < 0x3e300000) { /* when |x|<2**-28 */ 1391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(huge+x>one) return one+x;/* trigger inexact */ 1401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else k = 0; 1421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* x is now in primary range */ 1441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = x*x; 1451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 1461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(k==0) return one-((x*c)/(c-2.0)-x); 1471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else y = one-((lo-(x*c)/(2.0-c))-hi); 1481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(k >= -1021) { 1491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t hy; 1501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project GET_HIGH_WORD(hy,y); 1511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */ 1521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return y; 1531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 1541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t hy; 1551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project GET_HIGH_WORD(hy,y); 1561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */ 1571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return y*twom1000; 1581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 160