Skip to content

Commit

Permalink
Merge pull request #666 from cppalliance/ellint_coverage
Browse files Browse the repository at this point in the history
  • Loading branch information
@ckormanyos
ckormanyos authored Jun 20, 2024
2 parents d9dc532 + 43525f2 commit 7aceda2
Show file tree
Hide file tree
Showing 10 changed files with 460 additions and 120 deletions.
14 changes: 8 additions & 6 deletions include/boost/decimal/detail/cmath/cos.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
// Copyright 2023 Matt Borland
// Copyright 2023 Christopher Kormanyos
// Copyright 2023 - 2024 Matt Borland
// Copyright 2023 - 2024 Christopher Kormanyos
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

Expand Down Expand Up @@ -58,16 +58,18 @@ constexpr auto cos_impl(T x) noexcept
// | 2 | -sin(r) | -cos(r) | sin(r)/cos(r) |
// | 3 | -cos(r) | sin(r) | -cos(r)/sin(r) |

constexpr T my_pi_half { numbers::pi_v<T> / 2 };
const T two_x = x * 2;

const auto k = static_cast<unsigned>(x / my_pi_half);
const auto k = static_cast<unsigned>(two_x / numbers::pi_v<T>);
const auto n = static_cast<unsigned>(k % static_cast<unsigned>(UINT8_C(4)));

auto r = x - (my_pi_half * k);
const T two_r { two_x - (numbers::pi_v<T> * k) };

T r { two_r / 2 };

constexpr T half { 5, -1 };

const bool do_scaling { x > half };
const bool do_scaling { r > half };

if(do_scaling)
{
Expand Down
37 changes: 19 additions & 18 deletions include/boost/decimal/detail/cmath/ellint_1.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,26 +67,20 @@ constexpr auto ellint_1_impl(T m, T phi) noexcept
{
constexpr T small_phi_limit =
std::numeric_limits<T>::digits10 < 10 ? T { 1, -2 }
: std::numeric_limits<T>::digits10 < 20 ? T { 1, -2 }
: T { 1, -4 };
: std::numeric_limits<T>::digits10 < 20 ? T { 1, -3 }
: T { 1, -5 };

if (phi < small_phi_limit)
{
// PadeApproximant[EllipticF[phi, m2], {phi, 0, {6, 5}}]
// PadeApproximant[EllipticF[phi, m2], {phi, 0, {4, 3}}]
// FullSimplify[%]
// HornerForm[Numerator[Out[2]], phi]
// HornerForm[Denominator[Out[2]], phi]

// Then further collect factors of m2 in the numerator
// and denominator manually with HornerForm[poly, m],
// and finally get close to C-language form with CForm[].

const T phi_sq { phi * phi };

const T m2 { (m * m) };
const T m2 { m * m };

const T top { phi*(-2661120 + m2*(-4354560 + 8300880*m2) + phi_sq*(-349440 + m2*(-262080 + (6683040 - 6460020*m2)*m2) + (-19200 + m2*(-263296 + m2*(997488 + m2*(-1325988 + 621441*m2))))*phi_sq)) };
const T bot { -2661120 + m2*(-4354560 + 8300880*m2) + phi_sq*(-349440 + m2*(181440 + (7408800 - 7843500*m2)*m2) + (-19200 + m2*(-293760 + m2*(1021680 + m2*(-1957500 + 1306125*m2))))*phi_sq) };
const T top { phi * (-60 + (-12 + 17 * m2) * phi_sq) };
const T bot { -60 + 3 * (-4 + 9 * m2) * phi_sq };

result = top / bot;
}
Expand All @@ -98,9 +92,9 @@ constexpr auto ellint_1_impl(T m, T phi) noexcept

T phi_scaled = phi - (k_pi * numbers::pi_v<T>);

const bool b_neg { phi_scaled > my_pi_half };
const bool b_medium_phi_scale_and_negate { phi_scaled > my_pi_half };

if(b_neg)
if(b_medium_phi_scale_and_negate)
{
++k_pi;

Expand All @@ -109,14 +103,21 @@ constexpr auto ellint_1_impl(T m, T phi) noexcept

T Km { };

detail::ellint_detail::elliptic_series::agm(phi_scaled, m * m, result, Km);
const T m2 { m * m };

const bool m2_is_one { m2 == one };

detail::ellint_detail::elliptic_series::agm(phi_scaled, m2, result, Km, static_cast<T*>(nullptr), static_cast<T*>(nullptr));

if(b_neg)
if(b_medium_phi_scale_and_negate)
{
result = -result;
}

result += ((k_pi * Km) * 2);
if(!m2_is_one)
{
result += ((k_pi * Km) * 2);
}
}
}

Expand Down Expand Up @@ -151,7 +152,7 @@ constexpr auto comp_ellint_1_impl(T m) noexcept

T Fpm { };

detail::ellint_detail::elliptic_series::agm(zero, m * m, Fpm, result);
detail::ellint_detail::elliptic_series::agm(zero, m * m, Fpm, result, static_cast<T*>(nullptr), static_cast<T*>(nullptr));
}

