1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 /* 13 * See comments in log10.c. 14 */ 15 16 #include <math.h> 17 #include <stdint.h> 18 19 static const float 20 ivln10hi = 4.3432617188e-01, /* 0x3ede6000 */ 21 ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */ 22 log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ 23 log10_2lo = 7.9034151668e-07, /* 0x355427db */ 24 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ 25 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ 26 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ 27 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ 28 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ 29 30 float log10f(float x) 31 { 32 union {float f; uint32_t i;} u = {x}; 33 float_t hfsq,f,s,z,R,w,t1,t2,dk,hi,lo; 34 uint32_t ix; 35 int k; 36 37 ix = u.i; 38 k = 0; 39 if (ix < 0x00800000 || ix>>31) { /* x < 2**-126 */ 40 if (ix<<1 == 0) 41 return -1/(x*x); /* log(+-0)=-inf */ 42 if (ix>>31) 43 return (x-x)/0.0f; /* log(-#) = NaN */ 44 /* subnormal number, scale up x */ 45 k -= 25; 46 x *= 0x1p25f; 47 u.f = x; 48 ix = u.i; 49 } else if (ix >= 0x7f800000) { 50 return x; 51 } else if (ix == 0x3f800000) 52 return 0; 53 54 /* reduce x into [sqrt(2)/2, sqrt(2)] */ 55 ix += 0x3f800000 - 0x3f3504f3; 56 k += (int)(ix>>23) - 0x7f; 57 ix = (ix&0x007fffff) + 0x3f3504f3; 58 u.i = ix; 59 x = u.f; 60 61 f = x - 1.0f; 62 s = f/(2.0f + f); 63 z = s*s; 64 w = z*z; 65 t1= w*(Lg2+w*Lg4); 66 t2= z*(Lg1+w*Lg3); 67 R = t2 + t1; 68 hfsq = 0.5f*f*f; 69 70 hi = f - hfsq; 71 u.f = hi; 72 u.i &= 0xfffff000; 73 hi = u.f; 74 lo = f - hi - hfsq + s*(hfsq+R); 75 dk = k; 76 return dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + dk*log10_2hi; 77 } 78