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 @@ -77,12 +77,13 @@ HDRS ../sinf.h DEPENDS - .sincosf_utils + .common_constants libc.include.math libc.src.errno.errno libc.src.__support.FPUtil.fputil COMPILE_OPTIONS -O3 + -mfma ) add_entrypoint_object( 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,6 +31,18 @@ // > for i from 0 to 127 do { D(exp(i / 128)); }; extern const double EXP_M2[128]; +// Lookup table for sin(k * pi / 256) with k = 0, ..., 128. +// Table is generated with Sollya as follow: +// > display = hexadecimal; +// > for i from 0 to 128 do { D(sin(i * pi / 256)); }; +extern const double SIN_K_PI_OVER_256[129]; + +// Digits of 256/pi, generated by Sollya with: +// > a0 = D(256/pi); +// > a1 = D(256/pi - a0); +// > a2 = D(256/pi - a0 - a1); +// > a3 = D(256/pi - a0 - a1 - a2); +extern const double TWOFIFTYSIX_OVER_PI[4]; } // 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 @@ -226,4 +226,63 @@ 0x1.568bb722dd593p1, 0x1.593b7d72305bbp1, }; +// Lookup table for sin(k * pi / 256) with k = 0, ..., 128. +// Table is generated with Sollya as follow: +// > display = hexadecimal; +// > for i from 0 to 128 do { D(sin(i * pi / 256)); }; +const double SIN_K_PI_OVER_256[129] = { + 0x0.0000000000000p+0, 0x1.921d1fcdec784p-7, 0x1.92155f7a3667ep-6, + 0x1.2d865759455cdp-5, 0x1.91f65f10dd814p-5, 0x1.f656e79f820e0p-5, + 0x1.2d52092ce19f6p-4, 0x1.5f6d00a9aa419p-4, 0x1.917a6bc29b42cp-4, + 0x1.c3785c79ec2d5p-4, 0x1.f564e56a9730ep-4, 0x1.139f0cedaf577p-3, + 0x1.2c8106e8e613ap-3, 0x1.45576b1293e5ap-3, 0x1.5e214448b3fc6p-3, + 0x1.76dd9de50bf31p-3, 0x1.8f8b83c69a60bp-3, 0x1.a82a025b00451p-3, + 0x1.c0b826a7e4f63p-3, 0x1.d934fe5454311p-3, 0x1.f19f97b215f1bp-3, + 0x1.04fb80e37fdaep-2, 0x1.111d262b1f677p-2, 0x1.1d3443f4cdb3ep-2, + 0x1.294062ed59f06p-2, 0x1.35410c2e18152p-2, 0x1.4135c94176601p-2, + 0x1.4d1e24278e76ap-2, 0x1.58f9a75ab1fddp-2, 0x1.64c7ddd3f27c6p-2, + 0x1.7088530fa459fp-2, 0x1.7c3a9311dcce7p-2, 0x1.87de2a6aea963p-2, + 0x1.9372a63bc93d7p-2, 0x1.9ef7943a8ed8ap-2, 0x1.aa6c82b6d3fcap-2, + 0x1.b5d1009e15cc0p-2, 0x1.c1249d8011ee7p-2, 0x1.cc66e9931c45ep-2, + 0x1.d79775b86e389p-2, 0x1.e2b5d3806f63bp-2, 0x1.edc1952ef78d6p-2, + 0x1.f8ba4dbf89abap-2, 0x1.01cfc874c3eb7p-1, 0x1.073879922ffeep-1, + 0x1.0c9704d5d898fp-1, 0x1.11eb3541b4b23p-1, 0x1.1734d63dedb49p-1, + 0x1.1c73b39ae68c8p-1, 0x1.21a799933eb59p-1, 0x1.26d054cdd12dfp-1, + 0x1.2bedb25faf3eap-1, 0x1.30ff7fce17035p-1, 0x1.36058b10659f3p-1, + 0x1.3affa292050b9p-1, 0x1.3fed9534556d4p-1, 0x1.44cf325091dd6p-1, + 0x1.49a449b9b0939p-1, 0x1.4e6cabbe3e5e9p-1, 0x1.5328292a35596p-1, + 0x1.57d69348ceca0p-1, 0x1.5c77bbe65018cp-1, 0x1.610b7551d2cdfp-1, + 0x1.6591925f0783dp-1, 0x1.6a09e667f3bcdp-1, 0x1.6e74454eaa8afp-1, + 0x1.72d0837efff96p-1, 0x1.771e75f037261p-1, 0x1.7b5df226aafafp-1, + 0x1.7f8ece3571771p-1, 0x1.83b0e0bff976ep-1, 0x1.87c400fba2ebfp-1, + 0x1.8bc806b151741p-1, 0x1.8fbcca3ef940dp-1, 0x1.93a22499263fbp-1, + 0x1.9777ef4c7d742p-1, 0x1.9b3e047f38741p-1, 0x1.9ef43ef29af94p-1, + 0x1.a29a7a0462782p-1, 0x1.a63091b02fae2p-1, 0x1.a9b66290ea1a3p-1, + 0x1.ad2bc9e21d511p-1, 0x1.b090a58150200p-1, 0x1.b3e4d3ef55712p-1, + 0x1.b728345196e3ep-1, 0x1.ba5aa673590d2p-1, 0x1.bd7c0ac6f952ap-1, + 0x1.c08c426725549p-1, 0x1.c38b2f180bdb1p-1, 0x1.c678b3488739bp-1, + 0x1.c954b213411f5p-1, 0x1.cc1f0f3fcfc5cp-1, 0x1.ced7af43cc773p-1, + 0x1.d17e7743e35dcp-1, 0x1.d4134d14dc93ap-1, 0x1.d696173c9e68bp-1, + 0x1.d906bcf328d46p-1, 0x1.db6526238a09bp-1, 0x1.ddb13b6ccc23cp-1, + 0x1.dfeae622dbe2bp-1, 0x1.e212104f686e5p-1, 0x1.e426a4b2bc17ep-1, + 0x1.e6288ec48e112p-1, 0x1.e817bab4cd10dp-1, 0x1.e9f4156c62ddap-1, + 0x1.ebbd8c8df0b74p-1, 0x1.ed740e7684963p-1, 0x1.ef178a3e473c2p-1, + 0x1.f0a7efb9230d7p-1, 0x1.f2252f7763adap-1, 0x1.f38f3ac64e589p-1, + 0x1.f4e603b0b2f2dp-1, 0x1.f6297cff75cb0p-1, 0x1.f7599a3a12077p-1, + 0x1.f8764fa714ba9p-1, 0x1.f97f924c9099bp-1, 0x1.fa7557f08a517p-1, + 0x1.fb5797195d741p-1, 0x1.fc26470e19fd3p-1, 0x1.fce15fd6da67bp-1, + 0x1.fd88da3d12526p-1, 0x1.fe1cafcbd5b09p-1, 0x1.fe9cdad01883ap-1, + 0x1.ff095658e71adp-1, 0x1.ff621e3796d7ep-1, 0x1.ffa72effef75dp-1, + 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.0000000000000p+0, +}; + +// Digits of 256/pi, generated by Sollya with: +// > a0 = D(256/pi); +// > a1 = D(256/pi - a0); +// > a2 = D(256/pi - a0 - a1); +// > a3 = D(256/pi - a0 - a1 - a2); +const double TWOFIFTYSIX_OVER_PI[4] = { + 0x1.45f306dc9c883p6, -0x1.6b01ec5417056p-48, -0x1.6447e493ad4cep-102, + 0x1.e21c820ff28b2p-156}; + } // namespace __llvm_libc diff --git a/libc/src/math/generic/sinf.cpp b/libc/src/math/generic/sinf.cpp --- a/libc/src/math/generic/sinf.cpp +++ b/libc/src/math/generic/sinf.cpp @@ -7,63 +7,250 @@ //===----------------------------------------------------------------------===// #include "src/math/sinf.h" -#include "math_utils.h" -#include "sincosf_utils.h" - +#include "src/__support/FPUtil/BasicOperations.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FMA.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" #include "src/__support/common.h" -#include +#include "src/math/generic/common_constants.h" -#include +#include namespace __llvm_libc { -// Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative -// error is 0.5303 * 2^-23. A single-step range reduction is used for -// small values. Large inputs have their range reduced using fast integer -// arithmetic. -LLVM_LIBC_FUNCTION(float, sinf, (float y)) { - double x = y; - double s; - int n; - const sincos_t *p = &SINCOSF_TABLE[0]; - - if (abstop12(y) < abstop12(PIO4)) { - s = x * x; - - if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) { - if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) - // Force underflow for tiny y. - force_eval(s); - return y; - } +namespace { +static inline void get_sin_cos(int64_t k, double &sin_k, double &cos_k) { + int idx = k & 0x007f; + int quadrant = (k & 0x0180); + switch (quadrant >> 7) { + case 0: + sin_k = SIN_K_PI_OVER_256[idx]; + cos_k = SIN_K_PI_OVER_256[128 - idx]; + break; + case 1: + sin_k = SIN_K_PI_OVER_256[128 - idx]; + cos_k = -SIN_K_PI_OVER_256[idx]; + break; + case 2: + sin_k = -SIN_K_PI_OVER_256[idx]; + cos_k = -SIN_K_PI_OVER_256[128 - idx]; + break; + case 3: + sin_k = -SIN_K_PI_OVER_256[128 - idx]; + cos_k = SIN_K_PI_OVER_256[idx]; + break; + default: + __builtin_unreachable(); + } +} - return sinf_poly(x, s, p, 0); - } else if (likely(abstop12(y) < abstop12(120.0f))) { - x = reduce_fast(x, p, &n); +} // namespace - // Setup the signs for sin and cos. - s = p->sign[n & 3]; +INLINE_FMA +LLVM_LIBC_FUNCTION(float, sinf, (float x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); - if (n & 2) - p = &SINCOSF_TABLE[1]; + uint32_t x_u = xbits.uintval(); + uint32_t x_abs = x_u & 0x7fff'ffffU; - return sinf_poly(x * s, x * x, p, n); - } else if (abstop12(y) < abstop12(INFINITY)) { - uint32_t xi = as_uint32_bits(y); - int sign = xi >> 31; + double xd = static_cast(x); - x = reduce_large(xi, &n); + // Range reduction: + // For |x| > pi/16, we perform range reduction as follows: + // Find k and y such that: + // x = (k + y) * pi/256 + // k is an integer + // |y| < 1 + // This is done by performing the following computation: + // k = round(x * 256/pi) + // y = round(x * 256/pi - k) + // The digits of 256/pi are stored using 4 doubles. The last double stores + // digits ranging from 2^(-208) to 2^(-156) of 256/pi, so when multiplying + // by the largest values of single precision, the resulting output should be + // correct up to 2^(-208 + 128) ~ 2^-80. By the worst-case analysis of range + // reduction, |y| >= 2^-38, so this should give us more than 40 bits of + // accuracy. For the worst-case estimation of range reduction, see for + // instances: + // Elementary Functions by J-M. Muller, Chapter 11, + // Handbook of Floating-Point Arithmetic by J-M. Muller et. al., + // Chapter 10.2. + // + // Once k and y are computed, we then deduce the answer by the sine of sum + // formula: + // sin(x) = sin((k + y)*pi/256) + // = sin(y*pi/256) * cos(k*pi/256) + cos(y*pi/256) * sin(k*pi/256) + // The values of sin(k*pi/256) and cos(k*pi/256) are precomputed and stored + // using a vector of 129 doubles. Sin(y*pi/256) and cos(y*pi/256) are computed + // using degree-5 and degree-6 minimax polynomials generated by Sollya + // respectively, + double y, sin_k, cos_k; + int64_t k; - // Setup signs for sin and cos - include original sign. - s = p->sign[(n + sign) & 3]; + // |x| < 2^46 + if (x_abs < 0x5680'0000U) { + // |x| < 0x1.d12ed2p-12f + if (x_abs < 0x39e8'9769U) { + if (unlikely(x_abs == 0U)) { + // For signed zeros. + return x; + } + // When |x| < 2^-12, the relative error of the approximation sin(x) ~ x + // is: + // |sin(x) - x| / |sin(x)| < |x^3| / (6|x|) + // = x^2 / 6 + // < 2^-25 + // < epsilon(1)/2. + // So the correctly rounded values of sin(x) are: + // = x - sign(x)*eps(x) if rounding mode = FE_TOWARDZERO, + // or (rounding mode = FE_UPWARD and x is + // negative), + // = x otherwise. + // To simplify the rounding decision and make it more efficient, we use + // fma(x, -2^-25, x) instead. + // An exhaustive test shows that this formula work correctly for all + // rounding modes up to |x| < 0x1.c555dep-11f. + return fputil::fma(x, -0x1.0p-25f, x); + } - if ((n + sign) & 2) - p = &SINCOSF_TABLE[1]; + // |x| <= pi/16 + if (x_abs <= 0x3e49'0fdbU) { + double xsq = xd * xd; - return sinf_poly(x * s, x * x, p, n); + // Degree-9 polynomial approximation: + // sin(x) ~ x + a_3 x^3 + a_5 x^5 + a_7 x^7 + a_9 x^9 + // = x (1 + a_3 x^2 + ... + a_9 x^8) + // = x * P(x^2) + // generated by Sollya with the following commands: + // > display = hexadecimal; + // > Q = fpminimax(sin(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/16]); + double result = fputil::polyeval( + xsq, 1.0, -0x1.55555555554c6p-3, 0x1.1111111085e65p-7, + -0x1.a019f70fb4d4fp-13, 0x1.718d179815e74p-19); + return xd * result; + } + + // Exceptional cases. + switch (x_abs) { + case 0x4619'9998U: // |x| = 0x1.33333p+13f + if (xbits.get_sign()) { + if (fputil::get_round() == FE_UPWARD) + return 0x1.63f4bcp-2f; + return 0x1.63f4bap-2f; + } + if (fputil::get_round() == FE_DOWNWARD) + return -0x1.63f4bcp-2f; + return -0x1.63f4bap-2f; + case 0x4afd'ece4U: { // |x| = 0x1.fbd9c8p+22f; + int rounding = fputil::get_round(); + if (xbits.get_sign()) { + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return 0x1.ff6dcp-1f; + return 0x1.ff6dc2p-1f; + } + if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO) + return -0x1.ff6dcp-1f; + return -0x1.ff6dc2p-1f; + } + case 0x5239'47f6U: { // |x| = 0x1.728fecp+37f; + int rounding = fputil::get_round(); + if (xbits.get_sign()) { + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return 0x1.24f23ap-1f; + return 0x1.24f23cp-1f; + } + if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO) + return -0x1.24f23ap-1f; + return -0x1.24f23cp-1f; + } + } + + // Since casting from float to integer is round-toward-zero, we add + // sign(x)*0.5 before casting to make it round-to-nearest. + k = static_cast( + fputil::fma(xd, TWOFIFTYSIX_OVER_PI[0], xbits.get_sign() ? -0.5 : 0.5)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[0], -static_cast(k)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[1], y); + } else { + // Exceptional values + if (unlikely(x_abs == 0x6a19'76f1U || x_abs == 0x6f79'be45U)) { + // |x| = 0x1.32ede2p+85f or |x| = 0x1.f37c8ap+95f + int rounding = fputil::get_round(); + if (xbits.get_sign()) { + if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO) + return -0x1.fffffep-1f; + return -1.0f; + } + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return 0x1.fffffep-1f; + return 1.0f; + } + + // x is inf or nan. + if (unlikely(x_abs >= 0x7f80'0000U)) { + if (x_abs == 0x7f80'0000U) + errno = EDOM; + return x + + FPBits::build_nan(1 << (fputil::MantissaWidth::VALUE - 1)); + } + + // 2^46 <= |x| < 2^98 + if (x_abs < 0x7080'0000U) { + // - When x >= 2^46, double(xd*TWOFIFTYSIX_OVER_PI[0]) is an integer, so + // we just need to keep its lowest 9 bits, so that when adding with the + // lower part, it does not result in overflow when converting to integer. + // - When x >= 2^55, the lowest 9 bit of the product + // double(xd*TWOFIFTYSIX_OVER_PI[0]) can also be dropped. + fputil::FPBits prod(xd * TWOFIFTYSIX_OVER_PI[0]); + prod.bits &= (x_abs < 0x5b00'0000U) ? (~0x1ffULL) : (~0ULL); // |x| < 2^55 + double truncated_prod = + fputil::fma(xd, TWOFIFTYSIX_OVER_PI[0], -static_cast(prod)); + k = static_cast( + fputil::fma(xd, TWOFIFTYSIX_OVER_PI[1], truncated_prod) + + (xbits.get_sign() ? -0.5 : 0.5)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[1], + truncated_prod - static_cast(k)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[2], y); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[3], y); + } else { + // - When x >= 2^98, double(xd*TWOFIFTYSIX_OVER_PI[1])*4 is an integer, so + // we just need to keep its lowest 11 bits, so that when adding with the + // lower part, it does not result in overflow when converting to integer. + // - When x >= 2^109, the lowest 11 bit of the product + // double(xd*TWOFIFTYSIX_OVER_PI[1]) can also be dropped. + fputil::FPBits prod(xd * TWOFIFTYSIX_OVER_PI[1]); + prod.bits &= + (x_abs < 0x7600'0000U) ? (~0x0fffULL) : (~0ULL); // |x| < 2^109 + double truncated_prod = + fputil::fma(xd, TWOFIFTYSIX_OVER_PI[1], -static_cast(prod)); + k = static_cast( + fputil::fma(xd, TWOFIFTYSIX_OVER_PI[2], truncated_prod) + + (xbits.get_sign() ? -0.5 : 0.5)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[2], + truncated_prod - static_cast(k)); + y = fputil::fma(xd, TWOFIFTYSIX_OVER_PI[3], y); + } } - return invalid(y); + double ysq = y * y; + get_sin_cos(k, sin_k, cos_k); + + // Degree-5 minimax polynomial for sin(y * pi / 256) generated by Sollya with: + // > Q = fpminimax(sin(x*pi/256)/x, [|0, 2, 4|], [|D...|], [0, 1]); + double sin_y = + y * fputil::polyeval(ysq, 0x1.921fb54442d18p-7, -0x1.4abbce62424a9p-22, + 0x1.466b4cf57923bp-39); + // Degree-6 minimax polynomial for cos(y*pi/256) generated by Sollya with: + // > P = fpminimax(cos(x*pi/256), [|0, 2, 4, 6|], [|1, D...|], [0, 1]); + // Note that cosm1_y = cos(y * pi/256) - 1. + double cosm1_y = + ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be45dep-14, 0x1.03c1f0819cac8p-30, + -0x1.55d33e393617dp-48); + // Combine the results with the sine of sum formula: + // sin(x) = sin((k + y)*pi/256) + // = sin(y*pi/256) * cos(k*pi/256) + cos(y*pi/256) * sin(k*pi/256) + return fputil::fma(sin_y, cos_k, fputil::fma(cosm1_y, sin_k, sin_k)); } } // namespace __llvm_libc diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -23,15 +23,19 @@ add_fp_unittest( sinf_test + NO_RUN_POSTBUILD NEED_MPFR SUITE libc_math_exhaustive_tests SRCS sinf_test.cpp DEPENDS + .exhaustive_test libc.include.math libc.src.math.sinf libc.src.__support.FPUtil.fputil + LINK_OPTIONS + -lpthread ) add_fp_unittest( diff --git a/libc/test/src/math/exhaustive/sinf_test.cpp b/libc/test/src/math/exhaustive/sinf_test.cpp --- a/libc/test/src/math/exhaustive/sinf_test.cpp +++ b/libc/test/src/math/exhaustive/sinf_test.cpp @@ -6,20 +6,78 @@ // //===----------------------------------------------------------------------===// +#include "exhaustive_test.h" #include "src/__support/FPUtil/FPBits.h" #include "src/math/sinf.h" #include "utils/MPFRWrapper/MPFRUtils.h" -#include +#include "utils/UnitTest/FPMatcher.h" + +#include using FPBits = __llvm_libc::fputil::FPBits; namespace mpfr = __llvm_libc::testing::mpfr; -TEST(LlvmLibcsinffExhaustiveTest, AllValues) { - uint32_t bits = 0; - do { - FPBits xbits(bits); - float x = float(xbits); - ASSERT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), 1.0); - } while (bits++ < 0xffff'ffffU); +struct LlvmLibcSinfExhaustiveTest : public LlvmLibcExhaustiveTest { + bool check(uint32_t start, uint32_t stop, + mpfr::RoundingMode rounding) override { + mpfr::ForceRoundingMode r(rounding); + uint32_t bits = start; + bool result = true; + do { + FPBits xbits(bits); + float x = float(xbits); + result &= EXPECT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), + 0.5, rounding); + } while (++bits < stop); + return result; + } +}; + +static const int NUM_THREADS = std::thread::hardware_concurrency(); + +// Range: [0, +Inf); +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x7f80'0000U; + +TEST_F(LlvmLibcSinfExhaustiveTest, PostiveRangeRoundNearestTieToEven) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, PostiveRangeRoundUp) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, PostiveRangeRoundDown) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, PostiveRangeRoundTowardZero) { + test_full_range(POS_START, POS_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); +} + +// Range: (-Inf, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xff80'0000U; + +TEST_F(LlvmLibcSinfExhaustiveTest, NegativeRangeRoundNearestTieToEven) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Nearest); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, NegativeRangeRoundUp) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, mpfr::RoundingMode::Upward); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, NegativeRangeRoundDown) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::Downward); +} + +TEST_F(LlvmLibcSinfExhaustiveTest, NegativeRangeRoundTowardZero) { + test_full_range(NEG_START, NEG_STOP, NUM_THREADS, + mpfr::RoundingMode::TowardZero); } diff --git a/libc/test/src/math/sinf_test.cpp b/libc/test/src/math/sinf_test.cpp --- a/libc/test/src/math/sinf_test.cpp +++ b/libc/test/src/math/sinf_test.cpp @@ -51,26 +51,41 @@ float x = float(FPBits(v)); if (isnan(x) || isinf(x)) continue; - ASSERT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), 1.0); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, + __llvm_libc::sinf(x), 0.5); } } TEST(LlvmLibcSinfTest, SpecificBitPatterns) { - float x = float(FPBits(uint32_t(0xc70d39a1))); - EXPECT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), 1.0); + constexpr int N = 9; + constexpr uint32_t INPUTS[N] = { + 0x4049'0fdbU, // x = pi + 0x4619'9998U, // x = 0x1.33333p+13f + 0x4afd'ece4U, // x = 0x1.fbd9c8p+22f + 0x5239'47f6U, // x = 0x1.728fecp+37f + 0x55ca'fb2aU, // x = 0x1.95f654p+44f + 0x588e'f060U, // x = 0x1.1de0cp+50f + 0x6600'0001U, // x = 0x1.000002p+77f + 0x6a19'76f1U, // x = 0x1.32ede2p+85f + 0x6f79'be45U, // x = 0x1.f37c8ap+95f + }; + + for (int i = 0; i < N; ++i) { + float x = float(FPBits(INPUTS[i])); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, + __llvm_libc::sinf(x), 0.5); + } } // For small values, sin(x) is x. TEST(LlvmLibcSinfTest, SmallValues) { - float x = float(FPBits(uint32_t(0x17800000))); - float result = __llvm_libc::sinf(x); - EXPECT_MPFR_MATCH(mpfr::Operation::Sin, x, result, 1.0); - EXPECT_FP_EQ(x, result); - - x = float(FPBits(uint32_t(0x00400000))); - result = __llvm_libc::sinf(x); - EXPECT_MPFR_MATCH(mpfr::Operation::Sin, x, result, 1.0); - EXPECT_FP_EQ(x, result); + float x = float(FPBits(0x1780'0000U)); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), + 0.5); + + x = float(FPBits(0x0040'0000U)); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), + 0.5); } // SDCOMP-26094: check sinf in the cases for which the range reducer @@ -78,6 +93,7 @@ TEST(LlvmLibcSinfTest, SDCOMP_26094) { for (uint32_t v : SDCOMP26094_VALUES) { float x = float(FPBits((v))); - EXPECT_MPFR_MATCH(mpfr::Operation::Sin, x, __llvm_libc::sinf(x), 1.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, + __llvm_libc::sinf(x), 0.5); } }