return result;
Expand Down
30 changes: 15 additions & 15 deletions include/boost/decimal/detail/cmath/ellint_2.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,22 +65,19 @@ constexpr auto ellint_2_impl(T m, T phi) noexcept
}
else
{
constexpr int small_phi_order
{
std::numeric_limits<T>::digits10 < 10 ? 2
: std::numeric_limits<T>::digits10 < 20 ? 4
: 8
};
constexpr T small_phi_limit =
std::numeric_limits<T>::digits10 < 10 ? T { 1, -2 }
: std::numeric_limits<T>::digits10 < 20 ? T { 1, -3 }
: T { 1, -5 };

if (phi < T { 1, -small_phi_order })
if (phi < small_phi_limit)
{
// PadeApproximant[EllipticE[phi, m2], {phi, 0, {4, 3}}]
// FullSimplify[%]
// Then manually edit the interior field to regain HornerForm[poly, phi].

const T phi_sq { phi * phi };

const T m2 { (!signbit(m)) ? (m * m) : -(m * m) };
const T m2 { m * m };

const T top { phi * (-60 + (-12 + 19 * m2) * phi_sq) };
const T bot { -60 + 3 * (-4 + 3 * m2) * phi_sq };
Expand All @@ -91,12 +88,13 @@ constexpr auto ellint_2_impl(T m, T phi) noexcept
{
constexpr T my_pi_half { numbers::pi_v<T> / 2 };

T k_pi = static_cast<int>(phi / numbers::pi_v<T>);
int k_pi = static_cast<int>(phi / numbers::pi_v<T>);

T phi_scaled = phi - (k_pi * numbers::pi_v<T>);

const bool b_neg { phi_scaled > my_pi_half };
const bool b_medium_phi_scale_and_negate { phi_scaled > my_pi_half };

if(b_neg)
if(b_medium_phi_scale_and_negate)
{
++k_pi;

Expand All @@ -107,9 +105,11 @@ constexpr auto ellint_2_impl(T m, T phi) noexcept
T Km { };
T Em { };

detail::ellint_detail::elliptic_series::agm(phi_scaled, m * m, Fpm, Km, &Em, &result);
const T m2 { m * m };

detail::ellint_detail::elliptic_series::agm(phi_scaled, m2, Fpm, Km, &Em, &result);

if(b_neg)
if(b_medium_phi_scale_and_negate)
{
result = -result;
}
Expand Down Expand Up @@ -150,7 +150,7 @@ constexpr auto comp_ellint_2_impl(T m) noexcept
T Fpm { };
T Km { };

detail::ellint_detail::elliptic_series::agm(zero, m * m, Fpm, Km, &result);
detail::ellint_detail::elliptic_series::agm(zero, m * m, Fpm, Km, &result, static_cast<T*>(nullptr));
}

return result;
Expand Down
38 changes: 26 additions & 12 deletions include/boost/decimal/detail/cmath/impl/ellint_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,11 +28,17 @@ constexpr auto agm(T phi,
T m2,
T& Fpm,
T& Km,
T* const pEm = nullptr,
T* const pEpm = nullptr) noexcept
T* pEm,
T* pEpm) noexcept
BOOST_DECIMAL_REQUIRES_RETURN(detail::is_decimal_floating_point_v, T, void)
{
// Use the AGM algorithm implemented in e_float:
// See Chapter/Section 19.8(i) Elliptic Integrals: Quadratic Transformations:
// Arithmetic Geometric Mean (AGM), pp. 492-493 of:
// F.W.J. Olver et al., NIST Handbook of Mathematical Functions,
// Cambridge University Press. The same can also be found
// online at https://dlmf.nist.gov/19.8

// In particular, use the AGM algorithm implemented in e_float:
// C.M. Kormanyos, "Algorithm 910: A Portable C++ Multiple-Precision
// System for Special-Function Calculations", ACM TOMS (37) 4, February 2011.

Expand Down Expand Up @@ -63,18 +69,28 @@ constexpr auto agm(T phi,

constexpr T one { 1 };

const bool has_e { ((pEm != nullptr) || (pEpm != nullptr)) };

if(m2 == one)
{
if(pEm != nullptr) { *pEm = one; }

Km = std::numeric_limits<T>::quiet_NaN();

const T sp = sin(phi);
const T sp { sin(phi) };

Fpm = phi_is_pi_half ? std::numeric_limits<T>::quiet_NaN() : log((one + sp) / (one - sp)) / 2;

Fpm = phi_is_pi_half ? std::numeric_limits<T>::quiet_NaN()
: log((one + sp) / (one - sp)) / 2;
if(has_e)
{
if(pEm != nullptr)
{
*pEm = one;
}

if(pEpm != nullptr) { *pEpm = phi_is_pi_half ? one : sp; }
if(pEpm != nullptr)
{
*pEpm = phi_is_pi_half ? one : sp;
}
}
}
else
{
Expand All @@ -89,12 +105,10 @@ constexpr auto agm(T phi,

T an { };

const bool has_e { ((pEm != nullptr) || (pEpm != nullptr)) };

T cn_2ncn_inner_prod = (has_e ? m2 / 2 : zero);
T sin_phi_n_cn_inner_prod = zero;

const T break_check { 9, -1 - (std::numeric_limits<T>::digits / 2) };
const T break_check { b0 * T { 1, -1 - (std::numeric_limits<T>::digits / 2) } };

for(int n = 0; n < std::numeric_limits<std::uint32_t>::digits; ++n)
{
Expand Down
8 changes: 4 additions & 4 deletions include/boost/decimal/detail/cmath/impl/sin_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,8 +170,8 @@ constexpr auto sin_series_expansion<decimal128>(decimal128 x) noexcept

const decimal128 x2 { x * x };

const decimal128 top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) };
const decimal128 bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) };
const decimal128 top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) };
const decimal128 bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) };

