[libc-commits] [libc] [libc][stdfix] Add exp function for short _Accum and _Accum types. (PR #84391)
via libc-commits
libc-commits at lists.llvm.org
Thu Mar 7 14:23:03 PST 2024
llvmbot wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libc
Author: None (lntue)
<details>
<summary>Changes</summary>
---
Patch is 25.92 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/84391.diff
15 Files Affected:
- (modified) libc/config/baremetal/arm/entrypoints.txt (+2)
- (modified) libc/config/baremetal/riscv/entrypoints.txt (+2)
- (modified) libc/config/linux/x86_64/entrypoints.txt (+2)
- (modified) libc/docs/math/stdfix.rst (+1-1)
- (modified) libc/spec/llvm_libc_stdfix_ext.td (+3)
- (modified) libc/src/__support/fixed_point/fx_rep.h (+12-12)
- (modified) libc/src/stdfix/CMakeLists.txt (+26)
- (added) libc/src/stdfix/exphk.cpp (+92)
- (added) libc/src/stdfix/exphk.h (+20)
- (added) libc/src/stdfix/expk.cpp (+104)
- (added) libc/src/stdfix/expk.h (+20)
- (modified) libc/test/src/stdfix/CMakeLists.txt (+36)
- (added) libc/test/src/stdfix/ExpTest.h (+77)
- (added) libc/test/src/stdfix/exphk_test.cpp (+13)
- (added) libc/test/src/stdfix/expk_test.cpp (+13)
``````````diff
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 99796ad5edf5d5..6e4fdb03626436 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -288,6 +288,8 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.absr
libc.src.stdfix.abslk
libc.src.stdfix.abslr
+ libc.src.stdfix.exphk
+ libc.src.stdfix.expk
libc.src.stdfix.roundhk
libc.src.stdfix.roundhr
libc.src.stdfix.roundk
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index 99796ad5edf5d5..6e4fdb03626436 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -288,6 +288,8 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.absr
libc.src.stdfix.abslk
libc.src.stdfix.abslr
+ libc.src.stdfix.exphk
+ libc.src.stdfix.expk
libc.src.stdfix.roundhk
libc.src.stdfix.roundhr
libc.src.stdfix.roundk
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 1f36f127e3c473..0b77a9e170aae1 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -483,6 +483,8 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.absr
libc.src.stdfix.abslk
libc.src.stdfix.abslr
+ libc.src.stdfix.exphk
+ libc.src.stdfix.expk
libc.src.stdfix.roundhk
libc.src.stdfix.roundhr
libc.src.stdfix.roundk
diff --git a/libc/docs/math/stdfix.rst b/libc/docs/math/stdfix.rst
index 5e39d5c01d1e53..d8dcb0cfa4c521 100644
--- a/libc/docs/math/stdfix.rst
+++ b/libc/docs/math/stdfix.rst
@@ -110,7 +110,7 @@ floating point types, but are not part of the ISO/IEC TR 18037:2008 spec.
+===============+================+=============+===============+============+================+=============+================+=============+===============+============+================+=============+
| cos | | | | | | | | | | | | |
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
-| exp | | | | | | | | | | | | |
+| exp | | | | | | | | |check| | | |check| | | |
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
| log | | | | | | | | | | | | |
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
diff --git a/libc/spec/llvm_libc_stdfix_ext.td b/libc/spec/llvm_libc_stdfix_ext.td
index 75bde47810a6be..7bc7ec5464081b 100644
--- a/libc/spec/llvm_libc_stdfix_ext.td
+++ b/libc/spec/llvm_libc_stdfix_ext.td
@@ -5,6 +5,9 @@ def LLVMLibcStdfixExt : StandardSpec<"llvm_libc_stdfix_ext"> {
[], // types
[], // enums
[ // functions
+ GuardedFunctionSpec<"exphk", RetValSpec<ShortAccumType>, [ArgSpec<ShortAccumType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+ GuardedFunctionSpec<"expk", RetValSpec<AccumType>, [ArgSpec<AccumType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+
GuardedFunctionSpec<"sqrtuhr", RetValSpec<UnsignedShortFractType>, [ArgSpec<UnsignedShortFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
GuardedFunctionSpec<"sqrtur", RetValSpec<UnsignedFractType>, [ArgSpec<UnsignedFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
GuardedFunctionSpec<"sqrtulr", RetValSpec<UnsignedLongFractType>, [ArgSpec<UnsignedLongFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
diff --git a/libc/src/__support/fixed_point/fx_rep.h b/libc/src/__support/fixed_point/fx_rep.h
index 042cd2b20714c6..f13640a6c01918 100644
--- a/libc/src/__support/fixed_point/fx_rep.h
+++ b/libc/src/__support/fixed_point/fx_rep.h
@@ -45,7 +45,7 @@ template <> struct FXRep<short fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return SFRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return SFRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return SFRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0HR; }
LIBC_INLINE static constexpr Type EPS() { return SFRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5HR; }
@@ -65,7 +65,7 @@ template <> struct FXRep<unsigned short fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return USFRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return USFRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return USFRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0UHR; }
LIBC_INLINE static constexpr Type EPS() { return USFRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5UHR; }
@@ -85,7 +85,7 @@ template <> struct FXRep<fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return FRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return FRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return FRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0R; }
LIBC_INLINE static constexpr Type EPS() { return FRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5R; }
@@ -105,7 +105,7 @@ template <> struct FXRep<unsigned fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return UFRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return UFRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return UFRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0UR; }
LIBC_INLINE static constexpr Type EPS() { return UFRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5UR; }
@@ -125,7 +125,7 @@ template <> struct FXRep<long fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return LFRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return LFRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return LFRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0LR; }
LIBC_INLINE static constexpr Type EPS() { return LFRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5LR; }
@@ -145,7 +145,7 @@ template <> struct FXRep<unsigned long fract> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return ULFRACT_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return ULFRACT_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return ULFRACT_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0ULR; }
LIBC_INLINE static constexpr Type EPS() { return ULFRACT_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5ULR; }
@@ -165,7 +165,7 @@ template <> struct FXRep<short accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return SACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return SACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return SACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0HK; }
LIBC_INLINE static constexpr Type EPS() { return SACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5HK; }
@@ -185,7 +185,7 @@ template <> struct FXRep<unsigned short accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return USACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return USACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return USACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0UHK; }
LIBC_INLINE static constexpr Type EPS() { return USACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5UHK; }
@@ -205,7 +205,7 @@ template <> struct FXRep<accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return ACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return ACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return ACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0K; }
LIBC_INLINE static constexpr Type EPS() { return ACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5K; }
@@ -225,7 +225,7 @@ template <> struct FXRep<unsigned accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return UACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return UACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return UACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0UK; }
LIBC_INLINE static constexpr Type EPS() { return UACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5UK; }
@@ -245,7 +245,7 @@ template <> struct FXRep<long accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return LACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return LACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return LACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0LK; }
LIBC_INLINE static constexpr Type EPS() { return LACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5LK; }
@@ -265,7 +265,7 @@ template <> struct FXRep<unsigned long accum> {
SIGN_LEN + INTEGRAL_LEN + FRACTION_LEN;
LIBC_INLINE static constexpr Type MIN() { return ULACCUM_MIN; }
- LIBC_INLINE static constexpr Type MAX() { return ULACCUM_MIN; }
+ LIBC_INLINE static constexpr Type MAX() { return ULACCUM_MAX; }
LIBC_INLINE static constexpr Type ZERO() { return 0.0ULK; }
LIBC_INLINE static constexpr Type EPS() { return ULACCUM_EPSILON; }
LIBC_INLINE static constexpr Type ONE_HALF() { return 0.5ULK; }
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 3a1cb66b7abcaf..10d76ae31349f9 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -67,3 +67,29 @@ add_entrypoint_object(
DEPENDS
libc.src.__support.fixed_point.sqrt
)
+
+add_entrypoint_object(
+ exphk
+ HDRS
+ exphk.h
+ SRCS
+ exphk.cpp
+ COMPILE_OPTIONS
+ -O3
+ DEPENDS
+ libc.src.__support.fixed_point.fx_rep
+ libc.src.__support.CPP.bit
+)
+
+add_entrypoint_object(
+ expk
+ HDRS
+ expk.h
+ SRCS
+ expk.cpp
+ COMPILE_OPTIONS
+ -O3
+ DEPENDS
+ libc.src.__support.fixed_point.fx_rep
+ libc.src.__support.CPP.bit
+)
diff --git a/libc/src/stdfix/exphk.cpp b/libc/src/stdfix/exphk.cpp
new file mode 100644
index 00000000000000..19a972b390c71b
--- /dev/null
+++ b/libc/src/stdfix/exphk.cpp
@@ -0,0 +1,92 @@
+//===-- Implementation of exphk function ----------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "exphk.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/fx_bits.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace {
+
+// Look up tables for exp(hi) and exp(mid).
+// Generated with Sollya:
+// > for i from 0 to 89 do {
+// hi = floor(i/8) - 5;
+// m = i/8 - floor(i/8) - 0.5;
+// e_hi = nearestint(exp(hi) * 2^7) * 2^-7;
+// e_mid = nearestint(exp(m) * 2^7) * 2^-7;
+// print(hi, e_hi, m, e_mid);
+// };
+// Notice that when i = 88 and 89, e_hi will overflow short accum range.
+static constexpr short accum EXP_HI[12] = {
+ 0x1.0p-7hk, 0x1.0p-6hk, 0x1.8p-5hk, 0x1.1p-3hk, 0x1.78p-2hk, 0x1.0p0hk,
+ 0x1.5cp1hk, 0x1.d9p2hk, 0x1.416p4hk, 0x1.b4dp5hk, 0x1.28d4p7hk, SACCUM_MAX,
+};
+
+static constexpr short accum EXP_MID[8] = {
+ 0x1.38p-1hk, 0x1.6p-1hk, 0x1.9p-1hk, 0x1.c4p-1hk,
+ 0x1.0p0hk, 0x1.22p0hk, 0x1.48p0hk, 0x1.74p0hk,
+};
+
+} // anonymous namespace
+
+LLVM_LIBC_FUNCTION(short accum, exphk, (short accum x)) {
+ using FXRep = fixed_point::FXRep<short accum>;
+ using StorageType = typename FXRep::StorageType;
+ // Output overflow
+ if (LIBC_UNLIKELY(x >= 0x1.64p2hk))
+ return FXRep::MAX();
+ // Lower bound where exp(x) -> 0:
+ // floor(log(2^-8) * 2^7) * 2^-7
+ if (LIBC_UNLIKELY(x <= -0x1.63p2hk))
+ return FXRep::ZERO();
+
+ // Current range of x:
+ // -0x1.628p2 <= x <= 0x1.638p2
+ // Range reduction:
+ // x = hi + mid + lo,
+ // where:
+ // hi is an integer
+ // mid * 2^3 is an integer
+ // |lo| <= 2^-4.
+ // Then exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo)
+ // ~ exp(hi) * exp(mid) * (1 + lo)
+ // with relative errors < |lo|^2 <= 2^-8.
+ // exp(hi) and exp(mid) are extracted from small lookup tables.
+
+ // Round-to-nearest 1/8, tie-to-(+Int):
+ constexpr short accum ONE_SIXTEENTH = 0x1.0p-4hk;
+ // x_rounded = floor(x + 1/16).
+ short accum x_rounded = ((x + ONE_SIXTEENTH) >> (FXRep::FRACTION_LEN - 3))
+ << (FXRep::FRACTION_LEN - 3);
+ short accum lo = x - x_rounded;
+
+ // Range of x_rounded:
+ // x_rounded >= floor((-0x1.628p2 + 0x1.0p-4) * 2^3) * 2^-3
+ // = -0x1.6p2 = -5.5
+ // To get the indices, we shift the values so that it start with 0.
+ // Range of indices: 0 <= indices <= 89
+ StorageType indices = cpp::bit_cast<StorageType>((x_rounded + 0x1.6p2hk) >>
+ (FXRep::FRACTION_LEN - 3));
+ // So we have the following relation:
+ // indices = (hi + mid + 44/8) * 8
+ // That implies:
+ // hi + mid = indices/8 - 5.5
+ // So for lookup tables, we can use the upper 4 bits to get:
+ // exp( floor(indices / 8) - 5 )
+ // and lower 3 bits for:
+ // exp( (indices - floor(indices)) - 0.5 )
+ short accum exp_hi = EXP_HI[indices >> 3];
+ short accum exp_mid = EXP_MID[indices & 0x7];
+ // exp(x) ~ exp(hi) * exp(mid) * (1 + lo);
+ return (exp_hi * (exp_mid * (0x1.0p0hk + lo)));
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/exphk.h b/libc/src/stdfix/exphk.h
new file mode 100644
index 00000000000000..da03bb76d53f53
--- /dev/null
+++ b/libc/src/stdfix/exphk.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for exphk -------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_STDFIX_EXPHK_H
+#define LLVM_LIBC_SRC_STDFIX_EXPHK_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+short accum exphk(short accum x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_EXPHK_H
diff --git a/libc/src/stdfix/expk.cpp b/libc/src/stdfix/expk.cpp
new file mode 100644
index 00000000000000..57227fd27769cc
--- /dev/null
+++ b/libc/src/stdfix/expk.cpp
@@ -0,0 +1,104 @@
+//===-- Implementation of expk function ----------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "expk.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/fx_bits.h"
+
+namespace LIBC_NAMESPACE {
+
+namespace {
+
+// Look up tables for exp(hi) and exp(mid).
+// Generated with Sollya:
+// > for i from 0 to 23 do {
+// hi = i - 11;
+// e_hi = nearestint(exp(hi) * 2^15) * 2^-15;
+// print(e_hi, "k,");
+// };
+static constexpr accum EXP_HI[24] = {
+ 0x1p-15k, 0x1p-15k, 0x1p-13k, 0x1.6p-12k,
+ 0x1.ep-11k, 0x1.44p-9k, 0x1.bap-8k, 0x1.2cp-6k,
+ 0x1.97cp-5k, 0x1.153p-3k, 0x1.78b8p-2k, 0x1p0k,
+ 0x1.5bf1p1k, 0x1.d8e68p2k, 0x1.415e6p4k, 0x1.b4c9p5k,
+ 0x1.28d388p7k, 0x1.936dc6p8k, 0x1.1228858p10k, 0x1.749ea7cp11k,
+ 0x1.fa7157cp12k, 0x1.5829dcf8p14k, 0x1.d3c4489p15k, ACCUM_MAX,
+};
+
+// Generated with Sollya:
+// > for i from 0 to 15 do {
+// m = i/16 - 0.0625;
+// e_m = nearestint(exp(m) * 2^15) * 2^-15;
+// print(e_m, "k,");
+// };
+static constexpr accum EXP_MID[16] = {
+ 0x1.e0fcp-1k, 0x1p0k, 0x1.1082p0k, 0x1.2216p0k,
+ 0x1.34ccp0k, 0x1.48b6p0k, 0x1.5deap0k, 0x1.747ap0k,
+ 0x1.8c8p0k, 0x1.a612p0k, 0x1.c14cp0k, 0x1.de46p0k,
+ 0x1.fd1ep0k, 0x1.0efap1k, 0x1.2074p1k, 0x1.330ep1k,
+};
+
+} // anonymous namespace
+
+LLVM_LIBC_FUNCTION(accum, expk, (accum x)) {
+ using FXRep = fixed_point::FXRep<accum>;
+ using StorageType = typename FXRep::StorageType;
+ // Output overflow
+ // > floor(log(2^16) * 2^15) * 2^-15
+ if (LIBC_UNLIKELY(x >= 0x1.62e4p3k))
+ return FXRep::MAX();
+ // Lower bound where exp(x) -> 0:
+ // floor(log(2^-16) * 2^15) * 2^-15
+ if (LIBC_UNLIKELY(x <= -0x1.62e44p3k))
+ return FXRep::ZERO();
+
+ // Current range of x:
+ // -0x1.62e4p3 <= x <= 0x1.62e3cp3
+ // Range reduction:
+ // x = hi + mid + lo,
+ // where:
+ // hi is an integer
+ // mid * 2^4 is an integer
+ // |lo| <= 2^-5.
+ // Then exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo)
+ // ~ exp(hi) * exp(mid) * (1 + lo + lo^2 / 2)
+ // with relative errors < |lo|^3/2 <= 2^-16.
+ // exp(hi) and exp(mid) are extracted from small lookup tables.
+
+ // Round-to-nearest 1/16, tie-to-(+Int):
+ constexpr accum ONE_THIRTY_SECOND = 0x1.0p-5k;
+ // x_rounded = floor(x + 1/16).
+ accum x_rounded = ((x + ONE_THIRTY_SECOND) >> (FXRep::FRACTION_LEN - 4))
+ << (FXRep::FRACTION_LEN - 4);
+ accum lo = x - x_rounded;
+
+ // Range of x_rounded:
+ // x_rounded >= floor((-0x1.62e4p3 + 0x1.0p-5) * 2^4) * 2^-4
+ // = -0x1.62p3 = -11.0625
+ // To get the indices, we shift the values so that it start with 0.
+ // Range of indices: 0 <= indices <= 355.
+ StorageType indices = cpp::bit_cast<StorageType>((x_rounded + 0x1.62p3k) >>
+ (FXRep::FRACTION_LEN - 4));
+ // So we have the following relation:
+ // indices = (hi + mid + 177/16) * 16
+ // That implies:
+ // hi + mid = indices/16 - 11.0625
+ // So for lookup tables, we can use the upper 4 bits to get:
+ // exp( floor(indices / 16) - 11 )
+ // and lower 4 bits for:
+ // exp( (indices - floor(indices)) - 0.0625 )
+ accum exp_hi = EXP_HI[indices >> 4];
+ accum exp_mid = EXP_MID[indices & 0xf];
+ // exp(x) ~ exp(hi) * exp(mid) * (1 + lo);
+ accum l1 = 0x1.0p0k + (lo >> 1); // = 1 + lo / 2
+ accum l2 = 0x1.0p0k + lo * l1; // = 1 + lo * (1 + lo / 2) = 1 + lo + lo^2/2
+ return (exp_hi * (exp_mid * l2));
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/expk.h b/libc/src/stdfix/expk.h
new file mode 100644
index 00000000000000..4526686a200b47
--- /dev/null
+++ b/libc/src/stdfix/expk.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for expk --------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH ...
[truncated]
``````````
</details>
https://github.com/llvm/llvm-project/pull/84391
More information about the libc-commits
mailing list