1 /*
2 * Double-precision log2(x) function.
3 *
4 * Copyright (c) 2018-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include <float.h>
9 #include <math.h>
10 #include <stdint.h>
11 #include "math_config.h"
12 #include "test_defs.h"
13 #include "test_sig.h"
14
15 #define T __log2_data.tab
16 #define T2 __log2_data.tab2
17 #define B __log2_data.poly1
18 #define A __log2_data.poly
19 #define InvLn2hi __log2_data.invln2hi
20 #define InvLn2lo __log2_data.invln2lo
21 #define N (1 << LOG2_TABLE_BITS)
22 #define OFF 0x3fe6000000000000
23
24 /* Top 16 bits of a double. */
25 static inline uint32_t
top16(double x)26 top16 (double x)
27 {
28 return asuint64 (x) >> 48;
29 }
30
31 double
log2(double x)32 log2 (double x)
33 {
34 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
35 double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p;
36 uint64_t ix, iz, tmp;
37 uint32_t top;
38 int k, i;
39
40 ix = asuint64 (x);
41 top = top16 (x);
42
43 #if LOG2_POLY1_ORDER == 11
44 # define LO asuint64 (1.0 - 0x1.5b51p-5)
45 # define HI asuint64 (1.0 + 0x1.6ab2p-5)
46 #endif
47 if (unlikely (ix - LO < HI - LO))
48 {
49 /* Handle close to 1.0 inputs separately. */
50 /* Fix sign of zero with downward rounding when x==1. */
51 if (WANT_ROUNDING && unlikely (ix == asuint64 (1.0)))
52 return 0;
53 r = x - 1.0;
54 #if HAVE_FAST_FMA
55 hi = r * InvLn2hi;
56 lo = r * InvLn2lo + fma (r, InvLn2hi, -hi);
57 #else
58 double_t rhi, rlo;
59 rhi = asdouble (asuint64 (r) & -1ULL << 32);
60 rlo = r - rhi;
61 hi = rhi * InvLn2hi;
62 lo = rlo * InvLn2hi + r * InvLn2lo;
63 #endif
64 r2 = r * r; /* rounding error: 0x1p-62. */
65 r4 = r2 * r2;
66 #if LOG2_POLY1_ORDER == 11
67 /* Worst-case error is less than 0.54 ULP (0.55 ULP without fma). */
68 p = r2 * (B[0] + r * B[1]);
69 y = hi + p;
70 lo += hi - y + p;
71 lo += r4 * (B[2] + r * B[3] + r2 * (B[4] + r * B[5])
72 + r4 * (B[6] + r * B[7] + r2 * (B[8] + r * B[9])));
73 y += lo;
74 #endif
75 return eval_as_double (y);
76 }
77 if (unlikely (top - 0x0010 >= 0x7ff0 - 0x0010))
78 {
79 /* x < 0x1p-1022 or inf or nan. */
80 if (ix * 2 == 0)
81 return __math_divzero (1);
82 if (ix == asuint64 (INFINITY)) /* log(inf) == inf. */
83 return x;
84 if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0)
85 return __math_invalid (x);
86 /* x is subnormal, normalize it. */
87 ix = asuint64 (x * 0x1p52);
88 ix -= 52ULL << 52;
89 }
90
91 /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
92 The range is split into N subintervals.
93 The ith subinterval contains z and c is near its center. */
94 tmp = ix - OFF;
95 i = (tmp >> (52 - LOG2_TABLE_BITS)) % N;
96 k = (int64_t) tmp >> 52; /* arithmetic shift */
97 iz = ix - (tmp & 0xfffULL << 52);
98 invc = T[i].invc;
99 logc = T[i].logc;
100 z = asdouble (iz);
101 kd = (double_t) k;
102
103 /* log2(x) = log2(z/c) + log2(c) + k. */
104 /* r ~= z/c - 1, |r| < 1/(2*N). */
105 #if HAVE_FAST_FMA
106 /* rounding error: 0x1p-55/N. */
107 r = fma (z, invc, -1.0);
108 t1 = r * InvLn2hi;
109 t2 = r * InvLn2lo + fma (r, InvLn2hi, -t1);
110 #else
111 double_t rhi, rlo;
112 /* rounding error: 0x1p-55/N + 0x1p-65. */
113 r = (z - T2[i].chi - T2[i].clo) * invc;
114 rhi = asdouble (asuint64 (r) & -1ULL << 32);
115 rlo = r - rhi;
116 t1 = rhi * InvLn2hi;
117 t2 = rlo * InvLn2hi + r * InvLn2lo;
118 #endif
119
120 /* hi + lo = r/ln2 + log2(c) + k. */
121 t3 = kd + logc;
122 hi = t3 + t1;
123 lo = t3 - hi + t1 + t2;
124
125 /* log2(r+1) = r/ln2 + r^2*poly(r). */
126 /* Evaluation is optimized assuming superscalar pipelined execution. */
127 r2 = r * r; /* rounding error: 0x1p-54/N^2. */
128 r4 = r2 * r2;
129 #if LOG2_POLY_ORDER == 7
130 /* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma).
131 ~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma). */
132 p = A[0] + r * A[1] + r2 * (A[2] + r * A[3]) + r4 * (A[4] + r * A[5]);
133 y = lo + r2 * p + hi;
134 #endif
135 return eval_as_double (y);
136 }
137 #if USE_GLIBC_ABI
strong_alias(log2,__log2_finite)138 strong_alias (log2, __log2_finite)
139 hidden_alias (log2, __ieee754_log2)
140 # if LDBL_MANT_DIG == 53
141 long double log2l (long double x) { return log2 (x); }
142 # endif
143 #endif
144
145 TEST_SIG (S, D, 1, log2, 0.01, 11.1)
146 TEST_ULP (log2, 0.05)
147 TEST_ULP_NONNEAREST (log2, 0.5)
148 TEST_INTERVAL (log2, 0, 0xffff000000000000, 10000)
149 TEST_INTERVAL (log2, 0x1p-4, 0x1p4, 40000)
150 TEST_INTERVAL (log2, 0, inf, 40000)
151