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 Carlson elliptic integral coverage.

Browse Source
This commit is contained in:
jzmaddock
2024-02-25 19:40:10 +00:00
parent 7b4e2ddce0
commit 5435a91cc8
8 changed files with 76 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
include/boost/math/special_functions/ellint_d.hpp
View File

@@ -120,8 +120,7 @@ T ellint_d_imp(T k, const Policy& pol)
if (abs(k) >= 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)",
"Got k = %1%, function requires |k| <= 1", k, pol);
return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)", "Got k = %1%, function requires |k| <= 1", k, pol);
}
if(fabs(k) <= tools::root_epsilon<T>())
return constants::pi<T>() / 4;

6
include/boost/math/special_functions/ellint_rc.hpp
View File

@@ -40,13 +40,11 @@ T ellint_rc_imp(T x, T y, const Policy& pol)
if(x < 0)
{
return policies::raise_domain_error<T>(function,
"Argument x must be non-negative but got %1%", x, pol);
return policies::raise_domain_error<T>(function, "Argument x must be non-negative but got %1%", x, pol);
}
if(y == 0)
{
return policies::raise_domain_error<T>(function,
"Argument y must not be zero but got %1%", y, pol);
return policies::raise_domain_error<T>(function, "Argument y must not be zero but got %1%", y, pol);
}
// for y < 0, the integral is singular, return Cauchy principal value

16
include/boost/math/special_functions/ellint_rd.hpp
View File

@@ -38,23 +38,19 @@ T ellint_rd_imp(T x, T y, T z, const Policy& pol)
if(x < 0)
{
return policies::raise_domain_error<T>(function,
"Argument x must be >= 0, but got %1%", x, pol);
return policies::raise_domain_error<T>(function, "Argument x must be >= 0, but got %1%", x, pol);
}
if(y < 0)
{
return policies::raise_domain_error<T>(function,
"Argument y must be >= 0, but got %1%", y, pol);
return policies::raise_domain_error<T>(function, "Argument y must be >= 0, but got %1%", y, pol);
}
if(z <= 0)
{
return policies::raise_domain_error<T>(function,
"Argument z must be > 0, but got %1%", z, pol);
return policies::raise_domain_error<T>(function, "Argument z must be > 0, but got %1%", z, pol);
}
if(x + y == 0)
{
return policies::raise_domain_error<T>(function,
"At most one argument can be zero, but got, x + y = %1%", x + y, pol);
return policies::raise_domain_error<T>(function, "At most one argument can be zero, but got, x + y = %1%", x + y, pol);
}
//
// Special cases from http://dlmf.nist.gov/19.20#iv
@@ -74,14 +70,14 @@ T ellint_rd_imp(T x, T y, T z, const Policy& pol)
}
else
{
if((std::min)(x, y) / (std::max)(x, y) > T(1.3))
if((std::max)(x, y) / (std::min)(x, y) > T(1.3))
return 3 * (ellint_rc_imp(x, y, pol) - sqrt(x) / y) / (2 * (y - x));
// Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y)
}
}
if(x == y)
{
if((std::min)(x, z) / (std::max)(x, z) > T(1.3))
if((std::max)(x, z) / (std::min)(x, z) > T(1.3))
return 3 * (ellint_rc_imp(z, x, pol) - 1 / sqrt(z)) / (z - x);
// Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y)
}

10
include/boost/math/special_functions/ellint_rf.hpp
View File

