diff options
Diffstat (limited to 'boost/math/special_functions/ellint_rc.hpp')
-rw-r--r-- | boost/math/special_functions/ellint_rc.hpp | 63 |
1 files changed, 31 insertions, 32 deletions
diff --git a/boost/math/special_functions/ellint_rc.hpp b/boost/math/special_functions/ellint_rc.hpp index 5f6d5c64bc..846c752a14 100644 --- a/boost/math/special_functions/ellint_rc.hpp +++ b/boost/math/special_functions/ellint_rc.hpp @@ -1,4 +1,4 @@ -// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -8,6 +8,7 @@ // Summer of Code 2006. JM modified it to fit into the // Boost.Math conceptual framework better, and to correctly // handle the y < 0 case. +// Updated 2015 to use Carlson's latest methods. // #ifndef BOOST_MATH_ELLINT_RC_HPP @@ -20,6 +21,9 @@ #include <boost/math/policies/error_handling.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/log1p.hpp> +#include <boost/math/constants/constants.hpp> +#include <iostream> // Carlson's degenerate elliptic integral // R_C(x, y) = R_F(x, y, y) = 0.5 * \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1} dt @@ -30,11 +34,7 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> T ellint_rc_imp(T x, T y, const Policy& pol) { - T value, S, u, lambda, tolerance, prefix; - unsigned long k; - BOOST_MATH_STD_USING - using namespace boost::math::tools; static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)"; @@ -49,11 +49,9 @@ T ellint_rc_imp(T x, T y, const Policy& pol) "Argument y must not be zero but got %1%", y, pol); } - // error scales as the 6th power of tolerance - tolerance = pow(4 * tools::epsilon<T>(), T(1) / 6); - // for y < 0, the integral is singular, return Cauchy principal value - if (y < 0) + T prefix, result; + if(y < 0) { prefix = sqrt(x / (x - y)); x = x - y; @@ -62,30 +60,31 @@ T ellint_rc_imp(T x, T y, const Policy& pol) else prefix = 1; - // duplication: - k = 1; - do + if(x == 0) + { + result = constants::half_pi<T>() / sqrt(y); + } + else if(x == y) + { + result = 1 / sqrt(x); + } + else if(y > x) { - u = (x + y + y) / 3; - S = y / u - 1; // 1 - x / u = 2 * S - - if (2 * abs(S) < tolerance) - break; - - T sx = sqrt(x); - T sy = sqrt(y); - lambda = 2 * sx * sy + y; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - ++k; - }while(k < policies::get_max_series_iterations<Policy>()); - // Check to see if we gave up too soon: - policies::check_series_iterations<T>(function, k, pol); - - // Taylor series expansion to the 5th order - value = (1 + S * S * (T(3) / 10 + S * (T(1) / 7 + S * (T(3) / 8 + S * T(9) / 22)))) / sqrt(u); - - return value * prefix; + result = atan(sqrt((y - x) / x)) / sqrt(y - x); + } + else + { + if(y / x > 0.5) + { + T arg = sqrt((x - y) / x); + result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(x - y)); + } + else + { + result = log((sqrt(x) + sqrt(x - y)) / sqrt(y)) / sqrt(x - y); + } + } + return prefix * result; } } // namespace detail |