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(); +}