2 * Double-precision log(x) function.
4 * Copyright (c) 2018, Arm Limited.
5 * SPDX-License-Identifier: MIT
13 #define T __log_data.tab
14 #define T2 __log_data.tab2
15 #define B __log_data.poly1
16 #define A __log_data.poly
17 #define Ln2hi __log_data.ln2hi
18 #define Ln2lo __log_data.ln2lo
19 #define N (1 << LOG_TABLE_BITS)
20 #define OFF 0x3fe6000000000000
22 /* Top 16 bits of a double. */
23 static inline uint32_t top16(double x
)
25 return asuint64(x
) >> 48;
30 double_t w
, z
, r
, r2
, r3
, y
, invc
, logc
, kd
, hi
, lo
;
37 #define LO asuint64(1.0 - 0x1p-4)
38 #define HI asuint64(1.0 + 0x1.09p-4)
39 if (predict_false(ix
- LO
< HI
- LO
)) {
40 /* Handle close to 1.0 inputs separately. */
41 /* Fix sign of zero with downward rounding when x==1. */
42 if (WANT_ROUNDING
&& predict_false(ix
== asuint64(1.0)))
48 (B
[1] + r
* B
[2] + r2
* B
[3] +
49 r3
* (B
[4] + r
* B
[5] + r2
* B
[6] +
50 r3
* (B
[7] + r
* B
[8] + r2
* B
[9] + r3
* B
[10])));
51 /* Worst-case error is around 0.507 ULP. */
53 double_t rhi
= r
+ w
- w
;
54 double_t rlo
= r
- rhi
;
55 w
= rhi
* rhi
* B
[0]; /* B[0] == -0.5. */
58 lo
+= B
[0] * rlo
* (rhi
+ r
);
61 return eval_as_double(y
);
63 if (predict_false(top
- 0x0010 >= 0x7ff0 - 0x0010)) {
64 /* x < 0x1p-1022 or inf or nan. */
66 return __math_divzero(1);
67 if (ix
== asuint64(INFINITY
)) /* log(inf) == inf. */
69 if ((top
& 0x8000) || (top
& 0x7ff0) == 0x7ff0)
70 return __math_invalid(x
);
71 /* x is subnormal, normalize it. */
72 ix
= asuint64(x
* 0x1p
52);
76 /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
77 The range is split into N subintervals.
78 The ith subinterval contains z and c is near its center. */
80 i
= (tmp
>> (52 - LOG_TABLE_BITS
)) % N
;
81 k
= (int64_t)tmp
>> 52; /* arithmetic shift */
82 iz
= ix
- (tmp
& 0xfffULL
<< 52);
87 /* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */
88 /* r ~= z/c - 1, |r| < 1/(2*N). */
90 /* rounding error: 0x1p-55/N. */
91 r
= __builtin_fma(z
, invc
, -1.0);
93 /* rounding error: 0x1p-55/N + 0x1p-66. */
94 r
= (z
- T2
[i
].chi
- T2
[i
].clo
) * invc
;
98 /* hi + lo = r + log(c) + k*Ln2. */
99 w
= kd
* Ln2hi
+ logc
;
101 lo
= w
- hi
+ r
+ kd
* Ln2lo
;
103 /* log(x) = lo + (log1p(r) - r) + hi. */
104 r2
= r
* r
; /* rounding error: 0x1p-54/N^2. */
105 /* Worst case error if |y| > 0x1p-5:
106 0.5 + 4.13/N + abs-poly-error*2^57 ULP (+ 0.002 ULP without fma)
107 Worst case error if |y| > 0x1p-4:
108 0.5 + 2.06/N + abs-poly-error*2^56 ULP (+ 0.001 ULP without fma). */
110 r
* r2
* (A
[1] + r
* A
[2] + r2
* (A
[3] + r
* A
[4])) + hi
;
111 return eval_as_double(y
);