diff --git a/include/boost/decimal/detail/cmath/acos.hpp b/include/boost/decimal/detail/cmath/acos.hpp index a70df1e5d..a88fc4271 100644 --- a/include/boost/decimal/detail/cmath/acos.hpp +++ b/include/boost/decimal/detail/cmath/acos.hpp @@ -6,7 +6,6 @@ #define BOOST_DECIMAL_DETAIL_CMATH_ACOS_HPP #include -#include #include #include #include @@ -14,6 +13,7 @@ #include #include #include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE @@ -37,7 +37,6 @@ constexpr auto acos_impl(const T x) noexcept } #endif - constexpr auto half_pi {numbers::pi_v / 2}; const auto absx {fabs(static_cast(x))}; T result {}; @@ -46,21 +45,31 @@ constexpr auto acos_impl(const T x) noexcept { result = std::numeric_limits::quiet_NaN(); } - else if (x < T{-5, -1}) + else if (absx <= T{5, -1}) { - result = numbers::pi_v - 2 * detail::asin_series(sqrt((1 - absx) / 2)); - } - else if (x < -std::numeric_limits::epsilon()) - { - result = half_pi + detail::asin_series(absx); - } - else if (x < T{5, -1}) - { - result = half_pi - detail::asin_series(x); + result = detail::split_pi_values(detail::split_pi_detail::half_pi_hi) - + (detail::asin_series(x) - detail::split_pi_values(detail::split_pi_detail::half_pi_lo)); } else { - result = half_pi - (half_pi - 2 * detail::asin_series(sqrt((1 - x) / 2))); + // acos(|x|) = 2 asin(s) and acos(-|x|) = pi - 2 asin(s), with asin(s) = s + lo. + // Only the last operation rounds at the grid of the result, which keeps acos monotone. + T s {}; + const T lo {detail::asin_half_angle(absx, s)}; + + if (x > 0) + { + // s + s rounds for s in [0.05, 0.1); add the exact part it rounds off to the low part + const T two_s {s + s}; + const T two_s_lo {s - (two_s - s)}; + result = two_s + ((lo + lo) + two_s_lo); + } + else + { + // s + s can round here too, but by less than 1/10 ulp of a result above 2, so it is not added + result = detail::split_pi_values(detail::split_pi_detail::pi_hi) - + ((s + s) + ((lo + lo) - detail::split_pi_values(detail::split_pi_detail::pi_lo))); + } } return result; diff --git a/include/boost/decimal/detail/cmath/asin.hpp b/include/boost/decimal/detail/cmath/asin.hpp index eecaf665f..937c1f3c6 100644 --- a/include/boost/decimal/detail/cmath/asin.hpp +++ b/include/boost/decimal/detail/cmath/asin.hpp @@ -7,7 +7,6 @@ #define BOOST_DECIMAL_DETAIL_CMATH_ASIN_HPP #include -#include #include #include #include @@ -15,6 +14,7 @@ #include #include #include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -52,7 +52,7 @@ constexpr auto asin_impl(const T x) noexcept if (absx <= cbrt_eps) { - result = absx * (one + (absx / 6) * absx); + result = absx + (absx * absx) * (absx / 6); } else if (absx <= T { 5, -1 }) { @@ -60,11 +60,15 @@ constexpr auto asin_impl(const T x) noexcept } else { - constexpr T half_pi { numbers::pi_v / 2 }; - if (absx < one) { - result = half_pi - 2 * asin_series(sqrt((1 - absx) / 2)); + // asin(x) = 2 (pi/4_hi - s) - 2 (lo - pi/4_lo), as pi/2_hi - 2 s can round before lo is added; + // head is exact for s >= 0.1 (else off by < 0.1 ulp of the result), and so is head + head below 0.5 + T s { }; + const T lo { asin_half_angle(absx, s) - split_pi_values(split_pi_detail::quarter_pi_lo) }; + const T head { split_pi_values(split_pi_detail::quarter_pi_hi) - s }; + + result = head < T { 5, -1 } ? (head + head) - (lo + lo) : head + (head - (lo + lo)); } else if (absx > one) { @@ -76,7 +80,7 @@ constexpr auto asin_impl(const T x) noexcept } else { - result = half_pi; + result = split_pi_values(split_pi_detail::half_pi_hi); } } diff --git a/include/boost/decimal/detail/cmath/atan.hpp b/include/boost/decimal/detail/cmath/atan.hpp index f9581f948..2f5fc2656 100644 --- a/include/boost/decimal/detail/cmath/atan.hpp +++ b/include/boost/decimal/detail/cmath/atan.hpp @@ -8,10 +8,12 @@ #include // NOLINT(llvm-include-order) #include +#include #include #include #include -#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -23,16 +25,15 @@ namespace decimal { namespace detail { -template -constexpr auto atan_impl(const T x) noexcept +// atan(x + x_lo); with has_lo, x must be a positive normal number and x_lo a small correction +template +constexpr auto atan_impl(const T x, const T x_lo = T { }) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { const auto fpc { fpclassify(x) }; T result { }; - constexpr T my_pi_half { numbers::pi_v / 2 }; - if (fpc == FP_ZERO #ifndef BOOST_DECIMAL_FAST_MATH || fpc == FP_NAN @@ -48,51 +49,73 @@ constexpr auto atan_impl(const T x) noexcept #ifndef BOOST_DECIMAL_FAST_MATH else if (fpc == FP_INFINITE) { - result = my_pi_half; + result = detail::split_pi_values(detail::split_pi_detail::half_pi_hi); } #endif - else + // The breakpoints are those of fdlibm: 7/16, 11/16, 19/16 and 39/16 + else if (x < T { 4375, -4 }) { - constexpr T one { 1 }; - - if (x <= T { 48 }) + // x_lo changes atan by x_lo / (1 + x^2) + result = has_lo ? x + (detail::atan_tail(x) + x_lo / (T { 1 } + x * x)) : detail::atan_series(x); + } + else if (x < T { 24375, -4 }) + { + // atan(x) = atan(c) + atan(q) with q = (x - c) / (1 + x c) for c = 1/2, 1, 3/2; + // num is exact, and den_lo is the part of the denominator that den rounds off + T c_hi { }; + T c_lo { }; + T num { }; + T den { }; + T den_lo { }; + + if (x < T { 6875, -4 }) { - // Define small-ish arguments to be less than 39/16. - const bool is_smallish { x <= T { 24375, -4 } }; - - // The portion of the algorithm for arc-tangent regarding scaling large-valued - // argument is based on Chapter 11, page 194 of Cody and Waite, "Software Manual - // for the Elementary Functions", Prentice Hall, 1980. - - const T - fx_arg - { - (!is_smallish) - ? ((x * numbers::sqrt3_v) - one) / (numbers::sqrt3_v + x) - : x - }; - - constexpr T half { 5, -1 }; - constexpr T three_halves { 15, -1 }; - - result = (fx_arg <= std::numeric_limits::epsilon()) ? fx_arg - : (fx_arg <= T { 4375, -4 }) ? detail::atan_series (fx_arg) - : (fx_arg <= T { 6875, -4 }) ? detail::atan_values(0U) + detail::atan_series((fx_arg - half) / (one + fx_arg / 2)) - : (fx_arg <= T { 11875, -4 }) ? detail::atan_values(1U) + detail::atan_series((fx_arg - one) / (fx_arg + one)) - : detail::atan_values(2U) + detail::atan_series((fx_arg - three_halves) / (one + three_halves * fx_arg)) - ; - - if(!is_smallish) - { - constexpr T my_pi_over_six { numbers::pi_v / 6 }; - - result += my_pi_over_six; - } + c_hi = detail::atan_values(detail::atan_detail::atan_half_hi); + c_lo = detail::atan_values(detail::atan_detail::atan_half_lo); + const T h { x - T { 5, -1 } }; + num = h + h; + den = T { 2 } + x; + den_lo = x - (den - T { 2 }); + } + else if (x < T { 11875, -4 }) + { + c_hi = detail::split_pi_values(detail::split_pi_detail::quarter_pi_hi); + c_lo = detail::split_pi_values(detail::split_pi_detail::quarter_pi_lo); + num = x - T { 1 }; + den = x + T { 1 }; + den_lo = x - (den - T { 1 }); } else { - result = my_pi_half - detail::atan_series(one / x); + // 3 x and 2 + 3 x are exact here, since x >= 1 and 2 + 3 x < 10 + c_hi = detail::atan_values(detail::atan_detail::atan_three_halves_hi); + c_lo = detail::atan_values(detail::atan_detail::atan_three_halves_lo); + const T h { x - T { 15, -1 } }; + num = h + h; + den = T { 2 } + T { 3 } * x; } + + // q + q_lo = num / (den + den_lo); for |q| < 0.1 the rounding of q changes the result by + // less than 0.05 ulp, so only pay for the fma of the remainder when |q| >= 0.1 + const T q { num / den }; + const bool q_on_grid { fabs(q) >= T { 1, -1 } }; + const T q_lo { (q_on_grid ? detail::unchecked_fma(-q, den, num) - q * den_lo : -q * den_lo) / den }; + const T lo { detail::atan_tail(q) + ((has_lo ? c_lo + x_lo / (T { 1 } + x * x) : c_lo) + q_lo) }; + + // c_hi + q is exact only if q is on the grid of the result and the sum is below 1; + // else add q and the low part first, so only one rounding happens at the result grid + const T hi_q { c_hi + q }; + result = (q_on_grid && hi_q < T { 1 }) ? hi_q + lo : c_hi + (q + lo); + } + else + { + // atan(x) = pi/2 - atan(1 / x), with 1 / x <= 16/39 in the kernel's range + // with has_lo, the rounding of r = -1 / x and x_lo r^2 both change atan(r) by 1 / (1 + r^2) + const T r { T { -1 } / x }; + const T pi_lo { detail::split_pi_values(detail::split_pi_detail::half_pi_lo) }; + result = detail::split_pi_values(detail::split_pi_detail::half_pi_hi) + (has_lo ? + r + (detail::atan_tail(r) + (pi_lo + (detail::unchecked_fma(-r, x, T { -1 }) / x + x_lo * (r * r)) / (T { 1 } + r * r))) : + detail::atan_series(r) + pi_lo); } return result; diff --git a/include/boost/decimal/detail/cmath/atan2.hpp b/include/boost/decimal/detail/cmath/atan2.hpp index d1de750d7..e5e882d0e 100644 --- a/include/boost/decimal/detail/cmath/atan2.hpp +++ b/include/boost/decimal/detail/cmath/atan2.hpp @@ -9,6 +9,8 @@ #include #include #include +#include +#include #include #include #include @@ -24,20 +26,6 @@ namespace decimal { namespace detail { -namespace atan2_detail { - -template -struct pi_constants -{ - static constexpr T pi_over_two = numbers::pi_v / 2; - static constexpr T three_pi_over_four = 3 * numbers::pi_over_four_v; -}; - -template constexpr T pi_constants::pi_over_two; -template constexpr T pi_constants::three_pi_over_four; - -} // namespace atan2_detail - template constexpr auto atan2_impl(const T y, const T x) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) @@ -73,13 +61,13 @@ constexpr auto atan2_impl(const T y, const T x) noexcept #ifndef BOOST_DECIMAL_FAST_MATH else if (fpcy == FP_INFINITE && isfinitex) { - result = atan2_detail::pi_constants::pi_over_two; + result = split_pi_values(split_pi_detail::half_pi_hi); if (signy) { result = -result; } } else if (fpcy == FP_INFINITE && fpcx == FP_INFINITE && signx) { - result = atan2_detail::pi_constants::three_pi_over_four; + result = split_pi_values(split_pi_detail::three_quarter_pi); if (signy) { @@ -88,7 +76,7 @@ constexpr auto atan2_impl(const T y, const T x) noexcept } else if (fpcy == FP_INFINITE && fpcx == FP_INFINITE && !signx) { - result = numbers::pi_over_four_v; + result = split_pi_values(split_pi_detail::quarter_pi_hi); if (signy) { @@ -98,7 +86,7 @@ constexpr auto atan2_impl(const T y, const T x) noexcept #endif else if (fpcx == FP_ZERO) { - result = atan2_detail::pi_constants::pi_over_two; + result = split_pi_values(split_pi_detail::half_pi_hi); if (signy) { result = -result; } } @@ -122,7 +110,25 @@ constexpr auto atan2_impl(const T y, const T x) noexcept } else { - const auto ret_val {atan(fabs(y / x))}; + // For q in [10^k, tan(10^k)), q has a grid ten times coarser than atan(q); only there + // pass the remainder of y / x to atan as a low part (tan(1) < 1.5575) + const T ax {fabs(x)}; + const T q {fabs(y) / ax}; + T ret_val {}; + if (q < T {1}) + { + ret_val = atan(q); + int q_exp {}; + frexp10(q, &q_exp); + if (fpclassify(q) == FP_NORMAL && ret_val < T {1, q_exp + detail::precision_v - 1}) + { + ret_val = atan_impl(q, detail::unchecked_fma(-q, ax, fabs(y)) / ax); + } + } + else + { + ret_val = q < T {15575, -4} ? atan_impl(q, detail::unchecked_fma(-q, ax, fabs(y)) / ax) : atan(q); + } if (!signy && !signx) { diff --git a/include/boost/decimal/detail/cmath/impl/asin_impl.hpp b/include/boost/decimal/detail/cmath/impl/asin_impl.hpp index 8d120442c..5b850dc69 100644 --- a/include/boost/decimal/detail/cmath/impl/asin_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/asin_impl.hpp @@ -8,10 +8,13 @@ #include #include #include "../../int128.hpp" -#include -#include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include #include #endif @@ -22,214 +25,79 @@ namespace detail { namespace asin_detail { // Use a struct of arrays so that we can have static constexpr arrays of coefficients. -// If we don't do this we will end up constructing the values in the arrays each time the function is called -// which leads to massive slowdowns. -// // See https://github.com/boostorg/math/issues/923 for further information template struct asin_table_imp { -private: - using d32_coeffs_t = std::array; - using d64_coeffs_t = std::array; - using d128_coeffs_t = std::array; - - using d32_fast_coeffs_t = std::array; - using d64_fast_coeffs_t = std::array; - using d128_fast_coeffs_t = std::array; - -public: - - // 10th degree remez polynomial calculated from 0, 0.5 - // Estimated max error: 7.3651618860008751e-11 - static constexpr d32_coeffs_t d32_coeffs = {{ - decimal32_t {UINT64_C(263887099755925), -15}, - decimal32_t {UINT64_C(43491393212832818), -17, construction_sign::negative}, - decimal32_t {UINT64_C(38559884786102105), -17}, - decimal32_t {UINT64_C(13977130653211101), -17, construction_sign::negative}, - decimal32_t {UINT64_C(54573213517731915), -18}, - decimal32_t {UINT64_C(64851743877986187), -18}, - decimal32_t {UINT64_C(11606701725692841), -19}, - decimal32_t {UINT64_C(16658989049586517), -17}, - decimal32_t {UINT64_C(25906093603686159), -22}, - decimal32_t {UINT64_C(99999996600828589), -17}, - decimal32_t {UINT64_C(73651618860008751), -27} - }}; - - static constexpr d32_fast_coeffs_t d32_fast_coeffs = {{ - decimal_fast32_t {UINT64_C(263887099755925), -15}, - decimal_fast32_t {UINT64_C(43491393212832818), -17, construction_sign::negative}, - decimal_fast32_t {UINT64_C(38559884786102105), -17}, - decimal_fast32_t {UINT64_C(13977130653211101), -17, construction_sign::negative}, - decimal_fast32_t {UINT64_C(54573213517731915), -18}, - decimal_fast32_t {UINT64_C(64851743877986187), -18}, - decimal_fast32_t {UINT64_C(11606701725692841), -19}, - decimal_fast32_t {UINT64_C(16658989049586517), -17}, - decimal_fast32_t {UINT64_C(25906093603686159), -22}, - decimal_fast32_t {UINT64_C(99999996600828589), -17}, - decimal_fast32_t {UINT64_C(73651618860008751), -27} - }}; - - // 20th degree remez polynomial calculated from 0, 0.5 - // Estimated max error: 6.0872797932519911178133457751215133e-19 - static constexpr d64_coeffs_t d64_coeffs = {{ - decimal64_t {UINT64_C(2201841632531125594), -18}, - decimal64_t {UINT64_C(9319383818485265142), -18, construction_sign::negative}, - decimal64_t {UINT64_C(1876826158920611297), -17}, - decimal64_t {UINT64_C(2351630530022519158), -17, construction_sign::negative}, - decimal64_t {UINT64_C(2046603318375014621), -17}, - decimal64_t {UINT64_C(1304427904865204196), -17, construction_sign::negative}, - decimal64_t {UINT64_C(6308794339076719731), -18}, - decimal64_t {UINT64_C(2333806156857836980), -18, construction_sign::negative}, - decimal64_t {UINT64_C(6826985955727270693), -19}, - decimal64_t {UINT64_C(1326415745606167277), -19, construction_sign::negative}, - decimal64_t {UINT64_C(2747750823768175476), -20}, - decimal64_t {UINT64_C(2660509753516203115), -20}, - decimal64_t {UINT64_C(3977122944636320545), -22}, - decimal64_t {UINT64_C(4461135938842722307), -20}, - decimal64_t {UINT64_C(1826730778134521645), -24}, - decimal64_t {UINT64_C(7499992533825458566), -20}, - decimal64_t {UINT64_C(2034140780525051207), -27}, - decimal64_t {UINT64_C(1666666666327808185), -19}, - decimal64_t {UINT64_C(2987315928933390856), -31}, - decimal64_t {UINT64_C(9999999999999989542), -19}, - decimal64_t {UINT64_C(6087279793251991118), -37} - }}; - - static constexpr d64_fast_coeffs_t d64_fast_coeffs = {{ - decimal_fast64_t {UINT64_C(2201841632531125594), -18}, - decimal_fast64_t {UINT64_C(9319383818485265142), -18, construction_sign::negative}, - decimal_fast64_t {UINT64_C(1876826158920611297), -17}, - decimal_fast64_t {UINT64_C(2351630530022519158), -17, construction_sign::negative}, - decimal_fast64_t {UINT64_C(2046603318375014621), -17}, - decimal_fast64_t {UINT64_C(1304427904865204196), -17, construction_sign::negative}, - decimal_fast64_t {UINT64_C(6308794339076719731), -18}, - decimal_fast64_t {UINT64_C(2333806156857836980), -18, construction_sign::negative}, - decimal_fast64_t {UINT64_C(6826985955727270693), -19}, - decimal_fast64_t {UINT64_C(1326415745606167277), -19, construction_sign::negative}, - decimal_fast64_t {UINT64_C(2747750823768175476), -20}, - decimal_fast64_t {UINT64_C(2660509753516203115), -20}, - decimal_fast64_t {UINT64_C(3977122944636320545), -22}, - decimal_fast64_t {UINT64_C(4461135938842722307), -20}, - decimal_fast64_t {UINT64_C(1826730778134521645), -24}, - decimal_fast64_t {UINT64_C(7499992533825458566), -20}, - decimal_fast64_t {UINT64_C(2034140780525051207), -27}, - decimal_fast64_t {UINT64_C(1666666666327808185), -19}, - decimal_fast64_t {UINT64_C(2987315928933390856), -31}, - decimal_fast64_t {UINT64_C(9999999999999989542), -19}, - decimal_fast64_t {UINT64_C(6087279793251991118), -37} + // Chebyshev fits of P in asin(x) = x + x^3 P(x^2) for x in [0, 0.5], highest power first, as + // raw fixed point (see fixed_point_series.hpp): n = round(c * 2^62) or n = round(c * 2^126) + + // Q62, degree 5 in x^2, max kernel error 0.0479 ulp at 7 digits + static constexpr std::array d32_coeffs = {{ + INT64_C(155371621893934362), + INT64_C(79086910334655698), + INT64_C(143426478885863926), + INT64_C(205678429410250444), + INT64_C(345880832483771599), + INT64_C(768614490127431932) }}; - // 40th degree remez polynomial calculated from 0, 0.5 - // Estimated max error: 1.084502473818005718919720519483941e-34 - static constexpr d128_coeffs_t d128_coeffs = {{ - decimal128_t {int128::uint128_t{UINT64_C(236367828732266), UINT64_C(4865873281479238114)}, -31}, - decimal128_t {int128::uint128_t{UINT64_C(218966359248756), UINT64_C(1393338271545593644)}, -30, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(98104038983693), UINT64_C(4819646069944316372)}, -29}, - decimal128_t {int128::uint128_t{UINT64_C(282853615727310), UINT64_C(10104044375051504970)}, -29, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(58930987436658), UINT64_C(3829646337759276014)}, -28}, - decimal128_t {int128::uint128_t{UINT64_C(94467942291578), UINT64_C(14212526794757587650)}, -28, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(121156109355190), UINT64_C(6171523396929956760)}, -28}, - decimal128_t {int128::uint128_t{UINT64_C(127640043209581), UINT64_C(8369619306382995314)}, -28, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(112556984011870), UINT64_C(14401172681696800280)}, -28}, - decimal128_t {int128::uint128_t{UINT64_C(84240716950351), UINT64_C(10152945328926072964)}, -28, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(540724366020485), UINT64_C(8813105586620168570)}, -29}, - decimal128_t {int128::uint128_t{UINT64_C(300054630162323), UINT64_C(4862687399308912842)}, -29, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(144827005285082), UINT64_C(4790810090757542758)}, -29}, - decimal128_t {int128::uint128_t{UINT64_C(61085784025333), UINT64_C(3908625641731373429)}, -29, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(225929173229512), UINT64_C(18404095637827467688)}, -30}, - decimal128_t {int128::uint128_t{UINT64_C(73452862511516), UINT64_C(2655967943189644664)}, -30, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(210254502661653), UINT64_C(14174199201997297032)}, -31}, - decimal128_t {int128::uint128_t{UINT64_C(530269670900176), UINT64_C(3023877239296322874)}, -32, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(117870705400334), UINT64_C(8785618254907029456)}, -32}, - decimal128_t {int128::uint128_t{UINT64_C(230285265351731), UINT64_C(8107756519153341434)}, -33, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(397318429350031), UINT64_C(567549410172969484)}, -34}, - decimal128_t {int128::uint128_t{UINT64_C(54772616787306), UINT64_C(4168475956004989379)}, -34, construction_sign::negative}, - decimal128_t {int128::uint128_t{UINT64_C(79509164538790), UINT64_C(17928590725399689320)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(534376054761824), UINT64_C(1987644731805023176)}, -36}, - decimal128_t {int128::uint128_t{UINT64_C(92204817966183), UINT64_C(17576450582561384882)}, -37}, - decimal128_t {int128::uint128_t{UINT64_C(75623542590285), UINT64_C(990523592779300020)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(59680570668825), UINT64_C(14870623164911255928)}, -39}, - decimal128_t {int128::uint128_t{UINT64_C(94069144841714), UINT64_C(11353995396932754836)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(204081757431333), UINT64_C(1300964680833664202)}, -42}, - decimal128_t {int128::uint128_t{UINT64_C(121279716530202), UINT64_C(3054546075061258708)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(340541736068294), UINT64_C(674620373211314186)}, -45}, - decimal128_t {int128::uint128_t{UINT64_C(164700850853976), UINT64_C(1203142186405381614)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(246590930469756), UINT64_C(6088477928552847004)}, -48}, - decimal128_t {int128::uint128_t{UINT64_C(242009413501228), UINT64_C(3841246034215456962)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(64561634810301), UINT64_C(5259904364587721972)}, -51}, - decimal128_t {int128::uint128_t{UINT64_C(406575814682064), UINT64_C(3001055340328133406)}, -35}, - decimal128_t {int128::uint128_t{UINT64_C(447242814330412), UINT64_C(4234427805033793948)}, -56}, - decimal128_t {int128::uint128_t{UINT64_C(90350181040458), UINT64_C(12964998079628443792)}, -34}, - decimal128_t {int128::uint128_t{UINT64_C(430604756670586), UINT64_C(9888097447655546704)}, -61}, - decimal128_t {int128::uint128_t{UINT64_C(542101086242752), UINT64_C(4003012203950105568)}, -34}, - decimal128_t {int128::uint128_t{UINT64_C(58790996908969), UINT64_C(5250765973560640036)}, -67} + // Q62, degree 12 in x^2, max kernel error 0.0452 ulp at 16 digits + static constexpr std::array d64_coeffs = {{ + INT64_C(132622181071150947), + INT64_C(-68492239953733199), + INT64_C(80247392434212570), + INT64_C(25168307429827271), + INT64_C(47605578609573363), + INT64_C(52938361988655692), + INT64_C(64430847530605378), + INT64_C(80023786897097220), + INT64_C(103173437159152167), + INT64_C(140111986996306496), + INT64_C(205878840124498533), + INT64_C(345876451381981828), + INT64_C(768614336404564804) }}; - static constexpr d128_fast_coeffs_t d128_fast_coeffs = {{ - decimal_fast128_t {int128::uint128_t{UINT64_C(236367828732266), UINT64_C(4865873281479238114)}, -31}, - decimal_fast128_t {int128::uint128_t{UINT64_C(218966359248756), UINT64_C(1393338271545593644)}, -30, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(98104038983693), UINT64_C(4819646069944316372)}, -29}, - decimal_fast128_t {int128::uint128_t{UINT64_C(282853615727310), UINT64_C(10104044375051504970)}, -29, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(58930987436658), UINT64_C(3829646337759276014)}, -28}, - decimal_fast128_t {int128::uint128_t{UINT64_C(94467942291578), UINT64_C(14212526794757587650)}, -28, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(121156109355190), UINT64_C(6171523396929956760)}, -28}, - decimal_fast128_t {int128::uint128_t{UINT64_C(127640043209581), UINT64_C(8369619306382995314)}, -28, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(112556984011870), UINT64_C(14401172681696800280)}, -28}, - decimal_fast128_t {int128::uint128_t{UINT64_C(84240716950351), UINT64_C(10152945328926072964)}, -28, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(540724366020485), UINT64_C(8813105586620168570)}, -29}, - decimal_fast128_t {int128::uint128_t{UINT64_C(300054630162323), UINT64_C(4862687399308912842)}, -29, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(144827005285082), UINT64_C(4790810090757542758)}, -29}, - decimal_fast128_t {int128::uint128_t{UINT64_C(61085784025333), UINT64_C(3908625641731373429)}, -29, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(225929173229512), UINT64_C(18404095637827467688)}, -30}, - decimal_fast128_t {int128::uint128_t{UINT64_C(73452862511516), UINT64_C(2655967943189644664)}, -30, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(210254502661653), UINT64_C(14174199201997297032)}, -31}, - decimal_fast128_t {int128::uint128_t{UINT64_C(530269670900176), UINT64_C(3023877239296322874)}, -32, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(117870705400334), UINT64_C(8785618254907029456)}, -32}, - decimal_fast128_t {int128::uint128_t{UINT64_C(230285265351731), UINT64_C(8107756519153341434)}, -33, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(397318429350031), UINT64_C(567549410172969484)}, -34}, - decimal_fast128_t {int128::uint128_t{UINT64_C(54772616787306), UINT64_C(4168475956004989379)}, -34, construction_sign::negative}, - decimal_fast128_t {int128::uint128_t{UINT64_C(79509164538790), UINT64_C(17928590725399689320)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(534376054761824), UINT64_C(1987644731805023176)}, -36}, - decimal_fast128_t {int128::uint128_t{UINT64_C(92204817966183), UINT64_C(17576450582561384882)}, -37}, - decimal_fast128_t {int128::uint128_t{UINT64_C(75623542590285), UINT64_C(990523592779300020)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(59680570668825), UINT64_C(14870623164911255928)}, -39}, - decimal_fast128_t {int128::uint128_t{UINT64_C(94069144841714), UINT64_C(11353995396932754836)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(204081757431333), UINT64_C(1300964680833664202)}, -42}, - decimal_fast128_t {int128::uint128_t{UINT64_C(121279716530202), UINT64_C(3054546075061258708)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(340541736068294), UINT64_C(674620373211314186)}, -45}, - decimal_fast128_t {int128::uint128_t{UINT64_C(164700850853976), UINT64_C(1203142186405381614)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(246590930469756), UINT64_C(6088477928552847004)}, -48}, - decimal_fast128_t {int128::uint128_t{UINT64_C(242009413501228), UINT64_C(3841246034215456962)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(64561634810301), UINT64_C(5259904364587721972)}, -51}, - decimal_fast128_t {int128::uint128_t{UINT64_C(406575814682064), UINT64_C(3001055340328133406)}, -35}, - decimal_fast128_t {int128::uint128_t{UINT64_C(447242814330412), UINT64_C(4234427805033793948)}, -56}, - decimal_fast128_t {int128::uint128_t{UINT64_C(90350181040458), UINT64_C(12964998079628443792)}, -34}, - decimal_fast128_t {int128::uint128_t{UINT64_C(430604756670586), UINT64_C(9888097447655546704)}, -61}, - decimal_fast128_t {int128::uint128_t{UINT64_C(542101086242752), UINT64_C(4003012203950105568)}, -34}, - decimal_fast128_t {int128::uint128_t{UINT64_C(58790996908969), UINT64_C(5250765973560640036)}, -67} + // Q126, degree 28 in x^2, max kernel error 0.043 ulp at 34 digits + static constexpr std::array d128_coeffs = {{ + int128::int128_t {INT64_C(370876500756886001), UINT64_C(9221949471032297446)}, + int128::int128_t {INT64_C(-956307695304863090), UINT64_C(9148781051210321482)}, + int128::int128_t {INT64_C(1310549287927847162), UINT64_C(1958114557122675838)}, + int128::int128_t {INT64_C(-1134417282113220194), UINT64_C(17968718786407945997)}, + int128::int128_t {INT64_C(736228940689646156), UINT64_C(16663177741084984228)}, + int128::int128_t {INT64_C(-340829970488362816), UINT64_C(12459077857655578150)}, + int128::int128_t {INT64_C(145857634493384410), UINT64_C(17283691294030731877)}, + int128::int128_t {INT64_C(-29177674511655379), UINT64_C(5106089515100561328)}, + int128::int128_t {INT64_C(23589553566972586), UINT64_C(4940683816029200174)}, + int128::int128_t {INT64_C(11915337754018893), UINT64_C(6093997656988561099)}, + int128::int128_t {INT64_C(15585523597377591), UINT64_C(13114903993071182046)}, + int128::int128_t {INT64_C(16404669177712143), UINT64_C(5788416753019317356)}, + int128::int128_t {INT64_C(17904543839495541), UINT64_C(305658806213704697)}, + int128::int128_t {INT64_C(19557042090611134), UINT64_C(18251937558368413649)}, + int128::int128_t {INT64_C(21491177464483278), UINT64_C(13762596495847561352)}, + int128::int128_t {INT64_C(23765442023041221), UINT64_C(4709006491505993448)}, + int128::int128_t {INT64_C(26471251718294205), UINT64_C(7048697116100673247)}, + int128::int128_t {INT64_C(29732509641295605), UINT64_C(791966617446954024)}, + int128::int128_t {INT64_C(33723073331131660), UINT64_C(5494626099549260253)}, + int128::int128_t {INT64_C(38693594343105809), UINT64_C(14383093717251864696)}, + int128::int128_t {INT64_C(45017478183132254), UINT64_C(3729214876058472566)}, + int128::int128_t {INT64_C(53273278680384444), UINT64_C(15594621309481963446)}, + int128::int128_t {INT64_C(64401474671398092), UINT64_C(16567906768798270425)}, + int128::int128_t {INT64_C(80025501070968044), UINT64_C(5657558833186617601)}, + int128::int128_t {INT64_C(103173373281578635), UINT64_C(11738962649206475942)}, + int128::int128_t {INT64_C(140111988407082097), UINT64_C(14347467083897352460)}, + int128::int128_t {INT64_C(205878840108365531), UINT64_C(7905747461347776107)}, + int128::int128_t {INT64_C(345876451382054092), UINT64_C(14757395258966581313)}, + int128::int128_t {INT64_C(768614336404564650), UINT64_C(12297829382473037246)} }}; }; #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) -template -constexpr typename asin_table_imp::d32_coeffs_t asin_table_imp::d32_coeffs; - -template -constexpr typename asin_table_imp::d64_coeffs_t asin_table_imp::d64_coeffs; - -template -constexpr typename asin_table_imp::d128_coeffs_t asin_table_imp::d128_coeffs; - -template -constexpr typename asin_table_imp::d32_fast_coeffs_t asin_table_imp::d32_fast_coeffs; - -template -constexpr typename asin_table_imp::d64_fast_coeffs_t asin_table_imp::d64_fast_coeffs; - -template -constexpr typename asin_table_imp::d128_fast_coeffs_t asin_table_imp::d128_fast_coeffs; +template constexpr std::array asin_table_imp::d32_coeffs; +template constexpr std::array asin_table_imp::d64_coeffs; +template constexpr std::array asin_table_imp::d128_coeffs; #endif @@ -237,43 +105,43 @@ using asin_table = asin_table_imp; } //namespace asin_detail +// asin(x) - x = x^3 P(x^2) template -constexpr auto asin_series(T x) noexcept; - -template <> -constexpr auto asin_series(decimal32_t x) noexcept -{ - return remez_series_result(x, asin_detail::asin_table::d32_coeffs); -} - -template <> -constexpr auto asin_series(decimal_fast32_t x) noexcept -{ - return remez_series_result(x, asin_detail::asin_table::d32_fast_coeffs); -} - -template <> -constexpr auto asin_series(decimal64_t x) noexcept -{ - return remez_series_result(x, asin_detail::asin_table::d64_coeffs); -} +constexpr auto asin_tail(T x) noexcept -> T; -template <> -constexpr auto asin_series(decimal_fast64_t x) noexcept -{ - return remez_series_result(x, asin_detail::asin_table::d64_fast_coeffs); -} +template <> constexpr auto asin_tail(decimal32_t x) noexcept -> decimal32_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d32_coeffs); } +template <> constexpr auto asin_tail(decimal_fast32_t x) noexcept -> decimal_fast32_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d32_coeffs); } +template <> constexpr auto asin_tail(decimal64_t x) noexcept -> decimal64_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d64_coeffs); } +template <> constexpr auto asin_tail(decimal_fast64_t x) noexcept -> decimal_fast64_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d64_coeffs); } +template <> constexpr auto asin_tail(decimal128_t x) noexcept -> decimal128_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d128_coeffs); } +template <> constexpr auto asin_tail(decimal_fast128_t x) noexcept -> decimal_fast128_t { return fixed_point_odd_tail(x, asin_detail::asin_table::d128_coeffs); } -template <> -constexpr auto asin_series(decimal128_t x) noexcept +template +constexpr auto asin_series(T x) noexcept -> T { - return remez_series_result(x, asin_detail::asin_table::d128_coeffs); + return x + asin_tail(x); } -template <> -constexpr auto asin_series(decimal_fast128_t x) noexcept +// Half-angle step for a in (0.5, 1]: s = sqrt((1 - a) / 2), and asin(s) ~= s + the returned low +// part, which also covers the rounding error of s; asin(a) = pi/2 - 2 asin(s), acos(a) = 2 asin(s) +template +constexpr auto asin_half_angle(T a, T& s) noexcept -> T { - return remez_series_result(x, asin_detail::asin_table::d128_fast_coeffs); + const T w {1 - a}; + const T z {w / 2}; + s = sqrt(z); + + if (s == 0) + { + return s; + } + + // sqrt(w / 2) - s = (w - 2 s^2) / (4 s), times asin'(s) = 1 / sqrt(1 - z) ~= 1 / (1 - z / 2); + // s + s rounds for s in [0.05, 0.1), and two_s_lo is the exact part it rounds off + const T two_s {s + s}; + const T two_s_lo {s - (two_s - s)}; + const T c {(detail::unchecked_fma(-two_s, s, w) - two_s_lo * s) / (two_s * (T {2} - z))}; + return asin_tail(s) + c; } } //namespace detail diff --git a/include/boost/decimal/detail/cmath/impl/atan_impl.hpp b/include/boost/decimal/detail/cmath/impl/atan_impl.hpp index 2b31b426b..0e96111ce 100644 --- a/include/boost/decimal/detail/cmath/impl/atan_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/atan_impl.hpp @@ -6,9 +6,11 @@ #ifndef BOOST_DECIMAL_DETAIL_CMATH_IMPL_ATAN_IMPL_HPP #define BOOST_DECIMAL_DETAIL_CMATH_IMPL_ATAN_IMPL_HPP +#include #include -#include -#include +#include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -22,95 +24,131 @@ namespace detail { namespace atan_detail { +// Indices into atan_values +enum : std::size_t { atan_half_hi, atan_half_lo, atan_three_halves_hi, atan_three_halves_lo }; + +// Use a struct of arrays so that we can have static constexpr arrays of coefficients. +// See https://github.com/boostorg/math/issues/923 for further information template struct atan_table_imp { - static constexpr std::array<::boost::decimal::decimal32_t, 3> d32_atan_values = - {{ - ::boost::decimal::decimal32_t { UINT64_C(4636476090008061162), -19 }, // atan_half - ::boost::decimal::decimal32_t { UINT64_C(7853981633974483096), -19 }, // atan_one - ::boost::decimal::decimal32_t { UINT64_C(9827937232473290679), -19 }, // atan_three_halves + // Chebyshev fits of P in atan(x) = x + x^3 P(x^2) for x in [0, 0.4375], highest power first, as + // raw fixed point (see fixed_point_series.hpp): n = round(c * 2^62) or n = round(c * 2^126) + + // Q62, degree 4 in x^2, max kernel error 0.0185 ulp at 7 digits + static constexpr std::array d32_coeffs = {{ + INT64_C(-286653686818871604), + INT64_C(491992650021305608), + INT64_C(-657468245397313456), + INT64_C(922305383051950432), + INT64_C(-1537228519086262020) }}; - static constexpr std::array<::boost::decimal::decimal_fast32_t, 3> d32_fast_atan_values = - {{ - ::boost::decimal::decimal_fast32_t { UINT64_C(4636476090008061162), -19 }, // atan_half - ::boost::decimal::decimal_fast32_t { UINT64_C(7853981633974483096), -19 }, // atan_one - ::boost::decimal::decimal_fast32_t { UINT64_C(9827937232473290679), -19 }, // atan_three_halves - }}; + // Q62, degree 11 in x^2, max kernel error 0.027 ulp at 16 digits + static constexpr std::array d64_coeffs = {{ + INT64_C(68129288945790802), + INT64_C(-152776553565301478), + INT64_C(207168893735716832), + INT64_C(-240576890164062181), + INT64_C(271026318056636923), + INT64_C(-307426146734136447), + INT64_C(354744057432338816), + INT64_C(-419244149583371944), + INT64_C(512409556935913322), + INT64_C(-658812288339968584), + INT64_C(922337203685450372), + INT64_C(-1537228672809129148) + }}; - static constexpr std::array<::boost::decimal::decimal64_t, 3> d64_atan_values = - {{ - ::boost::decimal::decimal64_t { UINT64_C(4636476090008061162), -19 }, // atan_half - ::boost::decimal::decimal64_t { UINT64_C(7853981633974483096), -19 }, // atan_one - ::boost::decimal::decimal64_t { UINT64_C(9827937232473290679), -19 }, // atan_three_halves + // Q126, degree 24 in x^2, max kernel error 0.0314 ulp at 34 digits + static constexpr std::array d128_coeffs = {{ + int128::int128_t {INT64_C(-10485371323887896), UINT64_C(10517237177405548740)}, + int128::int128_t {INT64_C(36037310788783410), UINT64_C(17838482110476103184)}, + int128::int128_t {INT64_C(-65868101180795099), UINT64_C(11297532328258753760)}, + int128::int128_t {INT64_C(88709017543273301), UINT64_C(10070282870844457804)}, + int128::int128_t {INT64_C(-102651823325949239), UINT64_C(12908569194395175164)}, + int128::int128_t {INT64_C(111261874313837669), UINT64_C(1207799242924094493)}, + int128::int128_t {INT64_C(-117988633709135858), UINT64_C(11530285581856098480)}, + int128::int128_t {INT64_C(124595196758679357), UINT64_C(4918781308014821961)}, + int128::int128_t {INT64_C(-131756095136795816), UINT64_C(4568656608273312686)}, + int128::int128_t {INT64_C(139747322689771072), UINT64_C(2293705908972136111)}, + int128::int128_t {INT64_C(-148763994724226424), UINT64_C(13130196723226486758)}, + int128::int128_t {INT64_C(159023650305606894), UINT64_C(16383032007635900459)}, + int128::int128_t {INT64_C(-170803185516202007), UINT64_C(10439900632158261928)}, + int128::int128_t {INT64_C(184467440718861643), UINT64_C(10470388862571930467)}, + int128::int128_t {INT64_C(-200508087756951228), UINT64_C(6100221296306786481)}, + int128::int128_t {INT64_C(219604096115564635), UINT64_C(2688867774528155709)}, + int128::int128_t {INT64_C(-242720316759335550), UINT64_C(391790395505276335)}, + int128::int128_t {INT64_C(271275648142787510), UINT64_C(14028508071849608541)}, + int128::int128_t {INT64_C(-307445734561825861), UINT64_C(17059997906172122564)}, + int128::int128_t {INT64_C(354745078340568300), UINT64_C(5638858341042803538)}, + int128::int128_t {INT64_C(-419244183493398901), UINT64_C(11739098548382091930)}, + int128::int128_t {INT64_C(512409557603043100), UINT64_C(8198551786701958321)}, + int128::int128_t {INT64_C(-658812288346769701), UINT64_C(7905747462781355944)}, + int128::int128_t {INT64_C(922337203685477580), UINT64_C(14757395258965233795)}, + int128::int128_t {INT64_C(-1537228672809129302), UINT64_C(12297829382473037246)} }}; - static constexpr std::array<::boost::decimal::decimal_fast64_t, 3> d64_fast_atan_values = - {{ - ::boost::decimal::decimal_fast64_t { UINT64_C(4636476090008061162), -19 }, // atan_half - ::boost::decimal::decimal_fast64_t { UINT64_C(7853981633974483096), -19 }, // atan_one - ::boost::decimal::decimal_fast64_t { UINT64_C(9827937232473290679), -19 }, // atan_three_halves - }}; + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d32_values = {{ + decimal32_t {UINT64_C(4636476), -7}, + decimal32_t {UINT64_C(9000806), -15}, + decimal32_t {UINT64_C(9827937), -7}, + decimal32_t {UINT64_C(2324733), -14} + }}; + + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d32_fast_values = {{ + decimal_fast32_t {UINT64_C(4636476), -7}, + decimal_fast32_t {UINT64_C(9000806), -15}, + decimal_fast32_t {UINT64_C(9827937), -7}, + decimal_fast32_t {UINT64_C(2324733), -14} + }}; - static constexpr std::array<::boost::decimal::decimal128_t, 3> d128_atan_values = - {{ - ::boost::decimal::decimal128_t { boost::int128::uint128_t { UINT64_C(251343872473191), UINT64_C(15780610568723885484) }, -34 }, // atan_half - ::boost::decimal::decimal128_t { boost::int128::uint128_t { UINT64_C(425765197510819), UINT64_C(5970600460659265246) }, -34 }, // atan_one - ::boost::decimal::decimal128_t { boost::int128::uint128_t { UINT64_C(532773544924935), UINT64_C(16408933314882201700) }, -34 }, // atan_three_halves + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d64_values = {{ + decimal64_t {UINT64_C(4636476090008061), -16}, + decimal64_t {UINT64_C(1621425623146121), -32}, + decimal64_t {UINT64_C(9827937232473291), -16}, + decimal64_t {UINT64_C(3201428938898533), -32, construction_sign::negative} }}; - static constexpr std::array<::boost::decimal::decimal_fast128_t, 3> d128_fast_atan_values = - {{ - ::boost::decimal::decimal_fast128_t { boost::int128::uint128_t { UINT64_C(251343872473191), UINT64_C(15780610568723885484) }, -34 }, // atan_half - ::boost::decimal::decimal_fast128_t { boost::int128::uint128_t { UINT64_C(425765197510819), UINT64_C(5970600460659265246) }, -34 }, // atan_one - ::boost::decimal::decimal_fast128_t { boost::int128::uint128_t { UINT64_C(532773544924935), UINT64_C(16408933314882201700) }, -34 }, // atan_three_halves + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d64_fast_values = {{ + decimal_fast64_t {UINT64_C(4636476090008061), -16}, + decimal_fast64_t {UINT64_C(1621425623146121), -32}, + decimal_fast64_t {UINT64_C(9827937232473291), -16}, + decimal_fast64_t {UINT64_C(3201428938898533), -32, construction_sign::negative} }}; - // 10th degree remez polynomial calculated from 0, 0.4375 - // Estimated max error: 2.3032664387910605e-12 - static constexpr std::array<::boost::decimal::decimal32_t, 11> d32_coeffs = - {{ - ::boost::decimal::decimal32_t { UINT64_C(61037779951304161), -18, true }, - ::boost::decimal::decimal32_t { UINT64_C(10723099589331457), -17 }, - ::boost::decimal::decimal32_t { UINT64_C(22515613909953665), -18 }, - ::boost::decimal::decimal32_t { UINT64_C(15540713402718176), -17, true }, - ::boost::decimal::decimal32_t { UINT64_C(35999727706986597), -19 }, - ::boost::decimal::decimal32_t { UINT64_C(19938867353282852), -17 }, - ::boost::decimal::decimal32_t { UINT64_C(62252075283915644), -22 }, - ::boost::decimal::decimal32_t { UINT64_C(33333695504913247), -17, true }, - ::boost::decimal::decimal32_t { UINT64_C(10680927642397763), -24 }, - ::boost::decimal::decimal32_t { UINT64_C(99999999877886492), -17 }, - ::boost::decimal::decimal32_t { UINT64_C(23032664387910606), -29 }, + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d128_values = {{ + decimal128_t {int128::uint128_t{UINT64_C(251343872473191), UINT64_C(15780610568723885488)}, -34}, + decimal128_t {int128::uint128_t{UINT64_C(109967214061217), UINT64_C(15895449937643443526)}, -69}, + decimal128_t {int128::uint128_t{UINT64_C(53277354492493), UINT64_C(10864265368342995978)}, -33}, + decimal128_t {int128::uint128_t{UINT64_C(78587730147417), UINT64_C(12843742829858255927)}, -68} }}; - static constexpr std::array<::boost::decimal::decimal_fast32_t, 11> d32_fast_coeffs = - {{ - ::boost::decimal::decimal_fast32_t { UINT64_C(61037779951304161), -18, true }, - ::boost::decimal::decimal_fast32_t { UINT64_C(10723099589331457), -17 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(22515613909953665), -18 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(15540713402718176), -17, true }, - ::boost::decimal::decimal_fast32_t { UINT64_C(35999727706986597), -19 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(19938867353282852), -17 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(62252075283915644), -22 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(33333695504913247), -17, true }, - ::boost::decimal::decimal_fast32_t { UINT64_C(10680927642397763), -24 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(99999999877886492), -17 }, - ::boost::decimal::decimal_fast32_t { UINT64_C(23032664387910606), -29 }, - }}; + // atan(1/2) and atan(3/2), each as a rounded high part and the low part left over + static constexpr std::array d128_fast_values = {{ + decimal_fast128_t {int128::uint128_t{UINT64_C(251343872473191), UINT64_C(15780610568723885488)}, -34}, + decimal_fast128_t {int128::uint128_t{UINT64_C(109967214061217), UINT64_C(15895449937643443526)}, -69}, + decimal_fast128_t {int128::uint128_t{UINT64_C(53277354492493), UINT64_C(10864265368342995978)}, -33}, + decimal_fast128_t {int128::uint128_t{UINT64_C(78587730147417), UINT64_C(12843742829858255927)}, -68} + }}; }; #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) -template constexpr std::array atan_table_imp::d32_coeffs; -template constexpr std::array atan_table_imp::d32_fast_coeffs; - -template constexpr std::array atan_table_imp::d32_atan_values; -template constexpr std::array atan_table_imp::d32_fast_atan_values; -template constexpr std::array atan_table_imp::d64_atan_values; -template constexpr std::array atan_table_imp::d64_fast_atan_values; -template constexpr std::array atan_table_imp::d128_atan_values; -template constexpr std::array atan_table_imp::d128_fast_atan_values; +template constexpr std::array atan_table_imp::d32_coeffs; +template constexpr std::array atan_table_imp::d64_coeffs; +template constexpr std::array atan_table_imp::d128_coeffs; +template constexpr std::array atan_table_imp::d32_values; +template constexpr std::array atan_table_imp::d32_fast_values; +template constexpr std::array atan_table_imp::d64_values; +template constexpr std::array atan_table_imp::d64_fast_values; +template constexpr std::array atan_table_imp::d128_values; +template constexpr std::array atan_table_imp::d128_fast_values; #endif @@ -118,187 +156,35 @@ using atan_table = atan_table_imp; } //namespace atan_detail +// atan(x) - x = x^3 P(x^2) template -constexpr auto atan_series(T x) noexcept; +constexpr auto atan_tail(T x) noexcept -> T; + +template <> constexpr auto atan_tail(decimal32_t x) noexcept -> decimal32_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d32_coeffs); } +template <> constexpr auto atan_tail(decimal_fast32_t x) noexcept -> decimal_fast32_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d32_coeffs); } +template <> constexpr auto atan_tail(decimal64_t x) noexcept -> decimal64_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d64_coeffs); } +template <> constexpr auto atan_tail(decimal_fast64_t x) noexcept -> decimal_fast64_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d64_coeffs); } +template <> constexpr auto atan_tail(decimal128_t x) noexcept -> decimal128_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d128_coeffs); } +template <> constexpr auto atan_tail(decimal_fast128_t x) noexcept -> decimal_fast128_t { return fixed_point_odd_tail(x, atan_detail::atan_table::d128_coeffs); } template constexpr auto atan_values(std::size_t idx) noexcept -> T; -template <> constexpr auto atan_series (decimal32_t x) noexcept { return remez_series_result(x, atan_detail::atan_table::d32_coeffs); } - -template <> constexpr auto atan_series (decimal_fast32_t x) noexcept { return remez_series_result(x, atan_detail::atan_table::d32_fast_coeffs); } - -template <> -constexpr auto atan_series(decimal64_t x) noexcept -{ - // PadeApproximant[ArcTan[x]/x, {x, 0, {12, 12}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - const decimal64_t x2 { x * x }; - - const decimal64_t - top - { - decimal64_t { UINT64_C( 58561878375) } - + x2 * (decimal64_t { UINT64_C(163192434405) } - + x2 * (decimal64_t { UINT64_C(169269290190) } - + x2 * (decimal64_t { UINT64_C(80191217106) } - + x2 * (decimal64_t { UINT64_C(16979477515) } - + x2 * (decimal64_t { UINT32_C(1296036105) } - + x2 * decimal64_t { UINT32_C(15728640) }))))) - }; - - const decimal64_t - bot - { - decimal64_t { UINT64_C(58561878375) } - + x2 * (decimal64_t { UINT64_C(182713060530) } - + x2 * (decimal64_t { UINT64_C(218461268025) } - + x2 * (decimal64_t { UINT64_C(124835010300) } - + x2 * (decimal64_t { UINT64_C(34493884425) } - + x2 * (decimal64_t { UINT32_C(4058104050) } - + x2 * decimal64_t { UINT32_C(135270135) }))))) - }; - - return (x * top) / bot; -} - -template <> -constexpr auto atan_series(decimal_fast64_t x) noexcept -{ - // PadeApproximant[ArcTan[x]/x, {x, 0, {12, 12}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - const decimal_fast64_t x2 { x * x }; +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal32_t { return atan_detail::atan_table::d32_values[idx]; } +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal_fast32_t { return atan_detail::atan_table::d32_fast_values[idx]; } +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal64_t { return atan_detail::atan_table::d64_values[idx]; } +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal_fast64_t { return atan_detail::atan_table::d64_fast_values[idx]; } +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal128_t { return atan_detail::atan_table::d128_values[idx]; } +template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal_fast128_t { return atan_detail::atan_table::d128_fast_values[idx]; } - const decimal_fast64_t - top - { - decimal_fast64_t { UINT64_C( 58561878375) } - + x2 * (decimal_fast64_t { UINT64_C(163192434405) } - + x2 * (decimal_fast64_t { UINT64_C(169269290190) } - + x2 * (decimal_fast64_t { UINT64_C(80191217106) } - + x2 * (decimal_fast64_t { UINT64_C(16979477515) } - + x2 * (decimal_fast64_t { UINT32_C(1296036105) } - + x2 * decimal_fast64_t { UINT32_C(15728640) }))))) - }; - - const decimal_fast64_t - bot - { - decimal_fast64_t { UINT64_C(58561878375) } - + x2 * (decimal_fast64_t { UINT64_C(182713060530) } - + x2 * (decimal_fast64_t { UINT64_C(218461268025) } - + x2 * (decimal_fast64_t { UINT64_C(124835010300) } - + x2 * (decimal_fast64_t { UINT64_C(34493884425) } - + x2 * (decimal_fast64_t { UINT32_C(4058104050) } - + x2 * decimal_fast64_t { UINT32_C(135270135) }))))) - }; - - return (x * top) / bot; -} - -template <> -constexpr auto atan_series(decimal128_t x) noexcept -{ - // PadeApproximant[ArcTan[x]/x, {x, 0, {18, 18}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - const decimal128_t x2 { x * x }; - - const decimal128_t - top - { - decimal128_t { UINT64_C(21427381364263875) } - + x2 * (decimal128_t { UINT64_C(91886788553059500) } - + x2 * (decimal128_t { UINT64_C(163675410390191700) } - + x2 * (decimal128_t { UINT64_C(156671838074852100) } - + x2 * (decimal128_t { UINT64_C(87054123957610810) } - + x2 * (decimal128_t { UINT64_C(28283323008669300) } - + x2 * (decimal128_t { UINT64_C(5134145876036100) } - + x2 * (decimal128_t { UINT64_C(463911017673180) } - + x2 * (decimal128_t { UINT64_C(16016872057515) } - + x2 * decimal128_t { UINT64_C(90194313216) })))))))) - }; - - const decimal128_t - bot - { - decimal128_t { UINT64_C(21427381364263875) } - + x2 * (decimal128_t { UINT64_C(99029249007814125) } - + x2 * (decimal128_t { UINT64_C(192399683786610300) } - + x2 * (decimal128_t { UINT64_C(204060270682768500) } - + x2 * (decimal128_t { UINT64_C(128360492848838250) } - + x2 * (decimal128_t { UINT64_C(48688462804731750) } - + x2 * (decimal128_t { UINT64_C(10819658401051500) } - + x2 * (decimal128_t { UINT64_C(1298359008126180) } - + x2 * (decimal128_t { UINT64_C(70562989572075) } - + x2 * decimal128_t { UINT64_C(1120047453525) })))))))) - }; - - return (x * top) / bot; -} - -template <> -constexpr auto atan_series(decimal_fast128_t x) noexcept +template +constexpr auto atan_series(T x) noexcept -> T { - // PadeApproximant[ArcTan[x]/x, {x, 0, {18, 18}}] - // FullSimplify[%] - // HornerForm[Numerator[Out[2]]] - // HornerForm[Denominator[Out[2]]] - - const decimal_fast128_t x2 { x * x }; - - const decimal_fast128_t - top - { - decimal_fast128_t { UINT64_C(21427381364263875) } - + x2 * (decimal_fast128_t { UINT64_C(91886788553059500) } - + x2 * (decimal_fast128_t { UINT64_C(163675410390191700) } - + x2 * (decimal_fast128_t { UINT64_C(156671838074852100) } - + x2 * (decimal_fast128_t { UINT64_C(87054123957610810) } - + x2 * (decimal_fast128_t { UINT64_C(28283323008669300) } - + x2 * (decimal_fast128_t { UINT64_C(5134145876036100) } - + x2 * (decimal_fast128_t { UINT64_C(463911017673180) } - + x2 * (decimal_fast128_t { UINT64_C(16016872057515) } - + x2 * decimal_fast128_t { UINT64_C(90194313216) })))))))) - }; - - const decimal_fast128_t - bot - { - decimal_fast128_t { UINT64_C(21427381364263875) } - + x2 * (decimal_fast128_t { UINT64_C(99029249007814125) } - + x2 * (decimal_fast128_t { UINT64_C(192399683786610300) } - + x2 * (decimal_fast128_t { UINT64_C(204060270682768500) } - + x2 * (decimal_fast128_t { UINT64_C(128360492848838250) } - + x2 * (decimal_fast128_t { UINT64_C(48688462804731750) } - + x2 * (decimal_fast128_t { UINT64_C(10819658401051500) } - + x2 * (decimal_fast128_t { UINT64_C(1298359008126180) } - + x2 * (decimal_fast128_t { UINT64_C(70562989572075) } - + x2 * decimal_fast128_t { UINT64_C(1120047453525) })))))))) - }; - - return (x * top) / bot; + return x + atan_tail(x); } - -template <> constexpr auto atan_values (std::size_t idx) noexcept -> decimal32_t { return atan_detail::atan_table::d32_atan_values [idx]; } -template <> constexpr auto atan_values (std::size_t idx) noexcept -> decimal64_t { return atan_detail::atan_table::d64_atan_values [idx]; } -template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal128_t { return atan_detail::atan_table::d128_atan_values[idx]; } - -template <> constexpr auto atan_values (std::size_t idx) noexcept -> decimal_fast32_t { return atan_detail::atan_table::d32_fast_atan_values [idx]; } -template <> constexpr auto atan_values (std::size_t idx) noexcept -> decimal_fast64_t { return atan_detail::atan_table::d64_fast_atan_values [idx]; } -template <> constexpr auto atan_values(std::size_t idx) noexcept -> decimal_fast128_t { return atan_detail::atan_table::d128_fast_atan_values[idx]; } - } //namespace detail } //namespace decimal } //namespace boost -#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_ATAN_IMPL_HPP +#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_ATAN_IMPL_HPP diff --git a/include/boost/decimal/detail/cmath/impl/fixed_point_series.hpp b/include/boost/decimal/detail/cmath/impl/fixed_point_series.hpp new file mode 100644 index 000000000..fad2fcc0f --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/fixed_point_series.hpp @@ -0,0 +1,222 @@ +// 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_FIXED_POINT_SERIES_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_FIXED_POINT_SERIES_HPP + +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { + +namespace fixed_point_detail { + +template +struct fixed_point_table_imp +{ + // q62_recip[i] = floor(2^(62 + s) / 10^k) with k = i + 7 and s = q62_shift[i]; the Q62 raw + // value of t = sig * 10^-k is then (sig * q62_recip[i]) >> s, with no division + static constexpr std::array q62_recip = {{ + UINT64_C(15474250491067253436), + UINT64_C(12379400392853802748), + UINT64_C(9903520314283042199), + UINT64_C(15845632502852867518), + UINT64_C(12676506002282294014), + UINT64_C(10141204801825835211), + UINT64_C(16225927682921336339), + UINT64_C(12980742146337069071), + UINT64_C(10384593717069655257), + UINT64_C(16615349947311448411), + UINT64_C(13292279957849158729), + UINT64_C(10633823966279326983), + UINT64_C(17014118346046923173), + UINT64_C(13611294676837538538), + UINT64_C(10889035741470030830), + UINT64_C(17422457186352049329), + UINT64_C(13937965749081639463), + UINT64_C(11150372599265311570), + UINT64_C(17840596158824498513), + UINT64_C(14272476927059598810), + UINT64_C(11417981541647679048), + UINT64_C(18268770466636286477), + UINT64_C(14615016373309029182), + UINT64_C(11692013098647223345), + UINT64_C(9353610478917778676), + UINT64_C(14965776766268445882), + UINT64_C(11972621413014756705), + UINT64_C(9578097130411805364) + }}; + static constexpr std::array q62_shift = {{25, 28, 31, 35, 38, 41, 45, 48, 51, 55, 58, 61, 65, 68, 71, 75, 78, 81, 85, 88, 91, 95, 98, 101, 104, 108, 111, 114}}; + + // q126_recip[i] = floor(2^(126 + s) / 10^k) with k = i + 34 and s = q126_shift[i]; the Q126 + // raw value of t = sig * 10^-k is then (sig * q126_recip[i]) >> s, with no division + static constexpr std::array q126_recip = {{ + int128::uint128_t {UINT64_C(9578097130411805364), UINT64_C(13644483260788183358)}, + int128::uint128_t {UINT64_C(15324955408658888583), UINT64_C(10763126773035362404)}, + int128::uint128_t {UINT64_C(12259964326927110866), UINT64_C(15989199047912110569)}, + int128::uint128_t {UINT64_C(9807971461541688693), UINT64_C(9102010423587778132)}, + int128::uint128_t {UINT64_C(15692754338466701909), UINT64_C(10873867862998534689)}, + int128::uint128_t {UINT64_C(12554203470773361527), UINT64_C(12388443105140738074)}, + int128::uint128_t {UINT64_C(10043362776618689222), UINT64_C(2532056854628769813)}, + int128::uint128_t {UINT64_C(16069380442589902755), UINT64_C(7740639782147942024)}, + int128::uint128_t {UINT64_C(12855504354071922204), UINT64_C(6192511825718353619)}, + int128::uint128_t {UINT64_C(10284403483257537763), UINT64_C(8643358275316593218)}, + int128::uint128_t {UINT64_C(16455045573212060421), UINT64_C(10140024425764638826)}, + int128::uint128_t {UINT64_C(13164036458569648337), UINT64_C(4422670725869800738)}, + int128::uint128_t {UINT64_C(10531229166855718669), UINT64_C(14606183024921571560)}, + int128::uint128_t {UINT64_C(16849966666969149871), UINT64_C(12301846395648783526)}, + int128::uint128_t {UINT64_C(13479973333575319897), UINT64_C(6152128301777116498)}, + int128::uint128_t {UINT64_C(10783978666860255917), UINT64_C(15989749085647424168)}, + int128::uint128_t {UINT64_C(17254365866976409468), UINT64_C(10826203278068237376)}, + int128::uint128_t {UINT64_C(13803492693581127574), UINT64_C(16039660251938410547)}, + int128::uint128_t {UINT64_C(11042794154864902059), UINT64_C(16521077016292638761)}, + int128::uint128_t {UINT64_C(17668470647783843295), UINT64_C(15365676781842491048)}, + int128::uint128_t {UINT64_C(14134776518227074636), UINT64_C(12292541425473992838)}, + int128::uint128_t {UINT64_C(11307821214581659709), UINT64_C(6144684325637283947)}, + int128::uint128_t {UINT64_C(18092513943330655534), UINT64_C(17210192550503474962)}, + int128::uint128_t {UINT64_C(14474011154664524427), UINT64_C(17457502855144690293)}, + int128::uint128_t {UINT64_C(11579208923731619542), UINT64_C(6587304654631931588)}, + int128::uint128_t {UINT64_C(9263367138985295633), UINT64_C(16337890167931276240)}, + int128::uint128_t {UINT64_C(14821387422376473014), UINT64_C(4004531380238580045)}, + int128::uint128_t {UINT64_C(11857109937901178411), UINT64_C(6892973918932774359)}, + int128::uint128_t {UINT64_C(9485687950320942729), UINT64_C(1825030320404309164)}, + int128::uint128_t {UINT64_C(15177100720513508366), UINT64_C(10298746142130715309)}, + int128::uint128_t {UINT64_C(12141680576410806693), UINT64_C(4549648098962661924)}, + int128::uint128_t {UINT64_C(9713344461128645354), UINT64_C(11018416108653950185)}, + int128::uint128_t {UINT64_C(15541351137805832567), UINT64_C(6561419329620589327)}, + int128::uint128_t {UINT64_C(12433080910244666053), UINT64_C(16317181907922202431)}, + int128::uint128_t {UINT64_C(9946464728195732843), UINT64_C(1985699082112030975)}, + int128::uint128_t {UINT64_C(15914343565113172548), UINT64_C(17934513790346890853)}, + int128::uint128_t {UINT64_C(12731474852090538039), UINT64_C(3279564588051781713)}, + int128::uint128_t {UINT64_C(10185179881672430431), UINT64_C(6313000485183335694)} + }}; + static constexpr std::array q126_shift = {{114, 118, 121, 124, 128, 131, 134, 138, 141, 144, 148, 151, 154, 158, 161, 164, 168, 171, 174, 178, 181, 184, 188, 191, 194, 197, 201, 204, 207, 211, 214, 217, 221, 224, 227, 231, 234, 237}}; + + // 10^37, to turn a Q126 raw value n into the decimal (n * 10^37 >> 126) * 10^-37 + static constexpr int128::uint128_t pow10_37 {UINT64_C(542101086242752217), UINT64_C(68739955140067328)}; +}; + +#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) + +template constexpr std::array fixed_point_table_imp::q62_recip; +template constexpr std::array fixed_point_table_imp::q62_shift; +template constexpr std::array fixed_point_table_imp::q126_recip; +template constexpr std::array fixed_point_table_imp::q126_shift; +template constexpr int128::uint128_t fixed_point_table_imp::pow10_37; + +#endif + +using fixed_point_table = fixed_point_table_imp; + +// Qm fixed point: the raw integer n stands for n / 2^m. Q62 raw values are std::int64_t +// (step 2^-62), Q126 raw values are int128::int128_t (step 2^-126); both hold |v| < 2. +template +struct q_format; + +template <> +struct q_format +{ + // Raw Q62 value of t = sig * 10^e10, zero once t is below the Q62 step; t must be below 1, + // since e10 > -7 wraps idx past the table and gives zero + static constexpr auto from_decimal(const std::uint64_t sig, const int e10) noexcept -> std::uint64_t + { + const auto idx {static_cast(-e10 - 7)}; + if (sig == 0U || idx >= fixed_point_table::q62_recip.size()) + { + return 0U; + } + return static_cast((int128::uint128_t {sig} * fixed_point_table::q62_recip[idx]) >> fixed_point_table::q62_shift[idx]); + } + + // (a * b) >> 62 with the sign of a + static constexpr auto mul(const std::int64_t a, const std::uint64_t b) noexcept -> std::int64_t + { + const auto m {a < 0 ? UINT64_C(0) - static_cast(a) : static_cast(a)}; + const auto p {static_cast(static_cast((int128::uint128_t {m} * b) >> 62))}; + return a < 0 ? -p : p; + } + + // n / 2^62 as the decimal (n * 10^18 >> 62) * 10^-18 + template + static constexpr auto to_decimal(const std::int64_t n) noexcept -> T + { + const auto m {n < 0 ? UINT64_C(0) - static_cast(n) : static_cast(n)}; + const auto sig {static_cast((int128::uint128_t {m} * UINT64_C(1000000000000000000)) >> 62)}; + return T {sig, -18, n < 0 ? construction_sign::negative : construction_sign::positive}; + } +}; + +template <> +struct q_format +{ + // Raw Q126 value of t = sig * 10^e10, zero once t is below the Q126 step; t must be below 1, + // since e10 > -34 wraps idx past the table and gives zero + static constexpr auto from_decimal(const int128::uint128_t sig, const int e10) noexcept -> int128::uint128_t + { + const auto idx {static_cast(-e10 - 34)}; + if (sig == 0U || idx >= fixed_point_table::q126_recip.size()) + { + return 0U; + } + return static_cast(umul256(sig, fixed_point_table::q126_recip[idx]) >> fixed_point_table::q126_shift[idx]); + } + + // (a * b) >> 126 with the sign of a + static constexpr auto mul(const int128::int128_t a, const int128::uint128_t b) noexcept -> int128::int128_t + { + const int128::uint128_t m {a < 0 ? -a : a}; + const int128::int128_t p {static_cast(umul256(m, b) >> 126)}; + return a < 0 ? -p : p; + } + + // n / 2^126 as the decimal (n * 10^37 >> 126) * 10^-37 + template + static constexpr auto to_decimal(const int128::int128_t n) noexcept -> T + { + const int128::uint128_t m {n < 0 ? -n : n}; + const auto sig {static_cast(umul256(m, fixed_point_table::pow10_37) >> 126)}; + return T {sig, -37, n < 0 ? construction_sign::negative : construction_sign::positive}; + } +}; + +} // namespace fixed_point_detail + +// x^3 P(x^2) from raw fixed-point coefficients of P, highest power first: u = x^2 P(x^2) is +// formed in fixed point and rounded once to T, and the result is x * u +template +constexpr auto fixed_point_odd_tail(T x, const std::array& coeffs) noexcept -> T +{ + using format = fixed_point_detail::q_format; + + const T t {x * x}; + int e10 {}; + const auto sig {frexp10(t, &e10)}; + const auto tq {format::from_decimal(sig, e10)}; + + Raw r {coeffs[0]}; + for (std::size_t i {1}; i < N; ++i) + { + r = format::mul(r, tq) + coeffs[i]; + } + + return x * format::template to_decimal(format::mul(r, tq)); +} + +} // namespace detail +} // namespace decimal +} // namespace boost + +#endif // BOOST_DECIMAL_DETAIL_CMATH_IMPL_FIXED_POINT_SERIES_HPP diff --git a/include/boost/decimal/detail/cmath/impl/split_pi.hpp b/include/boost/decimal/detail/cmath/impl/split_pi.hpp new file mode 100644 index 000000000..d367ad484 --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/split_pi.hpp @@ -0,0 +1,130 @@ +// 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_SPLIT_PI_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_SPLIT_PI_HPP + +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { + +namespace split_pi_detail { + +// Indices into split_pi_values: each _hi value is the constant rounded to the type and its +// _lo value is the part left over; three_quarter_pi has no low part, since only atan2 uses it +enum : std::size_t { quarter_pi_hi, quarter_pi_lo, half_pi_hi, half_pi_lo, three_quarter_pi, pi_hi, pi_lo }; + +// Use a struct of arrays so that we can have static constexpr arrays of coefficients. +// See https://github.com/boostorg/math/issues/923 for further information +template +struct split_pi_table_imp +{ + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d32_values = {{ + decimal32_t {UINT64_C(7853982), -7}, + decimal32_t {UINT64_C(3660255), -14, construction_sign::negative}, + decimal32_t {UINT64_C(1570796), -6}, + decimal32_t {UINT64_C(3267949), -13}, + decimal32_t {UINT64_C(2356194), -6}, + decimal32_t {UINT64_C(3141593), -6}, + decimal32_t {UINT64_C(3464102), -13, construction_sign::negative} + }}; + + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d32_fast_values = {{ + decimal_fast32_t {UINT64_C(7853982), -7}, + decimal_fast32_t {UINT64_C(3660255), -14, construction_sign::negative}, + decimal_fast32_t {UINT64_C(1570796), -6}, + decimal_fast32_t {UINT64_C(3267949), -13}, + decimal_fast32_t {UINT64_C(2356194), -6}, + decimal_fast32_t {UINT64_C(3141593), -6}, + decimal_fast32_t {UINT64_C(3464102), -13, construction_sign::negative} + }}; + + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d64_values = {{ + decimal64_t {UINT64_C(7853981633974483), -16}, + decimal64_t {UINT64_C(9615660845819876), -33}, + decimal64_t {UINT64_C(1570796326794897), -15}, + decimal64_t {UINT64_C(3807686783083602), -31, construction_sign::negative}, + decimal64_t {UINT64_C(2356194490192345), -15}, + decimal64_t {UINT64_C(3141592653589793), -15}, + decimal64_t {UINT64_C(2384626433832795), -31} + }}; + + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d64_fast_values = {{ + decimal_fast64_t {UINT64_C(7853981633974483), -16}, + decimal_fast64_t {UINT64_C(9615660845819876), -33}, + decimal_fast64_t {UINT64_C(1570796326794897), -15}, + decimal_fast64_t {UINT64_C(3807686783083602), -31, construction_sign::negative}, + decimal_fast64_t {UINT64_C(2356194490192345), -15}, + decimal_fast64_t {UINT64_C(3141592653589793), -15}, + decimal_fast64_t {UINT64_C(2384626433832795), -31} + }}; + + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d128_values = {{ + decimal128_t {int128::uint128_t{UINT64_C(425765197510819), UINT64_C(5970600460659265253)}, -34}, + decimal128_t {int128::uint128_t{UINT64_C(114108442474915), UINT64_C(12088338259637095055)}, -68}, + decimal128_t {int128::uint128_t{UINT64_C(85153039502163), UINT64_C(15951515351099494343)}, -33}, + decimal128_t {int128::uint128_t{UINT64_C(239662122992084), UINT64_C(329523718765553795)}, -67}, + decimal128_t {int128::uint128_t{UINT64_C(127729559253245), UINT64_C(14703900989794465707)}, -33}, + decimal128_t {int128::uint128_t{UINT64_C(170306079004327), UINT64_C(13456286628489437071)}, -33}, + decimal128_t {int128::uint128_t{UINT64_C(62776840258584), UINT64_C(3343964766419005178)}, -67, construction_sign::negative} + }}; + + // pi/4, pi/2, 3 pi/4 and pi, each rounded to the type, and the rest of pi/4, pi/2 and pi + static constexpr std::array d128_fast_values = {{ + decimal_fast128_t {int128::uint128_t{UINT64_C(425765197510819), UINT64_C(5970600460659265253)}, -34}, + decimal_fast128_t {int128::uint128_t{UINT64_C(114108442474915), UINT64_C(12088338259637095055)}, -68}, + decimal_fast128_t {int128::uint128_t{UINT64_C(85153039502163), UINT64_C(15951515351099494343)}, -33}, + decimal_fast128_t {int128::uint128_t{UINT64_C(239662122992084), UINT64_C(329523718765553795)}, -67}, + decimal_fast128_t {int128::uint128_t{UINT64_C(127729559253245), UINT64_C(14703900989794465707)}, -33}, + decimal_fast128_t {int128::uint128_t{UINT64_C(170306079004327), UINT64_C(13456286628489437071)}, -33}, + decimal_fast128_t {int128::uint128_t{UINT64_C(62776840258584), UINT64_C(3343964766419005178)}, -67, construction_sign::negative} + }}; +}; + +#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) + +template constexpr std::array split_pi_table_imp::d32_values; +template constexpr std::array split_pi_table_imp::d32_fast_values; +template constexpr std::array split_pi_table_imp::d64_values; +template constexpr std::array split_pi_table_imp::d64_fast_values; +template constexpr std::array split_pi_table_imp::d128_values; +template constexpr std::array split_pi_table_imp::d128_fast_values; + +#endif + +using split_pi_table = split_pi_table_imp; + +} //namespace split_pi_detail + +template +constexpr auto split_pi_values(std::size_t idx) noexcept -> T; + +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal32_t { return split_pi_detail::split_pi_table::d32_values[idx]; } +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal_fast32_t { return split_pi_detail::split_pi_table::d32_fast_values[idx]; } +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal64_t { return split_pi_detail::split_pi_table::d64_values[idx]; } +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal_fast64_t { return split_pi_detail::split_pi_table::d64_fast_values[idx]; } +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal128_t { return split_pi_detail::split_pi_table::d128_values[idx]; } +template <> constexpr auto split_pi_values(std::size_t idx) noexcept -> decimal_fast128_t { return split_pi_detail::split_pi_table::d128_fast_values[idx]; } + +} //namespace detail +} //namespace decimal +} //namespace boost + +#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_SPLIT_PI_HPP diff --git a/test/Jamfile b/test/Jamfile index 3bd2e8116..beb1c7846 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -247,6 +247,7 @@ compile-fail test_illegal_fast_quantize.cpp ; compile-fail test_illegal_fast_samequantum.cpp ; run test_implicit_integral_conversion.cpp : : : off ; +run test_inverse_trig_rounding.cpp ; run test_laguerre.cpp ; run test_legal_implicit_conversions.cpp ; run test_legendre.cpp ; diff --git a/test/test_asin.cpp b/test/test_asin.cpp index 04ae7fb81..9d4fe0cc5 100644 --- a/test/test_asin.cpp +++ b/test/test_asin.cpp @@ -199,7 +199,9 @@ auto test_asin_edge() -> void BOOST_TEST_EQ(asin_tiny1 / nl::epsilon(), T(1)); BOOST_TEST_EQ(asin_tiny2 / nl::epsilon(), ctrl_tiny2); - constexpr T half_pi { numbers::pi_v / 2 }; + // pi/2 correctly rounded to T; numbers::pi_v / 2 rounds twice + using namespace boost::decimal::literals; + const T half_pi { static_cast(1.570796326794896619231321691639751_DL) }; std::random_device rd; std::mt19937_64 gen(rd()); @@ -255,9 +257,13 @@ void test_asin_1137() const T sqrt_tiny1 { sqrt(nl::epsilon()) }; const T sqrt_tiny2 { sqrt(nl::epsilon() * 1000/999) }; - BOOST_TEST_EQ(sqrt_tiny0, asin(sqrt_tiny0)); - BOOST_TEST_EQ(sqrt_tiny1, asin(sqrt_tiny1)); - BOOST_TEST_EQ(sqrt_tiny2, asin(sqrt_tiny2)); + // asin(x) = x (1 + eps/6 + ...) here, which can round one or two ulps above x + BOOST_TEST_LE(asin(sqrt_tiny0) - sqrt_tiny0, sqrt_tiny0 * nl::epsilon() / 2); + BOOST_TEST_LE(asin(sqrt_tiny1) - sqrt_tiny1, sqrt_tiny1 * nl::epsilon() / 2); + BOOST_TEST_LE(asin(sqrt_tiny2) - sqrt_tiny2, sqrt_tiny2 * nl::epsilon() / 2); + BOOST_TEST_GE(asin(sqrt_tiny0), sqrt_tiny0); + BOOST_TEST_GE(asin(sqrt_tiny1), sqrt_tiny1); + BOOST_TEST_GE(asin(sqrt_tiny2), sqrt_tiny2); const T cbrt_tiny0 { cbrt(nl::epsilon() * 999/1000) }; const T cbrt_tiny1 { cbrt(nl::epsilon()) }; diff --git a/test/test_atan.cpp b/test/test_atan.cpp index 3193540e1..d64352c1a 100644 --- a/test/test_atan.cpp +++ b/test/test_atan.cpp @@ -288,8 +288,11 @@ void test_atan() // Edge cases std::uniform_int_distribution one(1,1); - BOOST_TEST_EQ(atan(std::numeric_limits::infinity() * Dec(one(rng))), numbers::pi_v/2); - BOOST_TEST_EQ(atan(-std::numeric_limits::infinity() * Dec(one(rng))), -numbers::pi_v/2); + // pi/2 correctly rounded to Dec; numbers::pi_v / 2 rounds twice + using namespace boost::decimal::literals; + const Dec half_pi {static_cast(1.570796326794896619231321691639751_DL)}; + BOOST_TEST_EQ(atan(std::numeric_limits::infinity() * Dec(one(rng))), half_pi); + BOOST_TEST_EQ(atan(-std::numeric_limits::infinity() * Dec(one(rng))), -half_pi); BOOST_TEST(isnan(atan(std::numeric_limits::quiet_NaN() * Dec(one(rng))))); BOOST_TEST_EQ(atan(Dec(0) * Dec(one(rng))), Dec(0)); BOOST_TEST_EQ(atan(std::numeric_limits::epsilon() * Dec(one(rng))), std::numeric_limits::epsilon() * Dec(one(rng))); diff --git a/test/test_atan2.cpp b/test/test_atan2.cpp index 2bc8fa3d8..4611fdaba 100644 --- a/test/test_atan2.cpp +++ b/test/test_atan2.cpp @@ -82,19 +82,25 @@ void test() // Edge cases std::uniform_int_distribution one(1,1); + // pi/2, 3 pi/4 and pi/4 correctly rounded to Dec; numbers::pi_v / 2 rounds twice + using namespace boost::decimal::literals; + const Dec half_pi {static_cast(1.570796326794896619231321691639751_DL)}; + const Dec three_quarter_pi {static_cast(2.356194490192344928846982537459627_DL)}; + const Dec quarter_pi {static_cast(0.7853981633974483096156608458198757_DL)}; + BOOST_TEST(isnan(atan2(Dec{one(rng)}, std::numeric_limits::quiet_NaN()))); BOOST_TEST(isnan(atan2(std::numeric_limits::quiet_NaN(), Dec{one(rng)}))); BOOST_TEST_EQ(atan2(Dec{0 * one(rng)}, -Dec(1)), numbers::pi_v); BOOST_TEST_EQ(atan2(Dec{0 * -one(rng)}, -Dec(1)), numbers::pi_v); BOOST_TEST_EQ(atan2(Dec{0 * one(rng)}, Dec(1)), Dec{0 * one(rng)}); - BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), Dec{one(rng)}), numbers::pi_v / 2); - BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), Dec{one(rng)}), -numbers::pi_v / 2); - BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), -std::numeric_limits::infinity()), 3 * one(rng) * numbers::pi_v / 4); - BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), -std::numeric_limits::infinity()), -3 * one(rng) * numbers::pi_v / 4); - BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), std::numeric_limits::infinity()), one(rng) * numbers::pi_over_four_v); - BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), std::numeric_limits::infinity()), -one(rng) * numbers::pi_over_four_v); - BOOST_TEST_EQ(atan2(-Dec(1), Dec{0 * one(rng)}), -numbers::pi_v / 2); - BOOST_TEST_EQ(atan2(Dec(1), Dec{0 * one(rng)}), numbers::pi_v / 2); + BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), Dec{one(rng)}), half_pi); + BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), Dec{one(rng)}), -half_pi); + BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), -std::numeric_limits::infinity()), three_quarter_pi); + BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), -std::numeric_limits::infinity()), -three_quarter_pi); + BOOST_TEST_EQ(atan2(std::numeric_limits::infinity(), std::numeric_limits::infinity()), quarter_pi); + BOOST_TEST_EQ(atan2(-std::numeric_limits::infinity(), std::numeric_limits::infinity()), -quarter_pi); + BOOST_TEST_EQ(atan2(-Dec(1), Dec{0 * one(rng)}), -half_pi); + BOOST_TEST_EQ(atan2(Dec(1), Dec{0 * one(rng)}), half_pi); BOOST_TEST_EQ(atan2(-Dec{one(rng)}, -std::numeric_limits::infinity()), -numbers::pi_v); BOOST_TEST_EQ(atan2(Dec{one(rng)}, -std::numeric_limits::infinity()), numbers::pi_v); BOOST_TEST_EQ(atan2(-Dec{one(rng)}, std::numeric_limits::infinity()), -Dec{0 * one(rng)}); diff --git a/test/test_inverse_trig_rounding.cpp b/test/test_inverse_trig_rounding.cpp new file mode 100644 index 000000000..4d93bb499 --- /dev/null +++ b/test/test_inverse_trig_rounding.cpp @@ -0,0 +1,187 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include + +using namespace boost::decimal; +using namespace boost::decimal::literals; + +// hi + lo is the MPFR value to twice the precision of T. The old asin, acos, atan and atan2 were off +// by 1 to 1e18 ulps at these points; the error must now stay below one ulp. +template +void check(const T got, const T hi, const T lo) +{ + const T ulp {abs(nextafter(hi, hi * 2) - hi)}; + BOOST_TEST_LT(abs((got - hi) - lo), ulp); +} + +// Walk 400 consecutive values around each start; asin and atan must not decrease, acos must +// not increase. The starts are the points where the formulas change, and points between them. +template +void check_monotone() +{ + const T up {std::numeric_limits::max()}; + const T down {std::numeric_limits::lowest()}; + const T h {numbers::pi_v / 4 - T {5, -1}}; + const T starts[] { + // asin and acos: +-0.5, s = 0.1 and s = 0.05 where s + s starts to round, and pi/4 - s = 0.5 + T {-5, -1}, T {5, -1}, T {-98, -2}, T {98, -2}, T {-995, -3}, T {995, -3}, 1 - 2 * h * h, + // atan: the breakpoints 7/16, 11/16, 19/16 and 39/16, and the points where |q| = 0.1 + T {4375, -4}, T {6875, -4}, T {11875, -4}, T {24375, -4}, T {12} / 19, T {9} / 11, T {28} / 23, T {32} / 17, + T {-9, -1}, T {-75, -2}, T {-6, -1}, T {75, -2}, T {2, 0}, T {236, -2}}; + for (const T start : starts) + { + T x {start}; + for (int i {}; i < 200; ++i) + { + x = nextafter(x, down); + } + for (int i {}; i < 400; ++i) + { + const T next {nextafter(x, up)}; + if (abs(next) <= 1) + { + BOOST_TEST_LE(asin(x), asin(next)); + BOOST_TEST_GE(acos(x), acos(next)); + } + BOOST_TEST_LE(atan(x), atan(next)); + x = next; + } + } +} + +// atan2 must return pi/2, 3 pi/4 and pi/4 rounded once to T, as asin and atan do +template +void check_atan2() +{ + const T half_pi {static_cast(1.570796326794896619231321691639751_DL)}; + const T one {1}; + const T zero {0}; + BOOST_TEST_EQ(atan2(one, zero), half_pi); + BOOST_TEST_EQ(atan2(-one, zero), -half_pi); + BOOST_TEST_EQ(asin(one), half_pi); + #ifndef BOOST_DECIMAL_FAST_MATH + const T inf {std::numeric_limits::infinity()}; + BOOST_TEST_EQ(atan(inf), half_pi); + BOOST_TEST_EQ(atan2(inf, one), half_pi); + BOOST_TEST_EQ(atan2(inf, inf), static_cast(0.7853981633974483096156608458198757_DL)); + BOOST_TEST_EQ(atan2(inf, -inf), static_cast(2.356194490192344928846982537459627_DL)); + #endif +} + +int main() +{ + check_monotone(); + check_monotone(); + check_monotone(); + check_monotone(); + check_monotone(); + check_monotone(); + + check_atan2(); + check_atan2(); + check_atan2(); + check_atan2(); + check_atan2(); + check_atan2(); + + check(asin(5.065049e-01_DF), 5.311264e-01_DF, 1.602916e-08_DF); + check(asin(9.633317e-03_DF), 9.633466e-03_DF, 2.804126e-12_DF); + check(acos(9.999999e-01_DF), 4.472136e-04_DF, -7.732620e-13_DF); + check(acos(8.853146e-01_DF), 4.836262e-01_DF, -1.910526e-08_DF); + check(acos(-9.995746e-01_DF), 3.112423e+00_DF, 1.433340e-07_DF); + check(atan(4.364837e-01_DF), 4.115571e-01_DF, -2.236126e-09_DF); + check(atan(6.888915e-01_DF), 6.032316e-01_DF, 2.368609e-08_DF); + check(atan(3.206842e+01_DF), 1.539623e+00_DF, 1.023047e-07_DF); + check(asin(5.100377e-01_DF), 5.352286e-01_DF, 1.916637e-08_DF); + check(asin(5.492685e-01_DF), 5.814886e-01_DF, 1.502855e-08_DF); + check(acos(5.408637e-01_DF), 9.993327e-01_DF, -1.047165e-09_DF); + check(atan(6.839736e-01_DF), 5.998888e-01_DF, 1.746916e-08_DF); + check(atan(7.159895e-01_DF), 6.213768e-01_DF, -4.046992e-08_DF); + check(atan2(1.149130e+00_DF, 8.516064e-01_DF), 9.330234e-01_DF, -3.083453e-08_DF); + + check(asin(5.065049e-01_DFF), 5.311264e-01_DFF, 1.602916e-08_DFF); + check(asin(9.633317e-03_DFF), 9.633466e-03_DFF, 2.804126e-12_DFF); + check(acos(9.999999e-01_DFF), 4.472136e-04_DFF, -7.732620e-13_DFF); + check(acos(8.853146e-01_DFF), 4.836262e-01_DFF, -1.910526e-08_DFF); + check(acos(-9.995746e-01_DFF), 3.112423e+00_DFF, 1.433340e-07_DFF); + check(atan(4.364837e-01_DFF), 4.115571e-01_DFF, -2.236126e-09_DFF); + check(atan(6.888915e-01_DFF), 6.032316e-01_DFF, 2.368609e-08_DFF); + check(atan(3.206842e+01_DFF), 1.539623e+00_DFF, 1.023047e-07_DFF); + check(asin(5.100377e-01_DFF), 5.352286e-01_DFF, 1.916637e-08_DFF); + check(asin(5.492685e-01_DFF), 5.814886e-01_DFF, 1.502855e-08_DFF); + check(acos(5.408637e-01_DFF), 9.993327e-01_DFF, -1.047165e-09_DFF); + check(atan(6.839736e-01_DFF), 5.998888e-01_DFF, 1.746916e-08_DFF); + check(atan(7.159895e-01_DFF), 6.213768e-01_DFF, -4.046992e-08_DFF); + check(atan2(1.149130e+00_DFF, 8.516064e-01_DFF), 9.330234e-01_DFF, -3.083453e-08_DFF); + + check(asin(0.5_DD), 5.235987755982989e-01_DD, -2.692289276945342e-17_DD); + check(asin(5.004651001747461e-01_DD), 5.241359103330428e-01_DD, 2.990295215924292e-17_DD); + check(asin(1.191104442197430e-05_DD), 1.191104442225594e-05_DD, 2.559650313877577e-21_DD); + check(acos(9.999999999999999e-01_DD), 1.414213562373095e-08_DD, 6.058680174398549e-25_DD); + check(acos(9.999999999997268e-01_DD), 7.391887445030700e-07_DD, -1.334429551076057e-23_DD); + check(acos(5.408223492936285e-01_DD), 9.993818602282364e-01_DD, -2.634277016069807e-18_DD); + check(atan(4.364837490347062e-01_DD), 4.115571389515791e-01_DD, -3.417886906675521e-17_DD); + check(atan(7.136198652127231e-01_DD), 6.198084470537839e-01_DD, -1.771590134732764e-17_DD); + check(atan(7.812423720033109e+01_DD), 1.557996900819007e+00_DD, 1.214730049123154e-16_DD); + check(asin(5.513668750325045e-01_DD), 5.840017749177506e-01_DD, -1.421285587919738e-18_DD); + check(acos(5.616429838168425e-01_DD), 9.744260944481967e-01_DD, -1.956861786536074e-18_DD); + check(atan(6.862655641327955e-01_DD), 6.014486253684519e-01_DD, 1.662014588125191e-17_DD); + check(atan(7.051440412349395e-01_DD), 6.141700043207512e-01_DD, -3.855529142780040e-17_DD); + check(atan(6.949647737224145e-01_DD), 6.073386179069759e-01_DD, 4.042904847421836e-17_DD); + check(atan2(1.120258985683803e+01_DD, 1.091340798186341e+01_DD), 7.984731057062606e-01_DD, -4.169418992245100e-17_DD); + + check(asin(0.5_DDF), 5.235987755982989e-01_DDF, -2.692289276945342e-17_DDF); + check(asin(5.004651001747461e-01_DDF), 5.241359103330428e-01_DDF, 2.990295215924292e-17_DDF); + check(asin(1.191104442197430e-05_DDF), 1.191104442225594e-05_DDF, 2.559650313877577e-21_DDF); + check(acos(9.999999999999999e-01_DDF), 1.414213562373095e-08_DDF, 6.058680174398549e-25_DDF); + check(acos(9.999999999997268e-01_DDF), 7.391887445030700e-07_DDF, -1.334429551076057e-23_DDF); + check(acos(5.408223492936285e-01_DDF), 9.993818602282364e-01_DDF, -2.634277016069807e-18_DDF); + check(atan(4.364837490347062e-01_DDF), 4.115571389515791e-01_DDF, -3.417886906675521e-17_DDF); + check(atan(7.136198652127231e-01_DDF), 6.198084470537839e-01_DDF, -1.771590134732764e-17_DDF); + check(atan(7.812423720033109e+01_DDF), 1.557996900819007e+00_DDF, 1.214730049123154e-16_DDF); + check(asin(5.513668750325045e-01_DDF), 5.840017749177506e-01_DDF, -1.421285587919738e-18_DDF); + check(acos(5.616429838168425e-01_DDF), 9.744260944481967e-01_DDF, -1.956861786536074e-18_DDF); + check(atan(6.862655641327955e-01_DDF), 6.014486253684519e-01_DDF, 1.662014588125191e-17_DDF); + check(atan(7.051440412349395e-01_DDF), 6.141700043207512e-01_DDF, -3.855529142780040e-17_DDF); + check(atan(6.949647737224145e-01_DDF), 6.073386179069759e-01_DDF, 4.042904847421836e-17_DDF); + check(atan2(1.120258985683803e+01_DDF, 1.091340798186341e+01_DDF), 7.984731057062606e-01_DDF, -4.169418992245100e-17_DDF); + + check(asin(0.5_DL), 5.235987755982988730771072305465838e-01_DL, 1.403286156656251763682915743205130e-35_DL); + check(asin(5.003912494167833183647481136970553e-01_DL), 5.240506104602699416431335061998704e-01_DL, 3.828027545876172039940440784250984e-36_DL); + check(asin(1.021951104494350724501767305490306e-11_DL), 1.021951104494350724501785093981033e-11_DL, -1.663727041084267689046944494101754e-45_DL); + check(acos(9.999999999999999999999999999998349e-01_DL), 5.746303159423456662845786731898078e-16_DL, -4.395648457605540895340625419625079e-50_DL); + check(acos(5.404815354371389993198662570037106e-01_DL), 9.997869898888469193058044698685081e-01_DL, 8.558703671778992246137722745628842e-36_DL); + check(acos(-5.003912494167833183647481136970553e-01_DL), 2.094846937255166560874455197839622e+00_DL, -1.540733877544362750495720869195951e-34_DL); + check(atan(0.4_DL), 3.805063771123648863035879168104331e-01_DL, 4.497405713658100837576305622324200e-36_DL); + check(atan(0.5_DL), 4.636476090008061162142562314612144e-01_DL, 2.028537054286120263810933088720198e-36_DL); + check(atan(1_DL), 7.853981633974483096156608458198757e-01_DL, 2.104929234984377645524373614807695e-35_DL); + check(atan(7.068526585954127093318783992523534e-01_DL), 6.153102733138216952690804144016625e-01_DL, -2.355736704427937425845252354220856e-35_DL); + check(asin(5.357061056349391872935001409387989e-01_DL), 5.653437428133160123277855997530434e-01_DL, 1.024303028026974284782679076163778e-35_DL); + check(acos(5.419197605139804475954430930508049e-01_DL), 9.980766359735971232743255248684227e-01_DL, -5.902616948333705555725733542717919e-36_DL); + check(atan(6.660927254531924775945299356519535e-01_DL), 5.876051543687935417357987351933894e-01_DL, 2.002308346159484947226392029195896e-35_DL); + check(atan(6.968108909601345933169770898367005e-01_DL), 6.085824144794810607711470552450067e-01_DL, 3.380492439684669212621325350444062e-35_DL); + check(atan2(1.168271357189183012390725086719514e+01_DL, 9.291551526980655881838138847395345e+00_DL), 8.989126408521308262359466298659614e-01_DL, 4.007200880413205400346098160382424e-35_DL); + check(atan2(2.593355382902205084382369219298895e+00_DL, 2.587222617437111113066030265575797e+01_DL), 9.990334054463293937383915469668977e-02_DL, 3.668097627764618165721978443736755e-36_DL); + + check(asin(0.5_DLF), 5.235987755982988730771072305465838e-01_DLF, 1.403286156656251763682915743205130e-35_DLF); + check(asin(5.003912494167833183647481136970553e-01_DLF), 5.240506104602699416431335061998704e-01_DLF, 3.828027545876172039940440784250984e-36_DLF); + check(asin(1.021951104494350724501767305490306e-11_DLF), 1.021951104494350724501785093981033e-11_DLF, -1.663727041084267689046944494101754e-45_DLF); + check(acos(9.999999999999999999999999999998349e-01_DLF), 5.746303159423456662845786731898078e-16_DLF, -4.395648457605540895340625419625079e-50_DLF); + check(acos(5.404815354371389993198662570037106e-01_DLF), 9.997869898888469193058044698685081e-01_DLF, 8.558703671778992246137722745628842e-36_DLF); + check(acos(-5.003912494167833183647481136970553e-01_DLF), 2.094846937255166560874455197839622e+00_DLF, -1.540733877544362750495720869195951e-34_DLF); + check(atan(0.4_DLF), 3.805063771123648863035879168104331e-01_DLF, 4.497405713658100837576305622324200e-36_DLF); + check(atan(0.5_DLF), 4.636476090008061162142562314612144e-01_DLF, 2.028537054286120263810933088720198e-36_DLF); + check(atan(1_DLF), 7.853981633974483096156608458198757e-01_DLF, 2.104929234984377645524373614807695e-35_DLF); + check(atan(7.068526585954127093318783992523534e-01_DLF), 6.153102733138216952690804144016625e-01_DLF, -2.355736704427937425845252354220856e-35_DLF); + check(asin(5.357061056349391872935001409387989e-01_DLF), 5.653437428133160123277855997530434e-01_DLF, 1.024303028026974284782679076163778e-35_DLF); + check(acos(5.419197605139804475954430930508049e-01_DLF), 9.980766359735971232743255248684227e-01_DLF, -5.902616948333705555725733542717919e-36_DLF); + check(atan(6.660927254531924775945299356519535e-01_DLF), 5.876051543687935417357987351933894e-01_DLF, 2.002308346159484947226392029195896e-35_DLF); + check(atan(6.968108909601345933169770898367005e-01_DLF), 6.085824144794810607711470552450067e-01_DLF, 3.380492439684669212621325350444062e-35_DLF); + check(atan2(1.168271357189183012390725086719514e+01_DLF, 9.291551526980655881838138847395345e+00_DLF), 8.989126408521308262359466298659614e-01_DLF, 4.007200880413205400346098160382424e-35_DLF); + check(atan2(2.593355382902205084382369219298895e+00_DLF, 2.587222617437111113066030265575797e+01_DLF), 9.990334054463293937383915469668977e-02_DLF, 3.668097627764618165721978443736755e-36_DLF); + return boost::report_errors(); +}