OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.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 SunPro, a Sun Microsystems, Inc. business. | 6 * Developed at SunPro, 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 37 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
48 * | 48 * |
49 * u = 1+x; | 49 * u = 1+x; |
50 * if(u==1.0) return x ; else | 50 * if(u==1.0) return x ; else |
51 * return log(u)*(x/(u-1.0)); | 51 * return log(u)*(x/(u-1.0)); |
52 * | 52 * |
53 * See HP-15C Advanced Functions Handbook, p.193. | 53 * See HP-15C Advanced Functions Handbook, p.193. |
54 */ | 54 */ |
55 | 55 |
56 #include "libm.h" | 56 #include "libm.h" |
57 | 57 |
58 static const double | 58 static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ |
59 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ | 59 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ |
60 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ | 60 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ |
61 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ | 61 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ |
62 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ | 62 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ |
63 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ | 63 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ |
64 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ | 64 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ |
65 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ | 65 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ |
66 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ | 66 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ |
67 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ | |
68 | 67 |
69 double log1p(double x) | 68 double log1p(double x) { |
70 { | 69 union { |
71 » union {double f; uint64_t i;} u = {x}; | 70 double f; |
72 » double_t hfsq,f,c,s,z,R,w,t1,t2,dk; | 71 uint64_t i; |
73 » uint32_t hx,hu; | 72 } u = {x}; |
74 » int k; | 73 double_t hfsq, f, c, s, z, R, w, t1, t2, dk; |
| 74 uint32_t hx, hu; |
| 75 int k; |
75 | 76 |
76 » hx = u.i>>32; | 77 hx = u.i >> 32; |
77 » k = 1; | 78 k = 1; |
78 » if (hx < 0x3fda827a || hx>>31) { /* 1+x < sqrt(2)+ */ | 79 if (hx < 0x3fda827a || hx >> 31) { /* 1+x < sqrt(2)+ */ |
79 » » if (hx >= 0xbff00000) { /* x <= -1.0 */ | 80 if (hx >= 0xbff00000) { /* x <= -1.0 */ |
80 » » » if (x == -1) | 81 if (x == -1) |
81 » » » » return x/0.0; /* log1p(-1) = -inf */ | 82 return x / 0.0; /* log1p(-1) = -inf */ |
82 » » » return (x-x)/0.0; /* log1p(x<-1) = NaN */ | 83 return (x - x) / 0.0; /* log1p(x<-1) = NaN */ |
83 » » } | 84 } |
84 » » if (hx<<1 < 0x3ca00000<<1) { /* |x| < 2**-53 */ | 85 if (hx << 1 < 0x3ca00000 << 1) { /* |x| < 2**-53 */ |
85 » » » /* underflow if subnormal */ | 86 /* underflow if subnormal */ |
86 » » » if ((hx&0x7ff00000) == 0) | 87 if ((hx & 0x7ff00000) == 0) |
87 » » » » FORCE_EVAL((float)x); | 88 FORCE_EVAL((float)x); |
88 » » » return x; | 89 return x; |
89 » » } | 90 } |
90 » » if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ | 91 if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ |
91 » » » k = 0; | 92 k = 0; |
92 » » » c = 0; | 93 c = 0; |
93 » » » f = x; | 94 f = x; |
94 » » } | 95 } |
95 » } else if (hx >= 0x7ff00000) | 96 } else if (hx >= 0x7ff00000) |
96 » » return x; | 97 return x; |
97 » if (k) { | 98 if (k) { |
98 » » u.f = 1 + x; | 99 u.f = 1 + x; |
99 » » hu = u.i>>32; | 100 hu = u.i >> 32; |
100 » » hu += 0x3ff00000 - 0x3fe6a09e; | 101 hu += 0x3ff00000 - 0x3fe6a09e; |
101 » » k = (int)(hu>>20) - 0x3ff; | 102 k = (int)(hu >> 20) - 0x3ff; |
102 » » /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ | 103 /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ |
103 » » if (k < 54) { | 104 if (k < 54) { |
104 » » » c = k >= 2 ? 1-(u.f-x) : x-(u.f-1); | 105 c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1); |
105 » » » c /= u.f; | 106 c /= u.f; |
106 » » } else | 107 } else |
107 » » » c = 0; | 108 c = 0; |
108 » » /* reduce u into [sqrt(2)/2, sqrt(2)] */ | 109 /* reduce u into [sqrt(2)/2, sqrt(2)] */ |
109 » » hu = (hu&0x000fffff) + 0x3fe6a09e; | 110 hu = (hu & 0x000fffff) + 0x3fe6a09e; |
110 » » u.i = (uint64_t)hu<<32 | (u.i&0xffffffff); | 111 u.i = (uint64_t)hu << 32 | (u.i & 0xffffffff); |
111 » » f = u.f - 1; | 112 f = u.f - 1; |
112 » } | 113 } |
113 » hfsq = 0.5*f*f; | 114 hfsq = 0.5 * f * f; |
114 » s = f/(2.0+f); | 115 s = f / (2.0 + f); |
115 » z = s*s; | 116 z = s * s; |
116 » w = z*z; | 117 w = z * z; |
117 » t1 = w*(Lg2+w*(Lg4+w*Lg6)); | 118 t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); |
118 » t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); | 119 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); |
119 » R = t2 + t1; | 120 R = t2 + t1; |
120 » dk = k; | 121 dk = k; |
121 » return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; | 122 return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; |
122 } | 123 } |
OLD | NEW |