use #if LDBL_MANT_DIG == ... instead of custom LD80 etc macros
[libm] / src / math / logf.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
2 /*
3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
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 "libm.h"
17
18 static const float
19 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
20 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
21 two25  = 3.355443200e+07,  /* 0x4c000000 */
22 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
23 Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
24 Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
25 Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
26 Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
27
28 static const float zero = 0.0;
29
30 float logf(float x)
31 {
32         float hfsq,f,s,z,R,w,t1,t2,dk;
33         int32_t k,ix,i,j;
34
35         GET_FLOAT_WORD(ix, x);
36
37         k = 0;
38         if (ix < 0x00800000) {  /* x < 2**-126  */
39                 if ((ix & 0x7fffffff) == 0)
40                         return -two25/zero;  /* log(+-0)=-inf */
41                 if (ix < 0)
42                         return (x-x)/zero;   /* log(-#) = NaN */
43                 /* subnormal number, scale up x */
44                 k -= 25;
45                 x *= two25;
46                 GET_FLOAT_WORD(ix, x);
47         }
48         if (ix >= 0x7f800000)
49                 return x+x;
50         k += (ix>>23) - 127;
51         ix &= 0x007fffff;
52         i = (ix + (0x95f64<<3)) & 0x800000;
53         SET_FLOAT_WORD(x, ix|(i^0x3f800000));  /* normalize x or x/2 */
54         k += i>>23;
55         f = x - (float)1.0;
56         if ((0x007fffff & (0x8000 + ix)) < 0xc000) {  /* -2**-9 <= f < 2**-9 */
57                 if (f == zero) {
58                         if (k == 0)
59                                 return zero;
60                         dk = (float)k;
61                         return dk*ln2_hi + dk*ln2_lo;
62                 }
63                 R = f*f*((float)0.5 - (float)0.33333333333333333*f);
64                 if (k == 0)
65                         return f-R;
66                 dk = (float)k;
67                 return dk*ln2_hi - ((R-dk*ln2_lo)-f);
68         }
69         s = f/((float)2.0+f);
70         dk = (float)k;
71         z = s*s;
72         i = ix-(0x6147a<<3);
73         w = z*z;
74         j = (0x6b851<<3)-ix;
75         t1= w*(Lg2+w*Lg4);
76         t2= z*(Lg1+w*Lg3);
77         i |= j;
78         R = t2 + t1;
79         if (i > 0) {
80                 hfsq = (float)0.5*f*f;
81                 if (k == 0)
82                         return f - (hfsq-s*(hfsq+R));
83                 return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
84         } else {
85                 if (k == 0)
86                         return f - s*(f-R);
87                 return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
88         }
89 }