const decimal128 result { top / bot };

Expand Down Expand Up @@ -207,8 +207,8 @@ constexpr auto sin_series_expansion<decimal128_fast>(decimal128_fast x) noexcept

const decimal128_fast x2 { x * x };

const decimal128_fast top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) };
const decimal128_fast bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) };
const decimal128_fast top { x * (c0 + x2 * (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))))) };
const decimal128_fast bot { c0 + x2 * (d1 + x2 * (d2 + x2 * (d3 + x2 * (d4 + x2 * (d5 + x2 * d6))))) };

const decimal128_fast result { top / bot };

Expand Down
69 changes: 52 additions & 17 deletions include/boost/decimal/detail/cmath/sin.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,33 +58,68 @@ constexpr auto sin_impl(T x) noexcept
// | 2 | -sin(r) | -cos(r) | sin(r)/cos(r) |
// | 3 | -cos(r) | sin(r) | -cos(r)/sin(r) |

constexpr T my_pi_half { numbers::pi_v<T> / 2 };
const T two_x = x * 2;

const auto k = static_cast<unsigned>(x / my_pi_half);
const auto k = static_cast<unsigned>(two_x / numbers::pi_v<T>);
const auto n = static_cast<unsigned>(k % static_cast<unsigned>(UINT8_C(4)));

auto r = x - (my_pi_half * k);
const T two_r { two_x - (numbers::pi_v<T> * k) };

constexpr T half { 5, -1 };
T r { two_r / 2 };

const bool do_scaling { x > half };
constexpr T one { 1 };

if(do_scaling)
{
// Reduce the argument with factors of three.
r /= static_cast<unsigned>(UINT8_C(3));
}
bool do_scaling { two_r > one };

constexpr T sqrt_epsilon { sqrt(std::numeric_limits<T>::epsilon()) };

switch(n)
{
case static_cast<unsigned>(UINT8_C(1)):
case static_cast<unsigned>(UINT8_C(3)):
result = detail::cos_series_expansion(r);
case static_cast<unsigned>(UINT8_C(1)):
case static_cast<unsigned>(UINT8_C(3)):
{
const T d2r { numbers::pi_v<T> - two_r };

if (d2r < sqrt_epsilon)
{
result = d2r * (one - (d2r * d2r) / 12) / 2;

do_scaling = false;
}
else
{
if(do_scaling)
{
// Reduce the argument with one single factor of three.
r /= static_cast<unsigned>(UINT8_C(3));
}

result = detail::cos_series_expansion(r);
}
}
break;
case static_cast<unsigned>(UINT8_C(0)):
case static_cast<unsigned>(UINT8_C(2)):
default:
result = detail::sin_series_expansion(r);

case static_cast<unsigned>(UINT8_C(0)):
case static_cast<unsigned>(UINT8_C(2)):
default:
{
if (two_r < sqrt_epsilon)
{
// Normal[Series[Sin[x/2], {x, 0, 3}]]

result = (two_r * (one - (two_r * two_r) / 24)) / 2;
}
else
{
if(do_scaling)
{
// Reduce the argument with one single factor of three.
r /= static_cast<unsigned>(UINT8_C(3));
}

result = detail::sin_series_expansion(r);
}
}
break;
}

Expand Down
Loading

0 comments on commit 7aceda2

Please sign in to comment.