1 //===-- Single-precision 2^x function -------------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #include "src/math/exp2f.h"
10 #include "src/__support/FPUtil/BasicOperations.h"
11 #include "src/__support/FPUtil/FEnvImpl.h"
12 #include "src/__support/FPUtil/FMA.h"
13 #include "src/__support/FPUtil/FPBits.h"
14 #include "src/__support/FPUtil/PolyEval.h"
15 #include "src/__support/common.h"
16 
17 #include <errno.h>
18 
19 namespace __llvm_libc {
20 
21 // Lookup table for 2^(m * 2^(-6)) with m = 0, ..., 63.
22 // Table is generated with Sollya as follow:
23 // > display = hexadecimal;
24 // > for i from 0 to 63 do { D(2^(i / 64)); };
25 static constexpr double EXP_M[64] = {
26     0x1.0000000000000p0, 0x1.02c9a3e778061p0, 0x1.059b0d3158574p0,
27     0x1.0874518759bc8p0, 0x1.0b5586cf9890fp0, 0x1.0e3ec32d3d1a2p0,
28     0x1.11301d0125b51p0, 0x1.1429aaea92de0p0, 0x1.172b83c7d517bp0,
29     0x1.1a35beb6fcb75p0, 0x1.1d4873168b9aap0, 0x1.2063b88628cd6p0,
30     0x1.2387a6e756238p0, 0x1.26b4565e27cddp0, 0x1.29e9df51fdee1p0,
31     0x1.2d285a6e4030bp0, 0x1.306fe0a31b715p0, 0x1.33c08b26416ffp0,
32     0x1.371a7373aa9cbp0, 0x1.3a7db34e59ff7p0, 0x1.3dea64c123422p0,
33     0x1.4160a21f72e2ap0, 0x1.44e086061892dp0, 0x1.486a2b5c13cd0p0,
34     0x1.4bfdad5362a27p0, 0x1.4f9b2769d2ca7p0, 0x1.5342b569d4f82p0,
35     0x1.56f4736b527dap0, 0x1.5ab07dd485429p0, 0x1.5e76f15ad2148p0,
36     0x1.6247eb03a5585p0, 0x1.6623882552225p0, 0x1.6a09e667f3bcdp0,
37     0x1.6dfb23c651a2fp0, 0x1.71f75e8ec5f74p0, 0x1.75feb564267c9p0,
38     0x1.7a11473eb0187p0, 0x1.7e2f336cf4e62p0, 0x1.82589994cce13p0,
39     0x1.868d99b4492edp0, 0x1.8ace5422aa0dbp0, 0x1.8f1ae99157736p0,
40     0x1.93737b0cdc5e5p0, 0x1.97d829fde4e50p0, 0x1.9c49182a3f090p0,
41     0x1.a0c667b5de565p0, 0x1.a5503b23e255dp0, 0x1.a9e6b5579fdbfp0,
42     0x1.ae89f995ad3adp0, 0x1.b33a2b84f15fbp0, 0x1.b7f76f2fb5e47p0,
43     0x1.bcc1e904bc1d2p0, 0x1.c199bdd85529cp0, 0x1.c67f12e57d14bp0,
44     0x1.cb720dcef9069p0, 0x1.d072d4a07897cp0, 0x1.d5818dcfba487p0,
45     0x1.da9e603db3285p0, 0x1.dfc97337b9b5fp0, 0x1.e502ee78b3ff6p0,
46     0x1.ea4afa2a490dap0, 0x1.efa1bee615a27p0, 0x1.f50765b6e4540p0,
47     0x1.fa7c1819e90d8p0,
48 };
49 
50 INLINE_FMA
51 LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
52   using FPBits = typename fputil::FPBits<float>;
53   FPBits xbits(x);
54 
55   uint32_t x_u = xbits.uintval();
56   uint32_t x_abs = x_u & 0x7fff'ffffU;
57 
58   // Exceptional values.
59   switch (x_u) {
60   case 0x3b42'9d37U: // x = 0x1.853a6ep-9f
61     if (fputil::get_round() == FE_TONEAREST)
62       return 0x1.00870ap+0f;
63     break;
64   case 0x3c02'a9adU: // x = 0x1.05535ap-7f
65     if (fputil::get_round() == FE_TONEAREST)
66       return 0x1.016b46p+0f;
67     break;
68   case 0x3ca6'6e26U: { // x = 0x1.4cdc4cp-6f
69     int round_mode = fputil::get_round();
70     if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
71       return 0x1.03a16ap+0f;
72     return 0x1.03a168p+0f;
73   }
74   case 0x3d92'a282U: // x = 0x1.254504p-4f
75     if (fputil::get_round() == FE_UPWARD)
76       return 0x1.0d0688p+0f;
77     return 0x1.0d0686p+0f;
78   case 0xbcf3'a937U: // x = -0x1.e7526ep-6f
79     if (fputil::get_round() == FE_TONEAREST)
80       return 0x1.f58d62p-1f;
81     break;
82   case 0xb8d3'd026U: // x = -0x1.a7a04cp-14f
83     if (fputil::get_round() == FE_TONEAREST)
84       return 0x1.fff6d2p-1f;
85     break;
86   }
87 
88   // // When |x| >= 128, |x| < 2^-25, or x is nan
89   if (unlikely(x_abs >= 0x4300'0000U || x_abs <= 0x3280'0000U)) {
90     // |x| < 2^-25
91     if (x_abs <= 0x3280'0000U) {
92       return 1.0f + x;
93     }
94     // x >= 128
95     if (!xbits.get_sign()) {
96       // x is finite
97       if (x_u < 0x7f80'0000U) {
98         int rounding = fputil::get_round();
99         if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
100           return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
101 
102         errno = ERANGE;
103       }
104       // x is +inf or nan
105       return x + static_cast<float>(FPBits::inf());
106     }
107     // x < -150
108     if (x_u >= 0xc316'0000U) {
109       // exp(-Inf) = 0
110       if (xbits.is_inf())
111         return 0.0f;
112       // exp(nan) = nan
113       if (xbits.is_nan())
114         return x;
115       if (fputil::get_round() == FE_UPWARD)
116         return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL));
117       if (x != 0.0f)
118         errno = ERANGE;
119       return 0.0f;
120     }
121   }
122   // For -150 <= x < 128, to compute 2^x, we perform the following range
123   // reduction: find hi, mid, lo such that:
124   //   x = hi + mid + lo, in which
125   //     hi is an integer,
126   //     mid * 2^6 is an integer
127   //     -2^(-7) <= lo < 2^-7.
128   // In particular,
129   //   hi + mid = round(x * 2^6) * 2^(-6).
130   // Then,
131   //   2^(x) = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
132   // Multiply by 2^hi is simply adding hi to the exponent field.  We store
133   // exp(mid) in the lookup tables EXP_M.  exp(lo) is computed using a degree-4
134   // minimax polynomial generated by Sollya.
135 
136   // x_hi = round(hi + mid).
137   // The default rounding mode for float-to-int conversion in C++ is
138   // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5
139   // before conversion.
140   int x_hi =
141       static_cast<int>(x * 0x1.0p+6f + (xbits.get_sign() ? -0.5f : 0.5f));
142   // For 2-complement integers, arithmetic right shift is the same as dividing
143   // by a power of 2 and then round down (toward negative infinity).
144   int e_hi = (x_hi >> 6) +
145              static_cast<int>(fputil::FloatProperties<double>::EXPONENT_BIAS);
146   fputil::FPBits<double> exp_hi(
147       static_cast<uint64_t>(e_hi)
148       << fputil::FloatProperties<double>::MANTISSA_WIDTH);
149   // mid = x_hi & 0x0000'003fU;
150   double exp_hi_mid = static_cast<double>(exp_hi) * EXP_M[x_hi & 0x3f];
151   // Subtract (hi + mid) from x to get lo.
152   x -= static_cast<float>(x_hi) * 0x1.0p-6f;
153   double xd = static_cast<double>(x);
154   // Degree-4 minimax polynomial generated by Sollya with the following
155   // commands:
156   //   > display = hexadecimal;
157   //   > Q = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-7, 2^-7]);
158   //   > Q;
159   double exp_lo =
160       fputil::polyeval(xd, 0x1p0, 0x1.62e42fefa2417p-1, 0x1.ebfbdff82f809p-3,
161                        0x1.c6b0b92131c47p-5, 0x1.3b2ab6fb568a3p-7);
162   double result = exp_hi_mid * exp_lo;
163   return static_cast<float>(result);
164 }
165 
166 } // namespace __llvm_libc
167