Skip to content

Fix the accuracy and monotonicity of asin, acos, atan and atan2 #1478

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ckormanyos merged 1 commit into from
Oct 1, 2026
Merged

Fix the accuracy and monotonicity of asin, acos, atan and atan2 #1478

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 22 additions & 13 deletions include/boost/decimal/detail/cmath/acos.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,14 +6,14 @@
#define BOOST_DECIMAL_DETAIL_CMATH_ACOS_HPP

#include <boost/decimal/fwd.hpp>
#include <boost/decimal/numbers.hpp>
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/detail/promotion.hpp>
#include <boost/decimal/detail/config.hpp>
#include <boost/decimal/detail/cmath/fabs.hpp>
#include <boost/decimal/detail/cmath/sqrt.hpp>
#include <boost/decimal/detail/cmath/impl/asin_impl.hpp>
#include <boost/decimal/detail/cmath/impl/split_pi.hpp>


#ifndef BOOST_DECIMAL_BUILD_MODULE
Expand All @@ -37,7 +37,6 @@ constexpr auto acos_impl(const T x) noexcept
}
#endif

constexpr auto half_pi {numbers::pi_v<T> / 2};
const auto absx {fabs(static_cast<T>(x))};

T result {};
Expand All @@ -46,21 +45,31 @@ constexpr auto acos_impl(const T x) noexcept
{
result = std::numeric_limits<T>::quiet_NaN();
}
else if (x < T{-5, -1})
else if (absx <= T{5, -1})
{
result = numbers::pi_v<T> - 2 * detail::asin_series(sqrt((1 - absx) / 2));
}
else if (x < -std::numeric_limits<T>::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<T>(detail::split_pi_detail::half_pi_hi) -
(detail::asin_series(x) - detail::split_pi_values<T>(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<T>(detail::split_pi_detail::pi_hi) -
((s + s) + ((lo + lo) - detail::split_pi_values<T>(detail::split_pi_detail::pi_lo)));
}
}

return result;
Expand Down
16 changes: 10 additions & 6 deletions include/boost/decimal/detail/cmath/asin.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,14 +7,14 @@
#define BOOST_DECIMAL_DETAIL_CMATH_ASIN_HPP

#include <boost/decimal/fwd.hpp>
#include <boost/decimal/numbers.hpp>
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/detail/config.hpp>
#include <boost/decimal/detail/cmath/fpclassify.hpp>
#include <boost/decimal/detail/cmath/fabs.hpp>
#include <boost/decimal/detail/cmath/sqrt.hpp>
#include <boost/decimal/detail/cmath/impl/asin_impl.hpp>
#include <boost/decimal/detail/cmath/impl/split_pi.hpp>

#ifndef BOOST_DECIMAL_BUILD_MODULE
#include <type_traits>
Expand Down Expand Up @@ -52,19 +52,23 @@ 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 })
{
result = asin_series(absx);
}
else
{
constexpr T half_pi { numbers::pi_v<T> / 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<T>(split_pi_detail::quarter_pi_lo) };
const T head { split_pi_values<T>(split_pi_detail::quarter_pi_hi) - s };

result = head < T { 5, -1 } ? (head + head) - (lo + lo) : head + (head - (lo + lo));
}
else if (absx > one)
{
Expand All @@ -76,7 +80,7 @@ constexpr auto asin_impl(const T x) noexcept
}
else
{
result = half_pi;
result = split_pi_values<T>(split_pi_detail::half_pi_hi);
}
}

Expand Down
107 changes: 65 additions & 42 deletions include/boost/decimal/detail/cmath/atan.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,12 @@

#include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
#include <boost/decimal/detail/cmath/impl/atan_impl.hpp>
#include <boost/decimal/detail/cmath/impl/split_pi.hpp>
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/detail/config.hpp>
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/numbers.hpp>
#include <boost/decimal/detail/cmath/fabs.hpp>
#include <boost/decimal/detail/cmath/fma.hpp>

