summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/ellint_rc.hpp
diff options
context:
space:
mode:
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.hpp
index 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 <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