diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -536,7 +536,7 @@ HDRS ../exp2f.h DEPENDS - .common_constants + .explogxf libc.src.__support.FPUtil.fputil libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.nearest_integer @@ -1166,15 +1166,22 @@ -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 @@ -1190,9 +1197,8 @@ 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 @@ -1208,9 +1214,8 @@ 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 diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h --- a/libc/src/math/generic/common_constants.h +++ b/libc/src/math/generic/common_constants.h @@ -31,13 +31,6 @@ // > 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 diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp --- a/libc/src/math/generic/common_constants.cpp +++ b/libc/src/math/generic/common_constants.cpp @@ -225,19 +225,5 @@ 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 diff --git a/libc/src/math/generic/coshf.cpp b/libc/src/math/generic/coshf.cpp --- a/libc/src/math/generic/coshf.cpp +++ b/libc/src/math/generic/coshf.cpp @@ -9,7 +9,7 @@ #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 { diff --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp --- a/libc/src/math/generic/exp2f.cpp +++ b/libc/src/math/generic/exp2f.cpp @@ -8,6 +8,7 @@ #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" @@ -19,9 +20,6 @@ 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; @@ -78,24 +76,7 @@ } } - float kf = fputil::nearest_integer(x * mlp); - double dx = fputil::multiply_add(mmld, kf, x); - double mult_f, ml; - { - uint32_t ps = static_cast(kf) + (1 << (EXP_bits_p - 1)) + - (fputil::FPBits::EXPONENT_BIAS << EXP_bits_p); - fputil::FPBits 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 diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/explogxf.h @@ -0,0 +1,159 @@ +//===-- 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 + +#include + +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 +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(ps_dbl) + (1 << (EXP_bits_p - 1)) + + ((fputil::FPBits::EXPONENT_BIAS + MULT_POWER2) << EXP_bits_p); + int table_index = ps & (EXP_num_p - 1); + fputil::FPBits 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(kf) + (1 << (EXP_bits_p - 1)) + + (fputil::FPBits::EXPONENT_BIAS << EXP_bits_p); + fputil::FPBits 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; + 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 diff --git a/libc/src/math/generic/explogxf.cpp b/libc/src/math/generic/explogxf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/explogxf.cpp @@ -0,0 +1,89 @@ +//===-- 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 diff --git a/libc/src/math/generic/expxf.h b/libc/src/math/generic/expxf.h deleted file mode 100644 --- a/libc/src/math/generic/expxf.h +++ /dev/null @@ -1,81 +0,0 @@ -//===-- 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 - -#include - -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 -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(ps_dbl) + (1 << (EXP_bits_p - 1)) + - ((fputil::FPBits::EXPONENT_BIAS + MULT_POWER2) << EXP_bits_p); - int table_index = ps & (EXP_num_p - 1); - fputil::FPBits 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 diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp --- a/libc/src/math/generic/sinhf.cpp +++ b/libc/src/math/generic/sinhf.cpp @@ -8,7 +8,7 @@ #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 { diff --git a/libc/src/math/generic/tanhf.cpp b/libc/src/math/generic/tanhf.cpp --- a/libc/src/math/generic/tanhf.cpp +++ b/libc/src/math/generic/tanhf.cpp @@ -8,7 +8,7 @@ #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 { diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1319,14 +1319,16 @@ ) 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 ) diff --git a/libc/test/src/math/explogxf_test.cpp b/libc/test/src/math/explogxf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/explogxf_test.cpp @@ -0,0 +1,55 @@ +//===-- 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 + +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(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); +} diff --git a/libc/test/src/math/expxf_test.cpp b/libc/test/src/math/expxf_test.cpp deleted file mode 100644 --- a/libc/test/src/math/expxf_test.cpp +++ /dev/null @@ -1,38 +0,0 @@ -//===-- 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 - -#include -#include - -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(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); - } -} diff --git a/libc/test/src/math/in_float_range_test_helper.h b/libc/test/src/math/in_float_range_test_helper.h new file mode 100644 --- /dev/null +++ b/libc/test/src/math/in_float_range_test_helper.h @@ -0,0 +1,26 @@ +// +// Created by kirill on 8/30/22. +// + +#ifndef LLVM_IN_FLOAT_RANGE_TEST_HELPER_H +#define LLVM_IN_FLOAT_RANGE_TEST_HELPER_H + +#include + +#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((f)(x)), \ + (prec)); \ + } \ + } \ + } + +#endif // LLVM_IN_FLOAT_RANGE_TEST_HELPER_H