2  * Double-precision log2(x) function.
4  * Copyright (c) 2018-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
35   double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p;  in log2()  local
44 # define LO asuint64 (1.0 - 0x1.5b51p-5)  in log2()
45 # define HI asuint64 (1.0 + 0x1.6ab2p-5)  in log2()  macro
47   if (unlikely (ix - LO < HI - LO))  in log2()
53       r = x - 1.0;  in log2()
55       hi = r * InvLn2hi;  in log2()
56       lo = r * InvLn2lo + fma (r, InvLn2hi, -hi);  in log2()
59       rhi = asdouble (asuint64 (r) & -1ULL << 32);  in log2()
60       rlo = r - rhi;  in log2()
61       hi = rhi * InvLn2hi;  in log2()
64       r2 = r * r; /* rounding error: 0x1p-62.  */  in log2()
67       /* Worst-case error is less than 0.54 ULP (0.55 ULP without fma).  */  in log2()
69       y = hi + p;  in log2()
70       lo += hi - y + p;  in log2()
77   if (unlikely (top - 0x0010 >= 0x7ff0 - 0x0010))  in log2()
79       /* x < 0x1p-1022 or inf or nan.  */  in log2()
88       ix -= 52ULL << 52;  in log2()
91   /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.  in log2()
93      The ith subinterval contains z and c is near its center.  */  in log2()
94   tmp = ix - OFF;  in log2()
95   i = (tmp >> (52 - LOG2_TABLE_BITS)) % N;  in log2()
97   iz = ix - (tmp & 0xfffULL << 52);  in log2()
100   z = asdouble (iz);  in log2()
103   /* log2(x) = log2(z/c) + log2(c) + k.  */  in log2()
104   /* r ~= z/c - 1, |r| < 1/(2*N).  */  in log2()
106   /* rounding error: 0x1p-55/N.  */  in log2()
107   r = fma (z, invc, -1.0);  in log2()
109   t2 = r * InvLn2lo + fma (r, InvLn2hi, -t1);  in log2()
112   /* rounding error: 0x1p-55/N + 0x1p-65.  */  in log2()
113   r = (z - T2[i].chi - T2[i].clo) * invc;  in log2()
114   rhi = asdouble (asuint64 (r) & -1ULL << 32);  in log2()
115   rlo = r - rhi;  in log2()
120   /* hi + lo = r/ln2 + log2(c) + k.  */  in log2()
122   hi = t3 + t1;  in log2()
123   lo = t3 - hi + t1 + t2;  in log2()
127   r2 = r * r; /* rounding error: 0x1p-54/N^2.  */  in log2()
130   /* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma).  in log2()
131      ~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma).  */  in log2()
133   y = lo + r2 * p + hi;  in log2()
149 TEST_INTERVAL (log2, 0x1p-4, 0x1p4, 40000)