1/* origin: FreeBSD /usr/src/lib/msun/src/e_log.c */ 2/* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12/* log(x) 13 * Return the logarithm of x 14 * 15 * Method : 16 * 1. Argument Reduction: find k and f such that 17 * x = 2^k * (1+f), 18 * where sqrt(2)/2 < 1+f < sqrt(2) . 19 * 20 * 2. Approximation of log(1+f). 21 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 22 * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 23 * = 2s + s*R 24 * We use a special Remez algorithm on [0,0.1716] to generate 25 * a polynomial of degree 14 to approximate R The maximum error 26 * of this polynomial approximation is bounded by 2**-58.45. In 27 * other words, 28 * 2 4 6 8 10 12 14 29 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s 30 * (the values of Lg1 to Lg7 are listed in the program) 31 * and 32 * | 2 14 | -58.45 33 * | Lg1*s +...+Lg7*s - R(z) | <= 2 34 * | | 35 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 36 * In order to guarantee error in log below 1ulp, we compute log 37 * by 38 * log(1+f) = f - s*(f - R) (if f is not too large) 39 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 40 * 41 * 3. Finally, log(x) = k*ln2 + log(1+f). 42 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 43 * Here ln2 is split into two floating point number: 44 * ln2_hi + ln2_lo, 45 * where n*ln2_hi is always exact for |n| < 2000. 46 * 47 * Special cases: 48 * log(x) is NaN with signal if x < 0 (including -INF) ; 49 * log(+INF) is +INF; log(0) is -INF with signal; 50 * log(NaN) is that NaN with no signal. 51 * 52 * Accuracy: 53 * according to an error analysis, the error is always less than 54 * 1 ulp (unit in the last place). 55 * 56 * Constants: 57 * The hexadecimal values are the intended ones for the following 58 * constants. The decimal values may be used, provided that the 59 * compiler will convert from decimal to binary accurately enough 60 * to produce the hexadecimal values shown. 61 */ 62 63#include <math.h> 64#include <stdint.h> 65 66static const double 67ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 68ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 69Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 70Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 71Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 72Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 73Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 74Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 75Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 76 77double log(double x) 78{ 79 union {double f; uint64_t i;} u = {x}; 80 double_t hfsq,f,s,z,R,w,t1,t2,dk; 81 uint32_t hx; 82 int k; 83 84 hx = u.i>>32; 85 k = 0; 86 if (hx < 0x00100000 || hx>>31) { 87 if (u.i<<1 == 0) 88 return -1/(x*x); /* log(+-0)=-inf */ 89 if (hx>>31) 90 return (x-x)/0.0; /* log(-#) = NaN */ 91 /* subnormal number, scale x up */ 92 k -= 54; 93 x *= 0x1p54; 94 u.f = x; 95 hx = u.i>>32; 96 } else if (hx >= 0x7ff00000) { 97 return x; 98 } else if (hx == 0x3ff00000 && u.i<<32 == 0) 99 return 0; 100 101 /* reduce x into [sqrt(2)/2, sqrt(2)] */ 102 hx += 0x3ff00000 - 0x3fe6a09e; 103 k += (int)(hx>>20) - 0x3ff; 104 hx = (hx&0x000fffff) + 0x3fe6a09e; 105 u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); 106 x = u.f; 107 108 f = x - 1.0; 109 hfsq = 0.5*f*f; 110 s = f/(2.0+f); 111 z = s*s; 112 w = z*z; 113 t1 = w*(Lg2+w*(Lg4+w*Lg6)); 114 t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 115 R = t2 + t1; 116 dk = k; 117 return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; 118} 119