[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:22:34 PST 2024


https://github.com/lntue created https://github.com/llvm/llvm-project/pull/84391

None

>From 80ec4091500461b582d21ed383badaecadfe874d Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Thu, 7 Mar 2024 22:20:49 +0000
Subject: [PATCH] [libc][stdfix] Add exp function for short _Accum and _Accum
 types.

---
 libc/config/baremetal/arm/entrypoints.txt   |   2 +
 libc/config/baremetal/riscv/entrypoints.txt |   2 +
 libc/config/linux/x86_64/entrypoints.txt    |   2 +
 libc/docs/math/stdfix.rst                   |   2 +-
 libc/spec/llvm_libc_stdfix_ext.td           |   3 +
 libc/src/__support/fixed_point/fx_rep.h     |  24 ++---
 libc/src/stdfix/CMakeLists.txt              |  26 +++++
 libc/src/stdfix/exphk.cpp                   |  92 +++++++++++++++++
 libc/src/stdfix/exphk.h                     |  20 ++++
 libc/src/stdfix/expk.cpp                    | 104 ++++++++++++++++++++
 libc/src/stdfix/expk.h                      |  20 ++++
 libc/test/src/stdfix/CMakeLists.txt         |  36 +++++++
 libc/test/src/stdfix/ExpTest.h              |  77 +++++++++++++++
 libc/test/src/stdfix/exphk_test.cpp         |  13 +++
 libc/test/src/stdfix/expk_test.cpp          |  13 +++
 15 files changed, 423 insertions(+), 13 deletions(-)
 create mode 100644 libc/src/stdfix/exphk.cpp
 create mode 100644 libc/src/stdfix/exphk.h
 create mode 100644 libc/src/stdfix/expk.cpp
 create mode 100644 libc/src/stdfix/expk.h
 create mode 100644 libc/test/src/stdfix/ExpTest.h
 create mode 100644 libc/test/src/stdfix/exphk_test.cpp
 create mode 100644 libc/test/src/stdfix/expk_test.cpp

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 LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_STDFIX_EXPK_H
+#define LLVM_LIBC_SRC_STDFIX_EXPK_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+accum expk(accum x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_EXPK_H
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index d3e122884eb40e..74a1fb13127cc3 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -96,3 +96,39 @@ add_libc_test(
     libc.src.__support.FPUtil.basic_operations
     libc.src.__support.FPUtil.sqrt
 )
+
+add_libc_test(
+  exphk_test
+  SUITE
+    libc-stdfix-tests
+  HDRS
+    ExpTest.h
+  SRCS
+    exphk_test.cpp
+  COMPILE_OPTIONS
+    -O3
+  DEPENDS
+    libc.src.stdfix.exphk
+    libc.src.math.exp
+    libc.src.__support.CPP.bit
+    libc.src.__support.fixed_point.fx_rep
+    libc.src.__support.FPUtil.basic_operations
+)
+
+add_libc_test(
+  expk_test
+  SUITE
+    libc-stdfix-tests
+  HDRS
+    ExpTest.h
+  SRCS
+    expk_test.cpp
+  COMPILE_OPTIONS
+    -O3
+  DEPENDS
+    libc.src.stdfix.expk
+    libc.src.math.exp
+    libc.src.__support.CPP.bit
+    libc.src.__support.fixed_point.fx_rep
+    libc.src.__support.FPUtil.basic_operations
+)
diff --git a/libc/test/src/stdfix/ExpTest.h b/libc/test/src/stdfix/ExpTest.h
new file mode 100644
index 00000000000000..e588cebf621b90
--- /dev/null
+++ b/libc/test/src/stdfix/ExpTest.h
@@ -0,0 +1,77 @@
+//===-- Utility class to test integer sqrt ----------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/fixed_point/fx_rep.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+#include "src/math/exp.h"
+
+template <typename T> class ExpTest : public LIBC_NAMESPACE::testing::Test {
+
+  using FXRep = LIBC_NAMESPACE::fixed_point::FXRep<T>;
+  static constexpr T zero = FXRep::ZERO();
+  static constexpr T one = static_cast<T>(1);
+  static constexpr T eps = FXRep::EPS();
+
+public:
+  typedef T (*ExpFunc)(T);
+
+  void test_special_numbers(ExpFunc func) {
+    EXPECT_EQ(one, func(T(0)));
+    EXPECT_EQ(FXRep::MAX(), func(T(30)));
+    EXPECT_EQ(zero, func(T(-30)));
+  }
+
+  void test_range_with_step(ExpFunc func, T step, bool rel_error) {
+    constexpr int COUNT = 255;
+    constexpr double ERR = 3.0 * static_cast<double>(eps);
+    double x_d = 0.0;
+    T x = step;
+    for (int i = 0; i < COUNT; ++i) {
+      x += step;
+      x_d = static_cast<double>(x);
+      double y_d = static_cast<double>(func(x));
+      double result = LIBC_NAMESPACE::exp(x_d);
+      double errors = rel_error
+                          ? LIBC_NAMESPACE::fputil::abs((y_d / result) - 1.0)
+                          : LIBC_NAMESPACE::fputil::abs(y_d - result);
+      if (errors > ERR) {
+        // Print out the failure input and output.
+        EXPECT_EQ(x, T(0));
+        EXPECT_EQ(func(x), zero);
+      }
+      ASSERT_TRUE(errors <= ERR);
+    }
+  }
+
+  void test_positive_range(ExpFunc func) {
+    test_range_with_step(func, T(0x1.0p-6), /*rel_error*/ true);
+  }
+
+  void test_negative_range(ExpFunc func) {
+    test_range_with_step(func, T(-0x1.0p-6), /*rel_error*/ false);
+  }
+};
+
+#define LIST_EXP_TESTS(Name, T, func)                                          \
+  using LlvmLibcExp##Name##Test = ExpTest<T>;                                  \
+  TEST_F(LlvmLibcExp##Name##Test, SpecialNumbers) {                            \
+    test_special_numbers(&func);                                               \
+  }                                                                            \
+  TEST_F(LlvmLibcExp##Name##Test, PositiveRange) {                             \
+    test_positive_range(&func);                                                \
+  }                                                                            \
+  TEST_F(LlvmLibcExp##Name##Test, NegativeRange) {                             \
+    test_negative_range(&func);                                                \
+  }                                                                            \
+  static_assert(true, "Require semicolon.")
diff --git a/libc/test/src/stdfix/exphk_test.cpp b/libc/test/src/stdfix/exphk_test.cpp
new file mode 100644
index 00000000000000..24e92dc902faea
--- /dev/null
+++ b/libc/test/src/stdfix/exphk_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for exphk -----------------------------------------------===//
+//
+// 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 "ExpTest.h"
+
+#include "src/stdfix/exphk.h"
+
+LIST_EXP_TESTS(hk, short accum, LIBC_NAMESPACE::exphk);
diff --git a/libc/test/src/stdfix/expk_test.cpp b/libc/test/src/stdfix/expk_test.cpp
new file mode 100644
index 00000000000000..bc322037af04a7
--- /dev/null
+++ b/libc/test/src/stdfix/expk_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for expk ------------------------------------------------===//
+//
+// 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 "ExpTest.h"
+
+#include "src/stdfix/expk.h"
+
+LIST_EXP_TESTS(k, accum, LIBC_NAMESPACE::expk);



More information about the libc-commits mailing list