| 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 |