From 47c2f9254cb82b8cc64b570e077ad90ce8317cce Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 26 Dec 2014 17:21:10 +0000 Subject: [PATCH] [Ellint Pi] Add some more special case handling, plus tests. --- .../boost/math/special_functions/ellint_3.hpp | 25 ++++++++++- test/test_ellint_3.cpp | 17 +++----- test/test_ellint_3.hpp | 42 ++++++++++++++++--- 3 files changed, 65 insertions(+), 19 deletions(-) diff --git a/include/boost/math/special_functions/ellint_3.hpp b/include/boost/math/special_functions/ellint_3.hpp index 3a77f8ed1..31b93d665 100644 --- a/include/boost/math/special_functions/ellint_3.hpp +++ b/include/boost/math/special_functions/ellint_3.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -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()) { // 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()) || (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 m = boost::math::round((fabs(phi) - rphi) / constants::half_pi()); int sign = 1; + if((m != 0) && (k >= 1)) + { + return policies::raise_domain_error(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()); - BOOST_ASSERT(phi > 0); + BOOST_ASSERT(phi >= 0); T x, y, z, p, t; T cosp = cos(phi); x = cosp * cosp; diff --git a/test/test_ellint_3.cpp b/test/test_ellint_3.cpp index 892a3f938..f8af9cc22 100644 --- a/test/test_ellint_3.cpp +++ b/test/test_ellint_3.cpp @@ -43,11 +43,11 @@ void expected_results() #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS if(boost::math::policies::digits >() == boost::math::policies::digits >()) { - 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. diff --git a/test/test_ellint_3.hpp b/test/test_ellint_3.hpp index 8717ee1a7..3b4f7f437 100644 --- a/test/test_ellint_3.hpp +++ b/test/test_ellint_3.hpp @@ -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, 37> data1 = {{ + static const boost::array, 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() / 4, SC_(-1E-164), SC_(0.785398163397448309615660845819875721049292349843776455243736) } }, { { -ldexp(T(1.0), -52), boost::math::constants::pi() / 4, SC_(0.9375), SC_(0.866032844934895872810905364370384153285798081574191920571016) } }, { { -ldexp(T(1.0), -52), boost::math::constants::pi() / 4, SC_(-0.9375), SC_(0.866032844934895872810905364370384153285798081574191920571016) } }, + { { std::numeric_limits::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), ldexp(T(1), -510), std::numeric_limits::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } }, + { { std::numeric_limits::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), -ldexp(T(1), -510), std::numeric_limits::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::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), ldexp(T(1), -510), std::numeric_limits::max_exponent > 600 ? SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) : SC_(-2.28835573409367516299079046268930870596307631872422530813192e-246) } }, - { { std::numeric_limits::max_exponent > 600 ? -ldexp(T(1), 604) : T(0), -ldexp(T(1), -816), -ldexp(T(1), -510), std::numeric_limits::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::digits > 50 ? boost::math::constants::pi() : 0), SC_(0.5), SC_(0.000976568747692583207373162769739923352941882247128764119469) } }, + { { SC_(20.0), -ldexp(T(1), -10) - (std::numeric_limits::digits > 50 ? boost::math::constants::pi() : 0), SC_(0.5), SC_(-0.000976568747692583207373162769739923352941882247128764119469) } }, + { { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits::digits > 50 ? SC_(1.695796326794896619231321691639751442098584699687552910487472) : SC_(0.125)), SC_(0.5), (std::numeric_limits::digits > 50 ? SC_(-27.1556829547879272134963658442519620188811279238098426583417) : SC_(0.125747519536592933239312389845380877095789502237158621392763)) } }, + { { SC_(1.0) + ldexp(T(1), -6), (std::numeric_limits::digits > 50 ? SC_(-1.695796326794896619231321691639751442098584699687552910487472) : SC_(-0.125)), SC_(0.5), (std::numeric_limits::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(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); }