diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td --- a/libc/config/linux/api.td +++ b/libc/config/linux/api.td @@ -150,6 +150,8 @@ ]; let Functions = [ "cosf", + "expf", + "exp2f", "round", "sincosf", "sinf", diff --git a/libc/lib/CMakeLists.txt b/libc/lib/CMakeLists.txt --- a/libc/lib/CMakeLists.txt +++ b/libc/lib/CMakeLists.txt @@ -48,6 +48,8 @@ DEPENDS # math.h entrypoints libc.src.math.cosf + libc.src.math.expf + libc.src.math.exp2f libc.src.math.round libc.src.math.sincosf libc.src.math.sinf diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -185,6 +185,9 @@ FunctionSpec<"cosf", RetValSpec, [ArgSpec]>, FunctionSpec<"sinf", RetValSpec, [ArgSpec]>, + FunctionSpec<"expf", RetValSpec, [ArgSpec]>, + FunctionSpec<"exp2f", RetValSpec, [ArgSpec]>, + FunctionSpec<"round", RetValSpec, [ArgSpec]>, ] >; diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -1,5 +1,7 @@ -add_header_library( +add_object_library( math_utils + SRCS + math_utils.cpp HDRS math_utils.h DEPENDS @@ -68,3 +70,35 @@ libc.include.math libc.src.errno.__errno_location ) + +add_object_library( + exp_utils + HDRS + exp_utils.h + SRCS + exp_utils.cpp +) + +add_entrypoint_object( + expf + SRCS + expf.cpp + HDRS + expf.h + DEPENDS + .exp_utils + .math_utils + libc.include.math +) + +add_entrypoint_object( + exp2f + SRCS + exp2f.cpp + HDRS + exp2f.h + DEPENDS + .exp_utils + .math_utils + libc.include.math +) diff --git a/libc/src/math/cosf.cpp b/libc/src/math/cosf.cpp --- a/libc/src/math/cosf.cpp +++ b/libc/src/math/cosf.cpp @@ -58,7 +58,7 @@ return sinf_poly(x * s, x * x, p, n ^ 1); } - return invalidf(y); + return invalid(y); } } // namespace __llvm_libc diff --git a/libc/src/math/exp2f.h b/libc/src/math/exp2f.h new file mode 100644 --- /dev/null +++ b/libc/src/math/exp2f.h @@ -0,0 +1,18 @@ +//===-- Implementation header for exp2f -------------------------*- 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_MATH_EXP2F_H +#define LLVM_LIBC_SRC_MATH_EXP2F_H + +namespace __llvm_libc { + +float exp2f(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_EXP2F_H diff --git a/libc/src/math/exp2f.cpp b/libc/src/math/exp2f.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/exp2f.cpp @@ -0,0 +1,63 @@ +//===-- 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 "exp_utils.h" +#include "math_utils.h" + +#include "include/math.h" +#include "src/__support/common.h" + +#include + +#define T exp2f_data.tab +#define C exp2f_data.poly +#define SHIFT exp2f_data.shift_scaled + +namespace __llvm_libc { + +float LLVM_LIBC_ENTRYPOINT(exp2f)(float x) { + uint32_t abstop; + uint64_t ki, t; + // double_t for better performance on targets with FLT_EVAL_METHOD==2. + double_t kd, xd, z, r, r2, y, s; + + xd = static_cast(x); + abstop = top12_bits(x) & 0x7ff; + if (unlikely(abstop >= top12_bits(128.0f))) { + // |x| >= 128 or x is nan. + if (as_uint32_bits(x) == as_uint32_bits(-INFINITY)) + return 0.0f; + if (abstop >= top12_bits(INFINITY)) + return x + x; + if (x > 0.0f) + return overflow(0); + if (x <= -150.0f) + return underflow(0); + if (x < -149.0f) + return may_underflow(0); + } + + // x = k/N + r with r in [-1/(2N), 1/(2N)] and int k. + kd = static_cast(xd + SHIFT); + ki = as_uint64_bits(kd); + kd -= SHIFT; // k/N for int k. + r = xd - kd; + + // exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = as_double(t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return static_cast(y); +} + +} // namespace __llvm_libc diff --git a/libc/src/math/exp_utils.h b/libc/src/math/exp_utils.h new file mode 100644 --- /dev/null +++ b/libc/src/math/exp_utils.h @@ -0,0 +1,34 @@ +//===-- Collection of utils for exp and friends -----------------*- 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_MATH_EXP_UTILS_H +#define LLVM_LIBC_SRC_MATH_EXP_UTILS_H + +#include "math_utils.h" +#include + +#define EXP2F_TABLE_BITS 5 +#define EXP2F_POLY_ORDER 3 +#define N (1 << EXP2F_TABLE_BITS) + +namespace __llvm_libc { + +struct Exp2fDataTable { + uint64_t tab[1 << EXP2F_TABLE_BITS]; + double shift_scaled; + double poly[EXP2F_POLY_ORDER]; + double shift; + double invln2_scaled; + double poly_scaled[EXP2F_POLY_ORDER]; +}; + +extern const Exp2fDataTable exp2f_data; + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_EXP_UTILS_H diff --git a/libc/src/math/exp_utils.cpp b/libc/src/math/exp_utils.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/exp_utils.cpp @@ -0,0 +1,127 @@ +//===-- Implemention of exp and friends' utils ----------------------------===// +// +// 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 "exp_utils.h" + +namespace __llvm_libc { + +const Exp2fDataTable exp2f_data = { + // :tab[i] = uint(2^(i/N)) - (i << 52-BITS) + // used for computing 2^(k/N) for an int |k| < 150 N as + // double(tab[k%N] + (k << 52-BITS)) + { +// tab +#if N == 8 + 0x3ff0000000000000, + 0x3fef72b83c7d517b, + 0x3fef06fe0a31b715, + 0x3feebfdad5362a27, + 0x3feea09e667f3bcd, + 0x3feeace5422aa0db, + 0x3feee89f995ad3ad, + 0x3fef5818dcfba487, +#elif N == 16 + 0x3ff0000000000000, + 0x3fefb5586cf9890f, + 0x3fef72b83c7d517b, + 0x3fef387a6e756238, + 0x3fef06fe0a31b715, + 0x3feedea64c123422, + 0x3feebfdad5362a27, + 0x3feeab07dd485429, + 0x3feea09e667f3bcd, + 0x3feea11473eb0187, + 0x3feeace5422aa0db, + 0x3feec49182a3f090, + 0x3feee89f995ad3ad, + 0x3fef199bdd85529c, + 0x3fef5818dcfba487, + 0x3fefa4afa2a490da, +#elif N == 32 + 0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, + 0x3fef9301d0125b51, 0x3fef72b83c7d517b, 0x3fef54873168b9aa, + 0x3fef387a6e756238, 0x3fef1e9df51fdee1, 0x3fef06fe0a31b715, + 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d, + 0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429, + 0x3feea47eb03a5585, 0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, + 0x3feea11473eb0187, 0x3feea589994cce13, 0x3feeace5422aa0db, + 0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d, + 0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, + 0x3fef3720dcef9069, 0x3fef5818dcfba487, 0x3fef7c97337b9b5f, + 0x3fefa4afa2a490da, 0x3fefd0765b6e4540, +#elif N == 64 + 0x3ff0000000000000, 0x3fefec9a3e778061, 0x3fefd9b0d3158574, + 0x3fefc74518759bc8, 0x3fefb5586cf9890f, 0x3fefa3ec32d3d1a2, + 0x3fef9301d0125b51, 0x3fef829aaea92de0, 0x3fef72b83c7d517b, + 0x3fef635beb6fcb75, 0x3fef54873168b9aa, 0x3fef463b88628cd6, + 0x3fef387a6e756238, 0x3fef2b4565e27cdd, 0x3fef1e9df51fdee1, + 0x3fef1285a6e4030b, 0x3fef06fe0a31b715, 0x3feefc08b26416ff, + 0x3feef1a7373aa9cb, 0x3feee7db34e59ff7, 0x3feedea64c123422, + 0x3feed60a21f72e2a, 0x3feece086061892d, 0x3feec6a2b5c13cd0, + 0x3feebfdad5362a27, 0x3feeb9b2769d2ca7, 0x3feeb42b569d4f82, + 0x3feeaf4736b527da, 0x3feeab07dd485429, 0x3feea76f15ad2148, + 0x3feea47eb03a5585, 0x3feea23882552225, 0x3feea09e667f3bcd, + 0x3fee9fb23c651a2f, 0x3fee9f75e8ec5f74, 0x3fee9feb564267c9, + 0x3feea11473eb0187, 0x3feea2f336cf4e62, 0x3feea589994cce13, + 0x3feea8d99b4492ed, 0x3feeace5422aa0db, 0x3feeb1ae99157736, + 0x3feeb737b0cdc5e5, 0x3feebd829fde4e50, 0x3feec49182a3f090, + 0x3feecc667b5de565, 0x3feed503b23e255d, 0x3feede6b5579fdbf, + 0x3feee89f995ad3ad, 0x3feef3a2b84f15fb, 0x3feeff76f2fb5e47, + 0x3fef0c1e904bc1d2, 0x3fef199bdd85529c, 0x3fef27f12e57d14b, + 0x3fef3720dcef9069, 0x3fef472d4a07897c, 0x3fef5818dcfba487, + 0x3fef69e603db3285, 0x3fef7c97337b9b5f, 0x3fef902ee78b3ff6, + 0x3fefa4afa2a490da, 0x3fefba1bee615a27, 0x3fefd0765b6e4540, + 0x3fefe7c1819e90d8, +#endif + }, + as_double(0x4338000000000000) / N, // shift_scaled + { +// poly +#if N == 8 + as_double(0x3fac6a00335106e2), + as_double(0x3fcec0c313449f55), + as_double(0x3fe62e431111f69f), +#elif N == 16 + as_double(0x3fac6ac6aa313963), + as_double(0x3fcebfff4532d9ba), + as_double(0x3fe62e43001bc49f), +#elif N == 32 + as_double(0x3fac6af84b912394), + as_double(0x3fcebfce50fac4f3), + as_double(0x3fe62e42ff0c52d6), +#elif N == 64 + as_double(0x3fac6b04b4221b2a), + as_double(0x3fcebfc213e184d7), + as_double(0x3fe62e42fefb5b7f), +#endif + }, + as_double(0x4338000000000000), // shift + as_double(0x3ff71547652b82fe) * N, // invln2_scaled + { +// poly_scaled +#if N == 8 + as_double(0x3fac6a00335106e2) / N / N / N, + as_double(0x3fcec0c313449f55) / N / N, + as_double(0x3fe62e431111f69f) / N, +#elif N == 16 + as_double(0x3fac6ac6aa313963) / N / N / N, + as_double(0x3fcebfff4532d9ba) / N / N, + as_double(0x3fe62e43001bc49f) / N, +#elif N == 32 + as_double(0x3fac6af84b912394) / N / N / N, + as_double(0x3fcebfce50fac4f3) / N / N, + as_double(0x3fe62e42ff0c52d6) / N, +#elif N == 64 + as_double(0x3fac6b04b4221b2a) / N / N / N, + as_double(0x3fcebfc213e184d7) / N / N, + as_double(0x3fe62e42fefb5b7f) / N, +#endif + }, +}; + +} // namespace __llvm_libc diff --git a/libc/src/math/expf.h b/libc/src/math/expf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/expf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for expf --------------------------*- 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_MATH_EXPF_H +#define LLVM_LIBC_SRC_MATH_EXPF_H + +namespace __llvm_libc { + +float expf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_EXPF_H diff --git a/libc/src/math/expf.cpp b/libc/src/math/expf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/expf.cpp @@ -0,0 +1,69 @@ +//===-- Single-precision e^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 "exp_utils.h" +#include "math_utils.h" + +#include "include/math.h" +#include "src/__support/common.h" + +#include + +#define InvLn2N exp2f_data.invln2_scaled +#define T exp2f_data.tab +#define C exp2f_data.poly_scaled +#define SHIFT exp2f_data.shift + +namespace __llvm_libc { + +float LLVM_LIBC_ENTRYPOINT(expf)(float x) { + uint32_t abstop; + uint64_t ki, t; + // double_t for better performance on targets with FLT_EVAL_METHOD == 2. + double_t kd, xd, z, r, r2, y, s; + + xd = static_cast(x); + abstop = top12_bits(x) & 0x7ff; + if (unlikely(abstop >= top12_bits(88.0f))) { + // |x| >= 88 or x is nan. + if (as_uint32_bits(x) == as_uint32_bits(-INFINITY)) + return 0.0f; + if (abstop >= top12_bits(INFINITY)) + return x + x; + if (x > as_float(0x42b17217)) // x > log(0x1p128) ~= 88.72 + return overflow(0); + if (x < as_float(0xc2cff1b4)) // x < log(0x1p-150) ~= -103.97 + return underflow(0); + if (x < as_float(0xc2ce8ecf)) // x < log(0x1p-149) ~= -103.28 + return may_underflow(0); + } + + // x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. + z = InvLn2N * xd; + + // Round and convert z to int, the result is in [-150*N, 128*N] and + // ideally nearest int is used, otherwise the magnitude of r can be + // bigger which gives larger approximation error. + kd = static_cast(z + SHIFT); + ki = as_uint64_bits(kd); + kd -= SHIFT; + r = z - kd; + + // exp(x) = 2^(k/N) * 2^(r/N) ~= s *(C0*r^3 + C1*r^2 + C2*r + 1) + t = T[ki % N]; + t += ki << (52 - EXP2F_TABLE_BITS); + s = as_double(t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return static_cast(y); +} + +} // namespace __llvm_libc diff --git a/libc/src/math/math_utils.h b/libc/src/math/math_utils.h --- a/libc/src/math/math_utils.h +++ b/libc/src/math/math_utils.h @@ -11,24 +11,22 @@ #include "include/errno.h" #include "include/math.h" - #include "src/__support/common.h" #include "src/errno/llvmlibc_errno.h" +#include "utils/CPP/TypeTraits.h" #include namespace __llvm_libc { -static inline float with_errnof(float x, int err) { - if (math_errhandling & MATH_ERRNO) - llvmlibc_errno = err; - return x; -} - static inline uint32_t as_uint32_bits(float x) { return *reinterpret_cast(&x); } +static inline uint64_t as_uint64_bits(double x) { + return *reinterpret_cast(&x); +} + static inline float as_float(uint32_t x) { return *reinterpret_cast(&x); } @@ -37,12 +35,74 @@ return *reinterpret_cast(&x); } -static inline constexpr float invalidf(float x) { - float y = (x - x) / (x - x); - return isnan(x) ? y : with_errnof(y, EDOM); +static inline uint32_t top12_bits(float x) { return as_uint32_bits(x) >> 20; } + +static inline uint32_t top12_bits(double x) { return as_uint64_bits(x) >> 52; } + +// Values to trigger underflow and overflow. +template struct XFlowValues; + +template <> struct XFlowValues { + static const float overflow_value; + static const float underflow_value; + static const float may_underflow_value; +}; + +template <> struct XFlowValues { + static const double overflow_value; + static const double underflow_value; + static const double may_underflow_value; +}; + +template static inline T with_errno(T x, int err) { + if (math_errhandling & MATH_ERRNO) + llvmlibc_errno = err; + return x; } -static inline void force_eval_float(float x) { volatile float y UNUSED = x; } +template static inline void force_eval(T x) { + volatile T y UNUSED = x; +} + +template static inline T opt_barrier(T x) { + volatile T y = x; + return y; +} + +template struct IsFloatOrDouble { + static constexpr bool Value = + cpp::IsSame::Value || cpp::IsSame::Value; +}; + +template +using EnableIfFloatOrDouble = cpp::EnableIfType::Value, int>; + +template = 0> +T xflow(uint32_t sign, T y) { + // Underflow happens when two extremely small values are multiplied. + // Likewise, overflow happens when two large values are mulitplied. + y = opt_barrier(sign ? -y : y) * y; + return with_errno(y, ERANGE); +} + +template = 0> T overflow(uint32_t sign) { + return xflow(sign, XFlowValues::overflow_value); +} + +template = 0> T underflow(uint32_t sign) { + return xflow(sign, XFlowValues::underflow_value); +} + +template = 0> +T may_underflow(uint32_t sign) { + return xflow(sign, XFlowValues::may_underflow_value); +} + +template = 0> +static inline constexpr float invalid(T x) { + T y = (x - x) / (x - x); + return isnan(x) ? y : with_errno(y, EDOM); +} } // namespace __llvm_libc diff --git a/libc/src/math/math_utils.cpp b/libc/src/math/math_utils.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/math_utils.cpp @@ -0,0 +1,27 @@ +//===-- Implementation of math utils --------------------------------------===// +// +// 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 "math_utils.h" + +namespace __llvm_libc { + +const float XFlowValues::overflow_value = + as_float(0x70000000); // 0x1p97f +const float XFlowValues::underflow_value = + as_float(0x10000000); // 0x1p97f +const float XFlowValues::may_underflow_value = + as_float(0x1a200000); // 0x1.4p-75f + +const double XFlowValues::overflow_value = + as_double(0x7000000000000000); // 0x1p769 +const double XFlowValues::underflow_value = + as_double(0x1000000000000000); // 0x1p-767 +const double XFlowValues::may_underflow_value = + as_double(0x1e58000000000000); // 0x1.8p-538 + +} // namespace __llvm_libc diff --git a/libc/src/math/sincosf.cpp b/libc/src/math/sincosf.cpp --- a/libc/src/math/sincosf.cpp +++ b/libc/src/math/sincosf.cpp @@ -32,7 +32,7 @@ if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) { if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) // Force underflow for tiny y. - force_eval_float(x2); + force_eval(x2); *sinp = y; *cosp = 1.0f; return; @@ -69,7 +69,7 @@ // Needed to set errno for +-Inf, the add is a hack to work // around a gcc register allocation issue: just passing y // affects code generation in the fast path. - invalidf(y + y); + invalid(y + y); } } diff --git a/libc/src/math/sinf.cpp b/libc/src/math/sinf.cpp --- a/libc/src/math/sinf.cpp +++ b/libc/src/math/sinf.cpp @@ -32,7 +32,7 @@ if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) { if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) // Force underflow for tiny y. - force_eval_float(s); + force_eval(s); return y; } @@ -62,7 +62,7 @@ return sinf_poly(x * s, x * x, p, n); } - return invalidf(y); + return invalid(y); } } // 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 @@ -73,3 +73,30 @@ libc.src.math.sincosf libc.utils.CPP.standalone_cpp ) + +add_math_unittest( + expf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + expf_test.cpp + DEPENDS + .float_utils + libc.include.math + libc.src.math.expf +) + +add_math_unittest( + exp2f_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + exp2f_test.cpp + DEPENDS + .float_utils + libc.include.errno + libc.include.math + libc.src.math.exp2f +) diff --git a/libc/test/src/math/exp2f_test.cpp b/libc/test/src/math/exp2f_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exp2f_test.cpp @@ -0,0 +1,149 @@ +//===-- Unittests for exp2f -----------------------------------------------===// +// +// 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 "include/math.h" +#include "src/errno/llvmlibc_errno.h" +#include "src/math/exp2f.h" +#include "test/src/math/float.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/Test.h" + +#include + +using __llvm_libc::as_float; +using __llvm_libc::as_uint32_bits; + +using __llvm_libc::testing::FloatBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +// 12 additional bits of precision over the base precision of a |float| +// value. +static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12, + 0xFFF}; + +TEST(exp2fTest, SpecialNumbers) { + llvmlibc_errno = 0; + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::QNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::NegQNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::SNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::NegSNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::Inf)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::NegInf)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::Zero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::exp2f(as_float(FloatBits::NegZero)))); + EXPECT_EQ(llvmlibc_errno, 0); +} + +TEST(ExpfTest, Overflow) { + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::exp2f(as_float(0x7f7fffff)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::exp2f(as_float(0x43000000)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::exp2f(as_float(0x43000001)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); +} + +// Test with inputs which are the borders of underflow/overflow but still +// produce valid results without setting errno. +TEST(ExpfTest, Borderline) { + float x; + + llvmlibc_errno = 0; + x = as_float(0x42fa0001); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0x42ffffff); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc2fa0001); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc2fc0000); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc2fc0001); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc3150000); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); +} + +TEST(ExpfTest, Underflow) { + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::exp2f(as_float(0xff7fffff)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + float x = as_float(0xc3158000); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + x = as_float(0xc3165432); + EXPECT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + EXPECT_EQ(llvmlibc_errno, ERANGE); +} + +TEST(exp2fTest, InFloatRange) { + constexpr uint32_t count = 1000000; + constexpr uint32_t step = UINT32_MAX / count; + for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) { + float x = as_float(v); + if (isnan(x) || isinf(x)) + continue; + llvmlibc_errno = 0; + float result = __llvm_libc::exp2f(x); + + // If the computation resulted in an error or did not produce valid result + // in the single-precision floating point range, then ignore comparing with + // MPFR result as MPFR can still produce valid results because of its + // wider precision. + if (isnan(result) || isinf(result) || llvmlibc_errno != 0) + continue; + ASSERT_MPFR_MATCH(mpfr::OP_Exp2, x, __llvm_libc::exp2f(x), tolerance); + } +} diff --git a/libc/test/src/math/expf_test.cpp b/libc/test/src/math/expf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/expf_test.cpp @@ -0,0 +1,142 @@ +//===-- Unittests for expf ------------------------------------------------===// +// +// 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 "include/errno.h" +#include "include/math.h" +#include "src/errno/llvmlibc_errno.h" +#include "src/math/expf.h" +#include "test/src/math/float.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/Test.h" + +#include + +using __llvm_libc::as_float; +using __llvm_libc::as_uint32_bits; + +using __llvm_libc::testing::FloatBits; + +namespace mpfr = __llvm_libc::testing::mpfr; + +// 12 additional bits of precision over the base precision of a |float| +// value. +static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12, + 0xFFF}; + +TEST(ExpfTest, SpecialNumbers) { + llvmlibc_errno = 0; + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::QNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::NegQNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::SNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::NegSNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::Inf)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::NegInf)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::Zero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::expf(as_float(FloatBits::NegZero)))); + EXPECT_EQ(llvmlibc_errno, 0); +} + +TEST(ExpfTest, Overflow) { + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::expf(as_float(0x7f7fffff)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::expf(as_float(0x42cffff8)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Inf, + as_uint32_bits(__llvm_libc::expf(as_float(0x42d00008)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); +} + +TEST(ExpfTest, Underflow) { + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::expf(as_float(0xff7fffff)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::expf(as_float(0xc2cffff8)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); + + llvmlibc_errno = 0; + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::expf(as_float(0xc2d00008)))); + EXPECT_EQ(llvmlibc_errno, ERANGE); +} + +// Test with inputs which are the borders of underflow/overflow but still +// produce valid results without setting errno. +TEST(ExpfTest, Borderline) { + float x; + + llvmlibc_errno = 0; + x = as_float(0x42affff8); + ASSERT_MPFR_MATCH(mpfr::OP_Exp, x, __llvm_libc::expf(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0x42b00008); + ASSERT_MPFR_MATCH(mpfr::OP_Exp, x, __llvm_libc::expf(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc2affff8); + ASSERT_MPFR_MATCH(mpfr::OP_Exp, x, __llvm_libc::expf(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); + + x = as_float(0xc2b00008); + ASSERT_MPFR_MATCH(mpfr::OP_Exp, x, __llvm_libc::expf(x), tolerance); + EXPECT_EQ(llvmlibc_errno, 0); +} + +TEST(ExpfTest, InFloatRange) { + constexpr uint32_t count = 1000000; + constexpr uint32_t step = UINT32_MAX / count; + for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) { + float x = as_float(v); + if (isnan(x) || isinf(x)) + continue; + llvmlibc_errno = 0; + float result = __llvm_libc::expf(x); + + // If the computation resulted in an error or did not produce valid result + // in the single-precision floating point range, then ignore comparing with + // MPFR result as MPFR can still produce valid results because of its + // wider precision. + if (isnan(result) || isinf(result) || llvmlibc_errno != 0) + continue; + ASSERT_MPFR_MATCH(mpfr::OP_Exp, x, __llvm_libc::expf(x), tolerance); + } +} diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -39,10 +39,7 @@ uint32_t bits; }; -enum Operation { - OP_Cos, - OP_Sin, -}; +enum Operation { OP_Cos, OP_Sin, OP_Exp, OP_Exp2 }; namespace internal { diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -128,6 +128,12 @@ case OP_Sin: mpfr_sin(value, mpfrInput.value, MPFR_RNDN); break; + case OP_Exp: + mpfr_exp(value, mpfrInput.value, MPFR_RNDN); + break; + case OP_Exp2: + mpfr_exp2(value, mpfrInput.value, MPFR_RNDN); + break; } }