1 /* Double-precision 2^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 N (1 << EXP_TABLE_BITS)
37 #define Shift __exp_data.exp2_shift
38 #define T __exp_data.tab
39 #define C1 __exp_data.exp2_poly[0]
40 #define C2 __exp_data.exp2_poly[1]
41 #define C3 __exp_data.exp2_poly[2]
42 #define C4 __exp_data.exp2_poly[3]
43 #define C5 __exp_data.exp2_poly[4]
44 #define C6 __exp_data.exp2_poly[5]
45
46 /* Handle cases that may overflow or underflow when computing the result that
47 is scale*(1+TMP) without intermediate rounding. The bit representation of
48 scale is in SBITS, however it has a computed exponent that may have
49 overflown into the sign bit so that needs to be adjusted before using it as
50 a double. (int32_t)KI is the k used in the argument reduction and exponent
51 adjustment of scale, positive k here means the result may overflow and
52 negative k means the result may underflow. */
53 static inline double
specialcase(double_t tmp,uint64_t sbits,uint64_t ki)54 specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
55 {
56 double_t scale, y;
57
58 if ((ki & 0x80000000) == 0)
59 {
60 /* k > 0, the exponent of scale might have overflowed by 1. */
61 sbits -= 1ull << 52;
62 scale = asfloat64 (sbits);
63 y = 2 * (scale + scale * tmp);
64 return check_oflow (y);
65 }
66 /* k < 0, need special care in the subnormal range. */
67 sbits += 1022ull << 52;
68 scale = asfloat64 (sbits);
69 y = scale + scale * tmp;
70 if (y < 1.0)
71 {
72 /* Round y to the right precision before scaling it into the subnormal
73 range to avoid double rounding that can cause 0.5+E/2 ulp error where
74 E is the worst-case ulp error outside the subnormal range. So this
75 is only useful if the goal is better than 1 ulp worst-case error. */
76 double_t hi, lo;
77 lo = scale - y + scale * tmp;
78 hi = 1.0 + y;
79 lo = 1.0 - hi + y + lo;
80 y = eval_as_double (hi + lo) - 1.0;
81 /* Avoid -0.0 with downward rounding. */
82 if (WANT_ROUNDING && y == 0.0)
83 y = 0.0;
84 /* The underflow exception needs to be signaled explicitly. */
85 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
86 }
87 y = 0x1p-1022 * y;
88 return check_uflow (y);
89 }
90
91 /* Top 12 bits of a double (sign and exponent bits). */
92 static inline uint32_t
top12(double x)93 top12 (double x)
94 {
95 return asuint64 (x) >> 52;
96 }
97
98 double
exp2(double x)99 exp2 (double x)
100 {
101 uint32_t abstop;
102 uint64_t ki, idx, top, sbits;
103 /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
104 double_t kd, r, r2, scale, tail, tmp;
105
106 abstop = top12 (x) & 0x7ff;
107 if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
108 {
109 if (abstop - top12 (0x1p-54) >= 0x80000000)
110 /* Avoid spurious underflow for tiny x. */
111 /* Note: 0 is common input. */
112 return WANT_ROUNDING ? 1.0 + x : 1.0;
113 if (abstop >= top12 (1024.0))
114 {
115 if (asuint64 (x) == asuint64 ((double) -INFINITY))
116 return 0.0;
117 if (abstop >= top12 ((double) INFINITY))
118 return 1.0 + x;
119 if (!(asuint64 (x) >> 63))
120 return __math_oflow (0);
121 else if (asuint64 (x) >= asuint64 (-1075.0))
122 return __math_uflow (0);
123 }
124 if (2 * asuint64 (x) > 2 * asuint64 (928.0))
125 /* Large x is special cased below. */
126 abstop = 0;
127 }
128
129 /* exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)]. */
130 /* x = k/N + r, with int k and r in [-1/2N, 1/2N]. */
131 kd = eval_as_double (x + Shift);
132 ki = asuint64 (kd); /* k. */
133 kd -= Shift; /* k/N for int k. */
134 r = x - kd;
135 /* 2^(k/N) ~= scale * (1 + tail). */
136 idx = 2 * (ki % N);
137 top = ki << (52 - EXP_TABLE_BITS);
138 tail = asfloat64 (T[idx]);
139 /* This is only a valid scale when -1023*N < k < 1024*N. */
140 sbits = T[idx + 1] + top;
141 /* exp2(x) = 2^(k/N) * 2^r ~= scale + scale * (tail + 2^r - 1). */
142 /* Evaluation is optimized assuming superscalar pipelined execution. */
143 r2 = r * r;
144 /* Without fma the worst case error is 0.5/N ulp larger. */
145 /* Worst case error is less than 0.5+0.86/N+(abs poly error * 2^53) ulp. */
146 #if EXP2_POLY_ORDER == 4
147 tmp = tail + r * C1 + r2 * C2 + r * r2 * (C3 + r * C4);
148 #elif EXP2_POLY_ORDER == 5
149 tmp = tail + r * C1 + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
150 #elif EXP2_POLY_ORDER == 6
151 tmp = tail + r * C1 + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6);
152 #endif
153 if (unlikely (abstop == 0))
154 return specialcase (tmp, sbits, ki);
155 scale = asfloat64 (sbits);
156 /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-928, so there
157 is no spurious underflow here even without fma. */
158 return scale + scale * tmp;
159 }
160
161 _MATH_ALIAS_d_d(exp2)
162
163 #endif
164
