[libc][math] Added auxiliary function log2_eval for asinhf/acoshf/atanhf.

[orex]

Summary

1. double log2_eval(double) function added with better than float precision is added.
2. Some refactoring done to put all auxiliary functions and corresponding data

to one place to reuse the code.

3. Added tests for new functions.
4. Performance and precision tests of the function shows, that it more precise than exiting log2,

(no exceptional cases), but timing is ~5% higer that on current one.

[lntue]

libc/src/math/generic/expxf.h

0

explogxf.h is a better name for this file and cmake target.

15

Add comments about how the coefficients are generated.

42

Add comment about what the constant is and how it is generated.

49

Add comment about what the constant is and how it is generated.

libc/src/math/generic/supfuncf.cpp

28 ↗

(On Diff #456193)

Add comments about how these constants are defined and generated. The same for other tables below.

[orex]

Constants description added.

[orex]

Fix wrong upload.

[orex]

libc/src/math/generic/expxf.h

0

Please check also D132811 I've added to the file also atan_eval function, so therefore it is named such. If you have any other suggestion about the file name I'll appreciate your suggestion.

[lntue]

libc/src/math/generic/expxf.h

0

The constants and atan_eval function that you added in https://reviews.llvm.org/D132842 are independent with these exp and log utilities. Maybe you can put them into separate files / target called inv_trig_utils or inv_trig_helpers, as I expect asin and acos will also need some helpers, which could be put in there.

[orex]

Review fixes.

[orex]

The aux functions are split into two files according @lntue suggestion.

libc/src/math/generic/CMakeLists.txt

 
 add_entrypoint_object(
   exp2f
   SRCS
     exp2f.cpp
   HDRS
     ../exp2f.h
   DEPENDS
-    .common_constants
+    .explogxf
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
HDRS
     ../fmodf.h
   DEPENDS
     libc.include.math
     libc.src.__support.FPUtil.generic.fmod
   COMPILE_OPTIONS
     -O3
 )
 
+add_object_library(
+  explogxf
+  HDRS
+    explogxf.h
+  SRCS
+    explogxf.cpp
+)
+
 add_entrypoint_object(
   coshf
   SRCS
     coshf.cpp
   HDRS
     ../coshf.h
-    expxf.h
   DEPENDS
-    .common_constants
+    .explogxf
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   sinhf
   SRCS
     sinhf.cpp
   HDRS
     ../sinhf.h
-    expxf.h
   DEPENDS
-    .common_constants
+    .explogxf
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   tanhf
   SRCS
     tanhf.cpp
   HDRS
     ../tanhf.h
-    expxf.h 
   DEPENDS
-    .common_constants
+    .explogxf
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 

libc/src/math/generic/common_constants.h

 extern const double EXP_M1[195];
 
 // Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127.
 // Table is generated with Sollya as follow:
 // > display = hexadecimal;
 // > for i from 0 to 127 do { D(exp(i / 128)); };
 extern const double EXP_M2[128];
 
-static constexpr int EXP_bits_p = 5;
-static constexpr int EXP_num_p = 1 << EXP_bits_p;
-
-// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
-// printf("%.13a,\n", d[i]);
-extern const double EXP_2_POW[EXP_num_p];
-
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H

libc/src/math/generic/common_constants.cpp

const double EXP_M2[128] = {
     0x1.299be1f3e7f1cp1, 0x1.2bf16e7d2a38cp1, 0x1.2e4baacdb6614p1,
     0x1.30aaa04e80d05p1, 0x1.330e587b62b28p1, 0x1.3576dce33feadp1,
     0x1.37e437282d4eep1, 0x1.3a5670ff972edp1, 0x1.3ccd9432682b4p1,
     0x1.3f49aa9d30590p1, 0x1.41cabe304cb34p1, 0x1.4450d8f00edd4p1,
     0x1.46dc04f4e5338p1, 0x1.496c4c6b832dap1, 0x1.4c01b9950a111p1,
     0x1.4e9c56c731f5dp1, 0x1.513c2e6c731d7p1, 0x1.53e14b042f9cap1,
     0x1.568bb722dd593p1, 0x1.593b7d72305bbp1,
 };
-// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
-// printf("%.13a,\n", d[i]);
-const double EXP_2_POW[EXP_num_p] = {
-    -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2,
-    -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3,
-    -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3,
-    -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4,
-    -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
-    -0x1.5f134923757f3p-6, 0x0.0000000000000p+0,  0x1.66c34c5615d0fp-6,
-    0x1.6ab0d9f3121ecp-5,  0x1.1301d0125b50ap-4,  0x1.72b83c7d517aep-4,
-    0x1.d4873168b9aa8p-4,  0x1.1c3d373ab11c3p-3,  0x1.4f4efa8fef709p-3,
-    0x1.837f0518db8a9p-3,  0x1.b8d39b9d54e55p-3,  0x1.ef5326091a112p-3,
-    0x1.13821818624b4p-2,  0x1.2ff6b54d8a89cp-2,  0x1.4d0ad5a753e07p-2,
-    0x1.6ac1f752150a5p-2,  0x1.891fac0e95613p-2};
 
 } // namespace __llvm_libc

libc/src/math/generic/coshf.cpp

 //===-- Single-precision cosh 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 "src/math/coshf.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/multiply_add.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, coshf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
   xbits.set_sign(false);
   x = xbits.get_val();

libc/src/math/generic/exp2f.cpp

 //===-- Single-precision 2^x 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 "src/math/exp2f.h"
 #include "common_constants.h"
+#include "explogxf.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
-constexpr float mlp = EXP_num_p;
-constexpr float mmld = -1.0 / mlp;
-
 constexpr uint32_t exval1 = 0x3b42'9d37U;
 constexpr uint32_t exval2 = 0xbcf3'a937U;
 constexpr uint32_t exval_mask = exval1 & exval2;
 
 LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
if (unlikely(x_u == exval1)) { // x = 0x1.853a6ep-9f
       if (fputil::get_round() == FE_TONEAREST)
         return 0x1.00870ap+0f;
     } else if (unlikely(x_u == exval2)) { // x = -0x1.e7526ep-6f
       if (fputil::get_round() == FE_TONEAREST)
         return 0x1.f58d62p-1f;
     }
   }
 