@@ -40,17 +40,11 @@ namespace boost { namespace math { namespace detail{
if(x < 0 || y < 0 || z < 0)
{
return policies::raise_domain_error<T>(function,
"domain error, all arguments must be non-negative, "
"only sensible result is %1%.",
std::numeric_limits<T>::quiet_NaN(), pol);
return policies::raise_domain_error<T>(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
}
if(x + y == 0 || y + z == 0 || z + x == 0)
{
return policies::raise_domain_error<T>(function,
"domain error, at most one argument can be zero, "
"only sensible result is %1%.",
std::numeric_limits<T>::quiet_NaN(), pol);
return policies::raise_domain_error<T>(function, "domain error, at most one argument can be zero, only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
}
//
// Special cases from http://dlmf.nist.gov/19.20#i

11
include/boost/math/special_functions/ellint_rg.hpp
View File

@@ -28,10 +28,7 @@ namespace boost { namespace math { namespace detail{
if(x < 0 || y < 0 || z < 0)
{
return policies::raise_domain_error<T>(function,
"domain error, all arguments must be non-negative, "
"only sensible result is %1%.",
std::numeric_limits<T>::quiet_NaN(), pol);
return policies::raise_domain_error<T>(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
}
//
// Function is symmetric in x, y and z, but we require
@@ -73,10 +70,8 @@ namespace boost { namespace math { namespace detail{
}
else if(y == z)
{
if(x == 0)
return constants::pi<T>() * sqrt(y) / 4;
else
return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2);
BOOST_MATH_ASSERT(x > 0); // Ordering of x,y,z above takes care of x == 0 case.
return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2);
}
else if(y == 0)
{

22
include/boost/math/special_functions/ellint_rj.hpp
View File

@@ -40,11 +40,7 @@ T ellint_rc1p_imp(T y, const Policy& pol)
static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)";
if(y == -1)
{
return policies::raise_domain_error<T>(function,
"Argument y must not be zero but got %1%", y, pol);
}
BOOST_ASSERT(y != -1);
// for 1 + y < 0, the integral is singular, return Cauchy principal value
T result;
@@ -84,29 +80,23 @@ T ellint_rj_imp(T x, T y, T z, T p, const Policy& pol)
if(x < 0)
{
return policies::raise_domain_error<T>(function,
"Argument x must be non-negative, but got x = %1%", x, pol);
return policies::raise_domain_error<T>(function, "Argument x must be non-negative, but got x = %1%", x, pol);
}
if(y < 0)
{
return policies::raise_domain_error<T>(function,
"Argument y must be non-negative, but got y = %1%", y, pol);
return policies::raise_domain_error<T>(function, "Argument y must be non-negative, but got y = %1%", y, pol);
}
if(z < 0)
{
return policies::raise_domain_error<T>(function,
"Argument z must be non-negative, but got z = %1%", z, pol);
return policies::raise_domain_error<T>(function, "Argument z must be non-negative, but got z = %1%", z, pol);
}
if(p == 0)
{
return policies::raise_domain_error<T>(function,
"Argument p must not be zero, but got p = %1%", p, pol);
return policies::raise_domain_error<T>(function, "Argument p must not be zero, but got p = %1%", p, pol);
}
if(x + y == 0 || y + z == 0 || z + x == 0)
{
return policies::raise_domain_error<T>(function,
"At most one argument can be zero, "
"only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
return policies::raise_domain_error<T>(function, "At most one argument can be zero, only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
}
// for p < 0, the integral is singular, return Cauchy principal value

49
test/test_carlson.hpp
View File

@@ -9,6 +9,7 @@
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/ellint_rj.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/array.hpp>
#include <boost/random.hpp>
@@ -219,6 +220,12 @@ void t6(T, const char* type_name)
#include "ellint_rc_data.ipp"
do_test_ellint_rc<T>(ellint_rc_data, type_name, "RC: Random data");
//
// Error handling:
//
BOOST_CHECK_THROW(boost::math::ellint_rc(T(-1), T(1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_rc(T(1), T(0)), std::domain_error);
}
template <typename T>
@@ -270,6 +277,12 @@ void t12(T, const char* type_name)
#include "ellint_rd_data.ipp"
do_test_ellint_rd<T>(ellint_rd_data, type_name, "RD: Random data");
BOOST_CHECK_THROW(boost::math::ellint_rd(T(-1), T(1), T(2)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_rd(T(1), T(-1), T(2)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_rd(T(1), T(2), T(0)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_rd(T(0), T(0), T(2)), std::domain_error);
BOOST_CHECK_EQUAL(boost::math::ellint_rd(T(0.5), T(0), T(0.75)), boost::math::ellint_rd(T(0), T(0.5), T(0.75)));
}
template <typename T>
@@ -371,19 +384,55 @@ void test_spots(T val, const char* type_name)
BOOST_CHECK_CLOSE(ellint_rf(T(1), T(2), T(0)), T(1.3110287771461), tolerance);
BOOST_CHECK_CLOSE(ellint_rf(T(0.5), T(1), T(0)), T(1.8540746773014), tolerance);
BOOST_CHECK_CLOSE(ellint_rf(T(2), T(3), T(4)), T(0.58408284167715), tolerance);
BOOST_CHECK_THROW(ellint_rf(T(-1), T(1), T(1)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(1), T(-1), T(1)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(1), T(1), T(-1)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(0), T(0), T(1)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(1), T(0), T(0)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(0), T(1), T(0)), std::domain_error);
BOOST_CHECK_THROW(ellint_rf(T(0), T(0), T(0)), std::domain_error);
BOOST_CHECK_EQUAL(ellint_rf(T(0), T(2), T(2)), ellint_rf(T(2), T(0), T(2)));
BOOST_CHECK_EQUAL(ellint_rf(T(2), T(2), T(0)), ellint_rf(T(2), T(0), T(2)));
BOOST_CHECK_EQUAL(ellint_rf(T(0), T(2), T(3)), ellint_rf(T(2), T(3), T(0)));
BOOST_CHECK_EQUAL(ellint_rf(T(0), T(2), T(3)), ellint_rf(T(3), T(2), T(0)));
BOOST_CHECK_EQUAL(ellint_rf(T(0), T(2), T(3)), ellint_rf(T(3), T(0), T(2)));
BOOST_CHECK_EQUAL(ellint_rf(T(0), T(2), T(3)), ellint_rf(T(2), T(0), T(3)));
// RC:
BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(0), T(1)/4), boost::math::constants::pi<T>(), eps2);
BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(9)/4, T(2)), boost::math::constants::ln_two<T>(), eps2);
BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(1) / 4, T(-2)), boost::math::constants::ln_two<T>() / 3, eps2);
BOOST_CHECK_CLOSE_FRACTION(boost::math::detail::ellint_rc1p_imp(T(-2), boost::math::policies::policy<>()), ellint_rc(T(1), T(-1)), eps2);
BOOST_CHECK_CLOSE_FRACTION(boost::math::detail::ellint_rc1p_imp(T(-0.75), boost::math::policies::policy<>()), ellint_rc(T(1), T(0.25)), eps2);
// RJ:
BOOST_CHECK_CLOSE(ellint_rj(T(0), T(1), T(2), T(3)), T(0.77688623778582), tolerance);
BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(5)), T(0.14297579667157), tolerance);
BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-0.5)), T(0.24723819703052), tolerance);
BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-5)), T(-0.12711230042964), tolerance);
BOOST_CHECK_THROW(ellint_rj(T(-1), T(1), T(2), T(3)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(1), T(-1), T(2), T(3)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(1), T(2), T(-2), T(3)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(1), T(2), T(2), T(0)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(0), T(0), T(2), T(3)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(2), T(0), T(0), T(0)), std::domain_error);
BOOST_CHECK_THROW(ellint_rj(T(0), T(2), T(0), T(0)), std::domain_error);
// RD:
BOOST_CHECK_CLOSE(ellint_rd(T(0), T(2), T(1)), T(1.7972103521034), tolerance);
BOOST_CHECK_CLOSE(ellint_rd(T(2), T(3), T(4)), T(0.16510527294261), tolerance);
// RG
BOOST_CHECK_THROW(ellint_rg(T(-1), T(1), T(2)), std::domain_error);
BOOST_CHECK_THROW(ellint_rg(T(-1), T(-1), T(2)), std::domain_error);
BOOST_CHECK_THROW(ellint_rg(T(-1), T(2), T(-1)), std::domain_error);
BOOST_CHECK_EQUAL(ellint_rg(T(0), T(2), T(2)), ellint_rg(T(2), T(2), T(0)));
BOOST_CHECK_EQUAL(ellint_rg(T(0), T(2), T(2)), ellint_rg(T(2), T(0), T(2)));
// Sanity/consistency checks from Numerical Computation of Real or Complex
// Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
boost::mt19937 ran;

