[libc][math] Universal exp function for cosh/sinh calculation.

[orex]

Summary

Added a function and test, which can be used later for cosh/sinh
and possibly for expf/expm1f.

latest sinhf/coshf performance:

CORE_MATH_PERF_MODE=rdtsc ./perf.sh coshf
GNU libc version: 2.31
GNU libc release: stable
30.155
12.267
13.525

CORE_MATH_PERF_MODE=rdtsc ./perf.sh sinhf
GNU libc version: 2.31
GNU libc release: stable
13.390
38.197
13.507

[lntue]

libc/src/math/generic/expxf.h

33

Use hexadecimal float constant.

49

static inline since this is in a header file.

50–54

You might want to try to see if using fputil::nearest_integer can improve the throughput or not:

double ps = fputil::nearest_integer(fputil::multiply_add(LN2_INV, x, static_cast<double>(ADD_TO_POS)));
int64_t ps_shifted = static_cast<int64_t>(ps) + (1 << (EXP_bits_p - 1));
...

67–68

You might want to use a fix representation of 1.0/6.0, ... 1.0/720.0 with hexadecimal floating point literals to reduce an uncertainty due to compiler constant evaluation methods and rounding modes.

[lntue]

Can you add add_header_library for this in src/math/generic/CMakeLists.txt and use it as dependency for sinhf and coshf?

[lntue]

libc/src/math/generic/expxf.h

46

Can you add the range restriction for x in here, since it does not work for the whole double precision range?

58

Since in both sinhf and coshf, you have: 0.5 * ep, maybe you can fuse that * 0.5 here by:

bs.set_unbiased_exponent(fputil::FPBits<double>::EXPONENT_BIAS - 1 + dg)

If so, don't forget to update the comment about the returned value of mult_exp.

[orex]

Small fixes.

libc/src/math/generic/expxf.h

58

But, for tanh, I think that will be a problem. I introduced more general implementation.

[lntue]

libc/test/src/math/expxf_test.cpp

37

Does it work for all rounding modes?

[orex]

Small fixes in test.

[orex]

libc/test/src/math/expxf_test.cpp

37

This is an universal function, therefore it can have some exception and ranges, which should be treated separately in each math function, which use that one. Paul (@zimmermann6) independently confirm, that derivative functions gives perfect result for all rounding modes. My tests form others reviews shows the same.

I put all rounding modes to avoid such question in the future.

libc/src/math/generic/expxf.h

+//===-- Single-precision general exp 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
+
+#include "common_constants.h" // Lookup tables EXP_M
+#include "math_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include <src/__support/FPUtil/NearestIntegerOperations.h>
+
+#include <errno.h>
+
+namespace __llvm_libc {
+
+// The algorithm represents exp(x) as
+//   exp(x) = 2^(ln(2) * i) * 2^(ln(2) * j / NUM_P )) * exp(dx)
+// where i integer value, j integer in range [-NUM_P/2, NUM_P/2).
+// 2^(ln(2) * j / NUM_P )) is a table values: 1.0 + EXP_M
+// exp(dx) calculates by taylor expansion.
+
+// Inversion of ln(2). Multiplication by EXP_num_p due to sampling by 1 /
+// EXP_num_p Precise value of the constant is not needed.
+static constexpr double LN2_INV = 0x1.71547652b82fep+0 * EXP_num_p;
+
+// LN2_HIGH + LN2_LOW = ln(2) with precision higher than double(ln(2))
+// Minus sign is to use FMA directly.
+static constexpr double LN2_HIGH = -0x1.62e42fefa0000p-1 / EXP_num_p;
+static constexpr double LN2_LOW = -0x1.cf79abc9e3b3ap-40 / EXP_num_p;
+
+struct exe_eval_result_t {
+  // exp(x) = 2^MULT_POWER2 * mult_exp * (r + 1.0)
+  // where
+  //   MULT_POWER2 template parameter;
+  //   mult_exp = 2^e;
+  //   r in range [~-0.3, ~0.41]
+  double mult_exp;
+  double r;
+};
+
+// The function correctly calculates exp value with at least float precision
+// in range not narrow than [-log(2^-150), 90]
+template <int MULT_POWER2 = 0>
+inline static exe_eval_result_t exp_eval(double x) {
+  double ps_dbl = fputil::nearest_integer(LN2_INV * x);
+  // Negative sign due to multiply_add optimization
+  double mult_e1, ml;
+  {
+    int ps =
+        static_cast<int>(ps_dbl) + (1 << (EXP_bits_p - 1)) +
+        ((fputil::FPBits<double>::EXPONENT_BIAS + MULT_POWER2) << EXP_bits_p);
+    int table_index = ps & (EXP_num_p - 1);
+    fputil::FPBits<double> bs;
+    bs.set_unbiased_exponent(ps >> EXP_bits_p);
+    ml = EXP_2_POW[table_index];
+    mult_e1 = bs.get_val();
+  }
+  double dx = fputil::multiply_add(ps_dbl, LN2_LOW,
+                                   fputil::multiply_add(ps_dbl, LN2_HIGH, x));
+
+  // Taylor series coefficients
+  double pe = dx * fputil::polyeval(dx, 1.0, 0x1.0p-1, 0x1.5555555555555p-3,
+                                    0x1.5555555555555p-5, 0x1.1111111111111p-7,
+                                    0x1.6c16c16c16c17p-10);
+
+  double r = fputil::multiply_add(ml, pe, pe) + ml;
+  return {mult_e1, r};
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   HDRS
     FModTest.h
   DEPENDS
     libc.include.math
     libc.src.math.fmod
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  expxf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    expxf_test.cpp
+  DEPENDS
+    libc.src.math.generic.common_constants
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/expxf_test.cpp

+//===-- Unittests for expxf -----------------------------------------------===//
+//
+// 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 "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/expxf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcExpxfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
+  auto fx = [](float x) -> float {
+    auto result = __llvm_libc::exp_eval<-1>(x);
+    return static_cast<float>(2 * result.mult_exp * result.r +
+                              2 * result.mult_exp);
+  };
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    if (isnan(x) || isinf(x) || x < -70 || x > 70 || fabsf(x) < 0x1.0p-10)
+      continue;
+
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, fx(x), 0.5);
+  }
+}