#ifndef BOOST_DECIMAL_BUILD_MODULE
#include <type_traits>
Expand All @@ -23,16 +25,15 @@ namespace decimal {

namespace detail {

template <typename T>
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 <bool has_lo = false, typename T>
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<T> / 2 };

if (fpc == FP_ZERO
#ifndef BOOST_DECIMAL_FAST_MATH
|| fpc == FP_NAN
Expand All @@ -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<T>(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<T>) - one) / (numbers::sqrt3_v<T> + x)
: x
};

constexpr T half { 5, -1 };
constexpr T three_halves { 15, -1 };

result = (fx_arg <= std::numeric_limits<T>::epsilon()) ? fx_arg
: (fx_arg <= T { 4375, -4 }) ? detail::atan_series (fx_arg)
: (fx_arg <= T { 6875, -4 }) ? detail::atan_values<T>(0U) + detail::atan_series((fx_arg - half) / (one + fx_arg / 2))
: (fx_arg <= T { 11875, -4 }) ? detail::atan_values<T>(1U) + detail::atan_series((fx_arg - one) / (fx_arg + one))
: detail::atan_values<T>(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<T> / 6 };

result += my_pi_over_six;
}
c_hi = detail::atan_values<T>(detail::atan_detail::atan_half_hi);
c_lo = detail::atan_values<T>(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<T>(detail::split_pi_detail::quarter_pi_hi);
c_lo = detail::split_pi_values<T>(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<T>(detail::atan_detail::atan_three_halves_hi);
c_lo = detail::atan_values<T>(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<T>(detail::split_pi_detail::half_pi_lo) };
result = detail::split_pi_values<T>(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;
Expand Down
44 changes: 25 additions & 19 deletions include/boost/decimal/detail/cmath/atan2.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/detail/cmath/atan.hpp>
#include <boost/decimal/detail/cmath/fabs.hpp>
#include <boost/decimal/detail/cmath/frexp10.hpp>
#include <boost/decimal/detail/cmath/impl/split_pi.hpp>
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/detail/config.hpp>
#include <boost/decimal/numbers.hpp>
Expand All @@ -24,20 +26,6 @@ namespace decimal {

namespace detail {

namespace atan2_detail {

template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
struct pi_constants
{
static constexpr T pi_over_two = numbers::pi_v<T> / 2;
static constexpr T three_pi_over_four = 3 * numbers::pi_over_four_v<T>;
};

template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T> constexpr T pi_constants<T>::pi_over_two;
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T> constexpr T pi_constants<T>::three_pi_over_four;

} // namespace atan2_detail

template <typename T>
constexpr auto atan2_impl(const T y, const T x) noexcept
BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T)
Expand Down Expand Up @@ -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<T>::pi_over_two;
result = split_pi_values<T>(split_pi_detail::half_pi_hi);

if (signy) { result = -result; }
}
else if (fpcy == FP_INFINITE && fpcx == FP_INFINITE && signx)
{
result = atan2_detail::pi_constants<T>::three_pi_over_four;
result = split_pi_values<T>(split_pi_detail::three_quarter_pi);

if (signy)
{
Expand All @@ -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<T>;
result = split_pi_values<T>(split_pi_detail::quarter_pi_hi);

if (signy)
{
Expand All @@ -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<T>::pi_over_two;
result = split_pi_values<T>(split_pi_detail::half_pi_hi);

if (signy) { result = -result; }
}
Expand All @@ -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<T> - 1})
{
ret_val = atan_impl<true>(q, detail::unchecked_fma(-q, ax, fabs(y)) / ax);
}
}
else
{
ret_val = q < T {15575, -4} ? atan_impl<true>(q, detail::unchecked_fma(-q, ax, fabs(y)) / ax) : atan(q);
}

if (!signy && !signx)
{
Expand Down
Loading
Loading