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

Elliptic Integrals: extend range of ellint_1/2/3.

Browse Source 
See https://github.com/boostorg/math/issues/183.
This commit is contained in:
jzmaddock
2019-05-18 19:36:22 +01:00
parent 03cb3f04b5
commit a033166f7f
7 changed files with 58 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
doc/sf/ellint_legendre.qbk
View File

@@ -51,7 +51,7 @@ Returns the incomplete elliptic integral of the first kind ['F([phi], k)]:
[equation ellint2]
Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
Requires k[super 2]sin[super 2](phi) < 1, otherwise returns the result of __domain_error.
[optional_policy]
@@ -65,7 +65,7 @@ Returns the complete elliptic integral of the first kind ['K(k)]:
[equation ellint6]
Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
Requires |k| < 1, otherwise returns the result of __domain_error.
[optional_policy]
@@ -154,7 +154,7 @@ Returns the incomplete elliptic integral of the second kind ['E([phi], k)]:
[equation ellint3]
Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
Requires k[super 2]sin[super 2](phi) < 1, otherwise returns the result of __domain_error.
[optional_policy]
@@ -168,7 +168,7 @@ Returns the complete elliptic integral of the second kind ['E(k)]:
[equation ellint7]
Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
Requires |k| < 1, otherwise returns the result of __domain_error.
[optional_policy]
@@ -257,7 +257,7 @@ Returns the incomplete elliptic integral of the third kind ['[Pi](n, [phi], k)]:
[equation ellint4]
Requires ['-1 <= k <= 1] and ['n < 1/sin[super 2]([phi])], otherwise
Requires ['k[super 2]sin[super 2](phi) < 1] and ['n < 1/sin[super 2]([phi])], otherwise
returns the result of __domain_error (outside this range the result
would be complex).
@@ -273,7 +273,7 @@ Returns the complete elliptic integral of the first kind ['[Pi](n, k)]:
[equation ellint8]
Requires ['-1 <= k <= 1] and ['n < 1], otherwise returns the
Requires ['|k| < 1] and ['n < 1], otherwise returns the
result of __domain_error (outside this range the result would be complex).
[optional_policy]

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

@@ -51,12 +51,6 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
BOOST_MATH_INSTRUMENT_VARIABLE(k);
BOOST_MATH_INSTRUMENT_VARIABLE(function);
if (abs(k) > 1)
{
return policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, pol);
}
bool invert = false;
if(phi < 0)
{
@@ -104,6 +98,11 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
}
T sinp = sin(rphi);
sinp *= sinp;
if (sinp * k * k >= 1)
{
return policies::raise_domain_error<T>(function,
"Got k^2 * sin^2(phi) = %1%, but the function requires this < 1", sinp * k * k, pol);
}
T cosp = cos(rphi);
cosp *= cosp;
BOOST_MATH_INSTRUMENT_VARIABLE(sinp);

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

@@ -49,6 +49,9 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
using namespace boost::math::constants;
bool invert = false;
if (phi == 0)
return 0;
if(phi < 0)
{
phi = fabs(phi);
@@ -95,11 +98,7 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
rphi = constants::half_pi<T>() - rphi;
}
T k2 = k * k;
if(k2 > 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
}
else if(rphi < tools::root_epsilon<T>())
if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon<T>() * fabs(rphi))
{
// See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/
result = s * rphi;
@@ -108,6 +107,10 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
{
// http://dlmf.nist.gov/19.25#E10
T sinp = sin(rphi);
if (k2 * sinp * sinp >= 1)
{
return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
}
T cosp = cos(rphi);
T c = 1 / (sinp * sinp);
T cm1 = cosp * cosp / (sinp * sinp); // c - 1

12
include/boost/math/special_functions/ellint_3.hpp
View File

@@ -49,15 +49,15 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)";
if(abs(k) > 1)
{
return policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, pol);
}
T sphi = sin(fabs(phi));
T result = 0;
if (k * k * sphi * sphi > 1)
{
return policies::raise_domain_error<T>(function,
"Got k = %1%, function requires |k| <= 1", k, pol);
}
// Special cases first:
if(v == 0)
{
@@ -154,7 +154,7 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
}
}
if(v < 0)
if((v < 0) && fabs(k) <= 1)
{
//
// If we don't shift to 0 <= v <= 1 we get

16
test/test_ellint_1.hpp
View File

@@ -84,7 +84,7 @@ void test_spots(T, const char* type_name)
{
// Function values calculated on http://functions.wolfram.com/
// Note that Mathematica's EllipticF accepts k^2 as the second parameter.
static const boost::array<boost::array<typename table_type<T>::type, 3>, 19> data1 = {{
static const boost::array<boost::array<typename table_type<T>::type, 3>, 22> data1 = {{
{{ SC_(0.0), SC_(0.0), SC_(0.0) }},
{{ SC_(-10.0), SC_(0.0), SC_(-10.0) }},
{{ SC_(-1.0), SC_(-1.0), SC_(-1.2261911708835170708130609674719067527242483502207) }},
@@ -104,6 +104,10 @@ void test_spots(T, const char* type_name)
{{ SC_(-4.0), SC_(0.5), SC_(-4.2543274975235836861894752787874633017836785640477) }},
{{ SC_(-6.0), SC_(0.5), SC_(-6.4588766202317746302999080620490579800463614807916) }},
{{ SC_(-10.0), SC_(0.5), SC_(-10.697409951222544858346795279378531495869386960090) }},
// Some values where k is > 1:
{{ static_cast<typename table_type<T>::type>(2)/13, static_cast<typename table_type<T>::type>(15)/13, SC_(0.154661869446904722070471580919758948531148566762183486996920)}},
{{ static_cast<typename table_type<T>::type>(2)/13, static_cast<typename table_type<T>::type>(19)/13, SC_(0.155166467455029577314314021156113481657713115640002027219)}},
{{ static_cast<typename table_type<T>::type>(2)/13, static_cast<typename table_type<T>::type>(32)/13, SC_(0.15776272074094290829870142225970052217542486917945444918)}},
}};
do_test_ellint_f<T>(data1, type_name, "Elliptic Integral F: Mathworld Data");
@@ -131,5 +135,15 @@ void test_spots(T, const char* type_name)
#include "ellint_k_data.ipp"
do_test_ellint_k<T>(ellint_k_data, type_name, "Elliptic Integral K: Random Data");
//
// Test error handling:
//
BOOST_CHECK_GE(boost::math::ellint_1(T(1)), boost::math::tools::max_value<T>());
BOOST_CHECK_GE(boost::math::ellint_1(T(-1)), boost::math::tools::max_value<T>());
BOOST_CHECK_THROW(boost::math::ellint_1(T(1.0001)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_1(T(-1.0001)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_1(T(2.2), T(0.5)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_1(T(-2.2), T(0.5)), std::domain_error);
}

11
test/test_ellint_2.hpp
View File

File diff suppressed because one or more lines are too long

10
test/test_ellint_3.hpp
View File

@@ -87,7 +87,7 @@ void test_spots(T, const char* type_name)
{
BOOST_MATH_STD_USING
// function values calculated on http://functions.wolfram.com/
static const boost::array<boost::array<typename table_type<T>::type, 4>, 65> data1 = {{
static const boost::array<boost::array<typename table_type<T>::type, 4>, 70> data1 = {{
{{ SC_(1.0), SC_(-1.0), SC_(0.0), SC_(-1.557407724654902230506974807458360173087) }},
{{ SC_(0.0), SC_(-4.0), SC_(0.4), SC_(-4.153623371196831087495427530365430979011) }},
{{ SC_(0.0), SC_(8.0), SC_(-0.6), SC_(8.935930619078575123490612395578518914416) }},
@@ -166,6 +166,12 @@ void test_spots(T, const char* type_name)
{ { SC_(-1.0), SC_(0.0), SC_(0.5), SC_(0.0) } },
{ { SC_(100.0), SC_(0.0), SC_(0.5), SC_(0.0) } },
{ { SC_(-100.0), SC_(0.0), SC_(0.5), SC_(0.0) } },
// cases where |k| > 1:
{{ SC_(1.015625), SC_(0.125), SC_(2.5), SC_(0.12780840244950364924992513047014683502850891674019968726)}},
{{ SC_(1.015625), SC_(0.125), SC_(4.5), SC_(0.133468341825478826487248053944106425799790398297092314041)}},
{{ SC_(-1.015625), SC_(0.125), SC_(4.5), SC_(0.131998693459801470974303284674812630138938285546080214149)}},
{{ SC_(-1.015625), SC_(-0.125), SC_(4.5), SC_(-0.13199869345980147097430328467481263013893828554608021414)}},
{{ SC_(-1.015625), SC_(-0.125), SC_(1.5), SC_(-0.12508193549646497011359938978158598001227028251706394704)}},
} };
do_test_ellint_pi3<T>(data1, type_name, "Elliptic Integral PI: Mathworld Data");
@@ -207,7 +213,7 @@ void test_spots(T, const char* type_name)
do_test_ellint_pi2<T>(ellint_pi2_data, type_name, "Complete Elliptic Integral PI: Random Data");
// Special cases, exceptions etc:
BOOST_MATH_CHECK_THROW(boost::math::ellint_3(T(1.0001), T(-1), T(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_3(T(2.1), T(-1), T(0.5)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_3(T(0.5), T(20), T(1.5)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_3(T(1.0001), T(-1)), std::domain_error);
BOOST_MATH_CHECK_THROW(boost::math::ellint_3(T(0.5), T(1)), std::domain_error);