boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Improve Heuman Lambda precision:

Browse Source 
Make sure we pass 1 - k^2 down through the call stack otherwise we get cancellation errors.
This commit is contained in:
jzmaddock
2024-04-25 19:53:01 +01:00
parent 28a5ffb1b4
commit cfab531eb9
3 changed files with 39 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
include/boost/math/special_functions/ellint_1.hpp
View File

@@ -41,10 +41,12 @@ template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&);
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&);
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, T one_minus_k2);
// Elliptic integral (Legendre form) of the first kind
template <typename T, typename Policy>
T ellint_f_imp(T phi, T k, const Policy& pol)
T ellint_f_imp(T phi, T k, const Policy& pol, T one_minus_k2)
{
BOOST_MATH_STD_USING
using namespace boost::math::tools;
@@ -80,7 +82,7 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
> precision_tag_type;
result = 2 * phi * ellint_k_imp(k, pol, precision_tag_type()) / constants::pi<T>();
result = 2 * phi * ellint_k_imp(k, pol, one_minus_k2) / constants::pi<T>();
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
else
@@ -121,31 +123,40 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
BOOST_MATH_ASSERT(rphi != 0); // precondition, can't be true if sin(rphi) != 0.
//
// Use http://dlmf.nist.gov/19.25#E5, note that
// c-1 simplifies to cot^2(rphi) which avoid cancellation:
// c-1 simplifies to cot^2(rphi) which avoids cancellation.
// Likewise c - k^2 is the same as (c - 1) + (1 - k^2).
//
T c = 1 / sinp;
result = static_cast<T>(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol));
T c_minus_one = cosp / sinp;
T cross = fabs(c / (k * k));
T arg2;
if ((cross > 0.9) && (cross < 1.1))
arg2 = c_minus_one + one_minus_k2;
else
arg2 = c - k * k;
result = static_cast<T>(s * ellint_rf_imp(c_minus_one, arg2, c, pol));
}
else
result = s * sin(rphi);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(m != 0)
{
typedef std::integral_constant<int,
std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
> precision_tag_type;
result += m * ellint_k_imp(k, pol, precision_tag_type());
result += m * ellint_k_imp(k, pol, one_minus_k2);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
return invert ? T(-result) : result;
}
template <typename T, typename Policy>
inline T ellint_f_imp(T phi, T k, const Policy& pol)
{
return ellint_f_imp(phi, k, pol, T(1 - k * k));
}
// Complete elliptic integral (Legendre form) of the first kind
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
T ellint_k_imp(T k, const Policy& pol, T one_minus_k2)
{
BOOST_MATH_STD_USING
using namespace boost::math::tools;
@@ -162,12 +173,16 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
}
T x = 0;
T y = 1 - k * k;
T z = 1;
T value = ellint_rf_imp(x, y, z, pol);
T value = ellint_rf_imp(x, one_minus_k2, z, pol);
return value;
}
template <typename T, typename Policy>
inline T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
{
return ellint_k_imp(k, pol, T(1 - k * k));
}
//
// Special versions for double and 80-bit long double precision,

4
include/boost/math/special_functions/heuman_lambda.hpp
View File

@@ -63,8 +63,8 @@ T heuman_lambda_imp(T phi, T k, const Policy& pol)
return policies::raise_domain_error<T>(function, "When 1-k^2 == 1 then phi must be < Pi/2, but got phi = %1%", phi, pol);
}
else
ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol, precision_tag_type());
result = ratio + ellint_k_imp(k, pol, precision_tag_type()) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi<T>();
ratio = ellint_f_imp(phi, rkp, pol, k2) / ellint_k_imp(rkp, pol, k2);
result = ratio + ellint_k_imp(k, pol, precision_tag_type()) * jacobi_zeta_imp(phi, rkp, pol, k2) / constants::half_pi<T>();
}
return result;
}

13
include/boost/math/special_functions/jacobi_zeta.hpp
View File

@@ -27,7 +27,7 @@ namespace detail{
// Elliptic integral - Jacobi Zeta
template <typename T, typename Policy>
T jacobi_zeta_imp(T phi, T k, const Policy& pol)
T jacobi_zeta_imp(T phi, T k, const Policy& pol, T kp)
{
BOOST_MATH_STD_USING
using namespace boost::math::tools;
@@ -44,8 +44,9 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol)
T sinp = sin(phi);
T cosp = cos(phi);
T s2 = sinp * sinp;
T c2 = cosp * cosp;
T one_minus_ks2 = kp + c2 - kp * c2;
T k2 = k * k;
T kp = 1 - k2;
if(k == 1)
result = sinp * (boost::math::sign)(cosp); // We get here by simplifying JacobiZeta[w, 1] in Mathematica, and the fact that 0 <= phi.
else
@@ -54,11 +55,15 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol)
std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
> precision_tag_type;
result = k2 * sinp * cosp * sqrt(1 - k2 * s2) * ellint_rj_imp(T(0), kp, T(1), T(1 - k2 * s2), pol) / (3 * ellint_k_imp(k, pol, precision_tag_type()));
result = k2 * sinp * cosp * sqrt(one_minus_ks2) * ellint_rj_imp(T(0), kp, T(1), one_minus_ks2, pol) / (3 * ellint_k_imp(k, pol, kp));
}
return invert ? T(-result) : result;
}
template <typename T, typename Policy>
inline T jacobi_zeta_imp(T phi, T k, const Policy& pol)
{
return jacobi_zeta_imp(phi, k, pol, T(1 - k * k));
}
} // detail
template <class T1, class T2, class Policy>