2  * Double-precision log(x) function.
4  * Copyright (c) 2018-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
22 #define OFF 0x3fe6000000000000
44 # define LO asuint64 (1.0 - 0x1p-5)  in log()
45 # define HI asuint64 (1.0 + 0x1.1p-5)  in log()
47 # define LO asuint64 (1.0 - 0x1p-4)  in log()
48 # define HI asuint64 (1.0 + 0x1.09p-4)  in log()
50   if (unlikely (ix - LO < HI - LO))  in log()
55 	return 0;  in log()
56       r = x - 1.0;  in log()
60       /* Worst-case error is around 0.516 ULP.  */  in log()
63       w = B[0] * r2; /* B[0] == -0.5.  */  in log()
65       y += r - hi + w;  in log()
68       /* Worst-case error is around 0.516 ULP.  */  in log()
72       w = B[0] * r2; /* B[0] == -0.5.  */  in log()
74       y += r - hi + w;  in log()
81       /* Worst-case error is around 0.532 ULP.  */  in log()
82       w = B[0] * r2; /* B[0] == -0.5.  */  in log()
84       y += r - hi + w;  in log()
87       /* Worst-case error is around 0.507 ULP.  */  in log()
88       w = r * 0x1p27;  in log()
89       double_t rhi = r + w - w;  in log()
90       double_t rlo = r - rhi;  in log()
91       w = rhi * rhi * B[0]; /* B[0] == -0.5.  */  in log()
93       lo = r - hi + w;  in log()
94       lo += B[0] * rlo * (rhi + r);  in log()
101   if (unlikely (top - 0x0010 >= 0x7ff0 - 0x0010))  in log()
103       /* x < 0x1p-1022 or inf or nan.  */  in log()
104       if (ix * 2 == 0)  in log()
108       if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0)  in log()
111       ix = asuint64 (x * 0x1p52);  in log()
112       ix -= 52ULL << 52;  in log()
118   tmp = ix - OFF;  in log()
119   i = (tmp >> (52 - LOG_TABLE_BITS)) % N;  in log()
121   iz = ix - (tmp & 0xfffULL << 52);  in log()
126   /* log(x) = log1p(z/c-1) + log(c) + k*Ln2.  */  in log()
127   /* r ~= z/c - 1, |r| < 1/(2*N).  */  in log()
129   /* rounding error: 0x1p-55/N.  */  in log()
130   r = fma (z, invc, -1.0);  in log()
132   /* rounding error: 0x1p-55/N + 0x1p-66.  */  in log()
133   r = (z - T2[i].chi - T2[i].clo) * invc;  in log()
140   lo = w - hi + r + kd * Ln2lo;  in log()
142   /* log(x) = lo + (log1p(r) - r) + hi.  */  in log()
143   r2 = r * r; /* rounding error: 0x1p-54/N^2.  */  in log()
144   /* Worst case error if |y| > 0x1p-5:  in log()
145      0.5 + 4.13/N + abs-poly-error*2^57 ULP (+ 0.002 ULP without fma)  in log()
146      Worst case error if |y| > 0x1p-4:  in log()
147      0.5 + 2.06/N + abs-poly-error*2^56 ULP (+ 0.001 ULP without fma).  */  in log()
149   y = lo + r2 * A[0] + r * r2 * (A[1] + r * A[2] + r2 * (A[3] + r * A[4])) + hi;  in log()
152       + r2 * (A[0] + r * A[1] + r2 * (A[2] + r * A[3])  in log()
169 TEST_INTERVAL (log, 0, 0xffff000000000000, 10000)
170 TEST_INTERVAL (log, 0x1p-4, 0x1p4, 400000)
171 TEST_INTERVAL (log, 0, inf, 400000)