-
Notifications
You must be signed in to change notification settings - Fork 25
/
e_log.c
174 lines (160 loc) · 5.23 KB
/
e_log.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
/* @(#)e_log.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_log(x)
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#ifndef __FDLIBM_H__
#include "fdlibm.h"
#endif
#ifndef __have_fpu_log
double __ieee754_log(double x)
{
double hfsq, f, s, z, R, w, t1, t2, dk;
int32_t k, hx, i, j;
uint32_t lx;
static const double ln2_hi = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
static const double ln2_lo = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */
static const double Lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
static const double Lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
static const double Lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
static const double Lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
static const double Lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
static const double Lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
static const double Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
GET_DOUBLE_WORDS(hx, lx, x);
k = 0;
if (hx < IC(0x00100000))
{ /* x < 2**-1022 */
if (((hx & IC(0x7fffffff)) | lx) == 0)
return -two54 / zero; /* log(+-0)=-inf */
if (hx < 0)
return (x - x) / zero; /* log(-#) = NaN */
k -= 54;
x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx, x);
}
if (hx >= IC(0x7ff00000))
return x + x;
k += (hx >> 20) - 1023;
hx &= IC(0x000fffff);
i = (hx + IC(0x95f64)) & IC(0x100000);
SET_HIGH_WORD(x, hx | (i ^ IC(0x3ff00000))); /* normalize x or x/2 */
k += (i >> 20);
f = x - 1.0;
if ((IC(0x000fffff) & (2 + hx)) < 3)
{ /* |f| < 2**-20 */
if (f == zero)
{
if (k == 0)
return zero;
dk = (double) k;
return dk * ln2_hi + dk * ln2_lo;
}
R = f * f * (0.5 - 0.33333333333333333 * f);
if (k == 0)
return f - R;
dk = (double) k;
return dk * ln2_hi - ((R - dk * ln2_lo) - f);
}
s = f / (2.0 + f);
dk = (double) k;
z = s * s;
i = hx - IC(0x6147a);
w = z * z;
j = IC(0x6b851) - hx;
t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
i |= j;
R = t2 + t1;
if (i > 0)
{
hfsq = 0.5 * f * f;
if (k == 0)
return f - (hfsq - s * (hfsq + R));
return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f);
}
if (k == 0)
return f - s * (f - R);
return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
}
#endif
/* wrapper log(x) */
double __log(double x)
{
if (_LIB_VERSION != _IEEE_ && islessequal(x, 0.0))
{
if (x == 0.0)
{
feraiseexcept(FE_DIVBYZERO);
return __kernel_standard(x, x, -HUGE_VAL, KMATHERR_LOG_ZERO); /* log(0) */
} else
{
feraiseexcept(FE_INVALID);
return __kernel_standard(x, x, __builtin_nan(""), KMATHERR_LOG_MINUS); /* log(x<0) */
}
}
return __ieee754_log(x);
}
__typeof(__log) log __attribute__((weak, alias("__log")));
#ifdef __NO_LONG_DOUBLE_MATH
__typeof(__logl) __logl __attribute__((alias("__log")));
__typeof(__logl) logl __attribute__((weak, alias("__log")));
#endif