summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'boost/math/special_functions/ellint_rc.hpp')
-rw-r--r--boost/math/special_functions/ellint_rc.hpp63
1 files changed, 31 insertions, 32 deletions
 diff --git a/boost/math/special_functions/ellint_rc.hpp b/boost/math/special_functions/ellint_rc.hppindex 5f6d5c6..846c752 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 #include #include +#include +#include +#include // 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 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(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() / 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());- // Check to see if we gave up too soon:- policies::check_series_iterations(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