[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