From 2ab9315e179159b5761accb1706184fca7a744d2 Mon Sep 17 00:00:00 2001 From: Shen-Ta Hsieh Date: Tue, 29 Sep 2026 04:40:53 +0800 Subject: [PATCH] Reduce the argument of sin, cos and tan exactly and round once - The reduction used pi with only the digits of the type and an unsigned quotient, thus results lost digits from x = 1 and became garbage, inf or NaN above about 1e9. - The decimal32 and decimal64 sine polynomials were not accurate, and the kernels were chains of fma, which are slow. - sin and cos of an infinity gave an infinity, but they must give NaN. - The new impl/trig_reduce.hpp reduces with the digits of 2/pi in base-1e9 integer words. One table of 720 words covers the full range of all types. - For |x| < 10^19, a second reduction uses 64-bit words of 2/pi, which is 10 to 30 percent faster. It gives r in fixed point, thus the kernel does not convert it. Its tables have 352 bytes, and the guard words come from the worst cases below 2^64: 27, 58 and 117 leading zero bits. - The new impl/trig_fixed_point.hpp computes the kernels in binary fixed point: one word for decimal32, one or two words for decimal64 and up to three words for decimal128. The result is rounded once in the current rounding mode. - The checks for 0 and for |x| < pi/4 use the integer significand, not decimal compares. tan uses one reduction, and its reciprocal is in the new impl/tan_impl.hpp. - The results are correctly rounded for the 32-bit and 64-bit types, and all three functions are faster for all types. - The new test test_trig_rounding compares with exact values, also for the worst cases below 10^19 and the inputs on each side of 10^19. It walks chains of adjacent inputs in three rounding modes to test that the results do not go against the exact direction. It skips the checks that need constant evaluation or that fast math turns off. --- include/boost/decimal/detail/cmath/cos.hpp | 96 +-- .../decimal/detail/cmath/impl/cos_impl.hpp | 240 +++--- .../decimal/detail/cmath/impl/sin_impl.hpp | 254 +++---- .../decimal/detail/cmath/impl/tan_impl.hpp | 98 +++ .../detail/cmath/impl/trig_fixed_point.hpp | 315 ++++++++ .../decimal/detail/cmath/impl/trig_reduce.hpp | 682 ++++++++++++++++++ include/boost/decimal/detail/cmath/sin.hpp | 131 +--- include/boost/decimal/detail/cmath/tan.hpp | 109 +-- test/Jamfile | 1 + test/test_edges_and_behave.cpp | 4 +- test/test_sin_cos.cpp | 4 +- test/test_trig_rounding.cpp | 418 +++++++++++ 12 files changed, 1739 insertions(+), 613 deletions(-) create mode 100644 include/boost/decimal/detail/cmath/impl/tan_impl.hpp create mode 100644 include/boost/decimal/detail/cmath/impl/trig_fixed_point.hpp create mode 100644 include/boost/decimal/detail/cmath/impl/trig_reduce.hpp create mode 100644 test/test_trig_rounding.cpp diff --git a/include/boost/decimal/detail/cmath/cos.hpp b/include/boost/decimal/detail/cmath/cos.hpp index 40aad502b..b3b62a829 100644 --- a/include/boost/decimal/detail/cmath/cos.hpp +++ b/include/boost/decimal/detail/cmath/cos.hpp @@ -7,18 +7,17 @@ #define BOOST_DECIMAL_DETAIL_CMATH_COS_HPP #include -#include #include #include #include -#include -#include +#include +#include #include #include #ifndef BOOST_DECIMAL_BUILD_MODULE +#include #include -#include #endif namespace boost { @@ -30,90 +29,27 @@ template constexpr auto cos_impl(const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { - T result { }; - #ifndef BOOST_DECIMAL_FAST_MATH - const auto fpc = fpclassify(x); - - // First check non-finite values and small angles - - if ((fpc == FP_INFINITE) || (fpc == FP_NAN)) + if (isnan(x)) { - result = x; + return x; } - else - #endif - if (signbit(x)) + if (isinf(x)) { - result = cos(-x); + return std::numeric_limits::quiet_NaN(); } - else - { - if (x > std::numeric_limits::epsilon()) - { - // Perform argument reduction and subsequent scaling of the result. - - // Given x = k * (pi/2) + r, compute n = (k % 4). - - // | n | sin(x) | cos(x) | sin(x)/cos(x) | - // |----------------------------------------| - // | 0 | sin(r) | cos(r) | sin(r)/cos(r) | - // | 1 | cos(r) | -sin(r) | -cos(r)/sin(r) | - // | 2 | -sin(r) | -cos(r) | sin(r)/cos(r) | - // | 3 | -cos(r) | sin(r) | -cos(r)/sin(r) | - - const T two_x = x * 2; - - const auto k = static_cast(two_x / numbers::pi_v); - const auto n = k % static_cast(UINT8_C(4)); - - const T two_r { two_x - (numbers::pi_v * k) }; - - T r { two_r / 2 }; - - constexpr T half { 5, -1 }; - - const bool do_scaling { r > half }; - - if(do_scaling) - { - // Reduce the argument with factors of three. - r /= static_cast(UINT8_C(3)); - } - - switch(n) - { - case static_cast(UINT8_C(1)): - case static_cast(UINT8_C(3)): - result = detail::sin_series_expansion(r); - break; - case static_cast(UINT8_C(0)): - case static_cast(UINT8_C(2)): - default: - result = detail::cos_series_expansion(r); - break; - } - - if(do_scaling) - { - result *= (((result * result) * static_cast(UINT8_C(4))) - static_cast(UINT8_C(3))); - } - - if(signbit(result)) { result = -result; } - - const auto b_neg = ((n == static_cast(UINT8_C(1))) || (n == static_cast(UINT8_C(2)))); - - if(b_neg) { result = -result; } - } - else - { - constexpr T one { 1 }; + #endif - result = one; - } + // x = n*pi/2 + r: cos(x) is cos(r), -sin(r), -cos(r) or sin(r) for n = 0 to 3. + const auto r {trig::trig_prepare(x)}; + if (r.zero) + { + return T {1}; } - return result; + // The sign of sin(r) is r.neg, and cos(r) > 0. + return (r.n & 1U) == 0U ? trig::cos_of(r, r.n == 2U) + : trig::sin_of(r, r.neg != (r.n == 1U)); } } // namespace detail diff --git a/include/boost/decimal/detail/cmath/impl/cos_impl.hpp b/include/boost/decimal/detail/cmath/impl/cos_impl.hpp index cc68d1716..0bc226c93 100644 --- a/include/boost/decimal/detail/cmath/impl/cos_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/cos_impl.hpp @@ -6,204 +6,138 @@ #define BOOST_DECIMAL_DETAIL_CMATH_IMPL_COS_IMPL_HPP #include -#include -#include -#include #include -#include +#include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE -#include #include #endif namespace boost { namespace decimal { namespace detail { - -namespace cos_detail { +namespace trig { template struct cos_table_imp { - // 8th Degree Remez Polynomial from 0 to pi / 4 - // Estimated max error: 4.321978891364628e-14 - static constexpr std::array d32_coeffs = - {{ - decimal32_t {UINT64_C(22805960529562646), -21}, - decimal32_t {UINT64_C(39171880037888081), -22}, - decimal32_t {UINT64_C(1392392773950284), -18, construction_sign::negative}, - decimal32_t {UINT64_C(17339629614857501), -22}, - decimal32_t {UINT64_C(41666173896377827), -18}, - decimal32_t {UINT64_C(77764646000512304), -24}, - decimal32_t {UINT64_C(50000000610949535), -17, construction_sign::negative}, - decimal32_t {UINT64_C(18421494272283811), -26}, - decimal32_t {UINT64_C(99999999999908662), -17} - }}; - - static constexpr std::array d32_fast_coeffs = - {{ - decimal_fast32_t {UINT64_C(22805960529562646), -21}, - decimal_fast32_t {UINT64_C(39171880037888081), -22}, - decimal_fast32_t {UINT64_C(1392392773950284), -18, construction_sign::negative}, - decimal_fast32_t {UINT64_C(17339629614857501), -22}, - decimal_fast32_t {UINT64_C(41666173896377827), -18}, - decimal_fast32_t {UINT64_C(77764646000512304), -24}, - decimal_fast32_t {UINT64_C(50000000610949535), -17, construction_sign::negative}, - decimal_fast32_t {UINT64_C(18421494272283811), -26}, - decimal_fast32_t {UINT64_C(99999999999908662), -17} - }}; - - // 12th Degree Remez Polynomial from 0 to pi / 4 - // Estimated max error: 7.911867233315355155595617164843665e-20 - static constexpr std::array d64_coeffs = - {{ - decimal64_t {UINT64_C(1922641020040661424), -27}, - decimal64_t {UINT64_C(4960385936049718134), -28}, - decimal64_t {UINT64_C(2763064713566851512), -25, construction_sign::negative}, - decimal64_t {UINT64_C(6633276621376137827), -28}, - decimal64_t {UINT64_C(2480119161297283187), -23}, - decimal64_t {UINT64_C(1600210781837650114), -28}, - decimal64_t {UINT64_C(1388888932852646133), -21, construction_sign::negative}, - decimal64_t {UINT64_C(8054772849254568869), -30}, - decimal64_t {UINT64_C(4166666666572238908), -20}, - decimal64_t {UINT64_C(6574164404618517322), -32}, - decimal64_t {UINT64_C(5000000000000023748), -19, construction_sign::negative}, - decimal64_t {UINT64_C(3367952043014273196), -35}, - decimal64_t {UINT64_C(9999999999999999999), -19} - }}; - - static constexpr std::array d64_fast_coeffs = - {{ - decimal_fast64_t {UINT64_C(1922641020040661424), -27}, - decimal_fast64_t {UINT64_C(4960385936049718134), -28}, - decimal_fast64_t {UINT64_C(2763064713566851512), -25, construction_sign::negative}, - decimal_fast64_t {UINT64_C(6633276621376137827), -28}, - decimal_fast64_t {UINT64_C(2480119161297283187), -23}, - decimal_fast64_t {UINT64_C(1600210781837650114), -28}, - decimal_fast64_t {UINT64_C(1388888932852646133), -21, construction_sign::negative}, - decimal_fast64_t {UINT64_C(8054772849254568869), -30}, - decimal_fast64_t {UINT64_C(4166666666572238908), -20}, - decimal_fast64_t {UINT64_C(6574164404618517322), -32}, - decimal_fast64_t {UINT64_C(5000000000000023748), -19, construction_sign::negative}, - decimal_fast64_t {UINT64_C(3367952043014273196), -35}, - decimal_fast64_t {UINT64_C(9999999999999999999), -19} - }}; + // cos(r) = 1 - z * (1/2 - z * C(z)) with z = r^2 <= (pi/4)^2; the table has |coefficient| of C. + static constexpr fx<1> d32[5] = + { + // cos degree 4 on [0,0.61685] max rel err 9.2329e-17 + {{UINT64_C(0x02AAAAAAAAA5BB5B)}}, + {{UINT64_C(0x0016C16C16721120)}}, + {{UINT64_C(0x000068067E97FEFC)}}, + {{UINT64_C(0x00000127E00B90DA)}}, + {{UINT64_C(0x00000002377D1C31)}}, + }; + + // From z^2 on, one word scaled by 2^(62 + d64_shift) is enough for decimal64. + static constexpr int d64_shift {17}; + static constexpr fx<2> d64_head[2] = + { + {{UINT64_C(0x66A629652876004D), UINT64_C(0x02AAAAAAAAAAAAAA)}}, + {{UINT64_C(0x6CDA6F0ECA1D2ECE), UINT64_C(0x0016C16C16C16C0F)}}, + }; + + static constexpr std::uint64_t d64_tail[5] = { UINT64_C(0xD00D00D00C6653EB), UINT64_C(0x024FC9F6EBD7192B), UINT64_C(0x00047BB632432CE0), UINT64_C(0x0000064E4C907AB4), UINT64_C(0x00000006AAF461CB) }; + + // From z^2 on, 128 bits scaled by 2^(126 + d128_mid_shift), and from z^10 on, one word scaled by + // 2^(62 + d128_tail_shift). + static constexpr int d128_mid_shift {17}; + static constexpr int d128_tail_shift {81}; + static constexpr fx<3> d128_head[2] = + { + {{UINT64_C(0xEA45AD4A0D5AF6A0), UINT64_C(0xAAAAAAAAAAAAAA7A), UINT64_C(0x02AAAAAAAAAAAAAA)}}, + {{UINT64_C(0xB80DD449D4ED2C61), UINT64_C(0xC16C16C16C16B484), UINT64_C(0x0016C16C16C16C16)}}, + }; + + static constexpr boost::int128::uint128_t d128_mid[8] = + { + boost::int128::uint128_t {UINT64_C(0xD00D00D00D00D00D), UINT64_C(0x00D00CFE0B2E1E72)}, + boost::int128::uint128_t {UINT64_C(0x024FC9F6EF13EB8E), UINT64_C(0x5DE02D7E98F54548)}, + boost::int128::uint128_t {UINT64_C(0x00047BB63BFE3625), UINT64_C(0xED51352B37C62468)}, + boost::int128::uint128_t {UINT64_C(0x0000064E5D2A301F), UINT64_C(0x274825BC9C08DEA7)}, + boost::int128::uint128_t {UINT64_C(0x00000006B9FCF9CC), UINT64_C(0xEE079F20A4359338)}, + boost::int128::uint128_t {UINT64_C(0x0000000005A09E18), UINT64_C(0xEE5EF9C706842ACD)}, + boost::int128::uint128_t {UINT64_C(0x000000000003CA85), UINT64_C(0x747F6661142C8245)}, + boost::int128::uint128_t {UINT64_C(0x0000000000000219), UINT64_C(0xC72C8B0FDE1153D3)}, + }; + + static constexpr std::uint64_t d128_tail[2] = { UINT64_C(0xF96673E127BADABC), UINT64_C(0x0061AC811D606933) }; }; #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) template -constexpr std::array cos_table_imp::d32_coeffs; +constexpr fx<1> cos_table_imp::d32[5]; template -constexpr std::array cos_table_imp::d64_coeffs; +constexpr int cos_table_imp::d64_shift; template -constexpr std::array cos_table_imp::d32_fast_coeffs; +constexpr int cos_table_imp::d128_mid_shift; template -constexpr std::array cos_table_imp::d64_fast_coeffs; +constexpr int cos_table_imp::d128_tail_shift; -#endif +template +constexpr fx<2> cos_table_imp::d64_head[2]; -using cos_table = cos_table_imp; +template +constexpr std::uint64_t cos_table_imp::d64_tail[5]; + +template +constexpr fx<3> cos_table_imp::d128_head[2]; -} // namespace cos_detail +template +constexpr boost::int128::uint128_t cos_table_imp::d128_mid[8]; -template -constexpr auto cos_series_expansion(T x) noexcept; +template +constexpr std::uint64_t cos_table_imp::d128_tail[2]; -template <> -constexpr auto cos_series_expansion(decimal32_t x) noexcept -{ - return remez_series_result(x, cos_detail::cos_table::d32_coeffs); -} +#endif + +using cos_table = cos_table_imp; -template <> -constexpr auto cos_series_expansion(decimal_fast32_t x) noexcept +constexpr auto cos_poly(const fx<1>& z) noexcept -> fx<1> { - return remez_series_result(x, cos_detail::cos_table::d32_fast_coeffs); + return fx_alternating(z, cos_table::d32); } -template <> -constexpr auto cos_series_expansion(decimal64_t x) noexcept +constexpr auto cos_poly(const fx<2>& z) noexcept -> fx<2> { - return remez_series_result(x, cos_detail::cos_table::d64_coeffs); + return fx_alternating_split(z, cos_table::d64_head, cos_table::d64_tail); } -template <> -constexpr auto cos_series_expansion(decimal_fast64_t x) noexcept +constexpr auto cos_poly(const fx<3>& z) noexcept -> fx<3> { - return remez_series_result(x, cos_detail::cos_table::d64_fast_coeffs); + return fx_alternating_split(z, cos_table::d128_head, cos_table::d128_mid, cos_table::d128_tail); } -template <> -constexpr auto cos_series_expansion(decimal128_t x) noexcept +// cos(r) for z = r^2, which is below 1 for any r != 0. +template +constexpr auto cos_fx(const fx& z) noexcept -> fx { - // PadeApproximant[Cos[x], {x, 0, {14, 14}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - constexpr decimal128_t c0 { boost::int128::uint128_t { UINT64_C(307807346375396), UINT64_C(9191352932158695424) }, 3 }; - constexpr decimal128_t c1 { boost::int128::uint128_t { UINT64_C(149996550055690), UINT64_C(222763958071016960) }, 3, true }; - constexpr decimal128_t c2 { boost::int128::uint128_t { UINT64_C(108967212479807), UINT64_C(3937477076487471608) }, 2 }; - constexpr decimal128_t c3 { boost::int128::uint128_t { UINT64_C(277096228519262), UINT64_C(6277888927557284608) }, 0, true }; - constexpr decimal128_t c4 { boost::int128::uint128_t { UINT64_C(319580269604048), UINT64_C(10708241405247058432) }, -2 }; - constexpr decimal128_t c5 { boost::int128::uint128_t { UINT64_C(183739194803716), UINT64_C(9003931728965394944) }, -4, true }; - constexpr decimal128_t c6 { boost::int128::uint128_t { UINT64_C(518817586019902), UINT64_C(14598542072727738368) }, -7 }; - constexpr decimal128_t c7 { boost::int128::uint128_t { UINT64_C(58205916937364), UINT64_C(13388002334603019776) }, -9, true }; - - constexpr decimal128_t d1 { boost::int128::uint128_t { UINT64_C(390712313200823), UINT64_C(13016137105513388032) }, 1 }; - constexpr decimal128_t d2 { boost::int128::uint128_t { UINT64_C(249767150099857), UINT64_C(14534865724066009088) }, -1 }; - constexpr decimal128_t d3 { boost::int128::uint128_t { UINT64_C(105535117882474), UINT64_C(16245151810017622016) }, -3 }; - constexpr decimal128_t d4 { boost::int128::uint128_t { UINT64_C(322928599993793), UINT64_C(8055050913586880512) }, -6 }; - constexpr decimal128_t d5 { boost::int128::uint128_t { UINT64_C(72777849685460), UINT64_C(10172723920765296640) }, -8 }; - constexpr decimal128_t d6 { boost::int128::uint128_t { UINT64_C(114133059907344), UINT64_C(3036923607254532096) }, -11 }; - constexpr decimal128_t d7 { boost::int128::uint128_t { UINT64_C(98470690251347), UINT64_C(1521187190289973248) }, -14 }; - - const decimal128_t x2 { x * x }; - - const decimal128_t top { c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * (c6 + x2 * c7)))))) }; - const decimal128_t bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * (d6 + x2 * d7)))))) }; - - return decimal128_t { top / bot }; + fx half {}; + half.w[N - 1] = UINT64_C(1) << 61U; + const auto inner {fx_sub(half, fx_mul(z, cos_poly(z)))}; + return fx_clamp_below_one(fx_sub(fx_one(), fx_mul(z, inner))); } -template <> -constexpr auto cos_series_expansion(decimal_fast128_t x) noexcept +// cos(r) with the sign neg, for r != 0. +template +constexpr auto cos_of(const trig_arg::words>& r, const bool neg) noexcept -> T { - // PadeApproximant[Cos[x], {x, 0, {14, 14}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - constexpr decimal_fast128_t c0 { boost::int128::uint128_t { UINT64_C(307807346375396), UINT64_C(9191352932158695424) }, 3 }; - constexpr decimal_fast128_t c1 { boost::int128::uint128_t { UINT64_C(149996550055690), UINT64_C(222763958071016960) }, 3, true }; - constexpr decimal_fast128_t c2 { boost::int128::uint128_t { UINT64_C(108967212479807), UINT64_C(3937477076487471608) }, 2 }; - constexpr decimal_fast128_t c3 { boost::int128::uint128_t { UINT64_C(277096228519262), UINT64_C(6277888927557284608) }, 0, true }; - constexpr decimal_fast128_t c4 { boost::int128::uint128_t { UINT64_C(319580269604048), UINT64_C(10708241405247058432) }, -2 }; - constexpr decimal_fast128_t c5 { boost::int128::uint128_t { UINT64_C(183739194803716), UINT64_C(9003931728965394944) }, -4, true }; - constexpr decimal_fast128_t c6 { boost::int128::uint128_t { UINT64_C(518817586019902), UINT64_C(14598542072727738368) }, -7 }; - constexpr decimal_fast128_t c7 { boost::int128::uint128_t { UINT64_C(58205916937364), UINT64_C(13388002334603019776) }, -9, true }; - - constexpr decimal_fast128_t d1 { boost::int128::uint128_t { UINT64_C(390712313200823), UINT64_C(13016137105513388032) }, 1 }; - constexpr decimal_fast128_t d2 { boost::int128::uint128_t { UINT64_C(249767150099857), UINT64_C(14534865724066009088) }, -1 }; - constexpr decimal_fast128_t d3 { boost::int128::uint128_t { UINT64_C(105535117882474), UINT64_C(16245151810017622016) }, -3 }; - constexpr decimal_fast128_t d4 { boost::int128::uint128_t { UINT64_C(322928599993793), UINT64_C(8055050913586880512) }, -6 }; - constexpr decimal_fast128_t d5 { boost::int128::uint128_t { UINT64_C(72777849685460), UINT64_C(10172723920765296640) }, -8 }; - constexpr decimal_fast128_t d6 { boost::int128::uint128_t { UINT64_C(114133059907344), UINT64_C(3036923607254532096) }, -11 }; - constexpr decimal_fast128_t d7 { boost::int128::uint128_t { UINT64_C(98470690251347), UINT64_C(1521187190289973248) }, -14 }; - - const decimal_fast128_t x2 { x * x }; - - const decimal_fast128_t top { c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * (c6 + x2 * c7)))))) }; - const decimal_fast128_t bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * (d6 + x2 * d7)))))) }; - - return decimal_fast128_t { top / bot }; + const auto rf {fixed_r(r)}; + const auto f {cos_fx(fx_mul(rf, rf))}; + return trig_round(fx_scale(pow10(static_cast(38)), f), -38, neg); } +} // namespace trig } // namespace detail } // namespace decimal } // namespace boost diff --git a/include/boost/decimal/detail/cmath/impl/sin_impl.hpp b/include/boost/decimal/detail/cmath/impl/sin_impl.hpp index d90617753..48b01ad08 100644 --- a/include/boost/decimal/detail/cmath/impl/sin_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/sin_impl.hpp @@ -6,218 +6,136 @@ #define BOOST_DECIMAL_DETAIL_CMATH_IMPL_SIN_IMPL_HPP #include -#include -#include -#include #include -#include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE -#include #include #endif namespace boost { namespace decimal { namespace detail { - -namespace sin_detail { +namespace trig { template -struct sin_table_imp { - - // 5th Degree Remez Polynomial - // Estimated max error: 6.0855992690454531e-8 - static constexpr std::array d32_coeffs = - {{ - decimal32_t {UINT64_C(76426704684128569), -19}, - decimal32_t {UINT64_C(8163484279370784), -19}, - decimal32_t {UINT64_C(16704305092800237), -17, construction_sign::negative}, - decimal32_t {UINT64_C(74622903795259856), -21}, - decimal32_t {UINT64_C(9999946918542727), -16}, - decimal32_t {UINT64_C(60055992690454536), -24} - }}; - - static constexpr std::array d32_fast_coeffs = - {{ - decimal_fast32_t {UINT64_C(76426704684128569), -19}, - decimal_fast32_t {UINT64_C(8163484279370784), -19}, - decimal_fast32_t {UINT64_C(16704305092800237), -17, construction_sign::negative}, - decimal_fast32_t {UINT64_C(74622903795259856), -21}, - decimal_fast32_t {UINT64_C(9999946918542727), -16}, - decimal_fast32_t {UINT64_C(60055992690454536), -24} - }}; - - // 11th Degree Remez Polynomial - // Estimated max error: 5.2301715421592162270336342660001217e-18 - static constexpr std::array d64_coeffs = - {{ - decimal64_t {UINT64_C(2306518628003855678), -26, construction_sign::negative}, - decimal64_t {UINT64_C(5453073257634027470), -27, construction_sign::negative}, - decimal64_t {UINT64_C(2762996699568163845), -24}, - decimal64_t {UINT64_C(5023027013521532307), -27, construction_sign::negative}, - decimal64_t {UINT64_C(1984096861383546182), -22, construction_sign::negative}, - decimal64_t {UINT64_C(1026912296061211491), -27, construction_sign::negative}, - decimal64_t {UINT64_C(8333333562151404340), -21}, - decimal64_t {UINT64_C(3217043986646625014), -29, construction_sign::negative}, - decimal64_t {UINT64_C(1666666666640042905), -19, construction_sign::negative}, - decimal64_t {UINT64_C(1135995742940218051), -31, construction_sign::negative}, - decimal64_t {UINT64_C(1000000000000001896), -18}, - decimal64_t {UINT64_C(5230171542159216227), -36, construction_sign::negative} - }}; - - static constexpr std::array d64_fast_coeffs = - {{ - decimal_fast64_t {UINT64_C(2306518628003855678), -26, construction_sign::negative}, - decimal_fast64_t {UINT64_C(5453073257634027470), -27, construction_sign::negative}, - decimal_fast64_t {UINT64_C(2762996699568163845), -24}, - decimal_fast64_t {UINT64_C(5023027013521532307), -27, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1984096861383546182), -22, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1026912296061211491), -27, construction_sign::negative}, - decimal_fast64_t {UINT64_C(8333333562151404340), -21}, - decimal_fast64_t {UINT64_C(3217043986646625014), -29, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1666666666640042905), -19, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1135995742940218051), -31, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1000000000000001896), -18}, - decimal_fast64_t {UINT64_C(5230171542159216227), -36, construction_sign::negative} - }}; +struct sin_table_imp +{ + // sin(r) = r * (1 - z * S(z)) with z = r^2 <= (pi/4)^2; the table has |coefficient| of S. + static constexpr fx<1> d32[5] = + { + // sin degree 4 on [0,0.61685] max rel err 4.9991e-15 + {{UINT64_C(0x0AAAAAAAAA911C91)}}, + {{UINT64_C(0x0088888886439978)}}, + {{UINT64_C(0x00034033F272D0EC)}}, + {{UINT64_C(0x00000B8EBB68B293)}}, + {{UINT64_C(0x0000001A961A1733)}}, + }; + + // From z^2 on, one word scaled by 2^(62 + d64_shift) is enough for decimal64. + static constexpr int d64_shift {14}; + static constexpr fx<2> d64_head[2] = + { + {{UINT64_C(0xAA7D60D8D4E5483D), UINT64_C(0x0AAAAAAAAAAAAAAA)}}, + {{UINT64_C(0x7F5132D35AC08297), UINT64_C(0x0088888888888888)}}, + }; + + static constexpr std::uint64_t d64_tail[6] = { UINT64_C(0xD00D00D00D00A6F5), UINT64_C(0x02E3BC74AAD78393), UINT64_C(0x0006B99159F6B038), UINT64_C(0x00000B0922F75776), UINT64_C(0x0000000D73DC0530), UINT64_C(0x000000000C8F4EE3) }; + + // From z^3 on, 128 bits scaled by 2^(126 + d128_mid_shift), and from z^10 on, one word scaled by + // 2^(62 + d128_tail_shift). + static constexpr int d128_mid_shift {20}; + static constexpr int d128_tail_shift {76}; + static constexpr fx<3> d128_head[3] = + { + {{UINT64_C(0x0DD7D16E06FD9B5F), UINT64_C(0xAAAAAAAAAAAAAA05), UINT64_C(0x0AAAAAAAAAAAAAAA)}}, + {{UINT64_C(0x6725275D2EE91D9D), UINT64_C(0x888888888888419A), UINT64_C(0x0088888888888888)}}, + {{UINT64_C(0x717E4EB88E02AEE2), UINT64_C(0x4034034033F8A598), UINT64_C(0x0003403403403403)}}, + }; + + static constexpr boost::int128::uint128_t d128_mid[7] = + { + boost::int128::uint128_t {UINT64_C(0xB8EF1D2AB6399C7D), UINT64_C(0x560E37B7D30EF961)}, + boost::int128::uint128_t {UINT64_C(0x01AE64567F544E38), UINT64_C(0xFE73EEF2F07C2439)}, + boost::int128::uint128_t {UINT64_C(0x0002C248C2750DA1), UINT64_C(0x2F906C4EE1A467D8)}, + boost::int128::uint128_t {UINT64_C(0x0000035CFE7CE677), UINT64_C(0x03CF050EDDAE0588)}, + boost::int128::uint128_t {UINT64_C(0x000000032A58EE06), UINT64_C(0x156A97CEBE58298E)}, + boost::int128::uint128_t {UINT64_C(0x00000000025E9368), UINT64_C(0xCF9BEE177CFEF916)}, + boost::int128::uint128_t {UINT64_C(0x00000000000171B8), UINT64_C(0xEE9A64E2B76E0BF3)}, + }; + + static constexpr std::uint64_t d128_tail[2] = { UINT64_C(0xBB0CD32F8671EFEA), UINT64_C(0x004F5B056DD15B7F) }; }; #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) template -constexpr std::array sin_table_imp::d32_coeffs; +constexpr fx<1> sin_table_imp::d32[5]; template -constexpr std::array sin_table_imp::d64_coeffs; +constexpr int sin_table_imp::d64_shift; template -constexpr std::array sin_table_imp::d32_fast_coeffs; +constexpr int sin_table_imp::d128_mid_shift; template -constexpr std::array sin_table_imp::d64_fast_coeffs; +constexpr int sin_table_imp::d128_tail_shift; -#endif +template +constexpr fx<2> sin_table_imp::d64_head[2]; -using sin_table = sin_table_imp; +template +constexpr std::uint64_t sin_table_imp::d64_tail[6]; + +template +constexpr fx<3> sin_table_imp::d128_head[3]; -} //namespace sin_detail +template +constexpr boost::int128::uint128_t sin_table_imp::d128_mid[7]; -template -constexpr auto sin_series_expansion(T x) noexcept; +template +constexpr std::uint64_t sin_table_imp::d128_tail[2]; -template <> -constexpr auto sin_series_expansion(decimal32_t x) noexcept -{ - const auto b_neg = signbit(x); - x = abs(x); - auto result = remez_series_result(x, sin_detail::sin_table::d32_coeffs); - return b_neg ? -result : result; -} +#endif + +using sin_table = sin_table_imp; -template <> -constexpr auto sin_series_expansion(decimal_fast32_t x) noexcept +constexpr auto sin_poly(const fx<1>& z) noexcept -> fx<1> { - const auto b_neg = signbit(x); - x = abs(x); - auto result = remez_series_result(x, sin_detail::sin_table::d32_fast_coeffs); - return b_neg ? -result : result; + return fx_alternating(z, sin_table::d32); } -template <> -constexpr auto sin_series_expansion(decimal64_t x) noexcept +constexpr auto sin_poly(const fx<2>& z) noexcept -> fx<2> { - const auto b_neg = signbit(x); - x = abs(x); - auto result = remez_series_result(x, sin_detail::sin_table::d64_coeffs); - return b_neg ? -result : result; + return fx_alternating_split(z, sin_table::d64_head, sin_table::d64_tail); } -template <> -constexpr auto sin_series_expansion(decimal_fast64_t x) noexcept +constexpr auto sin_poly(const fx<3>& z) noexcept -> fx<3> { - const auto b_neg = signbit(x); - x = abs(x); - auto result = remez_series_result(x, sin_detail::sin_table::d64_fast_coeffs); - return b_neg ? -result : result; + return fx_alternating_split(z, sin_table::d128_head, sin_table::d128_mid, sin_table::d128_tail); } -template <> -constexpr auto sin_series_expansion(decimal128_t x) noexcept +// sin(r) / r for z = r^2, which is below 1 for any r != 0. +template +constexpr auto sinc_fx(const fx& z) noexcept -> fx { - const bool b_neg { signbit(x) }; - - x = abs(x); - - // PadeApproximant[Sin[x], {x, 0, {14, 13}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - constexpr decimal128_t c0 { boost::int128::uint128_t { UINT64_C(72470724512963), UINT64_C(12010094287581601792) }, -1 }; - constexpr decimal128_t c1 { boost::int128::uint128_t { UINT64_C(111100426260665), UINT64_C(12001293056709775360) }, -2, construction_sign::negative }; - constexpr decimal128_t c2 { boost::int128::uint128_t { UINT64_C(448976101608303), UINT64_C(8651619847551332352) }, -4 }; - constexpr decimal128_t c3 { boost::int128::uint128_t { UINT64_C(73569920121966), UINT64_C(7922026052315602944) }, -5, construction_sign::negative }; - constexpr decimal128_t c4 { boost::int128::uint128_t { UINT64_C(56791565109495), UINT64_C(18025512837605806080) }, -7 }; - constexpr decimal128_t c5 { boost::int128::uint128_t { UINT64_C(208944907042123), UINT64_C(1905626912845279232) }, -10, construction_sign::negative }; - constexpr decimal128_t c6 { boost::int128::uint128_t { UINT64_C(301324799882787), UINT64_C(8861120840873566208) }, -13 }; - - constexpr decimal128_t d1 { boost::int128::uint128_t { UINT64_C(96841145942737), UINT64_C(12517245955660587008) }, -3 }; - constexpr decimal128_t d2 { boost::int128::uint128_t { UINT64_C(64553072381691), UINT64_C(13718792646062137344) }, -5 }; - constexpr decimal128_t d3 { boost::int128::uint128_t { UINT64_C(279090388104865), UINT64_C(5072548100861788160) }, -8 }; - constexpr decimal128_t d4 { boost::int128::uint128_t { UINT64_C(84086452204639), UINT64_C(9046779044634853376) }, -10 }; - constexpr decimal128_t d5 { boost::int128::uint128_t { UINT64_C(171178955723736), UINT64_C(18053324302671642624) }, -13 }; - constexpr decimal128_t d6 { boost::int128::uint128_t { UINT64_C(189091057352841), UINT64_C(2258222749986258944) }, -16 }; - - const decimal128_t x2 { x * x }; - - const decimal128_t top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) }; - const decimal128_t bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) }; - - const decimal128_t result { top / bot }; - - return b_neg ? -result : result; + return fx_clamp_below_one(fx_sub(fx_one(), fx_mul(z, sin_poly(z)))); } -template <> -constexpr auto sin_series_expansion(decimal_fast128_t x) noexcept +// sin(|r|) with the sign neg, for r != 0. +template +constexpr auto sin_of(const trig_arg::words>& r, const bool neg) noexcept -> T { - const bool b_neg { signbit(x) }; - - x = abs(x); - - // PadeApproximant[Sin[x], {x, 0, {14, 13}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - constexpr decimal_fast128_t c0 { boost::int128::uint128_t { UINT64_C(72470724512963), UINT64_C(12010094287581601792) }, -1 }; - constexpr decimal_fast128_t c1 { boost::int128::uint128_t { UINT64_C(111100426260665), UINT64_C(12001293056709775360) }, -2, construction_sign::negative }; - constexpr decimal_fast128_t c2 { boost::int128::uint128_t { UINT64_C(448976101608303), UINT64_C(8651619847551332352) }, -4 }; - constexpr decimal_fast128_t c3 { boost::int128::uint128_t { UINT64_C(73569920121966), UINT64_C(7922026052315602944) }, -5, construction_sign::negative }; - constexpr decimal_fast128_t c4 { boost::int128::uint128_t { UINT64_C(56791565109495), UINT64_C(18025512837605806080) }, -7 }; - constexpr decimal_fast128_t c5 { boost::int128::uint128_t { UINT64_C(208944907042123), UINT64_C(1905626912845279232) }, -10, construction_sign::negative }; - constexpr decimal_fast128_t c6 { boost::int128::uint128_t { UINT64_C(301324799882787), UINT64_C(8861120840873566208) }, -13 }; - - constexpr decimal_fast128_t d1 { boost::int128::uint128_t { UINT64_C(96841145942737), UINT64_C(12517245955660587008) }, -3 }; - constexpr decimal_fast128_t d2 { boost::int128::uint128_t { UINT64_C(64553072381691), UINT64_C(13718792646062137344) }, -5 }; - constexpr decimal_fast128_t d3 { boost::int128::uint128_t { UINT64_C(279090388104865), UINT64_C(5072548100861788160) }, -8 }; - constexpr decimal_fast128_t d4 { boost::int128::uint128_t { UINT64_C(84086452204639), UINT64_C(9046779044634853376) }, -10 }; - constexpr decimal_fast128_t d5 { boost::int128::uint128_t { UINT64_C(171178955723736), UINT64_C(18053324302671642624) }, -13 }; - constexpr decimal_fast128_t d6 { boost::int128::uint128_t { UINT64_C(189091057352841), UINT64_C(2258222749986258944) }, -16 }; - - const decimal_fast128_t x2 { x * x }; - - const decimal_fast128_t top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) }; - const decimal_fast128_t bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) }; - - const decimal_fast128_t result { top / bot }; - - return b_neg ? -result : result; + const auto rf {fixed_r(r)}; + const auto f {sinc_fx(fx_mul(rf, rf))}; + return trig_round(fx_scale(r.sig, f), -r.k, neg); } +} // namespace trig } // namespace detail } // namespace decimal } // namespace boost -#endif +#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_SIN_IMPL_HPP diff --git a/include/boost/decimal/detail/cmath/impl/tan_impl.hpp b/include/boost/decimal/detail/cmath/impl/tan_impl.hpp new file mode 100644 index 000000000..c2edcb8ab --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/tan_impl.hpp @@ -0,0 +1,98 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#ifndef BOOST_DECIMAL_DETAIL_CMATH_IMPL_TAN_IMPL_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_TAN_IMPL_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { +namespace trig { + +// 1/c for c in [0.5, 1]: the seed 2^124 / (top word of c) is good to 61 bits, and two Newton steps +// give more than 190 bits. +template +constexpr auto fx_reciprocal(const fx& c) noexcept -> fx +{ + fx y {}; + y.w[N - 1] = static_cast(boost::int128::uint128_t {UINT64_C(1) << 60U, UINT64_C(0)} / c.w[N - 1]); + + fx two {}; + two.w[N - 1] = UINT64_C(1) << 63U; + for (int i {}; i < 2; ++i) + { + y = fx_mul(y, fx_sub(two, fx_mul(c, y))); + } + return y; +} + +// tan(r) for an even quadrant, else -cot(r), for r != 0. +template +constexpr auto tan_of(const trig_arg::words>& r) noexcept -> T +{ + constexpr int words {trig_traits::words}; + const auto rf {fixed_r(r)}; + const auto z {fx_mul(rf, rf)}; + const auto s {sinc_fx(z)}; + const auto c {cos_fx(z)}; + + if ((r.n & 1U) == 0U) + { + // tan(r) = r * (sin(r) / r) / cos(r), and the factor is above 1 for any r != 0. + auto m {fx_mul(s, fx_reciprocal(c))}; + if (m.w[words - 1] < (UINT64_C(1) << 62U) || fx_is_one(m)) + { + m = fx_one(); + m.w[0] += 1U; + } + return trig_round(fx_scale(r.sig, m), -r.k, r.neg); + } + + // -cot(r) = -(cos(r) / (sin(r) / r)) / r, and cos(r) / (sin(r) / r) is below 1 for any r != 0. + const auto m {fx_clamp_below_one(fx_mul(c, fx_reciprocal(s)))}; + + // floor(m * 10^75 / 2^(64N - 2)), then divide by sig: for m >= 0.7, the quotient has 37 or 38 digits. + constexpr std::uint64_t pow10_75[4] {UINT64_C(0), UINT64_C(10084168908774762496), + UINT64_C(12965995782233477362), UINT64_C(159309191113245227)}; + std::uint64_t p[static_cast(words + 4)] {}; + for (int i {}; i < words; ++i) + { + std::uint64_t carry {}; + for (int j {}; j < 4; ++j) + { + const auto t {static_cast(m.w[i]) * pow10_75[j] + p[i + j] + carry}; + p[i + j] = t.low; + carry = t.high; + } + p[i + 4] = carry; + } + std::uint64_t num[4] {}; + for (int i {}; i < 4; ++i) + { + num[i] = (p[words - 1 + i] >> 62U) | (p[words + i] << 2U); + } + const auto dm {impl::div_mod(u256 {num[3], num[2], num[1], num[0]}, r.sig)}; + return trig_round(boost::int128::uint128_t {dm.quotient.bytes[1], dm.quotient.bytes[0]}, r.k - 75, !r.neg); +} + +} // namespace trig +} // namespace detail +} // namespace decimal +} // namespace boost + +#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_TAN_IMPL_HPP diff --git a/include/boost/decimal/detail/cmath/impl/trig_fixed_point.hpp b/include/boost/decimal/detail/cmath/impl/trig_fixed_point.hpp new file mode 100644 index 000000000..9dd8e66e8 --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/trig_fixed_point.hpp @@ -0,0 +1,315 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#ifndef BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_FIXED_POINT_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_FIXED_POINT_HPP + +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { +namespace trig { + +#ifdef _MSC_VER +# pragma warning(push) +# pragma warning(disable : 4127) // Conditional expression is constant +#endif + +// Fixed point with N words of 64 bits, low word first: the value is the integer divided by +// 2^(64N - 2), in the range [0, 4). +template +struct fx +{ + std::uint64_t w[static_cast(N)]; +}; + +template +struct fixed_table_imp +{ + // floor(2^382 / 10^(46 + 9a)) for a = 0 to 7. Its top N + 1 words are floor(2^(64N + 190) / 10^(46 + 9a)). + static constexpr fx<4> pow10_neg[8] = + { + {{UINT64_C(0x51BA47FEE249CE5F), UINT64_C(0x329C3A0F675D7764), UINT64_C(0xAAC1C372ACE584C1), UINT64_C(0x00000024899C4858)}}, + {{UINT64_C(0x9CF16DDC1CC486D3), UINT64_C(0x464DD69685606BAB), UINT64_C(0xED737BB6C4183D55), UINT64_C(0x000000000000009C)}}, + {{UINT64_C(0xCD8890A87FB90364), UINT64_C(0xE7A694FC8E635D1E), UINT64_C(0x000002A1FFA89E94), UINT64_C(0x0000000000000000)}}, + {{UINT64_C(0x8AEB6360B1AF3451), UINT64_C(0xCD5F01A4AA8281E3), UINT64_C(0x0000000000000B4E), UINT64_C(0x0000000000000000)}}, + {{UINT64_C(0xC086FEFA16EB73A6), UINT64_C(0x0000309114B688A6), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000)}}, + {{UINT64_C(0xAD07A71F26B27E20), UINT64_C(0x000000000000D097), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000)}}, + {{UINT64_C(0x00037FE5DC91C0A5), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000)}}, + {{UINT64_C(0x00000000000F07DA), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000000)}}, + }; +}; + +#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) + +template +constexpr fx<4> fixed_table_imp::pow10_neg[8]; + +#endif + +using fixed_table = fixed_table_imp; + +template +constexpr auto fx_one() noexcept -> fx +{ + fx r {}; + r.w[N - 1] = UINT64_C(1) << 62U; + return r; +} + +template +constexpr auto fx_is_one(const fx& a) noexcept -> bool +{ + for (int i {}; i < N - 1; ++i) + { + if (a.w[i] != 0U) + { + return false; + } + } + return a.w[N - 1] == (UINT64_C(1) << 62U); +} + +// a * b, truncated. +template +constexpr auto fx_mul(const fx& a, const fx& b) noexcept -> fx +{ + std::uint64_t p[static_cast(2 * N)] {}; + for (int i {}; i < N; ++i) + { + std::uint64_t c {}; + for (int j {}; j < N; ++j) + { + const auto t {static_cast(a.w[i]) * b.w[j] + p[i + j] + c}; + p[i + j] = t.low; + c = t.high; + } + p[i + N] = c; + } + fx r {}; + for (int i {}; i < N; ++i) + { + r.w[i] = (p[N - 1 + i] >> 62U) | (p[N + i] << 2U); + } + return r; +} + +// a - b, for a >= b. +template +constexpr auto fx_sub(const fx& a, const fx& b) noexcept -> fx +{ + fx r {}; + std::uint64_t borrow {}; + for (int i {}; i < N; ++i) + { + const std::uint64_t t {a.w[i] - b.w[i] - borrow}; + borrow = (a.w[i] < b.w[i] || (a.w[i] == b.w[i] && borrow != 0U)) ? 1U : 0U; + r.w[i] = t; + } + return r; +} + +// min(a, 1 - 2^-(64N - 2)): a value that must be below one stays below one. +template +constexpr auto fx_clamp_below_one(const fx& a) noexcept -> fx +{ + if (a.w[N - 1] < (UINT64_C(1) << 62U)) + { + return a; + } + fx r {}; + for (int i {}; i < N - 1; ++i) + { + r.w[i] = ~UINT64_C(0); + } + r.w[N - 1] = (UINT64_C(1) << 62U) - 1U; + return r; +} + +// |c0| - z*(|c1| - z*(|c2| - ...)): the coefficients of the sin and cos kernels alternate in sign. +template +constexpr auto fx_alternating(const fx& z, const fx (&c)[M]) noexcept -> fx +{ + fx s {c[M - 1]}; + for (std::size_t i {M - 1}; i-- > 0U;) + { + s = fx_sub(c[i], fx_mul(z, s)); + } + return s; +} + +// The same, with the terms from z^M on in one word scaled by 2^(62 + Shift): their weight is small +// enough that 64 bits of them keep the error of the sum below 2^-(64N - 2). +template +constexpr auto fx_alternating_split(const fx& z, const fx (&head)[M], const std::uint64_t (&tail)[K]) noexcept -> fx +{ + static_assert(0 <= Shift && Shift < 64 * (N - 1), "the tail must fit in the N words"); + + const std::uint64_t z1 {z.w[N - 1]}; + std::uint64_t t {tail[K - 1]}; + for (std::size_t i {K - 1}; i-- > 0U;) + { + t = tail[i] - ((static_cast(z1) * t).high << 2U); + } + + // t * 2^-(62 + Shift) in N words + fx s {}; + constexpr int pos {64 * (N - 1) - Shift}; + constexpr int word {pos / 64}; + constexpr int bit {pos % 64}; + s.w[word] = t << bit; + if (bit != 0 && word + 1 < N) + { + s.w[word + 1] = t >> (64 - bit); + } + for (std::size_t i {M}; i-- > 0U;) + { + s = fx_sub(head[i], fx_mul(z, s)); + } + return s; +} + +// The high 128 bits of a * b. +constexpr auto mul_high(const boost::int128::uint128_t& a, const boost::int128::uint128_t& b) noexcept -> boost::int128::uint128_t +{ + const auto ll {static_cast(a.low) * b.low}; + const auto lh {static_cast(a.low) * b.high}; + const auto hl {static_cast(a.high) * b.low}; + const auto hh {static_cast(a.high) * b.high}; + const auto mid {static_cast(ll.high) + lh.low + hl.low}; + return hh + lh.high + hl.high + mid.high; +} + +// For three words: the terms from z^M on in 128 bits scaled by 2^(126 + MidShift), and the terms +// after those in one word scaled by 2^(62 + TailShift). +template +constexpr auto fx_alternating_split(const fx<3>& z, const fx<3> (&head)[M], const boost::int128::uint128_t (&mid)[K], + const std::uint64_t (&tail)[L]) noexcept -> fx<3> +{ + static_assert(0 < MidShift && MidShift < 64, "the middle terms must fit in the three words"); + static_assert(0 <= 64 + MidShift - TailShift && 64 + MidShift - TailShift < 64, "the tail must fit in 128 bits"); + + const std::uint64_t z1 {z.w[2]}; + std::uint64_t t {tail[L - 1]}; + for (std::size_t i {L - 1}; i-- > 0U;) + { + t = tail[i] - ((static_cast(z1) * t).high << 2U); + } + + const boost::int128::uint128_t z2 {z.w[2], z.w[1]}; + auto u {static_cast(t) << (64 + MidShift - TailShift)}; + for (std::size_t i {K}; i-- > 0U;) + { + u = mid[i] - (mul_high(z2, u) << 2U); + } + + // u * 2^-(126 + MidShift) in three words + constexpr int pos {64 - MidShift}; + fx<3> s {}; + s.w[0] = u.low << pos; + s.w[1] = (u.low >> (64 - pos)) | (u.high << pos); + s.w[2] = u.high >> (64 - pos); + for (std::size_t i {M}; i-- > 0U;) + { + s = fx_sub(head[i], fx_mul(z, s)); + } + return s; +} + +// sig * 10^-k as fixed point, for 38 <= k < 110, else 0. With k = 38 + 9a + b, sig * 10^(8 - b) is exact in +// three words; times 2^(64N + 190) / 10^(46 + 9a) it is sig * 10^-k * 2^(64N + 190): drop 192 bits. +template +constexpr auto fx_from(const boost::int128::uint128_t& sig, const int k) noexcept -> fx +{ + constexpr int rows {static_cast(sizeof(fixed_table::pow10_neg) / sizeof(fixed_table::pow10_neg[0]))}; + + fx r {}; + if (k < 38 || k >= 38 + 9 * rows) + { + return r; + } + const auto f {pow10(static_cast(8 - (k - 38) % 9))}; + const auto lo {static_cast(sig.low) * f}; + const auto hi {static_cast(sig.high) * f + lo.high}; + const std::uint64_t g[3] {lo.low, hi.low, hi.high}; + const auto& q {fixed_table::pow10_neg[(k - 38) / 9]}; + + std::uint64_t p[static_cast(N + 4)] {}; + for (int i {}; i < 3; ++i) + { + std::uint64_t c {}; + for (int j {}; j < N + 1; ++j) + { + const auto t {static_cast(g[i]) * q.w[j + 3 - N] + p[i + j] + c}; + p[i + j] = t.low; + c = t.high; + } + p[i + N + 1] = c; + } + for (int i {}; i < N; ++i) + { + r.w[i] = p[i + 3]; + } + return r; +} + +// floor(sig * f / 2^(64N - 2)), for f <= 1.5 so that the result has at most 39 digits. +template +constexpr auto fx_scale(const boost::int128::uint128_t& sig, const fx& f) noexcept -> boost::int128::uint128_t +{ + std::uint64_t p[static_cast(N + 2)] {}; + const std::uint64_t g[2] {sig.low, sig.high}; + for (int i {}; i < 2; ++i) + { + std::uint64_t c {}; + for (int j {}; j < N; ++j) + { + const auto t {static_cast(g[i]) * f.w[j] + p[i + j] + c}; + p[i + j] = t.low; + c = t.high; + } + p[i + N] = c; + } + return boost::int128::uint128_t {(p[N] >> 62U) | (p[N + 1] << 2U), (p[N - 1] >> 62U) | (p[N] << 2U)}; +} + +// sig * 10^e rounded once in the current mode, for a sig of 37 or 38 digits which truncates a value that +// is never exact. The value is in (sig, sig + 1), thus if the last digit is 0 or 5, sig + 1 rounds the same. +template +constexpr auto trig_round(boost::int128::uint128_t sig, int e, const bool neg) noexcept -> T +{ + if (sig < pow10(static_cast(37))) + { + sig *= 10U; + --e; + } + + // The constructor sees a 0 as exact and, where it drops one digit only, a 5 as a tie. + // 2^64 = 1 (mod 5), thus sig = high + low (mod 5). + if ((sig.high % 5U + sig.low % 5U) % 5U == 0U) + { + sig += 1U; + } + return T {sig, e, neg}; +} + +#ifdef _MSC_VER +# pragma warning(pop) +#endif + +} // namespace trig +} // namespace detail +} // namespace decimal +} // namespace boost + +#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_FIXED_POINT_HPP diff --git a/include/boost/decimal/detail/cmath/impl/trig_reduce.hpp b/include/boost/decimal/detail/cmath/impl/trig_reduce.hpp new file mode 100644 index 000000000..5455f31b2 --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/trig_reduce.hpp @@ -0,0 +1,682 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#ifndef BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_REDUCE_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_REDUCE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { +namespace trig { + +template +struct trig_table_imp +{ + // 2/pi = 0.636619772 367581343 ...: 6480 digits in words of 9 digits, which cover the window of + // the largest decimal128_t exponent. + static constexpr std::uint32_t two_over_pi[720] = + { + UINT32_C(636619772), UINT32_C(367581343), UINT32_C( 75535053), UINT32_C(490057448), UINT32_C(137838582), UINT32_C(961825794), UINT32_C(990669376), UINT32_C(235587190), + UINT32_C(536906140), UINT32_C(360455211), UINT32_C( 65012343), UINT32_C(824291370), UINT32_C(907031832), UINT32_C(147571647), UINT32_C(384458314), UINT32_C(611511869), + UINT32_C(642926799), UINT32_C(356916959), UINT32_C(867749636), UINT32_C(310292310), UINT32_C(985587701), UINT32_C(230754869), UINT32_C(571584869), UINT32_C(590646773), + UINT32_C(449560966), UINT32_C(894516047), UINT32_C(329520456), UINT32_C(890799022), UINT32_C(863761847), UINT32_C(560347610), UINT32_C(695824481), UINT32_C(957643747), + UINT32_C(751376342), UINT32_C(114892399), UINT32_C(785773600), UINT32_C(994689390), UINT32_C(957838443), UINT32_C(593292387), UINT32_C(132299624), UINT32_C(667945851), + UINT32_C(218797794), UINT32_C(608751526), UINT32_C(299146267), UINT32_C(856964155), UINT32_C(983496557), UINT32_C(394439935), UINT32_C(472396799), UINT32_C(849771502), + UINT32_C(340684715), UINT32_C(433724470), UINT32_C( 75068642), UINT32_C(186190147), UINT32_C(952038957), UINT32_C(841459037), UINT32_C(335072237), UINT32_C(209977986), + UINT32_C(541221308), UINT32_C(627102012), UINT32_C(881299111), UINT32_C(265588664), UINT32_C( 91786992), UINT32_C(478392663), UINT32_C(362424067), UINT32_C(212143992), + UINT32_C(535647949), UINT32_C(995331146), UINT32_C(617741119), UINT32_C( 20280064), UINT32_C(962710257), UINT32_C(555398285), UINT32_C(243520488), UINT32_C(797504590), + UINT32_C(725511058), UINT32_C(951562532), UINT32_C(272185831), UINT32_C(913927045), UINT32_C(249709256), UINT32_C(279843100), UINT32_C( 98001191), UINT32_C( 39428356), + UINT32_C(227611187), UINT32_C(140526100), UINT32_C(840065270), UINT32_C(984083699), UINT32_C(246424962), UINT32_C(245824812), UINT32_C(585936356), UINT32_C(993836765), + UINT32_C(740846301), UINT32_C(630224803), UINT32_C(486106427), UINT32_C(208868636), UINT32_C(563029898), UINT32_C(330890390), UINT32_C(985141599), UINT32_C(500621317), + UINT32_C(563255927), UINT32_C( 89637433), UINT32_C( 19188293), UINT32_C(314876162), UINT32_C(799903630), UINT32_C(630831397), UINT32_C(388157435), UINT32_C(931234869), + UINT32_C(370256146), UINT32_C(758046650), UINT32_C(182823773), UINT32_C(310525074), UINT32_C(600104490), UINT32_C(871884612), UINT32_C(845039801), UINT32_C(754671780), + UINT32_C(150502243), UINT32_C(345268467), UINT32_C(810390325), UINT32_C(128997664), UINT32_C(933372580), UINT32_C(424494147), UINT32_C(514252454), UINT32_C(546768668), + UINT32_C(568278987), UINT32_C(840517002), UINT32_C(313344212), UINT32_C(478434378), UINT32_C( 39358226), UINT32_C(874839818), UINT32_C(986041726), UINT32_C(495262070), + UINT32_C(323357771), UINT32_C(919883998), UINT32_C( 21017550), UINT32_C(264517783), UINT32_C(533227384), UINT32_C(203141166), UINT32_C( 60564161), UINT32_C(957195402), + UINT32_C(555264310), UINT32_C(478797229), UINT32_C(364155998), UINT32_C(314767562), UINT32_C(392374951), UINT32_C( 88247501), UINT32_C(728908757), UINT32_C(205465021), + UINT32_C( 44955121), UINT32_C(550155524), UINT32_C(427256270), UINT32_C(617363313), UINT32_C(114107733), UINT32_C(707198224), UINT32_C(283161544), UINT32_C(241410955), + UINT32_C(984980503), UINT32_C(982997105), UINT32_C(188094376), UINT32_C(382337204), UINT32_C(659318564), UINT32_C(742310849), UINT32_C(623017797), UINT32_C(828087159), + UINT32_C( 79169637), UINT32_C(961309179), UINT32_C( 80866598), UINT32_C(414261272), UINT32_C(614176015), UINT32_C(362759498), UINT32_C(870766355), UINT32_C( 52763866), + UINT32_C( 27857619), UINT32_C(107882750), UINT32_C(734627112), UINT32_C(419119181), UINT32_C(801413583), UINT32_C( 33207527), UINT32_C(354751751), UINT32_C( 64499259), + UINT32_C(812239862), UINT32_C(320876334), UINT32_C(395004140), UINT32_C(508516172), UINT32_C(926321994), UINT32_C(878747511), UINT32_C( 37862653), UINT32_C(848841368), + UINT32_C(177634219), UINT32_C(914015170), UINT32_C(954777174), UINT32_C(146477511), UINT32_C(317149437), UINT32_C(513738812), UINT32_C(920948583), UINT32_C(351694228), + UINT32_C(474545367), UINT32_C(717840732), UINT32_C(729167856), UINT32_C(660035132), UINT32_C(317325413), UINT32_C(991163989), UINT32_C(834597161), UINT32_C( 69802439), + UINT32_C(574756378), UINT32_C(353220134), UINT32_C(812215221), UINT32_C(892492863), UINT32_C(237727907), UINT32_C( 41291325), UINT32_C(256759238), UINT32_C(999289753), + UINT32_C(340697427), UINT32_C(959390004), UINT32_C(158002735), UINT32_C(520159146), UINT32_C(894398432), UINT32_C( 96010956), UINT32_C( 43499819), UINT32_C(419151694), + UINT32_C(273044559), UINT32_C(795613075), UINT32_C(989708333), UINT32_C(984459683), UINT32_C(315615107), UINT32_C(138972142), UINT32_C( 18273824), UINT32_C(334685917), + UINT32_C(233826893), UINT32_C(308141941), UINT32_C(570224808), UINT32_C(347357296), UINT32_C(398248847), UINT32_C( 13273576), UINT32_C( 83883174), UINT32_C(283099861), + UINT32_C(995234744), UINT32_C(265443874), UINT32_C(647868149), UINT32_C(898168411), UINT32_C(324877007), UINT32_C(384899339), UINT32_C(964644598), UINT32_C(266224151), + UINT32_C(878704559), UINT32_C(725131984), UINT32_C(310433111), UINT32_C(960403132), UINT32_C(144009353), UINT32_C( 91951634), UINT32_C(160955046), UINT32_C(229781723), + UINT32_C(704047640), UINT32_C(217351993), UINT32_C(556186196), UINT32_C(849931806), UINT32_C(428291412), UINT32_C( 20908840), UINT32_C(944070093), UINT32_C(252692719), + UINT32_C( 37244201), UINT32_C(312620437), UINT32_C(495654558), UINT32_C(581223170), UINT32_C(428720334), UINT32_C(471819506), UINT32_C(898583921), UINT32_C(895909169), + UINT32_C(792436803), UINT32_C(748503147), UINT32_C(673331583), UINT32_C(545135961), UINT32_C(743474666), UINT32_C(559026937), UINT32_C(805638014), UINT32_C(549308766), + UINT32_C(972455522), UINT32_C(655322903), UINT32_C(692110389), UINT32_C(380242192), UINT32_C(851112148), UINT32_C(261351132), UINT32_C(128683950), UINT32_C(939866273), + UINT32_C(963201307), UINT32_C(954026967), UINT32_C(165858734), UINT32_C( 33126467), UINT32_C(413257344), UINT32_C(642923980), UINT32_C(599412479), UINT32_C(278935033), + UINT32_C(776839366), UINT32_C(623816609), UINT32_C( 2573577), UINT32_C(251457761), UINT32_C(535534246), UINT32_C( 35190865), UINT32_C(800682588), UINT32_C(270075098), + UINT32_C(242366434), UINT32_C(867431431), UINT32_C(756904939), UINT32_C( 25326844), UINT32_C(531994623), UINT32_C(766387562), UINT32_C(879402754), UINT32_C(976920230), + UINT32_C( 76790822), UINT32_C(760152873), UINT32_C(570248813), UINT32_C(549694145), UINT32_C( 27233416), UINT32_C(626069188), UINT32_C(435246887), UINT32_C(183747330), + UINT32_C(259540749), UINT32_C(998994834), UINT32_C(212466393), UINT32_C(224405568), UINT32_C(578178406), UINT32_C(459538110), UINT32_C(810045644), UINT32_C(280994086), + UINT32_C(958980415), UINT32_C(466945615), UINT32_C(491440398), UINT32_C(699572694), UINT32_C(247248284), UINT32_C(696191559), UINT32_C(747554622), UINT32_C(769231394), + UINT32_C( 9222822), UINT32_C(857625455), UINT32_C(452809474), UINT32_C( 80429640), UINT32_C(229943691), UINT32_C(244628878), UINT32_C(720159129), UINT32_C(903812006), + UINT32_C(678340884), UINT32_C(921385675), UINT32_C( 94601741), UINT32_C(870585826), UINT32_C(263887604), UINT32_C(492339068), UINT32_C(397238834), UINT32_C(365134586), + UINT32_C(676767107), UINT32_C(755165733), UINT32_C(262266026), UINT32_C(792528656), UINT32_C(608403582), UINT32_C(846914495), UINT32_C(370428271), UINT32_C(380704044), + UINT32_C(538032027), UINT32_C(979073689), UINT32_C(427958499), UINT32_C(522063103), UINT32_C(923813588), UINT32_C(323419002), UINT32_C(390145062), UINT32_C(596137577), + UINT32_C(816823271), UINT32_C(545742732), UINT32_C(168001260), UINT32_C(382378973), UINT32_C(757010179), UINT32_C(402699657), UINT32_C(163459005), UINT32_C(769213285), + UINT32_C(329827804), UINT32_C(653978271), UINT32_C( 15757696), UINT32_C(144362175), UINT32_C(334211316), UINT32_C(973688139), UINT32_C(793746460), UINT32_C(586529144), + UINT32_C( 99106666), UINT32_C(419812562), UINT32_C(629374302), UINT32_C(120563633), UINT32_C(119523659), UINT32_C(146773739), UINT32_C(690950410), UINT32_C(539991319), + UINT32_C(828072647), UINT32_C(857284932), UINT32_C(561903051), UINT32_C(589936331), UINT32_C(564696389), UINT32_C(913055159), UINT32_C(672679975), UINT32_C(794999086), + UINT32_C( 79592749), UINT32_C( 66517840), UINT32_C(732215833), UINT32_C(310083694), UINT32_C(540274155), UINT32_C(569138729), UINT32_C(890398901), UINT32_C(132030674), + UINT32_C(277503346), UINT32_C(388916792), UINT32_C(977189896), UINT32_C(246552732), UINT32_C(455833226), UINT32_C(977394067), UINT32_C(714389532), UINT32_C(949570649), + UINT32_C(609738007), UINT32_C(991239761), UINT32_C(608758453), UINT32_C(933709445), UINT32_C(470579965), UINT32_C(530861666), UINT32_C(425369931), UINT32_C(745496740), + UINT32_C(244904434), UINT32_C(452847994), UINT32_C(533851388), UINT32_C(397673597), UINT32_C(709718236), UINT32_C(625133359), UINT32_C(619215284), UINT32_C(700046448), + UINT32_C(466688207), UINT32_C(650317214), UINT32_C(211716964), UINT32_C(537612464), UINT32_C(536449981), UINT32_C(273543707), UINT32_C(833961775), UINT32_C(387231396), + UINT32_C(389593123), UINT32_C(542118818), UINT32_C( 61221596), UINT32_C(560395479), UINT32_C(536353461), UINT32_C(934660889), UINT32_C(867449634), UINT32_C(901605616), + UINT32_C( 36471496), UINT32_C(848818092), UINT32_C(301338958), UINT32_C(901525976), UINT32_C(155367623), UINT32_C(473692463), UINT32_C(785290977), UINT32_C(356264500), + UINT32_C(649572425), UINT32_C(132781295), UINT32_C(533568526), UINT32_C(138225526), UINT32_C( 47008140), UINT32_C(434983823), UINT32_C(280449501), UINT32_C(743907262), + UINT32_C(136074962), UINT32_C(957736145), UINT32_C(359121552), UINT32_C(688401812), UINT32_C(676731807), UINT32_C(795183670), UINT32_C(695816711), UINT32_C(516974110), + UINT32_C(469628984), UINT32_C(237566410), UINT32_C(929131517), UINT32_C(872774596), UINT32_C(515798859), UINT32_C(813730210), UINT32_C(894366637), UINT32_C(192289919), + UINT32_C(943224507), UINT32_C(602932875), UINT32_C(378107177), UINT32_C(340182320), UINT32_C(780997026), UINT32_C(522481950), UINT32_C(646453746), UINT32_C(135968115), + UINT32_C( 18083422), UINT32_C(137657639), UINT32_C(620519309), UINT32_C( 98186364), UINT32_C(725288931), UINT32_C(362046664), UINT32_C(626028393), UINT32_C(502297349), + UINT32_C(181945248), UINT32_C(164486865), UINT32_C(523662424), UINT32_C(644662928), UINT32_C( 333224), UINT32_C(458424725), UINT32_C(121305034), UINT32_C(783806409), + UINT32_C(852866455), UINT32_C(430645921), UINT32_C(887973083), UINT32_C(108526576), UINT32_C(480637984), UINT32_C( 44253132), UINT32_C(208303833), UINT32_C(394012203), + UINT32_C(163823399), UINT32_C(319287469), UINT32_C(611593542), UINT32_C( 55329582), UINT32_C(808323055), UINT32_C(902017169), UINT32_C( 39390588), UINT32_C(284065707), + UINT32_C(897538017), UINT32_C(236663458), UINT32_C(113441299), UINT32_C(734417418), UINT32_C(628950231), UINT32_C(664546529), UINT32_C(648183123), UINT32_C(987886265), + UINT32_C(360886352), UINT32_C(218317725), UINT32_C(313112022), UINT32_C( 98452835), UINT32_C(560749684), UINT32_C(843697956), UINT32_C(416402086), UINT32_C(198723884), + UINT32_C(548830160), UINT32_C(228438536), UINT32_C(265725429), UINT32_C(817596639), UINT32_C( 77743155), UINT32_C(683173702), UINT32_C(471132088), UINT32_C(948045945), + UINT32_C(699700956), UINT32_C(994914852), UINT32_C(528087066), UINT32_C(944302658), UINT32_C(239309043), UINT32_C(829662640), UINT32_C(937514974), UINT32_C(516528438), + UINT32_C(994358860), UINT32_C(285229564), UINT32_C(162905741), UINT32_C(656718822), UINT32_C(889061919), UINT32_C(215260510), UINT32_C(383164960), UINT32_C(101378721), + UINT32_C(928810469), UINT32_C(369196004), UINT32_C( 81932249), UINT32_C(852135185), UINT32_C(898712762), UINT32_C( 7247321), UINT32_C(500615211), UINT32_C(518093733), + UINT32_C(678200854), UINT32_C(275908365), UINT32_C(162245727), UINT32_C(151516834), UINT32_C(482297999), UINT32_C(703159027), UINT32_C(607396841), UINT32_C(296825885), + UINT32_C(540764555), UINT32_C(259025608), UINT32_C(390422195), UINT32_C(831751405), UINT32_C(656165812), UINT32_C(206063358), UINT32_C(571293061), UINT32_C(624082413), + UINT32_C(247566346), UINT32_C(281088345), UINT32_C( 1079665), UINT32_C(575006111), UINT32_C(549442432), UINT32_C(458227793), UINT32_C(684128963), UINT32_C(109090968), + UINT32_C(660545693), UINT32_C(746797086), UINT32_C(536123762), UINT32_C(122992261), UINT32_C( 74037206), UINT32_C(635685476), UINT32_C(856572517), UINT32_C(485364246), + UINT32_C(286148562), UINT32_C(481591390), UINT32_C(473706011), UINT32_C(912314425), UINT32_C( 67879843), UINT32_C(236736893), UINT32_C(905340190), UINT32_C(986876069), + UINT32_C(801805784), UINT32_C(665531384), UINT32_C(832963469), UINT32_C(438040948), UINT32_C(521161777), UINT32_C(511763414), UINT32_C( 13781770), UINT32_C(533652250), + UINT32_C(522983805), UINT32_C(532124091), UINT32_C(725877378), UINT32_C(673314070), UINT32_C(653129660), UINT32_C(608407176), UINT32_C(905775828), UINT32_C(724868680), + UINT32_C(870259687), UINT32_C(857797586), UINT32_C(128888750), UINT32_C(633952978), UINT32_C( 47637605), UINT32_C(362017728), UINT32_C(559434514), UINT32_C(484332717), + UINT32_C(575843377), UINT32_C(559207659), UINT32_C(149559089), UINT32_C(324114524), UINT32_C( 52594782), UINT32_C( 85048207), UINT32_C(311225397), UINT32_C(828474651), + UINT32_C(113026395), UINT32_C(324021406), UINT32_C(209266639), UINT32_C(375763608), UINT32_C(872252578), UINT32_C(180848519), UINT32_C(158937885), UINT32_C(954965033), + UINT32_C( 72895440), UINT32_C(944108439), UINT32_C(924766082), UINT32_C(275293889), UINT32_C(593432053), UINT32_C(464273514), UINT32_C(531547171), UINT32_C(447892946), + UINT32_C(901442674), UINT32_C( 86742528), UINT32_C( 47795912), UINT32_C(293583367), UINT32_C(676266383), UINT32_C(354714117), UINT32_C(649674872), UINT32_C(869119500), + UINT32_C(244157842), UINT32_C(592783429), UINT32_C(824802435), UINT32_C(684913665), UINT32_C(577495386), UINT32_C(198359728), UINT32_C(113924945), UINT32_C(733864478), + UINT32_C(829297238), UINT32_C(183436293), UINT32_C(447514516), UINT32_C(252740066), UINT32_C( 42507030), UINT32_C(740486543), UINT32_C( 35478522), UINT32_C(980799688), + UINT32_C( 4310670), UINT32_C(732378792), UINT32_C(599024907), UINT32_C(297391746), UINT32_C(852433648), UINT32_C(408780835), UINT32_C(979276497), UINT32_C(761950046), + UINT32_C(842367376), UINT32_C(559631557), UINT32_C(823100738), UINT32_C(486476166), UINT32_C(123738175), UINT32_C(211235754), UINT32_C(512292950), UINT32_C(314461071), + UINT32_C(188457329), UINT32_C(296787943), UINT32_C(122255052), UINT32_C( 72353754), UINT32_C(656242870), UINT32_C(147328545), UINT32_C( 51868489), UINT32_C(704377141), + UINT32_C(604438528), UINT32_C(730510604), UINT32_C(804680902), UINT32_C(117171586), UINT32_C(223784328), UINT32_C(197536362), UINT32_C(763042768), UINT32_C( 15818584), + UINT32_C(766560086), UINT32_C(269344071), UINT32_C(638527491), UINT32_C(567994537), UINT32_C(364347612), UINT32_C(802318654), UINT32_C(841251444), UINT32_C(942795527), + UINT32_C( 56145701), UINT32_C( 16334839), UINT32_C(243259340), UINT32_C(761248527), UINT32_C(449889127), UINT32_C(242033804), UINT32_C(947607625), UINT32_C(865289437), + }; + + // pi/2 = 1.570796326 794896619 ... + static constexpr std::uint32_t pio2[6] = { UINT32_C(1), UINT32_C(570796326), UINT32_C(794896619), UINT32_C(231321691), UINT32_C(639751442), UINT32_C(98584699) }; + + // For |x| < 10^19 in base 2^64, high word first: 2/pi and ceil(10^(-9a) * 2^448) for a = 1 to 4, + // each word i with the weight 2^(-64(i + 1)), and pi/2 with word 0 as the integer part. + static constexpr std::uint64_t two_over_pi_bin[7] = + { + UINT64_C(0xA2F9836E4E441529), UINT64_C(0xFC2757D1F534DDC0), UINT64_C(0xDB6295993C439041), UINT64_C(0xFE5163ABDEBBC561), + UINT64_C(0xB7246E3A424DD2E0), UINT64_C(0x06492EEA09D1921C), UINT64_C(0xFE1DEB1CB129A73E) + }; + + static constexpr std::uint64_t pow10_neg9_bin[4][7] = + { + {UINT64_C(0x000000044B82FA09), UINT64_C(0xB5A52CB98B405447), UINT64_C(0xC4A98187EEBB22F0), UINT64_C(0x08D5D64F9C394AE9), + UINT64_C(0x213015356022EF32), UINT64_C(0x164179B6BF082CE3), UINT64_C(0xFD84BF5BB9D3E58A)}, + {UINT64_C(0x0000000000000012), UINT64_C(0x725DD1D243ABA0E7), UINT64_C(0x5FE645CC4873F9E6), UINT64_C(0x5AFE688C928E1F21), + UINT64_C(0x95818AE77F3C36A0), UINT64_C(0x8CCE4E0A36628033), UINT64_C(0xA40BE73647459D42)}, + {UINT64_C(0x0000000000000000), UINT64_C(0x0000004F3A68DBC8), UINT64_C(0xF03F243BAF513267), UINT64_C(0xAA9A3EE524F8E028), + UINT64_C(0x9064E3CFFA15AB8B), UINT64_C(0xB9CCC2933B76B4FA), UINT64_C(0x41402348EBC5909A)}, + {UINT64_C(0x0000000000000000), UINT64_C(0x0000000000000154), UINT64_C(0x484932D2E725A5BB), UINT64_C(0xCA17A3ABA173D3D5), + UINT64_C(0xFC130C23B7AA2DA1), UINT64_C(0x9B9A3CAB811D56FA), UINT64_C(0x9C85A535DF608EEE)}, + }; + + static constexpr std::uint64_t pio2_bin[5] = + { + UINT64_C(0x0000000000000001), UINT64_C(0x921FB54442D18469), UINT64_C(0x898CC51701B839A2), UINT64_C(0x52049C1114CF98E8), + UINT64_C(0x04177D4C76273644) + }; +}; + +#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) + +template +constexpr std::uint32_t trig_table_imp::two_over_pi[720]; + +template +constexpr std::uint32_t trig_table_imp::pio2[6]; + +template +constexpr std::uint64_t trig_table_imp::two_over_pi_bin[7]; + +template +constexpr std::uint64_t trig_table_imp::pow10_neg9_bin[4][7]; + +template +constexpr std::uint64_t trig_table_imp::pio2_bin[5]; + +#endif + +using trig_table = trig_table_imp; + +// Words of 2/pi in the product window, words of the fraction that go into r, and fixed-point words. +// The window keeps 9 * (window - 1 - frac) = 27, 45 and 81 digits for the leading zeros of the +// fraction: the worst cases have 9, 19 and 37. +// bin_frac: 64-bit fraction words for |x| < 10^19, where the fraction has at most 27, 58 and 117 leading zero bits. +template +struct trig_traits; + +template <> +struct trig_traits { static constexpr int window = 6; static constexpr int frac = 2; static constexpr int words = 1; static constexpr int bin_frac = 2; }; + +template <> +struct trig_traits { static constexpr int window = 6; static constexpr int frac = 2; static constexpr int words = 1; static constexpr int bin_frac = 2; }; + +template <> +struct trig_traits { static constexpr int window = 9; static constexpr int frac = 3; static constexpr int words = 2; static constexpr int bin_frac = 3; }; + +template <> +struct trig_traits { static constexpr int window = 9; static constexpr int frac = 3; static constexpr int words = 2; static constexpr int bin_frac = 3; }; + +template <> +struct trig_traits { static constexpr int window = 15; static constexpr int frac = 5; static constexpr int words = 3; static constexpr int bin_frac = 5; }; + +template <> +struct trig_traits { static constexpr int window = 15; static constexpr int frac = 5; static constexpr int words = 3; static constexpr int bin_frac = 5; }; + +constexpr std::uint64_t word_base {UINT64_C(1000000000)}; + +// x = n*pi/2 + r (mod 2pi), with |r| <= pi/4. +template +struct trig_arg +{ + boost::int128::uint128_t sig; // |r| = sig * 10^-k, with sig of 38 digits + int k; + fx rf; // |r| in fixed point, if fixed is set + bool fixed; + bool neg; // the sign of r, which includes the sign of x + unsigned n; // the quadrant 0 to 3, which includes the sign of x + bool zero; // x = 0: sig and k are not set +}; + +// |r| in fixed point: the binary reduction gives it, else it comes from sig and k. +template +constexpr auto fixed_r(const trig_arg& r) noexcept -> fx +{ + return r.fixed ? r.rf : fx_from(r.sig, r.k); +} + +// The significand in words of 9 digits, low word first; returns the count of words. +constexpr auto to_words(std::uint64_t m, std::uint64_t* out) noexcept -> int +{ + int n {}; + while (m != 0U) + { + out[n++] = m % word_base; + m /= word_base; + } + return n; +} + +constexpr auto to_words(std::uint32_t m, std::uint64_t* out) noexcept -> int +{ + return to_words(static_cast(m), out); +} + +constexpr auto to_words(const boost::int128::uint128_t& m, std::uint64_t* out) noexcept -> int +{ + constexpr std::uint64_t base2 {UINT64_C(1000000000000000000)}; + const auto hi {static_cast(m / base2)}; + const auto lo {static_cast(m % base2)}; + out[0] = lo % word_base; + out[1] = lo / word_base; + out[2] = hi % word_base; + out[3] = hi / word_base; + int n {4}; + while (n > 0 && out[n - 1] == 0U) + { + --n; + } + return n; +} + +// m * 10^shift in words of 9 digits, low word first, for 0 <= shift < 9; returns the count of words. +template +constexpr auto scaled_words(const Significand m, const int shift, std::uint64_t* out) noexcept -> int +{ + int count {to_words(m, out)}; + const auto scale {pow10(static_cast(shift))}; + std::uint64_t carry {}; + for (int i {}; i < count; ++i) + { + const std::uint64_t v {out[i] * scale + carry}; + out[i] = v % word_base; + carry = v / word_base; + } + if (carry != 0U) + { + out[count++] = carry; + } + return count; +} + +// The low Window words of m * 10^(9 * word_exp) * 2/pi, low word first. Word Window - 1 is the unit +// word, and 10^9 = 0 (mod 4), thus its low two bits are the quadrant and the higher words are not needed. +template +constexpr auto times_two_over_pi(const std::uint64_t* m, const int count, const int word_exp, std::uint64_t* out) noexcept -> void +{ + // two_over_pi[i] has the weight 10^(-9(i + 1)), thus word u of the window is table word Window - 2 - u + word_exp. + std::uint64_t w[static_cast(Window)] {}; + for (int u {}; u < Window; ++u) + { + const int i {Window - 2 - u + word_exp}; + w[u] = (i >= 0 && i < 720) ? static_cast(trig_table::two_over_pi[i]) : 0U; + } + + std::uint64_t carry {}; + for (int col {}; col < Window; ++col) + { + std::uint64_t acc {carry}; + for (int i {}; i < count && i <= col; ++i) + { + acc += m[i] * w[col - i]; + } + out[col] = acc % word_base; + carry = acc / word_base; + } +} + +// The fraction f in [0, 1) of Words words, high word first, to 1 - f when f >= 1/2; returns whether it did. +template +constexpr auto fold_half(std::uint64_t* f) noexcept -> bool +{ + if (f[0] < word_base / 2U) + { + return false; + } + for (int i {}; i < Words; ++i) + { + f[i] = word_base - 1U - f[i]; + } + for (int i {Words - 1}; i >= 0; --i) + { + if (++f[i] < word_base) + { + break; + } + f[i] = 0U; + } + return true; +} + +// The Frac words of the fraction from word lead on, high word first, times pi/2: 2 * Frac + 1 words, +// low word first. +template +constexpr auto times_half_pi(const std::uint64_t* f, const int lead, std::uint64_t* out) noexcept -> void +{ + std::uint64_t a[static_cast(Frac)] {}; + std::uint64_t b[static_cast(Frac + 1)] {}; + for (int u {}; u < Frac; ++u) + { + const int i {lead + Frac - 1 - u}; + a[u] = i < Words ? f[i] : 0U; + } + for (int v {}; v <= Frac; ++v) + { + b[v] = static_cast(trig_table::pio2[Frac - v]); + } + + std::uint64_t carry {}; + for (int col {}; col <= 2 * Frac; ++col) + { + std::uint64_t acc {carry}; + for (int u {col - Frac < 0 ? 0 : col - Frac}; u < Frac && u <= col; ++u) + { + acc += a[u] * b[col - u]; + } + out[col] = acc % word_base; + carry = acc / word_base; + } +} + +// The first 38 digits of the nonzero words r, low word first, as sig * 10^exp10 with word 0 at 10^0. +constexpr auto first_38_digits(const std::uint64_t* r, const int count, int& exp10) noexcept -> boost::int128::uint128_t +{ + int top {count - 1}; + while (r[top] == 0U) + { + --top; + } + boost::int128::uint128_t sig {r[top]}; + int digits {num_digits(r[top])}; + exp10 = 9 * top; + for (int i {top - 1}; i >= 0; --i) + { + if (digits + 9 <= 38) + { + sig = sig * word_base + r[i]; + digits += 9; + exp10 -= 9; + } + else + { + const int take {38 - digits}; + const auto low {pow10(static_cast(9 - take))}; + sig = sig * (word_base / low) + r[i] / low; + exp10 -= take; + break; + } + } + + const int sig_digits {num_digits(sig)}; + if (sig_digits < 38) + { + sig *= pow10(static_cast(38 - sig_digits)); + exp10 -= 38 - sig_digits; + } + return sig; +} + +// Payne-Hanek reduction in words of 9 decimal digits for x = m * 10^e: multiply m by the window of +// 2/pi which starts at the exponent, then the integer part mod 4 is the quadrant. +template +constexpr auto trig_reduce(const Significand m, const int e, const bool xneg) noexcept -> trig_arg::words> +{ + constexpr int words {trig_traits::words}; + constexpr int window {trig_traits::window}; + constexpr int frac {trig_traits::frac}; + static_assert(window - 2 + std::numeric_limits::max_exponent10 / 9 < 720, "2/pi must cover the largest exponent"); + + // m * 10^e = (m * 10^shift) * 10^(9 * word_exp) + const int shift {((e % 9) + 9) % 9}; + const int word_exp {(e - shift) / 9}; + std::uint64_t m_words[6] {}; + const int count {scaled_words(m, shift, m_words)}; + + std::uint64_t prod[static_cast(window)] {}; + times_two_over_pi(m_words, count, word_exp, prod); + + // The fraction, high word first. Above one half, go to the next quadrant with r < 0. + std::uint64_t f[static_cast(window - 1)] {}; + for (int i {}; i < window - 1; ++i) + { + f[i] = prod[window - 2 - i]; + } + unsigned n {static_cast(prod[window - 1] & 3U)}; + bool rneg {fold_half(f)}; + if (rneg) + { + n = (n + 1U) & 3U; + } + + // -x = (-n)*pi/2 + (-r) + if (xneg) + { + n = (4U - n) & 3U; + rneg = !rneg; + } + + trig_arg out {}; + out.n = n; + out.neg = rneg; + + int lead {}; + while (lead < window - 1 && f[lead] == 0U) + { + ++lead; + } + if (lead == window - 1) + { + // Only for an exact multiple of pi/2, which no nonzero x is. + out.zero = true; + return out; + } + + std::uint64_t r[static_cast(2 * frac + 1)] {}; + times_half_pi(f, lead, r); + + // Word 0 of r has the weight 10^-9(lead + 2 * frac). + int exp10 {}; + out.sig = first_38_digits(r, 2 * frac + 1, exp10); + out.k = 9 * (lead + 2 * frac) - exp10; + return out; +} + +// a * b + x + carry: returns the low word and puts the high word in carry. +constexpr auto mul_add(const std::uint64_t a, const std::uint64_t b, const std::uint64_t x, std::uint64_t& carry) noexcept -> std::uint64_t +{ + std::uint64_t hi {}; + std::uint64_t lo {boost::int128::detail::umul(a, b, hi)}; + lo += x; + hi += lo < x ? 1U : 0U; + lo += carry; + hi += lo < carry ? 1U : 0U; + carry = hi; + return lo; +} + +// a[0..count) * v, low word first, in place; returns the new count. +constexpr auto times_word(std::uint64_t* a, int count, const std::uint64_t v) noexcept -> int +{ + std::uint64_t carry {}; + for (int i {}; i < count; ++i) + { + a[i] = mul_add(a[i], v, 0U, carry); + } + if (carry != 0U) + { + a[count++] = carry; + } + return count; +} + +// x = m * 10^e < 10^19 as Frac + 1 words, low word first, with word Frac as the integer part. +// For e < 0, x = (m * 10^b) * 10^(-9a) with 9a = b - e, and the table error is below 2^-300. +template +constexpr auto binary_words(const Significand m, const int e, std::uint64_t* x) noexcept -> void +{ + constexpr int cw {Frac + 2}; + const auto m128 {static_cast(m)}; + if (e >= 0) + { + x[Frac] = m128.low * pow10(static_cast(e)); + return; + } + + const int a {(8 - e) / 9}; + std::uint64_t g[4] {m128.low, m128.high}; + const int count {times_word(g, m128.high != 0U ? 2 : 1, pow10(static_cast(9 * a + e)))}; + + const auto& c {trig_table::pow10_neg9_bin[a - 1]}; + std::uint64_t p[static_cast(cw + 3)] {}; + for (int i {}; i < count; ++i) + { + std::uint64_t carry {}; + for (int j {}; j < cw; ++j) + { + p[i + j] = mul_add(g[i], c[cw - 1 - j], p[i + j], carry); + } + p[i + cw] = carry; + } + for (int i {}; i <= Frac; ++i) + { + x[i] = p[cw - Frac + i]; + } +} + +// Payne-Hanek reduction in words of 64 bits for |x| < 10^19. The relative errors of r are below 2^-64, +// 2^-128 and 2^-190, less than those of trig_reduce. +template +constexpr auto trig_reduce_binary(const Significand m, const int e, const bool xneg) noexcept -> trig_arg::words> +{ + constexpr int words {trig_traits::words}; + constexpr int frac {trig_traits::bin_frac}; + constexpr int tw {frac + 2}; + constexpr int rw {words + 1}; + + std::uint64_t x[static_cast(frac + 1)] {}; + binary_words(m, e, x); + trig_arg out {}; + + // q = x * 2/pi with tw words of 2/pi: word 2 * frac + 2 is the integer part. The columns below + // frac only carry into words that are dropped. + std::uint64_t q[static_cast(2 * frac + 3)] {}; + for (int i {}; i <= frac; ++i) + { + std::uint64_t carry {}; + for (int j {frac - i < 0 ? 0 : frac - i}; j < tw; ++j) + { + q[i + j] = mul_add(x[i], trig_table::two_over_pi_bin[tw - 1 - j], q[i + j], carry); + } + q[i + tw] = carry; + } + + // The fraction, high word first. Above one half, go to the next quadrant with r = -(1 - f). + std::uint64_t f[static_cast(frac)] {}; + for (int i {}; i < frac; ++i) + { + f[i] = q[2 * frac + 1 - i]; + } + unsigned n {static_cast(q[2 * frac + 2] & 3U)}; + bool rneg {(f[0] >> 63U) != 0U}; + if (rneg) + { + std::uint64_t carry {1}; + for (int i {frac - 1}; i >= 0; --i) + { + f[i] = ~f[i] + carry; + carry = (carry != 0U && f[i] == 0U) ? 1U : 0U; + } + n = (n + 1U) & 3U; + } + if (xneg) + { + n = (4U - n) & 3U; + rneg = !rneg; + } + out.n = n; + out.neg = rneg; + + int lead {}; + while (lead < frac && f[lead] == 0U) + { + ++lead; + } + if (lead == frac) + { + // Only for an exact multiple of pi/2, which no nonzero x is. + out.zero = true; + return out; + } + + // r = (rw words of f from word lead) * pi/2 = (s from word rw) * 2^-64(lead + rw). The four spare + // words take the scale by 10^k below. + std::uint64_t s[static_cast(2 * rw + 6)] {}; + for (int u {}; u < rw; ++u) + { + const int i {lead + rw - 1 - u}; + const std::uint64_t a {i < frac ? f[i] : 0U}; + std::uint64_t carry {}; + for (int v {}; v <= rw; ++v) + { + s[u + v] = mul_add(a, trig_table::pio2_bin[rw - v], s[u + v], carry); + } + s[u + rw + 1] = carry; + } + std::uint64_t* const r {s + rw}; + + // rf = r * 2^(64N - 2): drop lead + 1 words and 2 bits. + for (int i {}; i < words; ++i) + { + const int w {i + lead + 1}; + const std::uint64_t lo {w <= rw ? r[w] : 0U}; + const std::uint64_t hi {w + 1 <= rw ? r[w + 1] : 0U}; + out.rf.w[i] = (lo >> 2U) | (hi << 62U); + } + out.fixed = true; + + // For r in [2^-(z + 1), 2^-z), sig = floor(r * 10^k) with k = 38 + floor(z * log10(2)) has 37 or 38 digits. + int count {rw + 1}; + while (r[count - 1] == 0U) + { + --count; + } + const int z {64 * (lead + rw - count) + countl_zero(r[count - 1])}; + int k {38 + ((z * 78913) >> 18)}; + int left {k}; + while (left >= 19) + { + count = times_word(r, count, pow10(static_cast(19))); + left -= 19; + } + count = times_word(r, count, pow10(static_cast(left))); + const int unit {lead + rw}; + out.sig = boost::int128::uint128_t {unit + 1 < count ? r[unit + 1] : 0U, unit < count ? r[unit] : 0U}; + if (out.sig < pow10(static_cast(37))) + { + out.sig *= 10U; + ++k; + } + out.k = k; + return out; +} + +// x = sig * 10^-k with sig of 38 digits; for |x| <= pi/4 this is r itself. +template +constexpr auto trig_prepare(const T x) noexcept -> trig_arg::words> +{ + constexpr int words {trig_traits::words}; + + // floor(pi/4 * 10^38) = 78539816339744830961566084581987572104 + constexpr boost::int128::uint128_t quarter_pi {UINT64_C(4257651975108193236), UINT64_C(12340784068543502728)}; + + trig_arg out {}; + const bool neg {signbit(x)}; + int e {}; + const auto m {frexp10(neg ? -x : x, &e)}; + out.neg = neg; + if (m == 0U) + { + out.zero = true; + return out; + } + + const auto digits {num_digits(m)}; + const auto sig {static_cast(m) * pow10(static_cast(38 - digits))}; + const int k {(38 - digits) - e}; + if (k > 38 || (k == 38 && sig <= quarter_pi)) + { + out.sig = sig; + out.k = k; + return out; + } + // |x| < 10^digits * 10^e + if (digits + e <= 19) + { + return trig_reduce_binary(m, e, neg); + } + return trig_reduce(m, e, neg); +} + +} // namespace trig +} // namespace detail +} // namespace decimal +} // namespace boost + +#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_TRIG_REDUCE_HPP diff --git a/include/boost/decimal/detail/cmath/sin.hpp b/include/boost/decimal/detail/cmath/sin.hpp index 829e8c952..7d35224cc 100644 --- a/include/boost/decimal/detail/cmath/sin.hpp +++ b/include/boost/decimal/detail/cmath/sin.hpp @@ -7,18 +7,17 @@ #define BOOST_DECIMAL_DETAIL_CMATH_SIN_HPP #include -#include #include #include #include -#include -#include +#include +#include #include #include #ifndef BOOST_DECIMAL_BUILD_MODULE +#include #include -#include #endif namespace boost { @@ -30,122 +29,27 @@ template constexpr auto sin_impl(const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { - T result { }; - - const auto fpc = fpclassify(x); - - // First check non-finite values and small angles. - if (fabs(x) < std::numeric_limits::epsilon() - #ifndef BOOST_DECIMAL_FAST_MATH - || (fpc == FP_INFINITE) || (fpc == FP_NAN) - #endif - ) + #ifndef BOOST_DECIMAL_FAST_MATH + if (isnan(x)) { - result = x; + return x; } - else if (signbit(x)) + if (isinf(x)) { - result = -sin(-x); + return std::numeric_limits::quiet_NaN(); } - else - { - if(x > 0) - { - // Perform argument reduction and subsequent scaling of the result. - - // Given x = k * (pi/2) + r, compute n = (k % 4). - - // | n | sin(x) | cos(x) | sin(x)/cos(x) | - // |----------------------------------------| - // | 0 | sin(r) | cos(r) | sin(r)/cos(r) | - // | 1 | cos(r) | -sin(r) | -cos(r)/sin(r) | - // | 2 | -sin(r) | -cos(r) | sin(r)/cos(r) | - // | 3 | -cos(r) | sin(r) | -cos(r)/sin(r) | - - const T two_x { x * 2 }; - - const unsigned k { static_cast(two_x / numbers::pi_v) }; - const auto n = k % static_cast(4); - - const T two_r { two_x - (numbers::pi_v * k) }; - - T r { two_r / 2 }; - - constexpr T one { 1 }; - - bool do_scaling { two_r > one }; - - constexpr T cbrt_epsilon { cbrt(std::numeric_limits::epsilon()) }; + #endif - switch(n) - { - case static_cast(UINT8_C(1)): - case static_cast(UINT8_C(3)): - { - const T d2r { numbers::pi_v - two_r }; - - if (d2r < cbrt_epsilon) - { - result = d2r * (one - (d2r * d2r) / 12) / 2; - - do_scaling = false; - } - else - { - if(do_scaling) - { - // Reduce the argument with one single factor of three. - r /= static_cast(UINT8_C(3)); - } - - result = detail::cos_series_expansion(r); - } - } - break; - - case static_cast(UINT8_C(0)): - case static_cast(UINT8_C(2)): - default: - { - if (two_r < cbrt_epsilon) - { - // Normal[Series[Sin[x/2], {x, 0, 3}]] - // FullSimplify[%] - // HornerForm[%] - - result = (two_r * (one - (two_r * two_r) / 24)) / 2; - } - else - { - if(do_scaling) - { - // Reduce the argument with one single factor of three. - r /= static_cast(UINT8_C(3)); - } - - result = detail::sin_series_expansion(r); - } - } - break; - } - - if(do_scaling) - { - result *= (static_cast(UINT8_C(3)) - ((result * result) * static_cast(UINT8_C(4)))); - } - - if(signbit(result)) - { - result = -result; - } - - const auto b_neg = (n > static_cast(UINT8_C(1))); - - if(b_neg) { result = -result; } - } + // x = n*pi/2 + r: sin(x) is sin(r), cos(r), -sin(r) or -cos(r) for n = 0 to 3. + const auto r {trig::trig_prepare(x)}; + if (r.zero) + { + return x; } - return result; + // The sign of sin(r) is r.neg, and cos(r) > 0. + return (r.n & 1U) == 0U ? trig::sin_of(r, r.neg != (r.n == 2U)) + : trig::cos_of(r, r.n == 3U); } } // namespace detail @@ -162,5 +66,4 @@ constexpr auto sin(const T x) noexcept } // namespace decimal } // namespace boost - #endif // BOOST_DECIMAL_DETAIL_CMATH_SIN_HPP diff --git a/include/boost/decimal/detail/cmath/tan.hpp b/include/boost/decimal/detail/cmath/tan.hpp index bcb32d657..32ecf543e 100644 --- a/include/boost/decimal/detail/cmath/tan.hpp +++ b/include/boost/decimal/detail/cmath/tan.hpp @@ -7,18 +7,16 @@ #define BOOST_DECIMAL_DETAIL_CMATH_TAN_HPP #include -#include #include #include #include -#include -#include -#include -#include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE +#include #include -#include #endif namespace boost { @@ -26,106 +24,29 @@ namespace decimal { namespace detail { -BOOST_DECIMAL_EXPORT template +template constexpr auto tan_impl(const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { - T result { }; - - const auto fpc = fpclassify(x); - - // First check non-finite values and small angles. #ifndef BOOST_DECIMAL_FAST_MATH - if (fabs(x) < std::numeric_limits::epsilon() || (fpc == FP_NAN)) - { - result = x; - } - else if (fpc == FP_INFINITE) - { - result = std::numeric_limits::quiet_NaN(); - } - else if (signbit(x)) + if (isnan(x)) { - result = -tan(-x); + return x; } - #else - if (fabs(x) < std::numeric_limits::epsilon()) + if (isinf(x)) { - result = x; + return std::numeric_limits::quiet_NaN(); } #endif - else - { - // Perform argument reduction. - - // Given x = k * (pi/2) + r, compute n = (k % 4). - - // | n | sin(x) | cos(x) | sin(x)/cos(x) | - // |----------------------------------------| - // | 0 | sin(r) | cos(r) | sin(r)/cos(r) | - // | 1 | cos(r) | -sin(r) | -cos(r)/sin(r) | - // | 2 | -sin(r) | -cos(r) | sin(r)/cos(r) | - // | 3 | -cos(r) | sin(r) | -cos(r)/sin(r) | - - const T two_x = x * 2; - - const auto k = static_cast(two_x / numbers::pi_v); - const auto n = k % static_cast(UINT8_C(4)); - const T two_r { two_x - (numbers::pi_v * k) }; - - const T r { two_r / 2 }; - - constexpr T cbrt_epsilon { cbrt(std::numeric_limits::epsilon()) }; - - constexpr T one { 1 }; - constexpr T two { 2 }; - - switch(n) - { - case static_cast(UINT8_C(1)): - case static_cast(UINT8_C(3)): - { - if (two_r < cbrt_epsilon) - { - // Normal[Series[Cos[x/2]/Sin[x/2], {x, 0, 3}]] - - result = (two / two_r) - (two_r * (one + (two_r * two_r) / 60) / 6); - } - else - { - result = cos(r) / sin(r); - } - - result = -result; - - break; - } - - case static_cast(UINT8_C(0)): - case static_cast(UINT8_C(2)): - default: - { - const T d2r { numbers::pi_v - two_r }; - - if (d2r < cbrt_epsilon) - { - // Use essentially the same series as shown above, but shifted via d2r. - - result = (two / d2r) - ((d2r * (one + (d2r * d2r) / 60)) / 6); - } - else - { - result = sin(r) / cos(r); - } - - break; - } - } - + // x = n*pi/2 + r: tan(x) is tan(r) for an even n, else -cot(r). + const auto r {trig::trig_prepare(x)}; + if (r.zero) + { + return x; } - return result; + return trig::tan_of(r); } } // namespace detail diff --git a/test/Jamfile b/test/Jamfile index 11eecd7c1..3bd2e8116 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -284,6 +284,7 @@ run test_to_string.cpp ; run test_total_ordering.cpp ; run test_toward_zero_rounding.cpp : : : off ; +run test_trig_rounding.cpp ; run test_upward_rounding.cpp : : : off ; run test_zeta.cpp ; diff --git a/test/test_edges_and_behave.cpp b/test/test_edges_and_behave.cpp index a8f931a92..088a07ac8 100644 --- a/test/test_edges_and_behave.cpp +++ b/test/test_edges_and_behave.cpp @@ -396,8 +396,8 @@ namespace local const auto result_sin_cos_non_normal_is_ok = ( - (isinf(sin_inf) && isnan(sin_nan)) - && (isinf(cos_inf) && isnan(cos_nan)) + (isnan(sin_inf) && isnan(sin_nan)) + && (isnan(cos_inf) && isnan(cos_nan)) ); BOOST_TEST(result_sin_cos_non_normal_is_ok); diff --git a/test/test_sin_cos.cpp b/test/test_sin_cos.cpp index 6c9381f5d..090cb2f8d 100644 --- a/test/test_sin_cos.cpp +++ b/test/test_sin_cos.cpp @@ -71,7 +71,7 @@ auto test_sin() -> void } } - BOOST_TEST(isinf(sin(std::numeric_limits::infinity() * Dec(dist(rng))))); + BOOST_TEST(isnan(sin(std::numeric_limits::infinity() * Dec(dist(rng))))); BOOST_TEST(isnan(sin(std::numeric_limits::quiet_NaN() * Dec(dist(rng))))); BOOST_TEST_EQ(abs(sin(Dec(0) * Dec(dist(rng)))), Dec(0)); @@ -119,7 +119,7 @@ auto test_cos() -> void } } - BOOST_TEST(isinf(cos(std::numeric_limits::infinity() * Dec(dist(rng))))); + BOOST_TEST(isnan(cos(std::numeric_limits::infinity() * Dec(dist(rng))))); BOOST_TEST(isnan(cos(std::numeric_limits::quiet_NaN() * Dec(dist(rng))))); BOOST_TEST_EQ(cos(Dec(0) * Dec(dist(rng))), Dec(1)); diff --git a/test/test_trig_rounding.cpp b/test/test_trig_rounding.cpp new file mode 100644 index 000000000..163e1ac76 --- /dev/null +++ b/test/test_trig_rounding.cpp @@ -0,0 +1,418 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// +// https://github.com/boostorg/decimal/issues/1471 +// +// sin, cos and tan must give the correctly rounded value, also for large arguments and in the +// directed rounding modes. The expected values come from MPFR at 600 bits. + +#include +#include +#include + +using namespace boost::decimal; +using namespace boost::decimal::literals; + +// Each type has inputs near pi/2, pi and 1000*pi/2, the value nearest to a multiple of pi/2 in its exponent +// range (4.327189e+48 for decimal32_t) and below 10^19 (76.96902), and the inputs on each side of 10^19. +template +void check(const T x, const T s, const T c, const T t) +{ + BOOST_TEST_EQ(sin(x), s); + BOOST_TEST_EQ(cos(x), c); + BOOST_TEST_EQ(tan(x), t); +} + +void test_nearest_decimal32_t() +{ + check(5.000000e-03_DF, 4.999979e-3_DF, 9.999875e-1_DF, 5.000042e-3_DF); + check(5.000000e-01_DF, 4.794255e-1_DF, 8.775826e-1_DF, 5.463025e-1_DF); + check(1.000000e+00_DF, 8.414710e-1_DF, 5.403023e-1_DF, 1.557408e0_DF); + check(-1.000000e+00_DF, -8.414710e-1_DF, 5.403023e-1_DF, -1.557408e0_DF); + check(3.000000e+00_DF, 1.411200e-1_DF, -9.899925e-1_DF, -1.425465e-1_DF); + check(1.000000e+02_DF, -5.063656e-1_DF, 8.623189e-1_DF, -5.872139e-1_DF); + check(1.570796e+00_DF, 1.000000e0_DF, 3.267949e-7_DF, 3.060023e6_DF); + check(3.141593e+00_DF, -3.464102e-7_DF, -1.000000e0_DF, 3.464102e-7_DF); + check(1.572367e+03_DF, 1.000000e0_DF, 1.231217e-4_DF, 8.122046e3_DF); + check(1.000000e+10_DF, -4.875060e-1_DF, 8.731196e-1_DF, -5.583496e-1_DF); + check(9.999999e+96_DF, 5.532571e-1_DF, 8.330106e-1_DF, 6.641658e-1_DF); + check(4.327189e+48_DF, -1.890807e-10_DF, 1.000000e0_DF, -1.890807e-10_DF); + check(7.696902e+01_DF, 1.000000e0_DF, 1.294993e-8_DF, 7.722047e7_DF); + check(9.999999e+18_DF, -9.628773e-1_DF, 2.699396e-1_DF, -3.567010e0_DF); + check(1.000000e+19_DF, -9.270632e-1_DF, -3.749052e-1_DF, 2.472794e0_DF); +} + +void test_nearest_decimal_fast32_t() +{ + check(5.000000e-03_DFF, 4.999979e-3_DFF, 9.999875e-1_DFF, 5.000042e-3_DFF); + check(5.000000e-01_DFF, 4.794255e-1_DFF, 8.775826e-1_DFF, 5.463025e-1_DFF); + check(1.000000e+00_DFF, 8.414710e-1_DFF, 5.403023e-1_DFF, 1.557408e0_DFF); + check(-1.000000e+00_DFF, -8.414710e-1_DFF, 5.403023e-1_DFF, -1.557408e0_DFF); + check(3.000000e+00_DFF, 1.411200e-1_DFF, -9.899925e-1_DFF, -1.425465e-1_DFF); + check(1.000000e+02_DFF, -5.063656e-1_DFF, 8.623189e-1_DFF, -5.872139e-1_DFF); + check(1.570796e+00_DFF, 1.000000e0_DFF, 3.267949e-7_DFF, 3.060023e6_DFF); + check(3.141593e+00_DFF, -3.464102e-7_DFF, -1.000000e0_DFF, 3.464102e-7_DFF); + check(1.572367e+03_DFF, 1.000000e0_DFF, 1.231217e-4_DFF, 8.122046e3_DFF); + check(1.000000e+10_DFF, -4.875060e-1_DFF, 8.731196e-1_DFF, -5.583496e-1_DFF); + check(9.999999e+96_DFF, 5.532571e-1_DFF, 8.330106e-1_DFF, 6.641658e-1_DFF); + check(4.327189e+48_DFF, -1.890807e-10_DFF, 1.000000e0_DFF, -1.890807e-10_DFF); + check(7.696902e+01_DFF, 1.000000e0_DFF, 1.294993e-8_DFF, 7.722047e7_DFF); + check(9.999999e+18_DFF, -9.628773e-1_DFF, 2.699396e-1_DFF, -3.567010e0_DFF); + check(1.000000e+19_DFF, -9.270632e-1_DFF, -3.749052e-1_DFF, 2.472794e0_DFF); +} + +void test_nearest_decimal64_t() +{ + check(5.000000000000000e-03_DD, 4.999979166692708e-3_DD, 9.999875000260416e-1_DD, 5.000041667083338e-3_DD); + check(5.000000000000000e-01_DD, 4.794255386042030e-1_DD, 8.775825618903727e-1_DD, 5.463024898437905e-1_DD); + check(1.000000000000000e+00_DD, 8.414709848078965e-1_DD, 5.403023058681397e-1_DD, 1.557407724654902e0_DD); + check(-1.000000000000000e+00_DD, -8.414709848078965e-1_DD, 5.403023058681397e-1_DD, -1.557407724654902e0_DD); + check(3.000000000000000e+00_DD, 1.411200080598672e-1_DD, -9.899924966004455e-1_DD, -1.425465430742778e-1_DD); + check(1.000000000000000e+02_DD, -5.063656411097588e-1_DD, 8.623188722876839e-1_DD, -5.872139151569291e-1_DD); + check(1.570796326794897e+00_DD, 1.000000000000000e0_DD, -3.807686783083602e-16_DD, -2.626266436731868e15_DD); + check(3.141592653589793e+00_DD, 2.384626433832795e-16_DD, -1.000000000000000e0_DD, -2.384626433832795e-16_DD); + check(1.572367123121692e+03_DD, 1.000000000000000e0_DD, -4.841494469866686e-13_DD, -2.065477935013599e12_DD); + check(1.000000000000000e+10_DD, -4.875060250875107e-1_DD, 8.731196226768560e-1_DD, -5.583496378112418e-1_DD); + check(9.999999999999999e+384_DD, 1.094503281143336e-1_DD, 9.939922664063663e-1_DD, 1.101118507793177e-1_DD); + check(8.919302781369317e+311_DD, -6.055274390996879e-20_DD, -1.000000000000000e0_DD, 6.055274390996879e-20_DD); + check(9.538513039607291e+06_DD, -6.771756623504492e-18_DD, -1.000000000000000e0_DD, 6.771756623504492e-18_DD); + check(9.999999999999999e+18_DD, -2.113595126444045e-1_DD, -9.774083877349937e-1_DD, 2.162448320442599e-1_DD); + check(1.000000000000000e+19_DD, -9.270631660486504e-1_DD, -3.749051695507178e-1_DD, 2.472793765846527e0_DD); +} + +void test_nearest_decimal_fast64_t() +{ + check(5.000000000000000e-03_DDF, 4.999979166692708e-3_DDF, 9.999875000260416e-1_DDF, 5.000041667083338e-3_DDF); + check(5.000000000000000e-01_DDF, 4.794255386042030e-1_DDF, 8.775825618903727e-1_DDF, 5.463024898437905e-1_DDF); + check(1.000000000000000e+00_DDF, 8.414709848078965e-1_DDF, 5.403023058681397e-1_DDF, 1.557407724654902e0_DDF); + check(-1.000000000000000e+00_DDF, -8.414709848078965e-1_DDF, 5.403023058681397e-1_DDF, -1.557407724654902e0_DDF); + check(3.000000000000000e+00_DDF, 1.411200080598672e-1_DDF, -9.899924966004455e-1_DDF, -1.425465430742778e-1_DDF); + check(1.000000000000000e+02_DDF, -5.063656411097588e-1_DDF, 8.623188722876839e-1_DDF, -5.872139151569291e-1_DDF); + check(1.570796326794897e+00_DDF, 1.000000000000000e0_DDF, -3.807686783083602e-16_DDF, -2.626266436731868e15_DDF); + check(3.141592653589793e+00_DDF, 2.384626433832795e-16_DDF, -1.000000000000000e0_DDF, -2.384626433832795e-16_DDF); + check(1.572367123121692e+03_DDF, 1.000000000000000e0_DDF, -4.841494469866686e-13_DDF, -2.065477935013599e12_DDF); + check(1.000000000000000e+10_DDF, -4.875060250875107e-1_DDF, 8.731196226768560e-1_DDF, -5.583496378112418e-1_DDF); + check(9.999999999999999e+384_DDF, 1.094503281143336e-1_DDF, 9.939922664063663e-1_DDF, 1.101118507793177e-1_DDF); + check(8.919302781369317e+311_DDF, -6.055274390996879e-20_DDF, -1.000000000000000e0_DDF, 6.055274390996879e-20_DDF); + check(9.538513039607291e+06_DDF, -6.771756623504492e-18_DDF, -1.000000000000000e0_DDF, 6.771756623504492e-18_DDF); + check(9.999999999999999e+18_DDF, -2.113595126444045e-1_DDF, -9.774083877349937e-1_DDF, 2.162448320442599e-1_DDF); + check(1.000000000000000e+19_DDF, -9.270631660486504e-1_DDF, -3.749051695507178e-1_DDF, 2.472793765846527e0_DDF); +} + +void test_nearest_decimal128_t() +{ + check(5.000000000000000000000000000000000e-03_DL, 4.999979166692708317832346652128958e-3_DL, 9.999875000260416449652874658951263e-1_DL, 5.000041667083337549645888880751339e-3_DL); + check(5.000000000000000000000000000000000e-01_DL, 4.794255386042030002732879352155714e-1_DL, 8.775825618903727161162815826038297e-1_DL, 5.463024898437905132551794657802854e-1_DL); + check(1.000000000000000000000000000000000e+00_DL, 8.414709848078965066525023216302990e-1_DL, 5.403023058681397174009366074429766e-1_DL, 1.557407724654902230506974807458360e0_DL); + check(-1.000000000000000000000000000000000e+00_DL, -8.414709848078965066525023216302990e-1_DL, 5.403023058681397174009366074429766e-1_DL, -1.557407724654902230506974807458360e0_DL); + check(3.000000000000000000000000000000000e+00_DL, 1.411200080598672221007448028081103e-1_DL, -9.899924966004454572715727947312613e-1_DL, -1.425465430742778052956354105339135e-1_DL); + check(1.000000000000000000000000000000000e+02_DL, -5.063656411097587936565576104597854e-1_DL, 8.623188722876839341019385139508425e-1_DL, -5.872139151569290766778096356445879e-1_DL); + check(1.570796326794896619231321691639751e+00_DL, 1.000000000000000000000000000000000e0_DL, 4.420985846996875529104874722961539e-34_DL, 2.261938930836633226244288822199802e33_DL); + check(3.141592653589793238462643383279503e+00_DL, -1.158028306006248941790250554076922e-34_DL, -1.000000000000000000000000000000000e0_DL, 1.158028306006248941790250554076922e-34_DL); + check(1.572367123121691515850553013331391e+03_DL, 1.000000000000000000000000000000000e0_DL, 1.935406832843872404633979597684501e-31_DL, 5.166872323844219553545949395867088e30_DL); + check(1.000000000000000000000000000000000e+10_DL, -4.875060250875106915277942943481060e-1_DL, 8.731196226768560011761913453076952e-1_DL, -5.583496378112418465618934073186368e-1_DL); + check(9.999999999999999999999999999999999e+6144_DL, 5.582907749092521238056875912416942e-1_DL, -8.296453523350967216289114223081673e-1_DL, -6.729270203682844056779140311680751e-1_DL); + check(2.344813655066356855719930664718056e+1414_DL, -1.000000000000000000000000000000000e0_DL, 1.030557387629248882465543827418861e-37_DL, -9.703486792719571064179898782826551e36_DL); + check(2.589477332551582209163675623249625e+08_DL, -1.000000000000000000000000000000000e0_DL, 1.532138639232766081410107492878833e-35_DL, -6.526824494816997196224318392388852e34_DL); + check(9.999999999999999999999999999999999e+18_DL, -9.270631660486500103289527389466423e-1_DL, -3.749051695507187572163865946598085e-1_DL, 2.472793765846520245629178386203776e0_DL); + check(1.000000000000000000000000000000000e+19_DL, -9.270631660486503852341222896649360e-1_DL, -3.749051695507178301532205460096107e-1_DL, 2.472793765846527360338186795636566e0_DL); +} + +void test_nearest_decimal_fast128_t() +{ + check(5.000000000000000000000000000000000e-03_DLF, 4.999979166692708317832346652128958e-3_DLF, 9.999875000260416449652874658951263e-1_DLF, 5.000041667083337549645888880751339e-3_DLF); + check(5.000000000000000000000000000000000e-01_DLF, 4.794255386042030002732879352155714e-1_DLF, 8.775825618903727161162815826038297e-1_DLF, 5.463024898437905132551794657802854e-1_DLF); + check(1.000000000000000000000000000000000e+00_DLF, 8.414709848078965066525023216302990e-1_DLF, 5.403023058681397174009366074429766e-1_DLF, 1.557407724654902230506974807458360e0_DLF); + check(-1.000000000000000000000000000000000e+00_DLF, -8.414709848078965066525023216302990e-1_DLF, 5.403023058681397174009366074429766e-1_DLF, -1.557407724654902230506974807458360e0_DLF); + check(3.000000000000000000000000000000000e+00_DLF, 1.411200080598672221007448028081103e-1_DLF, -9.899924966004454572715727947312613e-1_DLF, -1.425465430742778052956354105339135e-1_DLF); + check(1.000000000000000000000000000000000e+02_DLF, -5.063656411097587936565576104597854e-1_DLF, 8.623188722876839341019385139508425e-1_DLF, -5.872139151569290766778096356445879e-1_DLF); + check(1.570796326794896619231321691639751e+00_DLF, 1.000000000000000000000000000000000e0_DLF, 4.420985846996875529104874722961539e-34_DLF, 2.261938930836633226244288822199802e33_DLF); + check(3.141592653589793238462643383279503e+00_DLF, -1.158028306006248941790250554076922e-34_DLF, -1.000000000000000000000000000000000e0_DLF, 1.158028306006248941790250554076922e-34_DLF); + check(1.572367123121691515850553013331391e+03_DLF, 1.000000000000000000000000000000000e0_DLF, 1.935406832843872404633979597684501e-31_DLF, 5.166872323844219553545949395867088e30_DLF); + check(1.000000000000000000000000000000000e+10_DLF, -4.875060250875106915277942943481060e-1_DLF, 8.731196226768560011761913453076952e-1_DLF, -5.583496378112418465618934073186368e-1_DLF); + check(9.999999999999999999999999999999999e+6144_DLF, 5.582907749092521238056875912416942e-1_DLF, -8.296453523350967216289114223081673e-1_DLF, -6.729270203682844056779140311680751e-1_DLF); + check(2.344813655066356855719930664718056e+1414_DLF, -1.000000000000000000000000000000000e0_DLF, 1.030557387629248882465543827418861e-37_DLF, -9.703486792719571064179898782826551e36_DLF); + check(2.589477332551582209163675623249625e+08_DLF, -1.000000000000000000000000000000000e0_DLF, 1.532138639232766081410107492878833e-35_DLF, -6.526824494816997196224318392388852e34_DLF); + check(9.999999999999999999999999999999999e+18_DLF, -9.270631660486500103289527389466423e-1_DLF, -3.749051695507187572163865946598085e-1_DLF, 2.472793765846520245629178386203776e0_DLF); + check(1.000000000000000000000000000000000e+19_DLF, -9.270631660486503852341222896649360e-1_DLF, -3.749051695507178301532205460096107e-1_DLF, 2.472793765846527360338186795636566e0_DLF); +} + +void test_downward_decimal32_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000e-20_DF, 9.999999e-21_DF, 9.999999e-1_DF, 1.000000e-20_DF); + check(-1.000000e-20_DF, -1.000000e-20_DF, 9.999999e-1_DF, -1.000001e-20_DF); + check(-1.000000e+00_DF, -8.414710e-1_DF, 5.403023e-1_DF, -1.557408e0_DF); + check(1.570796e+00_DF, 9.999999e-1_DF, 3.267948e-7_DF, 3.060023e6_DF); + check(3.141593e+00_DF, -3.464103e-7_DF, -1.000000e0_DF, 3.464102e-7_DF); + check(4.327189e+48_DF, -1.890808e-10_DF, 9.999999e-1_DF, -1.890808e-10_DF); + check(7.696902e+01_DF, 9.999999e-1_DF, 1.294993e-8_DF, 7.722046e7_DF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_downward_decimal_fast32_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000e-20_DFF, 9.999999e-21_DFF, 9.999999e-1_DFF, 1.000000e-20_DFF); + check(-1.000000e-20_DFF, -1.000000e-20_DFF, 9.999999e-1_DFF, -1.000001e-20_DFF); + check(-1.000000e+00_DFF, -8.414710e-1_DFF, 5.403023e-1_DFF, -1.557408e0_DFF); + check(1.570796e+00_DFF, 9.999999e-1_DFF, 3.267948e-7_DFF, 3.060023e6_DFF); + check(3.141593e+00_DFF, -3.464103e-7_DFF, -1.000000e0_DFF, 3.464102e-7_DFF); + check(4.327189e+48_DFF, -1.890808e-10_DFF, 9.999999e-1_DFF, -1.890808e-10_DFF); + check(7.696902e+01_DFF, 9.999999e-1_DFF, 1.294993e-8_DFF, 7.722046e7_DFF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_downward_decimal64_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000000000000e-20_DD, 9.999999999999999e-21_DD, 9.999999999999999e-1_DD, 1.000000000000000e-20_DD); + check(-1.000000000000000e-20_DD, -1.000000000000000e-20_DD, 9.999999999999999e-1_DD, -1.000000000000001e-20_DD); + check(-1.000000000000000e+00_DD, -8.414709848078966e-1_DD, 5.403023058681397e-1_DD, -1.557407724654903e0_DD); + check(1.570796326794897e+00_DD, 9.999999999999999e-1_DD, -3.807686783083603e-16_DD, -2.626266436731868e15_DD); + check(3.141592653589793e+00_DD, 2.384626433832795e-16_DD, -1.000000000000000e0_DD, -2.384626433832796e-16_DD); + check(8.919302781369317e+311_DD, -6.055274390996880e-20_DD, -1.000000000000000e0_DD, 6.055274390996879e-20_DD); + check(9.538513039607291e+06_DD, -6.771756623504492e-18_DD, -1.000000000000000e0_DD, 6.771756623504491e-18_DD); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_downward_decimal_fast64_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000000000000e-20_DDF, 9.999999999999999e-21_DDF, 9.999999999999999e-1_DDF, 1.000000000000000e-20_DDF); + check(-1.000000000000000e-20_DDF, -1.000000000000000e-20_DDF, 9.999999999999999e-1_DDF, -1.000000000000001e-20_DDF); + check(-1.000000000000000e+00_DDF, -8.414709848078966e-1_DDF, 5.403023058681397e-1_DDF, -1.557407724654903e0_DDF); + check(1.570796326794897e+00_DDF, 9.999999999999999e-1_DDF, -3.807686783083603e-16_DDF, -2.626266436731868e15_DDF); + check(3.141592653589793e+00_DDF, 2.384626433832795e-16_DDF, -1.000000000000000e0_DDF, -2.384626433832796e-16_DDF); + check(8.919302781369317e+311_DDF, -6.055274390996880e-20_DDF, -1.000000000000000e0_DDF, 6.055274390996879e-20_DDF); + check(9.538513039607291e+06_DDF, -6.771756623504492e-18_DDF, -1.000000000000000e0_DDF, 6.771756623504491e-18_DDF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_downward_decimal128_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000000000000000000000000000000e-20_DL, 9.999999999999999999999999999999999e-21_DL, 9.999999999999999999999999999999999e-1_DL, 1.000000000000000000000000000000000e-20_DL); + check(-1.000000000000000000000000000000000e-20_DL, -1.000000000000000000000000000000000e-20_DL, 9.999999999999999999999999999999999e-1_DL, -1.000000000000000000000000000000001e-20_DL); + check(-1.000000000000000000000000000000000e+00_DL, -8.414709848078965066525023216302990e-1_DL, 5.403023058681397174009366074429766e-1_DL, -1.557407724654902230506974807458361e0_DL); + check(1.570796326794896619231321691639751e+00_DL, 9.999999999999999999999999999999999e-1_DL, 4.420985846996875529104874722961539e-34_DL, 2.261938930836633226244288822199802e33_DL); + check(3.141592653589793238462643383279503e+00_DL, -1.158028306006248941790250554076922e-34_DL, -1.000000000000000000000000000000000e0_DL, 1.158028306006248941790250554076921e-34_DL); + check(2.344813655066356855719930664718056e+1414_DL, -1.000000000000000000000000000000000e0_DL, 1.030557387629248882465543827418861e-37_DL, -9.703486792719571064179898782826551e36_DL); + check(2.589477332551582209163675623249625e+08_DL, -1.000000000000000000000000000000000e0_DL, 1.532138639232766081410107492878832e-35_DL, -6.526824494816997196224318392388852e34_DL); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_downward_decimal_fast128_t() +{ + fesetround(rounding_mode::fe_dec_downward); + check(1.000000000000000000000000000000000e-20_DLF, 9.999999999999999999999999999999999e-21_DLF, 9.999999999999999999999999999999999e-1_DLF, 1.000000000000000000000000000000000e-20_DLF); + check(-1.000000000000000000000000000000000e-20_DLF, -1.000000000000000000000000000000000e-20_DLF, 9.999999999999999999999999999999999e-1_DLF, -1.000000000000000000000000000000001e-20_DLF); + check(-1.000000000000000000000000000000000e+00_DLF, -8.414709848078965066525023216302990e-1_DLF, 5.403023058681397174009366074429766e-1_DLF, -1.557407724654902230506974807458361e0_DLF); + check(1.570796326794896619231321691639751e+00_DLF, 9.999999999999999999999999999999999e-1_DLF, 4.420985846996875529104874722961539e-34_DLF, 2.261938930836633226244288822199802e33_DLF); + check(3.141592653589793238462643383279503e+00_DLF, -1.158028306006248941790250554076922e-34_DLF, -1.000000000000000000000000000000000e0_DLF, 1.158028306006248941790250554076921e-34_DLF); + check(2.344813655066356855719930664718056e+1414_DLF, -1.000000000000000000000000000000000e0_DLF, 1.030557387629248882465543827418861e-37_DLF, -9.703486792719571064179898782826551e36_DLF); + check(2.589477332551582209163675623249625e+08_DLF, -1.000000000000000000000000000000000e0_DLF, 1.532138639232766081410107492878832e-35_DLF, -6.526824494816997196224318392388852e34_DLF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal32_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000e-20_DF, 1.000000e-20_DF, 1.000000e0_DF, 1.000001e-20_DF); + check(-1.000000e-20_DF, -9.999999e-21_DF, 1.000000e0_DF, -1.000000e-20_DF); + check(-1.000000e+00_DF, -8.414709e-1_DF, 5.403024e-1_DF, -1.557407e0_DF); + check(1.570796e+00_DF, 1.000000e0_DF, 3.267949e-7_DF, 3.060024e6_DF); + check(3.141593e+00_DF, -3.464102e-7_DF, -9.999999e-1_DF, 3.464103e-7_DF); + check(4.327189e+48_DF, -1.890807e-10_DF, 1.000000e0_DF, -1.890807e-10_DF); + check(7.696902e+01_DF, 1.000000e0_DF, 1.294994e-8_DF, 7.722047e7_DF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal_fast32_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000e-20_DFF, 1.000000e-20_DFF, 1.000000e0_DFF, 1.000001e-20_DFF); + check(-1.000000e-20_DFF, -9.999999e-21_DFF, 1.000000e0_DFF, -1.000000e-20_DFF); + check(-1.000000e+00_DFF, -8.414709e-1_DFF, 5.403024e-1_DFF, -1.557407e0_DFF); + check(1.570796e+00_DFF, 1.000000e0_DFF, 3.267949e-7_DFF, 3.060024e6_DFF); + check(3.141593e+00_DFF, -3.464102e-7_DFF, -9.999999e-1_DFF, 3.464103e-7_DFF); + check(4.327189e+48_DFF, -1.890807e-10_DFF, 1.000000e0_DFF, -1.890807e-10_DFF); + check(7.696902e+01_DFF, 1.000000e0_DFF, 1.294994e-8_DFF, 7.722047e7_DFF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal64_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000000000000e-20_DD, 1.000000000000000e-20_DD, 1.000000000000000e0_DD, 1.000000000000001e-20_DD); + check(-1.000000000000000e-20_DD, -9.999999999999999e-21_DD, 1.000000000000000e0_DD, -1.000000000000000e-20_DD); + check(-1.000000000000000e+00_DD, -8.414709848078965e-1_DD, 5.403023058681398e-1_DD, -1.557407724654902e0_DD); + check(1.570796326794897e+00_DD, 1.000000000000000e0_DD, -3.807686783083602e-16_DD, -2.626266436731867e15_DD); + check(3.141592653589793e+00_DD, 2.384626433832796e-16_DD, -9.999999999999999e-1_DD, -2.384626433832795e-16_DD); + check(8.919302781369317e+311_DD, -6.055274390996879e-20_DD, -9.999999999999999e-1_DD, 6.055274390996880e-20_DD); + check(9.538513039607291e+06_DD, -6.771756623504491e-18_DD, -9.999999999999999e-1_DD, 6.771756623504492e-18_DD); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal_fast64_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000000000000e-20_DDF, 1.000000000000000e-20_DDF, 1.000000000000000e0_DDF, 1.000000000000001e-20_DDF); + check(-1.000000000000000e-20_DDF, -9.999999999999999e-21_DDF, 1.000000000000000e0_DDF, -1.000000000000000e-20_DDF); + check(-1.000000000000000e+00_DDF, -8.414709848078965e-1_DDF, 5.403023058681398e-1_DDF, -1.557407724654902e0_DDF); + check(1.570796326794897e+00_DDF, 1.000000000000000e0_DDF, -3.807686783083602e-16_DDF, -2.626266436731867e15_DDF); + check(3.141592653589793e+00_DDF, 2.384626433832796e-16_DDF, -9.999999999999999e-1_DDF, -2.384626433832795e-16_DDF); + check(8.919302781369317e+311_DDF, -6.055274390996879e-20_DDF, -9.999999999999999e-1_DDF, 6.055274390996880e-20_DDF); + check(9.538513039607291e+06_DDF, -6.771756623504491e-18_DDF, -9.999999999999999e-1_DDF, 6.771756623504492e-18_DDF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal128_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000000000000000000000000000000e-20_DL, 1.000000000000000000000000000000000e-20_DL, 1.000000000000000000000000000000000e0_DL, 1.000000000000000000000000000000001e-20_DL); + check(-1.000000000000000000000000000000000e-20_DL, -9.999999999999999999999999999999999e-21_DL, 1.000000000000000000000000000000000e0_DL, -1.000000000000000000000000000000000e-20_DL); + check(-1.000000000000000000000000000000000e+00_DL, -8.414709848078965066525023216302989e-1_DL, 5.403023058681397174009366074429767e-1_DL, -1.557407724654902230506974807458360e0_DL); + check(1.570796326794896619231321691639751e+00_DL, 1.000000000000000000000000000000000e0_DL, 4.420985846996875529104874722961540e-34_DL, 2.261938930836633226244288822199803e33_DL); + check(3.141592653589793238462643383279503e+00_DL, -1.158028306006248941790250554076921e-34_DL, -9.999999999999999999999999999999999e-1_DL, 1.158028306006248941790250554076922e-34_DL); + check(2.344813655066356855719930664718056e+1414_DL, -9.999999999999999999999999999999999e-1_DL, 1.030557387629248882465543827418862e-37_DL, -9.703486792719571064179898782826550e36_DL); + check(2.589477332551582209163675623249625e+08_DL, -9.999999999999999999999999999999999e-1_DL, 1.532138639232766081410107492878833e-35_DL, -6.526824494816997196224318392388851e34_DL); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +void test_upward_decimal_fast128_t() +{ + fesetround(rounding_mode::fe_dec_upward); + check(1.000000000000000000000000000000000e-20_DLF, 1.000000000000000000000000000000000e-20_DLF, 1.000000000000000000000000000000000e0_DLF, 1.000000000000000000000000000000001e-20_DLF); + check(-1.000000000000000000000000000000000e-20_DLF, -9.999999999999999999999999999999999e-21_DLF, 1.000000000000000000000000000000000e0_DLF, -1.000000000000000000000000000000000e-20_DLF); + check(-1.000000000000000000000000000000000e+00_DLF, -8.414709848078965066525023216302989e-1_DLF, 5.403023058681397174009366074429767e-1_DLF, -1.557407724654902230506974807458360e0_DLF); + check(1.570796326794896619231321691639751e+00_DLF, 1.000000000000000000000000000000000e0_DLF, 4.420985846996875529104874722961540e-34_DLF, 2.261938930836633226244288822199803e33_DLF); + check(3.141592653589793238462643383279503e+00_DLF, -1.158028306006248941790250554076921e-34_DLF, -9.999999999999999999999999999999999e-1_DLF, 1.158028306006248941790250554076922e-34_DLF); + check(2.344813655066356855719930664718056e+1414_DLF, -9.999999999999999999999999999999999e-1_DLF, 1.030557387629248882465543827418862e-37_DLF, -9.703486792719571064179898782826550e36_DLF); + check(2.589477332551582209163675623249625e+08_DLF, -9.999999999999999999999999999999999e-1_DLF, 1.532138639232766081410107492878833e-35_DLF, -6.526824494816997196224318392388851e34_DLF); + fesetround(rounding_mode::fe_dec_to_nearest); +} + +// Over adjacent inputs where the exact sin and cos go in the directions sin_dir and cos_dir and tan +// goes up, the results must not go the other way. Equal results are correct rounding. +template +void check_chain(const T start, const int steps, const int sin_dir, const int cos_dir) +{ + const rounding_mode modes[] {rounding_mode::fe_dec_to_nearest, rounding_mode::fe_dec_upward, rounding_mode::fe_dec_downward}; + + for (const auto mode : modes) + { + fesetround(mode); + + T a {start}; + T sa {sin(a)}; + T ca {cos(a)}; + T ta {tan(a)}; + for (int i {}; i < steps; ++i) + { + const T b {nextafter(a, (std::numeric_limits::max)())}; + const T sb {sin(b)}; + const T cb {cos(b)}; + const T tb {tan(b)}; + + BOOST_TEST(sin_dir > 0 ? sb >= sa : sb <= sa); + BOOST_TEST(cos_dir > 0 ? cb >= ca : cb <= ca); + BOOST_TEST(tb >= ta); + + a = b; + sa = sb; + ca = cb; + ta = tb; + } + + fesetround(rounding_mode::fe_dec_to_nearest); + } +} + +// The chains cross the change of method at pi/4, and end 100 steps below pi/2, pi and 5*pi/2. +template +void test_monotonic(const T quarter_pi, const T half_pi, const T pi, const T five_half_pi) +{ + const auto back = [](T x, const int n) + { + for (int i {}; i < n; ++i) + { + x = nextafter(x, T{0}); + } + return x; + }; + + constexpr int steps {200}; + check_chain(back(quarter_pi, steps / 2), steps, 1, -1); + check_chain(back(half_pi, steps + 100), steps, 1, -1); + check_chain(back(pi, steps + 100), steps, -1, -1); + check_chain(back(five_half_pi, steps + 100), steps, 1, -1); +} + +// sin and cos of an infinity are NaN, as tan already is. +template +void test_special() +{ + #ifndef BOOST_DECIMAL_FAST_MATH + const T inf {std::numeric_limits::infinity()}; + + BOOST_TEST(isnan(sin(inf))); + BOOST_TEST(isnan(sin(-inf))); + BOOST_TEST(isnan(cos(inf))); + BOOST_TEST(isnan(cos(-inf))); + BOOST_TEST(isnan(tan(inf))); + BOOST_TEST(isnan(sin(std::numeric_limits::quiet_NaN()))); + #endif + BOOST_TEST(signbit(sin(-T{0}))); + BOOST_TEST(signbit(tan(-T{0}))); + BOOST_TEST_EQ(cos(-T{0}), T{1}); +} + +int main() +{ + test_nearest_decimal32_t(); + test_nearest_decimal_fast32_t(); + test_nearest_decimal64_t(); + test_nearest_decimal_fast64_t(); + test_nearest_decimal128_t(); + test_nearest_decimal_fast128_t(); + + #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION + // fesetround has no effect here, thus only the nearest mode is tested + test_downward_decimal32_t(); + test_downward_decimal_fast32_t(); + test_downward_decimal64_t(); + test_downward_decimal_fast64_t(); + test_downward_decimal128_t(); + test_downward_decimal_fast128_t(); + test_upward_decimal32_t(); + test_upward_decimal_fast32_t(); + test_upward_decimal64_t(); + test_upward_decimal_fast64_t(); + test_upward_decimal128_t(); + test_upward_decimal_fast128_t(); + #endif + + test_monotonic(0.7853981633974483096156608458198757_DF, 1.570796326794896619231321691639751_DF, + 3.141592653589793238462643383279503_DF, 7.853981633974483096156608458198757_DF); + test_monotonic(0.7853981633974483096156608458198757_DFF, 1.570796326794896619231321691639751_DFF, + 3.141592653589793238462643383279503_DFF, 7.853981633974483096156608458198757_DFF); + test_monotonic(0.7853981633974483096156608458198757_DD, 1.570796326794896619231321691639751_DD, + 3.141592653589793238462643383279503_DD, 7.853981633974483096156608458198757_DD); + test_monotonic(0.7853981633974483096156608458198757_DDF, 1.570796326794896619231321691639751_DDF, + 3.141592653589793238462643383279503_DDF, 7.853981633974483096156608458198757_DDF); + test_monotonic(0.7853981633974483096156608458198757_DL, 1.570796326794896619231321691639751_DL, + 3.141592653589793238462643383279503_DL, 7.853981633974483096156608458198757_DL); + test_monotonic(0.7853981633974483096156608458198757_DLF, 1.570796326794896619231321691639751_DLF, + 3.141592653589793238462643383279503_DLF, 7.853981633974483096156608458198757_DLF); + + test_special(); + test_special(); + test_special(); + test_special(); + test_special(); + test_special(); + + return boost::report_errors(); +}