OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log.c */ |
2 /* | 2 /* |
3 * ==================================================== | 3 * ==================================================== |
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 * | 5 * |
6 * Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 * Permission to use, copy, modify, and distribute this | 7 * Permission to use, copy, modify, and distribute this |
8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
9 * is preserved. | 9 * is preserved. |
10 * ==================================================== | 10 * ==================================================== |
(...skipping 45 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
56 * Constants: | 56 * Constants: |
57 * The hexadecimal values are the intended ones for the following | 57 * The hexadecimal values are the intended ones for the following |
58 * constants. The decimal values may be used, provided that the | 58 * constants. The decimal values may be used, provided that the |
59 * compiler will convert from decimal to binary accurately enough | 59 * compiler will convert from decimal to binary accurately enough |
60 * to produce the hexadecimal values shown. | 60 * to produce the hexadecimal values shown. |
61 */ | 61 */ |
62 | 62 |
63 #include <math.h> | 63 #include <math.h> |
64 #include <stdint.h> | 64 #include <stdint.h> |
65 | 65 |
66 static const double | 66 static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ |
67 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ | 67 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ |
68 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ | 68 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ |
69 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ | 69 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ |
70 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ | 70 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ |
71 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ | 71 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ |
72 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ | 72 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ |
73 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ | 73 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ |
74 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ | 74 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ |
75 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ | |
76 | 75 |
77 double log(double x) | 76 double log(double x) { |
78 { | 77 union { |
79 » union {double f; uint64_t i;} u = {x}; | 78 double f; |
80 » double_t hfsq,f,s,z,R,w,t1,t2,dk; | 79 uint64_t i; |
81 » uint32_t hx; | 80 } u = {x}; |
82 » int k; | 81 double_t hfsq, f, s, z, R, w, t1, t2, dk; |
| 82 uint32_t hx; |
| 83 int k; |
83 | 84 |
84 » hx = u.i>>32; | 85 hx = u.i >> 32; |
85 » k = 0; | 86 k = 0; |
86 » if (hx < 0x00100000 || hx>>31) { | 87 if (hx < 0x00100000 || hx >> 31) { |
87 » » if (u.i<<1 == 0) | 88 if (u.i << 1 == 0) |
88 » » » return -1/(x*x); /* log(+-0)=-inf */ | 89 return -1 / (x * x); /* log(+-0)=-inf */ |
89 » » if (hx>>31) | 90 if (hx >> 31) |
90 » » » return (x-x)/0.0; /* log(-#) = NaN */ | 91 return (x - x) / 0.0; /* log(-#) = NaN */ |
91 » » /* subnormal number, scale x up */ | 92 /* subnormal number, scale x up */ |
92 » » k -= 54; | 93 k -= 54; |
93 » » x *= 0x1p54; | 94 x *= 0x1p54; |
94 » » u.f = x; | 95 u.f = x; |
95 » » hx = u.i>>32; | 96 hx = u.i >> 32; |
96 » } else if (hx >= 0x7ff00000) { | 97 } else if (hx >= 0x7ff00000) { |
97 » » return x; | 98 return x; |
98 » } else if (hx == 0x3ff00000 && u.i<<32 == 0) | 99 } else if (hx == 0x3ff00000 && u.i << 32 == 0) |
99 » » return 0; | 100 return 0; |
100 | 101 |
101 » /* reduce x into [sqrt(2)/2, sqrt(2)] */ | 102 /* reduce x into [sqrt(2)/2, sqrt(2)] */ |
102 » hx += 0x3ff00000 - 0x3fe6a09e; | 103 hx += 0x3ff00000 - 0x3fe6a09e; |
103 » k += (int)(hx>>20) - 0x3ff; | 104 k += (int)(hx >> 20) - 0x3ff; |
104 » hx = (hx&0x000fffff) + 0x3fe6a09e; | 105 hx = (hx & 0x000fffff) + 0x3fe6a09e; |
105 » u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); | 106 u.i = (uint64_t)hx << 32 | (u.i & 0xffffffff); |
106 » x = u.f; | 107 x = u.f; |
107 | 108 |
108 » f = x - 1.0; | 109 f = x - 1.0; |
109 » hfsq = 0.5*f*f; | 110 hfsq = 0.5 * f * f; |
110 » s = f/(2.0+f); | 111 s = f / (2.0 + f); |
111 » z = s*s; | 112 z = s * s; |
112 » w = z*z; | 113 w = z * z; |
113 » t1 = w*(Lg2+w*(Lg4+w*Lg6)); | 114 t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); |
114 » t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); | 115 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); |
115 » R = t2 + t1; | 116 R = t2 + t1; |
116 » dk = k; | 117 dk = k; |
117 » return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; | 118 return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; |
118 } | 119 } |
OLD | NEW |