-  float kf = fputil::nearest_integer(x * mlp);
-  double dx = fputil::multiply_add(mmld, kf, x);
-  double mult_f, ml;
-  {
-    uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
-                  (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
-    fputil::FPBits<double> bs;
-    bs.set_unbiased_exponent(ps >> EXP_bits_p);
-    ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
-    mult_f = bs.get_val();
-  }
-
-  // N[Table[Ln[2]^n/n!,{n,1,6}],30]
-  double pe = fputil::polyeval(
-      dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
-      0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
-
-  return mult_f * ml * pe;
+  return exp2_eval(x);
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/explogxf.h

+//===-- Single-precision general exp/log functions ------------------------===//
+//
+// 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_EXPLOGXF_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_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 {
+
+static constexpr int EXP_bits_p = 5;
+static constexpr int EXP_num_p = 1 << EXP_bits_p;
+constexpr double mlp = EXP_num_p;
+constexpr double mmld = -1.0 / mlp;
+
+// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
+// printf("%.13a,\n", d[i]);
+extern const double EXP_2_POW[EXP_num_p];
+
+constexpr int LOG_P1_BITS = 6;
+constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS;
+
+// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_LOG2[LOG_P1_SIZE];
+
+// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_1_OVER[LOG_P1_SIZE];
+
+// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
+// K_LOG2_ODD starts from x^3
+extern const double K_LOG2_ODD[4];
+extern const double K_LOG2_EVEN[4];
+
+// 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};
+}
+
+inline static double exp2_eval(double x) {
+  double kf = fputil::nearest_integer(x * mlp);
+  double dx = fputil::multiply_add(mmld, kf, x);
+  double mult_f, ml;
+  {
+    uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
+                  (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
+    fputil::FPBits<double> bs;
+    bs.set_unbiased_exponent(ps >> EXP_bits_p);
+    ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
+    mult_f = bs.get_val();
+  }
+
+  // Taylor series coefficients for 2^x
+  double pe = fputil::polyeval(
+      dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
+      0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
+
+  return mult_f * ml * pe;
+}
+
+// x should be positive, normal finite value
+inline static double log2_eval(double x) {
+  using FPB = fputil::FPBits<double>;
+  FPB bs(x);
+
+  double result = 0;
+  result += bs.get_exponent();
+
+  int p1 =
+      (bs.get_mantissa() >> (FPB::FloatProp::MANTISSA_WIDTH - LOG_P1_BITS)) &
+      (LOG_P1_SIZE - 1);
+
+  bs.bits &= FPB::FloatProp::MANTISSA_MASK >> LOG_P1_BITS;
+  bs.set_unbiased_exponent(FPB::FloatProp::EXPONENT_BIAS);
+  double dx = (bs.get_val() - 1.0) * LOG_P1_1_OVER[p1];
+
+  // Taylor series for log(2,1+x)
+  double c1 = fputil::multiply_add(dx, K_LOG2_ODD[0], K_LOG2_EVEN[0]);
+  double c2 = fputil::multiply_add(dx, K_LOG2_ODD[1], K_LOG2_EVEN[1]);
+  double c3 = fputil::multiply_add(dx, K_LOG2_ODD[2], K_LOG2_EVEN[2]);
+  double c4 = fputil::multiply_add(dx, K_LOG2_ODD[3], K_LOG2_EVEN[3]);
+
+  // c0 = dx * (1.0 / ln(2)) + LOG_P1_LOG2[p1]
+  double c0 = fputil::multiply_add(dx, 0x1.71547652b82fep+0, LOG_P1_LOG2[p1]);
+  result += __llvm_libc::fputil::polyeval(dx * dx, c0, c1, c2, c3, c4);
+  return result;
+}
+
+// x should be positive, normal finite value
+inline static double log_eval(double x) {
+  // ln(x) = log[2,x] * ln(2)
+  return log2_eval(x) * 0x1.62e42fefa39efp-1;
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H

libc/src/math/generic/explogxf.cpp

+//===-- Single-precision general exp/log functions ------------------------===//
+//
+// 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 "explogxf.h"
+
+namespace __llvm_libc {
+
+// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
+// printf("%.13a,\n", d[i]);
+alignas(64) const double EXP_2_POW[EXP_num_p] = {
+    -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2,
+    -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3,
+    -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3,
+    -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4,
+    -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
+    -0x1.5f134923757f3p-6, 0x0.0000000000000p+0,  0x1.66c34c5615d0fp-6,
+    0x1.6ab0d9f3121ecp-5,  0x1.1301d0125b50ap-4,  0x1.72b83c7d517aep-4,
+    0x1.d4873168b9aa8p-4,  0x1.1c3d373ab11c3p-3,  0x1.4f4efa8fef709p-3,
+    0x1.837f0518db8a9p-3,  0x1.b8d39b9d54e55p-3,  0x1.ef5326091a112p-3,
+    0x1.13821818624b4p-2,  0x1.2ff6b54d8a89cp-2,  0x1.4d0ad5a753e07p-2,
+    0x1.6ac1f752150a5p-2,  0x1.891fac0e95613p-2};
+
+// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
+alignas(64) const double LOG_P1_LOG2[LOG_P1_SIZE] = {
+    0x0.0000000000000p+0, 0x1.6e79685c2d22ap-6, 0x1.6bad3758efd87p-5,
+    0x1.0eb389fa29f9bp-4, 0x1.663f6fac91316p-4, 0x1.bc84240adabbap-4,
+    0x1.08c588cda79e4p-3, 0x1.32ae9e278ae1ap-3, 0x1.5c01a39fbd688p-3,
+    0x1.84c2bd02f03b3p-3, 0x1.acf5e2db4ec94p-3, 0x1.d49ee4c325970p-3,
+    0x1.fbc16b902680ap-3, 0x1.11307dad30b76p-2, 0x1.24407ab0e073ap-2,
+    0x1.37124cea4cdedp-2, 0x1.49a784bcd1b8bp-2, 0x1.5c01a39fbd688p-2,
+    0x1.6e221cd9d0cdep-2, 0x1.800a563161c54p-2, 0x1.91bba891f1709p-2,
+    0x1.a33760a7f6051p-2, 0x1.b47ebf73882a1p-2, 0x1.c592fad295b56p-2,
+    0x1.d6753e032ea0fp-2, 0x1.e726aa1e754d2p-2, 0x1.f7a8568cb06cfp-2,
+    0x1.03fda8b97997fp-1, 0x1.0c10500d63aa6p-1, 0x1.140c9faa1e544p-1,
+    0x1.1bf311e95d00ep-1, 0x1.23c41d42727c8p-1, 0x1.2b803473f7ad1p-1,
+    0x1.3327c6ab49ca7p-1, 0x1.3abb3faa02167p-1, 0x1.423b07e986aa9p-1,
+    0x1.49a784bcd1b8bp-1, 0x1.510118708a8f9p-1, 0x1.5848226989d34p-1,
+    0x1.5f7cff41e09afp-1, 0x1.66a008e4788ccp-1, 0x1.6db196a76194ap-1,
+    0x1.74b1fd64e0754p-1, 0x1.7ba18f93502e4p-1, 0x1.82809d5be7073p-1,
+    0x1.894f74b06ef8bp-1, 0x1.900e6160002cdp-1, 0x1.96bdad2acb5f6p-1,
+    0x1.9d5d9fd5010b3p-1, 0x1.a3ee7f38e181fp-1, 0x1.aa708f58014d3p-1,
+    0x1.b0e4126bcc86cp-1, 0x1.b74948f5532dap-1, 0x1.bda071cc67e6ep-1,
+    0x1.c3e9ca2e1a055p-1, 0x1.ca258dca93316p-1, 0x1.d053f6d260896p-1,
+    0x1.d6753e032ea0fp-1, 0x1.dc899ab3ff56cp-1, 0x1.e29142e0e0140p-1,
+    0x1.e88c6b3626a73p-1, 0x1.ee7b471b3a950p-1, 0x1.f45e08bcf0655p-1,
+    0x1.fa34e1177c233p-1,
+};
+
+// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
+alignas(64) const double LOG_P1_1_OVER[LOG_P1_SIZE] = {
+    0x1.0000000000000p+0, 0x1.f81f81f81f820p-1, 0x1.f07c1f07c1f08p-1,
+    0x1.e9131abf0b767p-1, 0x1.e1e1e1e1e1e1ep-1, 0x1.dae6076b981dbp-1,
+    0x1.d41d41d41d41dp-1, 0x1.cd85689039b0bp-1, 0x1.c71c71c71c71cp-1,
+    0x1.c0e070381c0e0p-1, 0x1.bacf914c1bad0p-1, 0x1.b4e81b4e81b4fp-1,
+    0x1.af286bca1af28p-1, 0x1.a98ef606a63bep-1, 0x1.a41a41a41a41ap-1,
+    0x1.9ec8e951033d9p-1, 0x1.999999999999ap-1, 0x1.948b0fcd6e9e0p-1,
+    0x1.8f9c18f9c18fap-1, 0x1.8acb90f6bf3aap-1, 0x1.8618618618618p-1,
+    0x1.8181818181818p-1, 0x1.7d05f417d05f4p-1, 0x1.78a4c8178a4c8p-1,
+    0x1.745d1745d1746p-1, 0x1.702e05c0b8170p-1, 0x1.6c16c16c16c17p-1,
+    0x1.6816816816817p-1, 0x1.642c8590b2164p-1, 0x1.6058160581606p-1,
+    0x1.5c9882b931057p-1, 0x1.58ed2308158edp-1, 0x1.5555555555555p-1,
+    0x1.51d07eae2f815p-1, 0x1.4e5e0a72f0539p-1, 0x1.4afd6a052bf5bp-1,
+    0x1.47ae147ae147bp-1, 0x1.446f86562d9fbp-1, 0x1.4141414141414p-1,
+    0x1.3e22cbce4a902p-1, 0x1.3b13b13b13b14p-1, 0x1.3813813813814p-1,
+    0x1.3521cfb2b78c1p-1, 0x1.323e34a2b10bfp-1, 0x1.2f684bda12f68p-1,
+    0x1.2c9fb4d812ca0p-1, 0x1.29e4129e4129ep-1, 0x1.27350b8812735p-1,
+    0x1.2492492492492p-1, 0x1.21fb78121fb78p-1, 0x1.1f7047dc11f70p-1,
+    0x1.1cf06ada2811dp-1, 0x1.1a7b9611a7b96p-1, 0x1.1811811811812p-1,
+    0x1.15b1e5f75270dp-1, 0x1.135c81135c811p-1, 0x1.1111111111111p-1,
+    0x1.0ecf56be69c90p-1, 0x1.0c9714fbcda3bp-1, 0x1.0a6810a6810a7p-1,
+    0x1.0842108421084p-1, 0x1.0624dd2f1a9fcp-1, 0x1.0410410410410p-1,
+    0x1.0204081020408p-1};
+
+// Taylos series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
+// K_LOG2_ODD starts from x^3
+alignas(64) const
+    double K_LOG2_ODD[4] = {0x1.ec709dc3a03fdp-2, 0x1.2776c50ef9bfep-2,
+                            0x1.a61762a7aded9p-3, 0x1.484b13d7c02a9p-3};
+
+alignas(64) const
+    double K_LOG2_EVEN[4] = {-0x1.71547652b82fep-1, -0x1.71547652b82fep-2,
+                             -0x1.ec709dc3a03fdp-3, -0x1.2776c50ef9bfep-3};
+
+} // namespace __llvm_libc

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/src/math/generic/sinhf.cpp

 //===-- Single-precision sinh 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 "src/math/sinhf.h"
 #include "src/__support/FPUtil/FPBits.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, sinhf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
   bool sign = xbits.get_sign();
   uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK;

libc/src/math/generic/tanhf.cpp

 //===-- Single-precision tanh 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 "src/math/tanhf.h"
 #include "src/__support/FPUtil/FPBits.h"
-#include "src/math/generic/expxf.h"
+#include "src/math/generic/explogxf.h"
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, tanhf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
   bool sign = xbits.get_sign();
   uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK;

libc/test/src/CMakeLists.txt

function(add_fp_unittest name)
   )
 endfunction(add_fp_unittest)
 
 add_subdirectory(__support)
 add_subdirectory(ctype)
 add_subdirectory(errno)
 add_subdirectory(fenv)
 add_subdirectory(inttypes)
-# add_subdirectory(math)
+add_subdirectory(math)
 add_subdirectory(string)
 add_subdirectory(stdlib)
 add_subdirectory(stdio)
 
 if(${LIBC_TARGET_OS} STREQUAL "linux")
   add_subdirectory(fcntl)
   add_subdirectory(sys)
   add_subdirectory(unistd)

libc/test/src/math/CMakeLists.txt

HDRS
     FModTest.h
   DEPENDS
     libc.include.math
     libc.src.math.fmod
     libc.src.__support.FPUtil.fputil
 )
 
 add_fp_unittest(
-  expxf_test
+  explogxf_test
   NEED_MPFR
   SUITE
     libc_math_unittests
+  HDRS
+    in_float_range_test_helper.h
   SRCS
-    expxf_test.cpp
+    explogxf_test.cpp
   DEPENDS
-    libc.src.math.generic.common_constants
+    libc.src.math.generic.explogxf
     libc.src.__support.FPUtil.fputil
 )
 
 add_fp_unittest(
   coshf_test
   NEED_MPFR
   SUITE
     libc_math_unittests

libc/test/src/math/explogxf_test.cpp

+//===-- Unittests for explogxf --------------------------------------------===//
+//
+// 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 "in_float_range_test_helper.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/generic/explogxf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+constexpr int def_count = 100003;
+constexpr float def_prec = 0.500001f;
+
+TEST(LlvmLibcExpxfTest, InFloatRange) {
+  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);
+  };
+  auto f_check = [](float x) -> bool {
+    return !(
+        (isnan(x) || isinf(x) || x < -70 || x > 70 || fabsf(x) < 0x1.0p-10));
+  };
+  CHECK_DATA(0.0f, neg_inf, mpfr::Operation::Exp, fx, f_check, def_count,
+             def_prec);
+}
+
+TEST(LlvmLibcExp2xfTest, InFloatRange) {
+  auto f_check = [](float x) -> bool {
+    return !(
+        (isnan(x) || isinf(x) || x < -130 || x > 130 || fabsf(x) < 0x1.0p-10));
+  };
+  CHECK_DATA(0.0f, neg_inf, mpfr::Operation::Exp2, __llvm_libc::exp2_eval,
+             f_check, def_count, def_prec);
+}
+
+TEST(LlvmLibcLog2xfTest, InFloatRange) {
+  CHECK_DATA(0.0f, inf, mpfr::Operation::Log2, __llvm_libc::log2_eval, isnormal,
+             def_count, def_prec);
+}
+
+TEST(LlvmLibcLogxfTest, InFloatRange) {
+  CHECK_DATA(0.0f, inf, mpfr::Operation::Log, __llvm_libc::log_eval, isnormal,
+             def_count, def_prec);
+}

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);
-  }
-}

libc/test/src/math/in_float_range_test_helper.h

+//
+// Created by kirill on 8/30/22.
+//
+
+#ifndef LLVM_IN_FLOAT_RANGE_TEST_HELPER_H
+#define LLVM_IN_FLOAT_RANGE_TEST_HELPER_H
+
+#include <stdint.h>
+
+#define CHECK_DATA(start, stop, mfp_op, f, f_check, count, prec)               \
+  {                                                                            \
+    uint64_t ustart = FPBits(start).uintval();                                 \
+    uint64_t ustop = FPBits(stop).uintval();                                   \
+    for (uint64_t i = 0;; ++i) {                                               \
+      uint64_t v = ustart + (ustop - ustart) * i / count;                      \
+      if (v > ustop)                                                           \
+        break;                                                                 \
+      float x = FPBits(uint32_t(v)).get_val();                                 \
+      if ((f_check)(x)) {                                                      \
+        EXPECT_MPFR_MATCH_ALL_ROUNDING(mfp_op, x, static_cast<float>((f)(x)),  \
+                                       (prec));                                \
+      }                                                                        \
+    }                                                                          \
+  }
+
+#endif // LLVM_IN_FLOAT_RANGE_TEST_HELPER_H