1//===-- Common utils for exp10f ---------------------------------*- C++ -*-===//
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#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H
11
12#include "src/__support/FPUtil/FPBits.h"
13#include "src/__support/FPUtil/PolyEval.h"
14#include "src/__support/FPUtil/nearest_integer.h"
15#include "src/__support/macros/config.h"
16
17namespace LIBC_NAMESPACE_DECL {
18
19struct ExpBase {
20 // Base = e
21 static constexpr int MID_BITS = 5;
22 static constexpr int MID_MASK = (1 << MID_BITS) - 1;
23 // log2(e) * 2^5
24 static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS);
25 // High and low parts of -log(2) * 2^(-5)
26 static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS);
27 static constexpr double M_LOGB_2_LO =
28 -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS);
29 // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya
30 // with:
31 // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN));
32 static constexpr int64_t EXP_2_MID[1 << MID_BITS] = {
33 0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f,
34 0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa,
35 0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715,
36 0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d,
37 0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429,
38 0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74,
39 0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db,
40 0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d,
41 0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c,
42 0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f,
43 0x3ffea4afa2a490da, 0x3fff50765b6e4540,
44 };
45
46 // Approximating e^dx with degree-5 minimax polynomial generated by Sollya:
47 // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]);
48 // Then:
49 // e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5.
50 static constexpr double COEFFS[4] = {
51 0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5,
52 0x1.11112a0e34bdbp-7};
53
54 LIBC_INLINE static double powb_lo(double dx) {
55 using fputil::multiply_add;
56 double dx2 = dx * dx;
57 double c0 = 1.0 + dx;
58 // c1 = COEFFS[0] + COEFFS[1] * dx
59 double c1 = multiply_add(x: dx, y: ExpBase::COEFFS[1], z: ExpBase::COEFFS[0]);
60 // c2 = COEFFS[2] + COEFFS[3] * dx
61 double c2 = multiply_add(x: dx, y: ExpBase::COEFFS[3], z: ExpBase::COEFFS[2]);
62 // r = c4 + c5 * dx^4
63 // = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7
64 return fputil::polyeval(x: dx2, a0: c0, a: c1, a: c2);
65 }
66};
67
68struct Exp10Base : public ExpBase {
69 // log2(10) * 2^5
70 static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS);
71 // High and low parts of -log10(2) * 2^(-5).
72 // Notice that since |x * log2(10)| < 150:
73 // |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13
74 // So when the FMA instructions are not available, in order for the product
75 // k * M_LOGB_2_HI
76 // to be exact, we only store the high part of log10(2) up to 38 bits
77 // (= 53 - 15) of precision.
78 // It is generated by Sollya with:
79 // > round(log10(2), 44, RN);
80 static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS);
81 // > round(log10(2) - 0x1.34413509f8p-2, D, RN);
82 static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS);
83
84 // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya:
85 // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]);
86 // Then:
87 // 10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
88 static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1,
89 0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0,
90 0x1.1429e74a98f43p-1};
91
92 LIBC_INLINE static double powb_lo(double dx) {
93 using fputil::multiply_add;
94 double dx2 = dx * dx;
95 // c0 = 1 + COEFFS[0] * dx
96 double c0 = multiply_add(x: dx, y: Exp10Base::COEFFS[0], z: 1.0);
97 // c1 = COEFFS[1] + COEFFS[2] * dx
98 double c1 = multiply_add(x: dx, y: Exp10Base::COEFFS[2], z: Exp10Base::COEFFS[1]);
99 // c2 = COEFFS[3] + COEFFS[4] * dx
100 double c2 = multiply_add(x: dx, y: Exp10Base::COEFFS[4], z: Exp10Base::COEFFS[3]);
101 // r = c0 + dx^2 * (c1 + c2 * dx^2)
102 // = c0 + c1 * dx^2 + c2 * dx^4
103 // = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
104 return fputil::polyeval(x: dx2, a0: c0, a: c1, a: c2);
105 }
106};
107
108// Output of range reduction for exp_b: (2^(mid + hi), lo)
109// where:
110// b^x = 2^(mid + hi) * b^lo
111struct exp_b_reduc_t {
112 double mh; // 2^(mid + hi)
113 double lo;
114};
115
116// The function correctly calculates b^x value with at least float precision
117// in a limited range.
118// Range reduction:
119// b^x = 2^(hi + mid) * b^lo
120// where:
121// x = (hi + mid) * log_b(2) + lo
122// hi is an integer,
123// 0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer
124// -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1)
125// Base class needs to provide the following constants:
126// - MID_BITS : number of bits after decimal points used for mid
127// - MID_MASK : 2^MID_BITS - 1, mask to extract mid bits
128// - LOG2_B : log2(b) * 2^MID_BITS for scaling
129// - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS)
130// - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS)
131// - EXP_2_MID : look up table for bit fields of 2^mid
132// Return:
133// { 2^(hi + mid), lo }
134template <class Base>
135LIBC_INLINE constexpr exp_b_reduc_t exp_b_range_reduc(float x) {
136 double xd = static_cast<double>(x);
137 // kd = round((hi + mid) * log2(b) * 2^MID_BITS)
138 double kd = fputil::nearest_integer(Base::LOG2_B * xd);
139 // k = round((hi + mid) * log2(b) * 2^MID_BITS)
140 int k = static_cast<int>(kd);
141 // hi = floor(kd * 2^(-MID_BITS))
142 // exp_hi = shift hi to the exponent field of double precision.
143 uint64_t exp_hi = static_cast<uint64_t>(k >> Base::MID_BITS)
144 << fputil::FPBits<double>::FRACTION_LEN;
145 // mh = 2^hi * 2^mid
146 // mh_bits = bit field of mh
147 uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi;
148 double mh = fputil::FPBits<double>(mh_bits).get_val();
149 // dx = lo = x - (hi + mid) * log(2)
150 double dx = fputil::multiply_add(
151 kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd));
152 return {.mh: mh, .lo: dx};
153}
154
155} // namespace LIBC_NAMESPACE_DECL
156
157#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_UTILS_H
158