OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.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 * ==================================================== |
11 */ | 11 */ |
12 | 12 |
13 #include "libm.h" | 13 #include "libm.h" |
14 | 14 |
15 static const float | 15 static const float ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ |
16 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ | 16 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ |
17 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ | 17 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ |
18 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ | 18 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ |
19 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ | 19 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ |
20 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ | 20 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ |
21 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ | 21 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ |
22 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ | |
23 | 22 |
24 float log1pf(float x) | 23 float log1pf(float x) { |
25 { | 24 union { |
26 » union {float f; uint32_t i;} u = {x}; | 25 float f; |
27 » float_t hfsq,f,c,s,z,R,w,t1,t2,dk; | 26 uint32_t i; |
28 » uint32_t ix,iu; | 27 } u = {x}; |
29 » int k; | 28 float_t hfsq, f, c, s, z, R, w, t1, t2, dk; |
| 29 uint32_t ix, iu; |
| 30 int k; |
30 | 31 |
31 » ix = u.i; | 32 ix = u.i; |
32 » k = 1; | 33 k = 1; |
33 » if (ix < 0x3ed413d0 || ix>>31) { /* 1+x < sqrt(2)+ */ | 34 if (ix < 0x3ed413d0 || ix >> 31) { /* 1+x < sqrt(2)+ */ |
34 » » if (ix >= 0xbf800000) { /* x <= -1.0 */ | 35 if (ix >= 0xbf800000) { /* x <= -1.0 */ |
35 » » » if (x == -1) | 36 if (x == -1) |
36 » » » » return x/0.0f; /* log1p(-1)=+inf */ | 37 return x / 0.0f; /* log1p(-1)=+inf */ |
37 » » » return (x-x)/0.0f; /* log1p(x<-1)=NaN */ | 38 return (x - x) / 0.0f; /* log1p(x<-1)=NaN */ |
38 » » } | 39 } |
39 » » if (ix<<1 < 0x33800000<<1) { /* |x| < 2**-24 */ | 40 if (ix << 1 < 0x33800000 << 1) { /* |x| < 2**-24 */ |
40 » » » /* underflow if subnormal */ | 41 /* underflow if subnormal */ |
41 » » » if ((ix&0x7f800000) == 0) | 42 if ((ix & 0x7f800000) == 0) |
42 » » » » FORCE_EVAL(x*x); | 43 FORCE_EVAL(x * x); |
43 » » » return x; | 44 return x; |
44 » » } | 45 } |
45 » » if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ | 46 if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ |
46 » » » k = 0; | 47 k = 0; |
47 » » » c = 0; | 48 c = 0; |
48 » » » f = x; | 49 f = x; |
49 » » } | 50 } |
50 » } else if (ix >= 0x7f800000) | 51 } else if (ix >= 0x7f800000) |
51 » » return x; | 52 return x; |
52 » if (k) { | 53 if (k) { |
53 » » u.f = 1 + x; | 54 u.f = 1 + x; |
54 » » iu = u.i; | 55 iu = u.i; |
55 » » iu += 0x3f800000 - 0x3f3504f3; | 56 iu += 0x3f800000 - 0x3f3504f3; |
56 » » k = (int)(iu>>23) - 0x7f; | 57 k = (int)(iu >> 23) - 0x7f; |
57 » » /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ | 58 /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ |
58 » » if (k < 25) { | 59 if (k < 25) { |
59 » » » c = k >= 2 ? 1-(u.f-x) : x-(u.f-1); | 60 c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1); |
60 » » » c /= u.f; | 61 c /= u.f; |
61 » » } else | 62 } else |
62 » » » c = 0; | 63 c = 0; |
63 » » /* reduce u into [sqrt(2)/2, sqrt(2)] */ | 64 /* reduce u into [sqrt(2)/2, sqrt(2)] */ |
64 » » iu = (iu&0x007fffff) + 0x3f3504f3; | 65 iu = (iu & 0x007fffff) + 0x3f3504f3; |
65 » » u.i = iu; | 66 u.i = iu; |
66 » » f = u.f - 1; | 67 f = u.f - 1; |
67 » } | 68 } |
68 » s = f/(2.0f + f); | 69 s = f / (2.0f + f); |
69 » z = s*s; | 70 z = s * s; |
70 » w = z*z; | 71 w = z * z; |
71 » t1= w*(Lg2+w*Lg4); | 72 t1 = w * (Lg2 + w * Lg4); |
72 » t2= z*(Lg1+w*Lg3); | 73 t2 = z * (Lg1 + w * Lg3); |
73 » R = t2 + t1; | 74 R = t2 + t1; |
74 » hfsq = 0.5f*f*f; | 75 hfsq = 0.5f * f * f; |
75 » dk = k; | 76 dk = k; |
76 » return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; | 77 return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; |
77 } | 78 } |
OLD | NEW |