| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */ |
| 2 /* | 2 /* |
| 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
| 4 */ | 4 */ |
| 5 /* | 5 /* |
| 6 * ==================================================== | 6 * ==================================================== |
| 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 8 * | 8 * |
| 9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
| 11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
| 12 * is preserved. | 12 * is preserved. |
| 13 * ==================================================== | 13 * ==================================================== |
| 14 */ | 14 */ |
| 15 | 15 |
| 16 #include <math.h> | 16 #include <math.h> |
| 17 #include <stdint.h> | 17 #include <stdint.h> |
| 18 | 18 |
| 19 static const float | 19 static const float ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ |
| 20 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ | 20 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ |
| 21 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ | 21 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ |
| 22 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ | 22 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ |
| 23 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ | 23 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ |
| 24 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ | 24 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ |
| 25 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ | 25 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ |
| 26 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ | |
| 27 | 26 |
| 28 float logf(float x) | 27 float logf(float x) { |
| 29 { | 28 union { |
| 30 » union {float f; uint32_t i;} u = {x}; | 29 float f; |
| 31 » float_t hfsq,f,s,z,R,w,t1,t2,dk; | 30 uint32_t i; |
| 32 » uint32_t ix; | 31 } u = {x}; |
| 33 » int k; | 32 float_t hfsq, f, s, z, R, w, t1, t2, dk; |
| 33 uint32_t ix; |
| 34 int k; |
| 34 | 35 |
| 35 » ix = u.i; | 36 ix = u.i; |
| 36 » k = 0; | 37 k = 0; |
| 37 » if (ix < 0x00800000 || ix>>31) { /* x < 2**-126 */ | 38 if (ix < 0x00800000 || ix >> 31) { /* x < 2**-126 */ |
| 38 » » if (ix<<1 == 0) | 39 if (ix << 1 == 0) |
| 39 » » » return -1/(x*x); /* log(+-0)=-inf */ | 40 return -1 / (x * x); /* log(+-0)=-inf */ |
| 40 » » if (ix>>31) | 41 if (ix >> 31) |
| 41 » » » return (x-x)/0.0f; /* log(-#) = NaN */ | 42 return (x - x) / 0.0f; /* log(-#) = NaN */ |
| 42 » » /* subnormal number, scale up x */ | 43 /* subnormal number, scale up x */ |
| 43 » » k -= 25; | 44 k -= 25; |
| 44 » » x *= 0x1p25f; | 45 x *= 0x1p25f; |
| 45 » » u.f = x; | 46 u.f = x; |
| 46 » » ix = u.i; | 47 ix = u.i; |
| 47 » } else if (ix >= 0x7f800000) { | 48 } else if (ix >= 0x7f800000) { |
| 48 » » return x; | 49 return x; |
| 49 » } else if (ix == 0x3f800000) | 50 } else if (ix == 0x3f800000) |
| 50 » » return 0; | 51 return 0; |
| 51 | 52 |
| 52 » /* reduce x into [sqrt(2)/2, sqrt(2)] */ | 53 /* reduce x into [sqrt(2)/2, sqrt(2)] */ |
| 53 » ix += 0x3f800000 - 0x3f3504f3; | 54 ix += 0x3f800000 - 0x3f3504f3; |
| 54 » k += (int)(ix>>23) - 0x7f; | 55 k += (int)(ix >> 23) - 0x7f; |
| 55 » ix = (ix&0x007fffff) + 0x3f3504f3; | 56 ix = (ix & 0x007fffff) + 0x3f3504f3; |
| 56 » u.i = ix; | 57 u.i = ix; |
| 57 » x = u.f; | 58 x = u.f; |
| 58 | 59 |
| 59 » f = x - 1.0f; | 60 f = x - 1.0f; |
| 60 » s = f/(2.0f + f); | 61 s = f / (2.0f + f); |
| 61 » z = s*s; | 62 z = s * s; |
| 62 » w = z*z; | 63 w = z * z; |
| 63 » t1= w*(Lg2+w*Lg4); | 64 t1 = w * (Lg2 + w * Lg4); |
| 64 » t2= z*(Lg1+w*Lg3); | 65 t2 = z * (Lg1 + w * Lg3); |
| 65 » R = t2 + t1; | 66 R = t2 + t1; |
| 66 » hfsq = 0.5f*f*f; | 67 hfsq = 0.5f * f * f; |
| 67 » dk = k; | 68 dk = k; |
| 68 » return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; | 69 return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; |
| 69 } | 70 } |
| OLD | NEW |