1 /* ef_log.c -- float version of e_log.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 #include "fdlibm.h"
17
18 #if __OBSOLETE_MATH_FLOAT
19 static const float ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
20 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
21 two25 = 3.355443200e+07, /* 0x4c000000 */
22 Lg1 = 6.6666668653e-01, /* 3F2AAAAB */
23 Lg2 = 4.0000000596e-01, /* 3ECCCCCD */
24 Lg3 = 2.8571429849e-01, /* 3E924925 */
25 Lg4 = 2.2222198546e-01, /* 3E638E29 */
26 Lg5 = 1.8183572590e-01, /* 3E3A3325 */
27 Lg6 = 1.5313838422e-01, /* 3E1CD04F */
28 Lg7 = 1.4798198640e-01; /* 3E178897 */
29
30 static const float zero = 0.0;
31
32 float
logf(float x)33 logf(float x)
34 {
35 float hfsq, f, s, z, R, w, t1, t2, dk;
36 __int32_t k, ix, i, j;
37
38 GET_FLOAT_WORD(ix, x);
39
40 k = 0;
41 if (FLT_UWORD_IS_ZERO(ix & 0x7fffffff))
42 return __math_divzerof(1); /* log(+-0)=-inf */
43 if (ix < 0)
44 return __math_invalidf(x); /* log(-#) = NaN */
45 if (!FLT_UWORD_IS_FINITE(ix))
46 return x + x;
47 if (FLT_UWORD_IS_SUBNORMAL(ix)) {
48 k -= 25;
49 x *= two25; /* subnormal number, scale up x */
50 GET_FLOAT_WORD(ix, x);
51 }
52 k += (ix >> 23) - 127;
53 ix &= 0x007fffff;
54 i = (ix + (0x95f64 << 3)) & 0x800000;
55 SET_FLOAT_WORD(x, ix | (i ^ 0x3f800000)); /* normalize x or x/2 */
56 k += (i >> 23);
57 f = x - (float)1.0;
58 if ((0x007fffff & (15 + ix)) < 16) { /* |f| < 2**-20 */
59 if (f == zero) {
60 if (k == 0)
61 return zero;
62 else {
63 dk = (float)k;
64 return dk * ln2_hi + dk * ln2_lo;
65 }
66 }
67 R = f * f * ((float)0.5 - (float)0.33333333333333333 * f);
68 if (k == 0)
69 return f - R;
70 else {
71 dk = (float)k;
72 return dk * ln2_hi - ((R - dk * ln2_lo) - f);
73 }
74 }
75 s = f / ((float)2.0 + f);
76 dk = (float)k;
77 z = s * s;
78 i = ix - (0x6147a << 3);
79 w = z * z;
80 j = (0x6b851 << 3) - ix;
81 t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
82 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
83 i |= j;
84 R = t2 + t1;
85 if (i > 0) {
86 hfsq = (float)0.5 * f * f;
87 if (k == 0)
88 return f - (hfsq - s * (hfsq + R));
89 else
90 return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f);
91 } else {
92 if (k == 0)
93 return f - s * (f - R);
94 else
95 return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
96 }
97 }
98
99 _MATH_ALIAS_f_f(log)
100
101 #else
102 #include "../common/sf_log.c"
103 #endif /* __OBSOLETE_MATH_FLOAT */
104
