k_log.h revision 226376
118334Speter 218334Speter/* @(#)e_log.c 1.3 95/01/18 */ 350397Sobrien/* 418334Speter * ==================================================== 518334Speter * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 618334Speter * 718334Speter * Developed at SunSoft, a Sun Microsystems, Inc. business. 818334Speter * Permission to use, copy, modify, and distribute this 918334Speter * software is freely granted, provided that this notice 1018334Speter * is preserved. 1118334Speter * ==================================================== 1218334Speter */ 1318334Speter 1418334Speter#include <sys/cdefs.h> 1518334Speter__FBSDID("$FreeBSD: head/lib/msun/src/k_log.h 226376 2011-10-15 05:23:28Z das $"); 1618334Speter 1718334Speter/* 1818334Speter * k_log1p(f): 1918334Speter * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. 2018334Speter * 2118334Speter * The following describes the overall strategy for computing 2218334Speter * logarithms in base e. The argument reduction and adding the final 2318334Speter * term of the polynomial are done by the caller for increased accuracy 2418334Speter * when different bases are used. 2518334Speter * 2618334Speter * Method : 2718334Speter * 1. Argument Reduction: find k and f such that 2818334Speter * x = 2^k * (1+f), 2918334Speter * where sqrt(2)/2 < 1+f < sqrt(2) . 3018334Speter * 3118334Speter * 2. Approximation of log(1+f). 3218334Speter * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 3318334Speter * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 3418334Speter * = 2s + s*R 3518334Speter * We use a special Reme algorithm on [0,0.1716] to generate 3618334Speter * a polynomial of degree 14 to approximate R The maximum error 3718334Speter * of this polynomial approximation is bounded by 2**-58.45. In 3818334Speter * other words, 3918334Speter * 2 4 6 8 10 12 14 4018334Speter * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s 4118334Speter * (the values of Lg1 to Lg7 are listed in the program) 4218334Speter * and 4318334Speter * | 2 14 | -58.45 4418334Speter * | Lg1*s +...+Lg7*s - R(z) | <= 2 4518334Speter * | | 4618334Speter * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 4718334Speter * In order to guarantee error in log below 1ulp, we compute log 4818334Speter * by 4918334Speter * log(1+f) = f - s*(f - R) (if f is not too large) 5018334Speter * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 5118334Speter * 5218334Speter * 3. Finally, log(x) = k*ln2 + log(1+f). 5318334Speter * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 5418334Speter * Here ln2 is split into two floating point number: 5518334Speter * ln2_hi + ln2_lo, 5618334Speter * where n*ln2_hi is always exact for |n| < 2000. 5718334Speter * 5818334Speter * Special cases: 5918334Speter * log(x) is NaN with signal if x < 0 (including -INF) ; 6018334Speter * log(+INF) is +INF; log(0) is -INF with signal; 6118334Speter * log(NaN) is that NaN with no signal. 6218334Speter * 6318334Speter * Accuracy: 6418334Speter * according to an error analysis, the error is always less than 6518334Speter * 1 ulp (unit in the last place). 6618334Speter * 6718334Speter * Constants: 6818334Speter * The hexadecimal values are the intended ones for the following 6918334Speter * constants. The decimal values may be used, provided that the 7018334Speter * compiler will convert from decimal to binary accurately enough 7118334Speter * to produce the hexadecimal values shown. 7218334Speter */ 7318334Speter 7418334Speterstatic const double 7518334SpeterLg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 7618334SpeterLg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 7718334SpeterLg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 7818334SpeterLg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 7918334SpeterLg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 8018334SpeterLg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 8118334SpeterLg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 8218334Speter 8318334Speter/* 8418334Speter * We always inline k_log1p(), since doing so produces a 8518334Speter * substantial performance improvement (~40% on amd64). 8618334Speter */ 8718334Speterstatic inline double 8818334Speterk_log1p(double f) 8918334Speter{ 9018334Speter double hfsq,s,z,R,w,t1,t2; 9118334Speter 9218334Speter s = f/(2.0+f); 9318334Speter z = s*s; 9418334Speter w = z*z; 9518334Speter t1= w*(Lg2+w*(Lg4+w*Lg6)); 9618334Speter t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 9718334Speter R = t2+t1; 9818334Speter hfsq=0.5*f*f; 9918334Speter return s*(hfsq+R); 10018334Speter} 10118334Speter