1//===-- Common utils for expf16 functions -----------------------*- 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_EXPF16_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "src/__support/CPP/array.h"
17#include "src/__support/FPUtil/PolyEval.h"
18#include "src/__support/FPUtil/nearest_integer.h"
19#include "src/__support/macros/properties/types.h"
20
21namespace LIBC_NAMESPACE_DECL {
22
23// Generated by Sollya with the following commands:
24// > display = hexadecimal;
25// > for i from -18 to 12 do print(round(exp(i), SG, RN));
26LIBC_INLINE_VAR constexpr cpp::array<float, 31> EXP_HI = {
27 0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f,
28 0x1.be6c7p-21f, 0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f,
29 0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f, 0x1.de16bap-11f,
30 0x1.44e52p-9f, 0x1.b993fep-8f, 0x1.2c155cp-6f, 0x1.97db0cp-5f,
31 0x1.152aaap-3f, 0x1.78b564p-2f, 0x1p+0f, 0x1.5bf0a8p+1f,
32 0x1.d8e64cp+2f, 0x1.415e5cp+4f, 0x1.b4c902p+5f, 0x1.28d38ap+7f,
33 0x1.936dc6p+8f, 0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f,
34 0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f,
35};
36
37// Generated by Sollya with the following commands:
38// > display = hexadecimal;
39// > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN));
40LIBC_INLINE_VAR constexpr cpp::array<float, 8> EXP_MID = {
41 0x1p+0f, 0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f,
42 0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f,
43};
44
45struct ExpRangeReduction {
46 float exp_hi_mid;
47 float exp_lo;
48};
49
50[[maybe_unused]] LIBC_INLINE ExpRangeReduction exp_range_reduction(float16 x) {
51 // For -18 < x < 12, to compute exp(x), we perform the following range
52 // reduction: find hi, mid, lo, such that:
53 // x = hi + mid + lo, in which
54 // hi is an integer,
55 // mid * 2^3 is an integer,
56 // -2^(-4) <= lo < 2^(-4).
57 // In particular,
58 // hi + mid = round(x * 2^3) * 2^(-3).
59 // Then,
60 // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
61 // We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID
62 // respectively. exp(lo) is computed using a degree-3 minimax polynomial
63 // generated by Sollya.
64
65 float xf = x;
66 float kf = fputil::nearest_integer(x: xf * 0x1.0p+3f);
67 int x_hi_mid = static_cast<int>(kf);
68 int x_hi = x_hi_mid >> 3;
69 int x_mid = x_hi_mid & 0x7;
70 // lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
71 float lo = fputil::multiply_add(x: kf, y: -0x1.0p-3f, z: xf);
72
73 float exp_hi = EXP_HI[x_hi + 18];
74 float exp_mid = EXP_MID[x_mid];
75 // Degree-3 minimax polynomial generated by Sollya with the following
76 // commands:
77 // > display = hexadecimal;
78 // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]);
79 // > 1 + x * P;
80 float exp_lo =
81 fputil::polyeval(x: lo, a0: 0x1p+0f, a: 0x1p+0f, a: 0x1.001p-1f, a: 0x1.555ddep-3f);
82 return {.exp_hi_mid: exp_hi * exp_mid, .exp_lo: exp_lo};
83}
84
85} // namespace LIBC_NAMESPACE_DECL
86
87#endif // LIBC_TYPES_HAS_FLOAT16
88
89#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H
90