diff --git a/include/boost/decimal/detail/cmath/asin.hpp b/include/boost/decimal/detail/cmath/asin.hpp index 4b6febdeb..eecaf665f 100644 --- a/include/boost/decimal/detail/cmath/asin.hpp +++ b/include/boost/decimal/detail/cmath/asin.hpp @@ -1,4 +1,5 @@ -// Copyright 2024 Matt Borland +// Copyright 2024 - 2025 Matt Borland +// Copyright 2024 - 2025 Christopher Kormanyos // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt @@ -41,29 +42,42 @@ constexpr auto asin_impl(const T x) noexcept return x; } - const auto absx {fabs(x)}; - T result {}; + const auto absx { fabs(x) }; - if (absx <= std::numeric_limits::epsilon()) + T result { }; + + constexpr T cbrt_eps { cbrt(std::numeric_limits::epsilon()) }; + + constexpr T one { 1 }; + + if (absx <= cbrt_eps) { - result = absx; + result = absx * (one + (absx / 6) * absx); } - else if (absx <= T{5, -1}) + else if (absx <= T { 5, -1 }) { result = asin_series(absx); } - else if (absx <= T{1, 0}) - { - constexpr T half_pi {numbers::pi_v / 2}; - result = half_pi - 2 * asin_series(sqrt((1 - absx) / 2)); - } else { - #ifndef BOOST_DECIMAL_FAST_MATH - result = std::numeric_limits::quiet_NaN(); - #else - result = T{0}; - #endif + constexpr T half_pi { numbers::pi_v / 2 }; + + if (absx < one) + { + result = half_pi - 2 * asin_series(sqrt((1 - absx) / 2)); + } + else if (absx > one) + { + #ifndef BOOST_DECIMAL_FAST_MATH + result = std::numeric_limits::quiet_NaN(); + #else + result = T{0}; + #endif + } + else + { + result = half_pi; + } } // arcsin(-x) == -arcsin(x) diff --git a/test/test_asin.cpp b/test/test_asin.cpp index 9d3ae4427..3f9f0181e 100644 --- a/test/test_asin.cpp +++ b/test/test_asin.cpp @@ -1,4 +1,5 @@ -// Copyright 2024 Matt Borland +// Copyright 2024 - 2025 Matt Borland +// Copyright 2024 - 2025 Christopher Kormanyos // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt @@ -14,7 +15,7 @@ # pragma clang diagnostic ignored "-Wconversion" # pragma clang diagnostic ignored "-Wsign-conversion" # pragma clang diagnostic ignored "-Wfloat-equal" -# if __clang_major__ >= 20 +# if (__clang_major__ >= 20) # pragma clang diagnostic ignored "-Wfortify-source" # endif #elif defined(__GNUC__) @@ -28,9 +29,9 @@ #include #include -#include -#include + #include +#include #if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH) static constexpr auto N = static_cast(128U); // Number of trials @@ -42,6 +43,9 @@ static std::mt19937_64 rng(42); using namespace boost::decimal; +template auto my_zero() -> T; +template auto my_one () -> T; + template void test_asin() { @@ -137,6 +141,100 @@ void print_value(T value, const char* str) << "\nExp: " << ptr << "\n" << std::endl; } +template +auto test_asin_edge() -> void +{ + using nl = std::numeric_limits; + + const T tiny0 { nl::epsilon() * 999 / 1000 }; + const T tiny1 { nl::epsilon() }; + const T tiny2 { nl::epsilon() * 1000 / 999 }; + + const T asin_tiny0 { asin(tiny0) }; + const T asin_tiny1 { asin(tiny1) }; + const T asin_tiny2 { asin(tiny2) }; + + // tiny1: 1 + // tiny2: 1.001001 + // tiny0: 0.999 + // tiny1: 1 + // tiny2: 1.001001001001001 + // tiny0: 0.999 + // tiny1: 1 + // tiny2: 1.001001001001001001001001001001001 + + constexpr T ctrl_tiny2 + { + std::numeric_limits::digits10 < 10 ? T("1.001001") + : std::numeric_limits::digits10 < 20 ? T("1.001001001001001") + : T("1.001001001001001001001001001001001") + }; + + BOOST_TEST_EQ(asin_tiny0 / nl::epsilon(), T(999, -3)); + BOOST_TEST_EQ(asin_tiny1 / nl::epsilon(), T(1)); + BOOST_TEST_EQ(asin_tiny2 / nl::epsilon(), ctrl_tiny2); + + constexpr T half_pi { numbers::pi_v / 2 }; + + BOOST_TEST_EQ(asin(my_zero() + my_one()), half_pi); + BOOST_TEST_EQ(asin(my_zero() - my_one()), -half_pi); +} + +template +void test_asin_1137() +{ + using nl = std::numeric_limits; + + const T tiny0 { nl::epsilon() * 999/1000 }; + const T tiny1 { nl::epsilon() }; + const T tiny2 { nl::epsilon() * 1000/999 }; + + BOOST_TEST(tiny0 != tiny1); + BOOST_TEST(tiny1 != tiny2); + + std::stringstream strm { }; + + BOOST_TEST_EQ(tiny0, asin(tiny0)); + BOOST_TEST_EQ(tiny1, asin(tiny1)); + BOOST_TEST_EQ(tiny2, asin(tiny2)); + + const T sqrt_tiny0 { sqrt(nl::epsilon() * 999/1000) }; + const T sqrt_tiny1 { sqrt(nl::epsilon()) }; + const T sqrt_tiny2 { sqrt(nl::epsilon() * 1000/999) }; + + BOOST_TEST_EQ(sqrt_tiny0, asin(sqrt_tiny0)); + BOOST_TEST_EQ(sqrt_tiny1, asin(sqrt_tiny1)); + BOOST_TEST_EQ(sqrt_tiny2, asin(sqrt_tiny2)); + + const T cbrt_tiny0 { cbrt(nl::epsilon() * 999/1000) }; + const T cbrt_tiny1 { cbrt(nl::epsilon()) }; + const T cbrt_tiny2 { cbrt(nl::epsilon() * 1000/999) }; + const T cbrt_tiny3 { cbrt(nl::epsilon() * 1004/999) }; + + auto mini_series + { + [](const T eps) + { + return eps * (1 + (eps / 6) * eps); + } + }; + + auto is_close + { + [](const T a, const T b) + { + const T delta { fabs(a - b) }; + + return (delta < (std::numeric_limits::epsilon() * 4)); + } + }; + + BOOST_TEST(is_close(asin(cbrt_tiny0), mini_series(cbrt_tiny0))); + BOOST_TEST(is_close(asin(cbrt_tiny1), mini_series(cbrt_tiny1))); + BOOST_TEST(is_close(asin(cbrt_tiny2), mini_series(cbrt_tiny2))); + BOOST_TEST(is_close(asin(cbrt_tiny3), mini_series(cbrt_tiny3))); +} + int main() { #ifdef BOOST_DECIMAL_GENERATE_CONSTANT_SIGS @@ -187,12 +285,22 @@ int main() test_asin(); test_asin(); - #if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH) test_asin(); #endif test_asin(); + test_asin_edge(); + test_asin_edge(); + test_asin_edge(); + + test_asin_1137(); + test_asin_1137(); + test_asin_1137(); + return boost::report_errors(); } + +template auto my_zero() -> T { T zero { 0 }; return zero; } +template auto my_one () -> T { T one { 1 }; return one; }