e_log.c revision 251292
1141296Sdas 2141296Sdas/* @(#)e_log.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 14176451Sdas#include <sys/cdefs.h> 15176451Sdas__FBSDID("$FreeBSD: head/lib/msun/src/e_log.c 251292 2013-06-03 09:14:31Z das $"); 162116Sjkh 172116Sjkh/* __ieee754_log(x) 182116Sjkh * Return the logrithm of x 192116Sjkh * 20141296Sdas * Method : 21141296Sdas * 1. Argument Reduction: find k and f such that 22141296Sdas * 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 29141296Sdas * We use a special Reme algorithm on [0,0.1716] to generate 30141296Sdas * 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 38141296Sdas * | 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) 45141296Sdas * 46141296Sdas * 3. Finally, log(x) = k*ln2 + log(1+f). 472116Sjkh * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 48141296Sdas * 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: 53141296Sdas * 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: 62141296Sdas * The hexadecimal values are the intended ones for the following 63141296Sdas * constants. The decimal values may be used, provided that the 64141296Sdas * compiler will convert from decimal to binary accurately enough 652116Sjkh * to produce the hexadecimal values shown. 662116Sjkh */ 672116Sjkh 68251292Sdas#include <float.h> 69251292Sdas 702116Sjkh#include "math.h" 712116Sjkh#include "math_private.h" 722116Sjkh 732116Sjkhstatic const double 742116Sjkhln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 752116Sjkhln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 762116Sjkhtwo54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ 772116SjkhLg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 782116SjkhLg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 792116SjkhLg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 802116SjkhLg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 812116SjkhLg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 822116SjkhLg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 832116SjkhLg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 842116Sjkh 852116Sjkhstatic const double zero = 0.0; 86251024Sdasstatic volatile double vzero = 0.0; 872116Sjkh 8897413Salfreddouble 89117912Speter__ieee754_log(double x) 902116Sjkh{ 912116Sjkh double hfsq,f,s,z,R,w,t1,t2,dk; 922116Sjkh int32_t k,hx,i,j; 932116Sjkh u_int32_t lx; 942116Sjkh 952116Sjkh EXTRACT_WORDS(hx,lx,x); 962116Sjkh 972116Sjkh k=0; 982116Sjkh if (hx < 0x00100000) { /* x < 2**-1022 */ 99141296Sdas if (((hx&0x7fffffff)|lx)==0) 100251024Sdas return -two54/vzero; /* log(+-0)=-inf */ 1012116Sjkh if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 1022116Sjkh k -= 54; x *= two54; /* subnormal number, scale up x */ 1032116Sjkh GET_HIGH_WORD(hx,x); 104141296Sdas } 1052116Sjkh if (hx >= 0x7ff00000) return x+x; 1062116Sjkh k += (hx>>20)-1023; 1072116Sjkh hx &= 0x000fffff; 1082116Sjkh i = (hx+0x95f64)&0x100000; 1092116Sjkh SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ 1102116Sjkh k += (i>>20); 1112116Sjkh f = x-1.0; 112170707Sbde if((0x000fffff&(2+hx))<3) { /* -2**-20 <= f < 2**-20 */ 113177711Sdas if(f==zero) { 114177711Sdas if(k==0) { 115177711Sdas return zero; 116177711Sdas } else { 117177711Sdas dk=(double)k; 118177711Sdas return dk*ln2_hi+dk*ln2_lo; 119177711Sdas } 120177711Sdas } 1212116Sjkh R = f*f*(0.5-0.33333333333333333*f); 1222116Sjkh if(k==0) return f-R; else {dk=(double)k; 1232116Sjkh return dk*ln2_hi-((R-dk*ln2_lo)-f);} 1242116Sjkh } 125141296Sdas s = f/(2.0+f); 1262116Sjkh dk = (double)k; 1272116Sjkh z = s*s; 1282116Sjkh i = hx-0x6147a; 1292116Sjkh w = z*z; 1302116Sjkh j = 0x6b851-hx; 131141296Sdas t1= w*(Lg2+w*(Lg4+w*Lg6)); 132141296Sdas t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 1332116Sjkh i |= j; 1342116Sjkh R = t2+t1; 1352116Sjkh if(i>0) { 1362116Sjkh hfsq=0.5*f*f; 1372116Sjkh if(k==0) return f-(hfsq-s*(hfsq+R)); else 1382116Sjkh return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); 1392116Sjkh } else { 1402116Sjkh if(k==0) return f-s*(f-R); else 1412116Sjkh return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); 1422116Sjkh } 1432116Sjkh} 144251292Sdas 145251292Sdas#if (LDBL_MANT_DIG == 53) 146251292Sdas__weak_reference(log, logl); 147251292Sdas#endif 148