Skip to content

[libc][stdfix] Add exp function for short _Accum and _Accum types. #84391

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Mar 7, 2024

Conversation

lntue
Copy link
Contributor

@lntue lntue commented Mar 7, 2024

No description provided.

@llvmbot
Copy link
Member

llvmbot commented Mar 7, 2024

@llvm/pr-subscribers-libc

Author: None (lntue)

Changes

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 --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]

Copy link
Contributor

@michaelrj-google michaelrj-google left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

overall LGTM


} // anonymous namespace

LLVM_LIBC_FUNCTION(short accum, exphk, (short accum x)) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

in future it would probably be better to have a shared backing function for exp_fixed but given there are only two functions for the moment this seems fine.

@lntue lntue merged commit 14171b8 into llvm:main Mar 7, 2024
@nickdesaulniers
Copy link
Member

Hi @lntue it looks like this is failing in post submit: https://lab.llvm.org/buildbot/#/builders/225/builds/32368/steps/7/logs/stdio

@nickdesaulniers
Copy link
Member

Specifically, it's just the libc-x86_64-debian-dbg-runtimes-build build. Some of the actual values are off by one from the expected.

@lntue
Copy link
Contributor Author

lntue commented Mar 8, 2024

Hi @lntue it looks like this is failing in post submit: https://lab.llvm.org/buildbot/#/builders/225/builds/32368/steps/7/logs/stdio

Fix with #84411

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants