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

[Ellint Pi] Add some more special case handling, plus tests.

Browse Source
This commit is contained in:
jzmaddock
2014-12-26 17:21:10 +00:00
parent 18dd27295d
commit 47c2f9254c
3 changed files with 65 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -24,6 +24,7 @@
#include <boost/math/special_functions/ellint_1.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/atanh.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/workaround.hpp>
@@ -55,6 +56,7 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
}
T sphi = sin(fabs(phi));
T result = 0;
if(v > 1 / (sphi * sphi))
{
@@ -69,6 +71,15 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
// A&S 17.7.18 & 19
return (k == 0) ? phi : ellint_f_imp(phi, k, pol);
}
if(v == 1)
{
// http://functions.wolfram.com/08.06.03.0008.01
T m = k * k;
result = sqrt(1 - m * sphi * sphi) * tan(phi) - ellint_e_imp(phi, k, pol);
result /= 1 - m;
result += ellint_f_imp(phi, k, pol);
return result;
}
if(phi == constants::half_pi<T>())
{
// Have to filter this case out before the next
@@ -80,7 +91,6 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
//
return ellint_pi_imp(v, k, vc, pol);
}
T result = 0;
if((phi > constants::half_pi<T>()) || (phi < 0))
{
// Carlson's algorithm works only for |phi| <= pi/2,
@@ -109,6 +119,10 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
T rphi = boost::math::tools::fmod_workaround(T(fabs(phi)), T(constants::half_pi<T>()));
T m = boost::math::round((fabs(phi) - rphi) / constants::half_pi<T>());
int sign = 1;
if((m != 0) && (k >= 1))
{
return policies::raise_domain_error<T>(function, "Got k=1 and phi=%1% but the result is complex in that domain", phi, pol);
}
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
{
m += 1;
@@ -193,6 +207,13 @@ 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(k == 1)
{
// See http://functions.wolfram.com/08.06.03.0013.01
result = sqrt(v) * atanh(sqrt(v) * sin(phi)) - log(1 / cos(phi) + tan(phi));
result /= v - 1;
return result;
}
#if 0 // disabled but retained for future reference: see below.
if(v > 1)
{
@@ -231,7 +252,7 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
// normalised above.
//
BOOST_ASSERT(fabs(phi) < constants::half_pi<T>());
BOOST_ASSERT(phi > 0);
BOOST_ASSERT(phi >= 0);
T x, y, z, p, t;
T cosp = cos(phi);
x = cosp * cosp;

17
test/test_ellint_3.cpp
View File

@@ -43,11 +43,11 @@ void expected_results()
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
{
largest_type = "(long\\s+)?double";
largest_type = "(long\\s+)?double|real_concept";
}
else
{
largest_type = "long double";
largest_type = "long double|real_concept";
}
#else
largest_type = "(long\\s+)?double";
@@ -67,9 +67,9 @@ void expected_results()
".*", // compiler
".*", // stdlib
".*", // platform
"real_concept", // test type(s)
".*Large.*", // test data group
".*", 50, 20); // test function
largest_type, // test type(s)
".*Mathworld.*", // test data group
".*", 100, 50); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
@@ -77,13 +77,6 @@ void expected_results()
largest_type, // test type(s)
".*", // test data group
".*", 15, 8); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
"real_concept", // test type(s)
".*", // test data group
".*", 15, 8); // test function
//
// Finish off by printing out the compiler/stdlib/platform names,
// we do this to make it easier to mark up expected error rates.

42
test/test_ellint_3.hpp
View File

@@ -81,7 +81,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<T, 4>, 37> data1 = {{
static const boost::array<boost::array<T, 4>, 61> 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) }},
@@ -112,15 +112,44 @@ void test_spots(T, const char* type_name)
{ { SC_(-1E-170), boost::math::constants::pi<T>() / 4, SC_(-1E-164), SC_(0.785398163397448309615660845819875721049292349843776455243736) } },
{ { -ldexp(T(1.0), -52), boost::math::constants::pi<T>() / 4, SC_(0.9375), SC_(0.866032844934895872810905364370384153285798081574191920571016) } },
{ { -ldexp(T(1.0), -52), boost::math::constants::pi<T>() / 4, SC_(-0.9375), SC_(0.866032844934895872810905364370384153285798081574191920571016) } },
{ { std::numeric_limits<T>::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), ldexp(T(1), -510), std::numeric_limits<T>::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } },
{ { std::numeric_limits<T>::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), -ldexp(T(1), -510), std::numeric_limits<T>::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), -ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
// Special cases where k = 0:
{ { SC_(0.5), SC_(1.0), SC_(0.0), SC_(1.17881507892743738986863357869566288974084658835353613038547) } },
{ { SC_(-0.5), SC_(1.0), SC_(0.0), SC_(0.888286691263535380266337576823783210424994266596287990733270) } },
{ { SC_(0.5), SC_(-1.0), SC_(0.0), SC_(-1.17881507892743738986863357869566288974084658835353613038547) } },
{ { SC_(-0.5), SC_(-1.0), SC_(0.0), SC_(-0.888286691263535380266337576823783210424994266596287990733270) } },
// k == 0 and phi > pi/2:
{ { SC_(0.5), SC_(2.0), SC_(0.0), SC_(3.03379730757207227653600089552126882582809860566558143254794) } },
{ { SC_(-0.5), SC_(2.0), SC_(0.0), SC_(1.57453655812023739911111328195028658229986230310938753315640) } },
{ { SC_(0.5), SC_(-2.0), SC_(0.0), SC_(-3.03379730757207227653600089552126882582809860566558143254794) } },
{ { SC_(-0.5), SC_(-2.0), SC_(0.0), SC_(-1.57453655812023739911111328195028658229986230310938753315640) } },
{ { std::numeric_limits<T>::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), ldexp(T(1), -510), std::numeric_limits<T>::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } },
{ { std::numeric_limits<T>::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), -ldexp(T(1), -510), std::numeric_limits<T>::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
{ { -ldexp(T(1), -622), -ldexp(T(1), -800), -ldexp(T(1), -132), SC_(-1.49969681389563095481764443762806535353996169623910829793733e-241) } },
// Special cases where k = 1:
{ { SC_(0.5), SC_(1.0), SC_(1.0), SC_(1.4830998734200773326887632776553375078936815318419194718912351) } },
{ { SC_(-0.5), SC_(1.0), SC_(1.0), SC_(1.07048347329000030842347009377117215811122412769516781788253) } },
{ { SC_(0.5), SC_(-1.0), SC_(1.0), SC_(-1.4830998734200773326887632776553375078936815318419194718912) } },
{ { SC_(-0.5), SC_(-1.0), SC_(1.0), SC_(-1.07048347329000030842347009377117215811122412769516781788253) } },
// special cases where v = 1:
{ { SC_(1.0), SC_(0.5), SC_(0.5), SC_(0.55225234291197632914658859230278152249148960801635386133501) } },
{ { SC_(1.0), SC_(-0.5), SC_(0.5), SC_(-0.55225234291197632914658859230278152249148960801635386133501) } },
{ { SC_(1.0), SC_(2.0), SC_(0.5), SC_(-2.87534521505997989921579168327307068134740792740155171368532) } },
{ { SC_(1.0), SC_(-2.0), SC_(0.5), SC_(2.87534521505997989921579168327307068134740792740155171368532) } },
{ { SC_(1.0), SC_(2.0), ldexp(T(1), -200), SC_(-2.18503986326151899164330610231368254343201774622766316456295) } },
{ { SC_(1.0), SC_(-2.0), ldexp(T(1), -200), SC_(2.18503986326151899164330610231368254343201774622766316456295) } },
{ { SC_(1.0), ldexp(T(1.0), -150), ldexp(T(1), -200), SC_(7.006492321624085354618647916449580656401309709382578858e-46) } },
{ { SC_(1.0), -ldexp(T(1.0), -150), ldexp(T(1), -200), SC_(-7.006492321624085354618647916449580656401309709382578858e-46) } },
// Previosuly unsupported region with v > 1 and |phi| > PI/2:
{ { SC_(20.0), ldexp(T(1), -10) + (std::numeric_limits<T>::digits > 50 ? boost::math::constants::pi<T>() : 0), SC_(0.5), SC_(0.000976568747692583207373162769739923352941882247128764119469) } },
{ { SC_(20.0), -ldexp(T(1), -10) - (std::numeric_limits<T>::digits > 50 ? boost::math::constants::pi<T>() : 0), SC_(0.5), SC_(-0.000976568747692583207373162769739923352941882247128764119469) } },
{ { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits<T>::digits > 50 ? SC_(1.695796326794896619231321691639751442098584699687552910487472) : SC_(0.125)), SC_(0.5), (std::numeric_limits<T>::digits > 50 ? SC_(-27.1556829547879272134963658442519620188811279238098426583417) : SC_(0.125747519536592933239312389845380877095789502237158621392763)) } },
{ { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits<T>::digits > 50 ? SC_(-1.695796326794896619231321691639751442098584699687552910487472) : SC_(-0.125)), SC_(0.5), (std::numeric_limits<T>::digits > 50 ? SC_(27.1556829547879272134963658442519620188811279238098426583417) : SC_(-0.125747519536592933239312389845380877095789502237158621392763)) } },
// Phi = 0:
{ { SC_(1.0), SC_(0.0), SC_(0.5), SC_(0.0) } },
{ { 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) } },
} };
do_test_ellint_pi3<T>(data1, type_name, "Elliptic Integral PI: Mathworld Data");
@@ -167,4 +196,7 @@ void test_spots(T, const char* type_name)
BOOST_CHECK_THROW(boost::math::ellint_3(T(1.0001), T(-1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_3(T(0.5), T(1)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_3(T(0.5), T(2)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_3(T(1), T(0.5), T(2)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_3(T(1), T(-0.5), T(2)), std::domain_error);
BOOST_CHECK_THROW(boost::math::ellint_3(T(1), T(-0.5), T(-2)), std::domain_error);
}