e_log.c revision 117912
12116Sjkh/* @(#)e_log.c 5.1 93/09/24 */ 22116Sjkh/* 32116Sjkh * ==================================================== 42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52116Sjkh * 62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 72116Sjkh * Permission to use, copy, modify, and distribute this 88870Srgrimes * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 132116Sjkh#ifndef lint 1450476Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/e_log.c 117912 2003-07-23 04:53:47Z peter $"; 152116Sjkh#endif 162116Sjkh 172116Sjkh/* __ieee754_log(x) 182116Sjkh * Return the logrithm of x 192116Sjkh * 208870Srgrimes * Method : 218870Srgrimes * 1. Argument Reduction: find k and f such that 228870Srgrimes * x = 2^k * (1+f), 232116Sjkh * where sqrt(2)/2 < 1+f < sqrt(2) . 242116Sjkh * 252116Sjkh * 2. Approximation of log(1+f). 262116Sjkh * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 272116Sjkh * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 282116Sjkh * = 2s + s*R 298870Srgrimes * We use a special Reme algorithm on [0,0.1716] to generate 308870Srgrimes * a polynomial of degree 14 to approximate R The maximum error 312116Sjkh * of this polynomial approximation is bounded by 2**-58.45. In 322116Sjkh * other words, 332116Sjkh * 2 4 6 8 10 12 14 342116Sjkh * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s 352116Sjkh * (the values of Lg1 to Lg7 are listed in the program) 362116Sjkh * and 372116Sjkh * | 2 14 | -58.45 388870Srgrimes * | Lg1*s +...+Lg7*s - R(z) | <= 2 392116Sjkh * | | 402116Sjkh * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 412116Sjkh * In order to guarantee error in log below 1ulp, we compute log 422116Sjkh * by 432116Sjkh * log(1+f) = f - s*(f - R) (if f is not too large) 442116Sjkh * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 458870Srgrimes * 468870Srgrimes * 3. Finally, log(x) = k*ln2 + log(1+f). 472116Sjkh * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 488870Srgrimes * Here ln2 is split into two floating point number: 492116Sjkh * ln2_hi + ln2_lo, 502116Sjkh * where n*ln2_hi is always exact for |n| < 2000. 512116Sjkh * 522116Sjkh * Special cases: 538870Srgrimes * log(x) is NaN with signal if x < 0 (including -INF) ; 542116Sjkh * log(+INF) is +INF; log(0) is -INF with signal; 552116Sjkh * log(NaN) is that NaN with no signal. 562116Sjkh * 572116Sjkh * Accuracy: 582116Sjkh * according to an error analysis, the error is always less than 592116Sjkh * 1 ulp (unit in the last place). 602116Sjkh * 612116Sjkh * Constants: 628870Srgrimes * The hexadecimal values are the intended ones for the following 638870Srgrimes * constants. The decimal values may be used, provided that the 648870Srgrimes * compiler will convert from decimal to binary accurately enough 652116Sjkh * to produce the hexadecimal values shown. 662116Sjkh */ 672116Sjkh 682116Sjkh#include "math.h" 692116Sjkh#include "math_private.h" 702116Sjkh 712116Sjkhstatic const double 722116Sjkhln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 732116Sjkhln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 742116Sjkhtwo54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ 752116SjkhLg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 762116SjkhLg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 772116SjkhLg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 782116SjkhLg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 792116SjkhLg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 802116SjkhLg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 812116SjkhLg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 822116Sjkh 832116Sjkhstatic const double zero = 0.0; 842116Sjkh 8597413Salfreddouble 86117912Speter__ieee754_log(double x) 872116Sjkh{ 882116Sjkh double hfsq,f,s,z,R,w,t1,t2,dk; 892116Sjkh int32_t k,hx,i,j; 902116Sjkh u_int32_t lx; 912116Sjkh 922116Sjkh EXTRACT_WORDS(hx,lx,x); 932116Sjkh 942116Sjkh k=0; 952116Sjkh if (hx < 0x00100000) { /* x < 2**-1022 */ 968870Srgrimes if (((hx&0x7fffffff)|lx)==0) 972116Sjkh return -two54/zero; /* log(+-0)=-inf */ 982116Sjkh if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 992116Sjkh k -= 54; x *= two54; /* subnormal number, scale up x */ 1002116Sjkh GET_HIGH_WORD(hx,x); 1018870Srgrimes } 1022116Sjkh if (hx >= 0x7ff00000) return x+x; 1032116Sjkh k += (hx>>20)-1023; 1042116Sjkh hx &= 0x000fffff; 1052116Sjkh i = (hx+0x95f64)&0x100000; 1062116Sjkh SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ 1072116Sjkh k += (i>>20); 1082116Sjkh f = x-1.0; 1092116Sjkh if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ 1102116Sjkh if(f==zero) if(k==0) return zero; else {dk=(double)k; 1112116Sjkh return dk*ln2_hi+dk*ln2_lo;} 1122116Sjkh R = f*f*(0.5-0.33333333333333333*f); 1132116Sjkh if(k==0) return f-R; else {dk=(double)k; 1142116Sjkh return dk*ln2_hi-((R-dk*ln2_lo)-f);} 1152116Sjkh } 1168870Srgrimes s = f/(2.0+f); 1172116Sjkh dk = (double)k; 1182116Sjkh z = s*s; 1192116Sjkh i = hx-0x6147a; 1202116Sjkh w = z*z; 1212116Sjkh j = 0x6b851-hx; 1228870Srgrimes t1= w*(Lg2+w*(Lg4+w*Lg6)); 1238870Srgrimes t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 1242116Sjkh i |= j; 1252116Sjkh R = t2+t1; 1262116Sjkh if(i>0) { 1272116Sjkh hfsq=0.5*f*f; 1282116Sjkh if(k==0) return f-(hfsq-s*(hfsq+R)); else 1292116Sjkh return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); 1302116Sjkh } else { 1312116Sjkh if(k==0) return f-s*(f-R); else 1322116Sjkh return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); 1332116Sjkh } 1342116Sjkh} 135