diff --git a/include/boost/math/special_functions/ellint_d.hpp b/include/boost/math/special_functions/ellint_d.hpp index 67c3702c1..da1e87ba3 100644 --- a/include/boost/math/special_functions/ellint_d.hpp +++ b/include/boost/math/special_functions/ellint_d.hpp @@ -120,8 +120,7 @@ T ellint_d_imp(T k, const Policy& pol) if (abs(k) >= 1) { - return policies::raise_domain_error("boost::math::ellint_d<%1%>(%1%)", - "Got k = %1%, function requires |k| <= 1", k, pol); + return policies::raise_domain_error("boost::math::ellint_d<%1%>(%1%)", "Got k = %1%, function requires |k| <= 1", k, pol); } if(fabs(k) <= tools::root_epsilon()) return constants::pi() / 4; diff --git a/include/boost/math/special_functions/ellint_rc.hpp b/include/boost/math/special_functions/ellint_rc.hpp index c25f00cc5..2f9a1f8cf 100644 --- a/include/boost/math/special_functions/ellint_rc.hpp +++ b/include/boost/math/special_functions/ellint_rc.hpp @@ -40,13 +40,11 @@ T ellint_rc_imp(T x, T y, const Policy& pol) if(x < 0) { - return policies::raise_domain_error(function, - "Argument x must be non-negative but got %1%", x, pol); + return policies::raise_domain_error(function, "Argument x must be non-negative but got %1%", x, pol); } if(y == 0) { - return policies::raise_domain_error(function, - "Argument y must not be zero but got %1%", y, pol); + return policies::raise_domain_error(function, "Argument y must not be zero but got %1%", y, pol); } // for y < 0, the integral is singular, return Cauchy principal value diff --git a/include/boost/math/special_functions/ellint_rd.hpp b/include/boost/math/special_functions/ellint_rd.hpp index 08cfecb04..2a79e54ca 100644 --- a/include/boost/math/special_functions/ellint_rd.hpp +++ b/include/boost/math/special_functions/ellint_rd.hpp @@ -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(function, - "Argument x must be >= 0, but got %1%", x, pol); + return policies::raise_domain_error(function, "Argument x must be >= 0, but got %1%", x, pol); } if(y < 0) { - return policies::raise_domain_error(function, - "Argument y must be >= 0, but got %1%", y, pol); + return policies::raise_domain_error(function, "Argument y must be >= 0, but got %1%", y, pol); } if(z <= 0) { - return policies::raise_domain_error(function, - "Argument z must be > 0, but got %1%", z, pol); + return policies::raise_domain_error(function, "Argument z must be > 0, but got %1%", z, pol); } if(x + y == 0) { - return policies::raise_domain_error(function, - "At most one argument can be zero, but got, x + y = %1%", x + y, pol); + return policies::raise_domain_error(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) } diff --git a/include/boost/math/special_functions/ellint_rf.hpp b/include/boost/math/special_functions/ellint_rf.hpp index 5d5d73c25..c781ac035 100644 --- a/include/boost/math/special_functions/ellint_rf.hpp +++ b/include/boost/math/special_functions/ellint_rf.hpp @@ -40,17 +40,11 @@ namespace boost { namespace math { namespace detail{ if(x < 0 || y < 0 || z < 0) { - return policies::raise_domain_error(function, - "domain error, all arguments must be non-negative, " - "only sensible result is %1%.", - std::numeric_limits::quiet_NaN(), pol); + return policies::raise_domain_error(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", std::numeric_limits::quiet_NaN(), pol); } if(x + y == 0 || y + z == 0 || z + x == 0) { - return policies::raise_domain_error(function, - "domain error, at most one argument can be zero, " - "only sensible result is %1%.", - std::numeric_limits::quiet_NaN(), pol); + return policies::raise_domain_error(function, "domain error, at most one argument can be zero, only sensible result is %1%.", std::numeric_limits::quiet_NaN(), pol); } // // Special cases from http://dlmf.nist.gov/19.20#i diff --git a/include/boost/math/special_functions/ellint_rg.hpp b/include/boost/math/special_functions/ellint_rg.hpp index 96accdac0..051c104bc 100644 --- a/include/boost/math/special_functions/ellint_rg.hpp +++ b/include/boost/math/special_functions/ellint_rg.hpp @@ -28,10 +28,7 @@ namespace boost { namespace math { namespace detail{ if(x < 0 || y < 0 || z < 0) { - return policies::raise_domain_error(function, - "domain error, all arguments must be non-negative, " - "only sensible result is %1%.", - std::numeric_limits::quiet_NaN(), pol); + return policies::raise_domain_error(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", std::numeric_limits::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() * 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) { diff --git a/include/boost/math/special_functions/ellint_rj.hpp b/include/boost/math/special_functions/ellint_rj.hpp index 714223321..8665e7bd5 100644 --- a/include/boost/math/special_functions/ellint_rj.hpp +++ b/include/boost/math/special_functions/ellint_rj.hpp @@ -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(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(function, - "Argument x must be non-negative, but got x = %1%", x, pol); + return policies::raise_domain_error(function, "Argument x must be non-negative, but got x = %1%", x, pol); } if(y < 0) { - return policies::raise_domain_error(function, - "Argument y must be non-negative, but got y = %1%", y, pol); + return policies::raise_domain_error(function, "Argument y must be non-negative, but got y = %1%", y, pol); } if(z < 0) { - return policies::raise_domain_error(function, - "Argument z must be non-negative, but got z = %1%", z, pol); + return policies::raise_domain_error(function, "Argument z must be non-negative, but got z = %1%", z, pol); } if(p == 0) { - return policies::raise_domain_error(function, - "Argument p must not be zero, but got p = %1%", p, pol); + return policies::raise_domain_error(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(function, - "At most one argument can be zero, " - "only possible result is %1%.", std::numeric_limits::quiet_NaN(), pol); + return policies::raise_domain_error(function, "At most one argument can be zero, only possible result is %1%.", std::numeric_limits::quiet_NaN(), pol); } // for p < 0, the integral is singular, return Cauchy principal value diff --git a/test/test_carlson.hpp b/test/test_carlson.hpp index 97412096a..02a3a7df3 100644 --- a/test/test_carlson.hpp +++ b/test/test_carlson.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -219,6 +220,12 @@ void t6(T, const char* type_name) #include "ellint_rc_data.ipp" do_test_ellint_rc(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 @@ -270,6 +277,12 @@ void t12(T, const char* type_name) #include "ellint_rd_data.ipp" do_test_ellint_rd(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 @@ -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(), eps2); BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(9)/4, T(2)), boost::math::constants::ln_two(), eps2); BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(1) / 4, T(-2)), boost::math::constants::ln_two() / 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; diff --git a/test/test_ellint_d.hpp b/test/test_ellint_d.hpp index a573c92ff..de53936f1 100644 --- a/test/test_ellint_d.hpp +++ b/test/test_ellint_d.hpp @@ -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, 11> data1 = {{ + static const std::array, 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(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::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::ellint_d(T(0.5), std::numeric_limits::infinity()), std::numeric_limits::infinity()); + } + BOOST_MATH_CHECK_THROW(boost::math::ellint_d(T(1.5), T(1.0)), std::domain_error); }