diff --git a/include/boost/math/special_functions/ellint_rf.hpp b/include/boost/math/special_functions/ellint_rf.hpp index ac5725783..5f8a48836 100644 --- a/include/boost/math/special_functions/ellint_rf.hpp +++ b/include/boost/math/special_functions/ellint_rf.hpp @@ -27,82 +27,68 @@ namespace boost { namespace math { namespace detail{ -template -T ellint_rf_imp(T x, T y, T z, const Policy& pol) -{ - T value, X, Y, Z, E2, E3, u, lambda, tolerance; - unsigned long k; + template + T ellint_rf_imp(T x, T y, T z, const Policy& pol) + { + BOOST_MATH_STD_USING + using namespace boost::math; - BOOST_MATH_STD_USING - using namespace boost::math::tools; + static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; - static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; - - if (x < 0 || y < 0 || z < 0) - { - return policies::raise_domain_error(function, + 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); - } - if (x + y == 0 || y + z == 0 || z + x == 0) - { - return policies::raise_domain_error(function, + } + 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); - } + } - // Carlson scales error as the 6th power of tolerance, - // but this seems not to work for types larger than - // 80-bit reals, this heuristic seems to work OK: - if(policies::digits() > 64) - { - tolerance = pow(tools::epsilon(), T(1)/4.25f); - BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); - } - else - { - tolerance = pow(4*tools::epsilon(), T(1)/6); - BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); - } + T xn = x; + T yn = y; + T zn = z; + T An = (x + y + z) / 3; + T A0 = An; + T Q = pow(3 * boost::math::tools::epsilon(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn)); + T fn = 1; - // duplication - k = 1; - do - { - u = (x + y + z) / 3; - X = (u - x) / u; - Y = (u - y) / u; - Z = (u - z) / u; - // Termination condition: - if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) - break; + // duplication + unsigned k = 1; + for(; k < boost::math::policies::get_max_series_iterations(); ++k) + { + T root_x = sqrt(xn); + T root_y = sqrt(yn); + T root_z = sqrt(zn); + T lambda = root_x * root_y + root_x * root_z + root_y * root_z; + An = (An + lambda) / 4; + xn = (xn + lambda) / 4; + yn = (yn + lambda) / 4; + zn = (zn + lambda) / 4; + Q /= 4; + fn *= 4; + if(Q < fabs(An)) + break; + } + // Check to see if we gave up too soon: + policies::check_series_iterations(function, k, pol); + BOOST_MATH_INSTRUMENT_VARIABLE(k); - T sx = sqrt(x); - T sy = sqrt(y); - T sz = sqrt(z); - lambda = sy * (sx + sz) + sz * sx; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - z = (z + lambda) / 4; - ++k; - } - while(k < policies::get_max_series_iterations()); + T X = (A0 - x) / (An * fn); + T Y = (A0 - y) / (An * fn); + T Z = -X - Y; - // Check to see if we gave up too soon: - policies::check_series_iterations(function, k, pol); - BOOST_MATH_INSTRUMENT_VARIABLE(k); - - // Taylor series expansion to the 5th order - E2 = X * Y - Z * Z; - E3 = X * Y * Z; - value = (1 + E2*(E2/24 - E3*T(3)/44 - T(0.1)) + E3/14) / sqrt(u); - BOOST_MATH_INSTRUMENT_VARIABLE(value); - - return value; -} + // Taylor series expansion to the 7th order + T E2 = X * Y - Z * Z; + T E3 = X * Y * Z; + return (1 + E3 * (T(1) / 14 + 3 * E3 / 104) + E2 * (T(-1) / 10 + E2 / 24 - (3 * E3) / 44 - 5 * E2 * E2 / 208 + E2 * E3 / 16)) / sqrt(An); + } } // namespace detail diff --git a/test/test_carlson.cpp b/test/test_carlson.cpp index 1cf840883..de69de452 100644 --- a/test/test_carlson.cpp +++ b/test/test_carlson.cpp @@ -81,14 +81,14 @@ void expected_results() ".*", // platform largest_type, // test type(s) ".*RJ.*", // test data group - ".*", 180, 50); // test function + ".*", 250, 50); // test function add_expected_result( ".*", // compiler ".*", // stdlib ".*", // platform "real_concept", // test type(s) ".*RJ.*", // test data group - ".*", 180, 50); // test function + ".*", 250, 50); // test function add_expected_result( ".*", // compiler ".*", // stdlib diff --git a/test/test_carlson.hpp b/test/test_carlson.hpp index 7f40a6098..57e212c6a 100644 --- a/test/test_carlson.hpp +++ b/test/test_carlson.hpp @@ -131,14 +131,14 @@ void test_spots(T, const char* type_name) // Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227 // RF: T tolerance = (std::max)(T(1e-13f), tools::epsilon() * 5) * 100; // Note 5eps expressed as a persentage!!! - T eps2 = 2 * tools::epsilon(); + T eps2 = 5 * tools::epsilon(); 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); // 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)), log(T(2)), eps2); - BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(1)/4, T(-2)), log(T(2))/3, 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); // 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);