1 /* Double-precision log2(x) function.
2 Copyright (c) 2018 Arm Ltd. All rights reserved.
3
4 SPDX-License-Identifier: BSD-3-Clause
5
6 Redistribution and use in source and binary forms, with or without
7 modification, are permitted provided that the following conditions
8 are met:
9 1. Redistributions of source code must retain the above copyright
10 notice, this list of conditions and the following disclaimer.
11 2. Redistributions in binary form must reproduce the above copyright
12 notice, this list of conditions and the following disclaimer in the
13 documentation and/or other materials provided with the distribution.
14 3. The name of the company may not be used to endorse or promote
15 products derived from this software without specific prior written
16 permission.
17
18 THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS'' AND ANY EXPRESS OR IMPLIED
19 WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
20 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
22 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
23 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
24 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
25 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
26 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
27 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
28
29 #include "fdlibm.h"
30 #if !__OBSOLETE_MATH_DOUBLE
31
32 #include <math.h>
33 #include <stdint.h>
34 #include "math_config.h"
35
36 #define T __log2_data.tab
37 #define T2 __log2_data.tab2
38 #define B __log2_data.poly1
39 #define A __log2_data.poly
40 #define InvLn2hi __log2_data.invln2hi
41 #define InvLn2lo __log2_data.invln2lo
42 #define N (1 << LOG2_TABLE_BITS)
43 #define OFF 0x3fe6000000000000
44
45 /* Top 16 bits of a double. */
46 static inline uint32_t
top16(double x)47 top16 (double x)
48 {
49 return asuint64 (x) >> 48;
50 }
51
52 double
53 (log2) (double x)
54 {
55 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
56 double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p;
57 uint64_t ix, iz, tmp;
58 uint32_t top;
59 int k, i;
60
61 ix = asuint64 (x);
62 top = top16 (x);
63
64 #if LOG2_POLY1_ORDER == 11
65 # define LO asuint64 (1.0 - 0x1.5b51p-5)
66 # define HI asuint64 (1.0 + 0x1.6ab2p-5)
67 #endif
68 if (unlikely (ix - LO < HI - LO))
69 {
70 /* Handle close to 1.0 inputs separately. */
71 /* Fix sign of zero with downward rounding when x==1. */
72 if (WANT_ROUNDING && unlikely (ix == asuint64 (1.0)))
73 return 0;
74 r = x - 1.0;
75 #if _HAVE_FAST_FMA
76 hi = r * InvLn2hi;
77 lo = r * InvLn2lo + fma (r, InvLn2hi, -hi);
78 #else
79 double_t rhi, rlo;
80 rhi = asfloat64 (asuint64 (r) & -1ULL << 32);
81 rlo = r - rhi;
82 hi = rhi * InvLn2hi;
83 lo = rlo * InvLn2hi + r * InvLn2lo;
84 #endif
85 r2 = r * r; /* rounding error: 0x1p-62. */
86 r4 = r2 * r2;
87 #if LOG2_POLY1_ORDER == 11
88 /* Worst-case error is less than 0.54 ULP (0.55 ULP without fma). */
89 p = r2 * (B[0] + r * B[1]);
90 y = hi + p;
91 lo += hi - y + p;
92 lo += r4 * (B[2] + r * B[3] + r2 * (B[4] + r * B[5])
93 + r4 * (B[6] + r * B[7] + r2 * (B[8] + r * B[9])));
94 y += lo;
95 #endif
96 return y;
97 }
98 if (unlikely (top - 0x0010 >= 0x7ff0 - 0x0010))
99 {
100 /* x < 0x1p-1022 or inf or nan. */
101 if (ix * 2 == 0)
102 return __math_divzero (1);
103 if (ix == asuint64 ((double) INFINITY)) /* log(inf) == inf. */
104 return x;
105 if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0)
106 return __math_invalid (x);
107 /* x is subnormal, normalize it. */
108 ix = asuint64 (x * 0x1p52);
109 ix -= 52ULL << 52;
110 }
111
112 /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
113 The range is split into N subintervals.
114 The ith subinterval contains z and c is near its center. */
115 tmp = ix - OFF;
116 i = (tmp >> (52 - LOG2_TABLE_BITS)) % N;
117 k = (int64_t) tmp >> 52; /* arithmetic shift */
118 iz = ix - (tmp & 0xfffULL << 52);
119 invc = T[i].invc;
120 logc = T[i].logc;
121 z = asfloat64 (iz);
122 kd = (double_t) k;
123
124 /* log2(x) = log2(z/c) + log2(c) + k. */
125 /* r ~= z/c - 1, |r| < 1/(2*N). */
126 #if _HAVE_FAST_FMA
127 /* rounding error: 0x1p-55/N. */
128 r = fma (z, invc, -1.0);
129 t1 = r * InvLn2hi;
130 t2 = r * InvLn2lo + fma (r, InvLn2hi, -t1);
131 #else
132 double_t rhi, rlo;
133 /* rounding error: 0x1p-55/N + 0x1p-65. */
134 r = (z - T2[i].chi - T2[i].clo) * invc;
135 rhi = asfloat64 (asuint64 (r) & -1ULL << 32);
136 rlo = r - rhi;
137 t1 = rhi * InvLn2hi;
138 t2 = rlo * InvLn2hi + r * InvLn2lo;
139 #endif
140
141 /* hi + lo = r/ln2 + log2(c) + k. */
142 t3 = kd + logc;
143 hi = t3 + t1;
144 lo = t3 - hi + t1 + t2;
145
146 /* log2(r+1) = r/ln2 + r^2*poly(r). */
147 /* Evaluation is optimized assuming superscalar pipelined execution. */
148 r2 = r * r; /* rounding error: 0x1p-54/N^2. */
149 r4 = r2 * r2;
150 #if LOG2_POLY_ORDER == 7
151 /* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma).
152 ~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma). */
153 p = A[0] + r * A[1] + r2 * (A[2] + r * A[3]) + r4 * (A[4] + r * A[5]);
154 y = lo + r2 * p + hi;
155 #endif
156 return y;
157 }
158
159 _MATH_ALIAS_d_d(log2)
160
161 #endif
162
