From a429de81874ab490acdddf5502e1a50d84751927 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Mon, 5 Feb 2024 03:16:08 +0100 Subject: [PATCH] Implement and test asin() --- math/softfloat/soft_double.h | 110 ++++++++++++++++++++++++++++++++--- test/test_soft_double.cpp | 64 +++++++++++++++++++- 2 files changed, 166 insertions(+), 8 deletions(-) diff --git a/math/softfloat/soft_double.h b/math/softfloat/soft_double.h index 2d24c2a..aa24fa1 100644 --- a/math/softfloat/soft_double.h +++ b/math/softfloat/soft_double.h @@ -1,5 +1,5 @@ /////////////////////////////////////////////////////////////////// -// Copyright Christopher Kormanyos 2012 - 2023. // +// Copyright Christopher Kormanyos 2012 - 2024. // // Distributed under the Boost Software License, // // Version 1.0. (See accompanying file LICENSE_1_0.txt // // or copy at http://www.boost.org/LICENSE_1_0.txt) // @@ -592,6 +592,9 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. constexpr auto sin (soft_double x) -> soft_double; constexpr auto cos (soft_double x) -> soft_double; constexpr auto tan (soft_double x) -> soft_double; + constexpr auto asin (soft_double x) -> soft_double; + constexpr auto acos (soft_double x) -> soft_double; + constexpr auto atan (soft_double x) -> soft_double; constexpr auto sinh (soft_double x) -> soft_double; constexpr auto cosh (soft_double x) -> soft_double; constexpr auto tanh (soft_double x) -> soft_double; @@ -614,7 +617,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. using representation_type = std::uint64_t; - constexpr soft_double() noexcept = default; + constexpr soft_double() noexcept : my_value(static_cast(UINT8_C(0))) { } template::value @@ -2409,6 +2412,57 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. return soft_double::my_value_one() + ((x2 * top) / (bot * static_cast(INT8_C(24)))); } + constexpr auto asin_pade(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param) + { + // Usa a Pade approximation of (asin(x) / x). + // Then simplify, extract the coefficients of numerator and denominator. + // Scale both numeroator and denominator by a factor of 10^29. + // Subsequently extract the integer parts of the coefficients. + + // PadeApproximant[ArcSin[x]/x, {x, 0, {12, 12}}] + // FullSimplify[%] + // IntegerPart[N[CoefficientList[Numerator[Out[2]]/10^29, x^2], 24]] + // IntegerPart[N[CoefficientList[Denominator[Out[2]]/10^29, x^2], 24]] + + // Numerator: + // 140095773199074572 + // -388690947968573359 + // 402271633106633823 + // -190922593544635779 + // 40839234741969911 + // -3218334247465525 + // 45053024618672 + + // Denominator: + // 140095773199074572 + // -412040243501752455 + // 460437824033661973 + // -243013488210192408 + // 60946917703998365 + // -6321063389564933 + // 174545740275468 + + const auto x2 = x * x; + + const soft_double top = ((((((( + INT64_C(45053024618672)) + * x2 - INT64_C(3218334247465525)) + * x2 + INT64_C(40839234741969911)) + * x2 - INT64_C(190922593544635779)) + * x2 + INT64_C(402271633106633823)) + * x2 - INT64_C(388690947968573359)) + * x2 + INT64_C(140095773199074572)); + + const soft_double bot = ((((((( + INT64_C(174545740275468)) + * x2 - INT64_C(6321063389564933)) + * x2 + INT64_C(60946917703998365)) + * x2 - INT64_C(243013488210192408)) + * x2 + INT64_C(460437824033661973)) + * x2 - INT64_C(412040243501752455)) + * x2 + INT64_C(140095773199074572)); + + return (x * top) / bot; + } + } // namespace detail constexpr auto sin(soft_double x) -> soft_double // NOLINT(misc-no-recursion) @@ -2607,6 +2661,53 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. return c; } + constexpr auto tan(soft_double x) -> soft_double + { + return sin(x) / cos(x); + } + + constexpr auto asin(soft_double x) -> soft_double // NOLINT(misc-no-recursion) + { + soft_double result { }; + + if(x < static_cast(INT8_C(0))) + { + result = -asin(-x); + } + else if(x > static_cast(INT8_C(0))) + { + if(x > static_cast(INT8_C(1))) + { + result = soft_double::my_value_quiet_NaN(); + } + else if(x < static_cast(INT8_C(1))) + { + if(x < soft_double::my_value_epsilon()) + { + result = x; + } + else if(x < soft_double::my_value_half()) + { + result = detail::asin_pade(x); + } + else + { + result = soft_double::my_value_pi_half() - (detail::asin_pade(sqrt((soft_double::my_value_one() - x) / 2)) * 2); + } + } + else + { + result = soft_double::my_value_pi_half(); + } + } + else + { + // x == 0 and result = 0. + } + + return result; + } + constexpr auto floor(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param) { auto result = soft_double { }; @@ -2781,11 +2882,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. return exp(a * log(x)); // NOLINT(performance-unnecessary-value-param) } - constexpr auto tan(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param) - { - return sin(x) / cos(x); - } - constexpr auto sinh(soft_double x) -> soft_double // NOLINT(performance-unnecessary-value-param) { const soft_double ep = exp(x); // NOLINT(performance-unnecessary-value-param) diff --git a/test/test_soft_double.cpp b/test/test_soft_double.cpp index 7ba6fbb..7dbf3c3 100644 --- a/test/test_soft_double.cpp +++ b/test/test_soft_double.cpp @@ -1,5 +1,5 @@ /////////////////////////////////////////////////////////////////// -// Copyright Christopher Kormanyos 2012 - 2023. // +// Copyright Christopher Kormanyos 2012 - 2024. // // Distributed under the Boost Software License, // // Version 1.0. (See accompanying file LICENSE_1_0.txt // // or copy at http://www.boost.org/LICENSE_1_0.txt) // @@ -454,6 +454,66 @@ auto test_cos() -> bool return result_is_ok; } +auto test_asin() -> bool +{ + bool result_is_ok = true; + + for(int i = 0; i < 50; ++i) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + { + const math::softfloat::float64_t x = math::softfloat::float64_t(i) / 100; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const auto d = static_cast(x); + + using std::asin; + + const math::softfloat::float64_t as_x = asin(x); + const double as_d = asin(d); + + if(i == 0) + { + const auto result_asin_zero_is_ok = (as_x.crepresentation() == 0U); + + result_is_ok = (result_asin_zero_is_ok && result_is_ok); + } + else + { + const double closeness = std::fabs(1.0 - (static_cast(as_x) / as_d)); + + const auto result_asin_is_ok = (closeness < std::numeric_limits::epsilon() * 32.0); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + + result_is_ok = (result_asin_is_ok && result_is_ok); + } + } + + { + const math::softfloat::float64_t x = math::softfloat::float64_t(-33) / 100; + const auto d = static_cast(x); + + using std::asin; + + const math::softfloat::float64_t as_x = asin(x); + const double as_d = asin(d); + + const double closeness = std::fabs(1.0 - (static_cast(as_x) / as_d)); + + const auto result_asin_is_ok = ((as_x < 0) && (closeness < std::numeric_limits::epsilon() * 32.0)); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + + result_is_ok = (result_asin_is_ok && result_is_ok); + } + + for(int i = 121; i < 124; ++i) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + { + const math::softfloat::float64_t x = math::softfloat::float64_t(i) / 100; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + + const math::softfloat::float64_t as_x = asin(x); + + const auto result_asin_is_ok = isnan(as_x); + + result_is_ok = (result_asin_is_ok && result_is_ok); + } + + return result_is_ok; +} + auto test_floor(const std::uint32_t n) -> bool { bool result_is_ok = true; @@ -751,6 +811,7 @@ auto test_soft_double() -> bool std::cout << "testing log____... "; const auto result_log____is_ok = test_log ( ); std::cout << std::boolalpha << result_log____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp) std::cout << "testing sin____... "; const auto result_sin____is_ok = test_sin ( ); std::cout << std::boolalpha << result_sin____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp) std::cout << "testing cos____... "; const auto result_cos____is_ok = test_cos ( ); std::cout << std::boolalpha << result_cos____is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp) + std::cout << "testing asin___... "; const auto result_asin___is_ok = test_asin ( ); std::cout << std::boolalpha << result_asin___is_ok << std::endl; // NOLINT(bugprone-reserved-identifier,cert-dcl37-c,cert-dcl51-cpp) local::eng32.seed(::util::util_pseudorandom_time_point_seed::value()); local::eng64.seed(::util::util_pseudorandom_time_point_seed::value()); @@ -798,6 +859,7 @@ auto test_soft_double() -> bool && result_log____is_ok && result_sin____is_ok && result_cos____is_ok + && result_asin___is_ok && result_floor__is_ok && result_ceil___is_ok && result_add____is_ok