8
test/test_ellint_d.hpp
View File

@@ -84,7 +84,7 @@ void test_spots(T, const char* type_name)
BOOST_MATH_STD_USING
// Function values calculated on http://functions.wolfram.com/
// Note that Mathematica's EllipticE accepts k^2 as the second parameter.
static const std::array<std::array<T, 3>, 11> data1 = {{
static const std::array<std::array<T, 3>, 12> data1 = {{
{ { SC_(0.5), SC_(0.5), SC_(0.040348098248931543984282958654503585) } },
{{ SC_(0), SC_(0.5), SC_(0) }},
{ { SC_(1), SC_(0.5), SC_(0.28991866293419922467977188008516755) } },
@@ -94,6 +94,7 @@ void test_spots(T, const char* type_name)
{ { SC_(-10), T(0.5), SC_(-5.2996914501577855803123384771117708) } },
{ { SC_(10), SC_(-0.5), SC_(5.2996914501577855803123384771117708) } },
{ { SC_(0.125), SC_(1.5), SC_(0.000655956467603362564458676111698495009248974444516843) } },
{ { SC_(1.208925819614629174706176e24) /* 2^80 */, SC_(0.5), SC_(672000998924580555450487.42418840712)}},
}};
do_test_ellint_d2<T>(data1, type_name, "Elliptic Integral E: Mathworld Data");
@@ -120,5 +121,10 @@ void test_spots(T, const char* type_name)
BOOST_MATH_CHECK_THROW(boost::math::ellint_d(T(-1)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_d(T(1.5)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_d(T(-1.5)), std::domain_error);
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(boost::math::ellint_d(T(0.5), std::numeric_limits<T>::infinity()), std::numeric_limits<T>::infinity());
}
BOOST_MATH_CHECK_THROW(boost::math::ellint_d(T(1.5), T(1.0)), std::domain_error);
}