summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/ellint_rd.hpp
diff options
context:
space:
mode:
authorDongHun Kwak <dh0128.kwak@samsung.com>2016-03-21 06:45:20 (GMT)
committerDongHun Kwak <dh0128.kwak@samsung.com>2016-03-21 06:46:37 (GMT)
commit733b5d5ae2c5d625211e2985ac25728ac3f54883 (patch)
treea5b214744b256f07e1dc2bd7273035a7808c659f /boost/math/special_functions/ellint_rd.hpp
parent08c1e93fa36a49f49325a07fe91ff92c964c2b6c (diff)
downloadboost-733b5d5ae2c5d625211e2985ac25728ac3f54883.zip
boost-733b5d5ae2c5d625211e2985ac25728ac3f54883.tar.gz
boost-733b5d5ae2c5d625211e2985ac25728ac3f54883.tar.bz2
Imported Upstream version 1.58.0upstream/1.58.0refs/changes/33/62933/1
Change-Id: If0072143aa26874812e0db6872e1efb10a3e5e94 Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
Diffstat (limited to 'boost/math/special_functions/ellint_rd.hpp')
-rw-r--r--boost/math/special_functions/ellint_rd.hpp217
1 files changed, 144 insertions, 73 deletions
diff --git a/boost/math/special_functions/ellint_rd.hpp b/boost/math/special_functions/ellint_rd.hpp
index 61014d3..03b73b1 100644
--- a/boost/math/special_functions/ellint_rd.hpp
+++ b/boost/math/special_functions/ellint_rd.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)
@@ -7,6 +7,7 @@
// XZ wrote the original of this file as part of the Google
// Summer of Code 2006. JM modified it slightly to fit into the
// Boost.Math conceptual framework better.
+// Updated 2015 to use Carlson's latest methods.
#ifndef BOOST_MATH_ELLINT_RD_HPP
#define BOOST_MATH_ELLINT_RD_HPP
@@ -16,6 +17,8 @@
#endif
#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/special_functions/ellint_rc.hpp>
+#include <boost/math/special_functions/pow.hpp>
#include <boost/math/tools/config.hpp>
#include <boost/math/policies/error_handling.hpp>
@@ -28,78 +31,146 @@ namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
T ellint_rd_imp(T x, T y, T z, const Policy& pol)
{
- T value, u, lambda, sigma, factor, tolerance;
- T X, Y, Z, EA, EB, EC, ED, EE, S1, S2;
- unsigned long k;
-
- BOOST_MATH_STD_USING
- using namespace boost::math::tools;
-
- static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)";
-
- if (x < 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument x must be >= 0, but got %1%", x, pol);
- }
- if (y < 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument y must be >= 0, but got %1%", y, pol);
- }
- if (z <= 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument z must be > 0, but got %1%", z, pol);
- }
- if (x + y == 0)
- {
- return policies::raise_domain_error<T>(function,
- "At most one argument can be zero, but got, x + y = %1%", x+y, pol);
- }
-
- // error scales as the 6th power of tolerance
- tolerance = pow(tools::epsilon<T>() / 3, T(1)/6);
-
- // duplication
- sigma = 0;
- factor = 1;
- k = 1;
- do
- {
- u = (x + y + z + z + z) / 5;
- X = (u - x) / u;
- Y = (u - y) / u;
- Z = (u - z) / u;
- if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance)
- break;
- T sx = sqrt(x);
- T sy = sqrt(y);
- T sz = sqrt(z);
- lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x);
- sigma += factor / (sz * (z + lambda));
- factor /= 4;
- x = (x + lambda) / 4;
- y = (y + lambda) / 4;
- z = (z + 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
- EA = X * Y;
- EB = Z * Z;
- EC = EA - EB;
- ED = EA - 6 * EB;
- EE = ED + EC + EC;
- S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14);
- S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26));
- value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u));
-
- return value;
+ BOOST_MATH_STD_USING
+ using std::swap;
+
+ static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)";
+
+ if(x < 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument x must be >= 0, but got %1%", x, pol);
+ }
+ if(y < 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument y must be >= 0, but got %1%", y, pol);
+ }
+ if(z <= 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument z must be > 0, but got %1%", z, pol);
+ }
+ if(x + y == 0)
+ {
+ return policies::raise_domain_error<T>(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
+ //
+ using std::swap;
+ if(x == z)
+ swap(x, y);
+ if(y == z)
+ {
+ if(x == y)
+ {
+ return 1 / (x * sqrt(x));
+ }
+ else if(x == 0)
+ {
+ return 3 * constants::pi<T>() / (4 * y * sqrt(y));
+ }
+ else
+ {
+ if((std::min)(x, y) / (std::max)(x, y) > 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) > 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)
+ }
+ if(y == 0)
+ swap(x, y);
+ if(x == 0)
+ {
+ //
+ // Special handling for common case, from
+ // Numerical Computation of Real or Complex Elliptic Integrals, eq.47
+ //
+ T xn = sqrt(y);
+ T yn = sqrt(z);
+ T x0 = xn;
+ T y0 = yn;
+ T sum = 0;
+ T sum_pow = 0.25f;
+
+ while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn))
+ {
+ T t = sqrt(xn * yn);
+ xn = (xn + yn) / 2;
+ yn = t;
+ sum_pow *= 2;
+ sum += sum_pow * boost::math::pow<2>(xn - yn);
+ }
+ T RF = constants::pi<T>() / (xn + yn);
+ //
+ // This following calculation suffers from serious cancellation when y ~ z
+ // unless we combine terms. We have:
+ //
+ // ( ((x0 + y0)/2)^2 - z ) / (z(y-z))
+ //
+ // Substituting y = x0^2 and z = y0^2 and simplifying we get the following:
+ //
+ T pt = (x0 + 3 * y0) / (4 * z * (x0 + y0));
+ //
+ // Since we've moved the demoninator from eq.47 inside the expression, we
+ // need to also scale "sum" by the same value:
+ //
+ pt -= sum / (z * (y - z));
+ return pt * RF * 3;
+ }
+
+ T xn = x;
+ T yn = y;
+ T zn = z;
+ T An = (x + y + 3 * z) / 5;
+ T A0 = An;
+ // This has an extra 1.2 fudge factor which is really only needed when x, y and z are close in magnitude:
+ T Q = pow(tools::epsilon<T>() / 4, -T(1) / 8) * (std::max)((std::max)(An - x, An - y), An - z) * 1.2f;
+ T lambda, rx, ry, rz;
+ unsigned k = 0;
+ T fn = 1;
+ T RD_sum = 0;
+
+ for(; k < policies::get_max_series_iterations<Policy>(); ++k)
+ {
+ rx = sqrt(xn);
+ ry = sqrt(yn);
+ rz = sqrt(zn);
+ lambda = rx * ry + rx * rz + ry * rz;
+ RD_sum += fn / (rz * (zn + lambda));
+ An = (An + lambda) / 4;
+ xn = (xn + lambda) / 4;
+ yn = (yn + lambda) / 4;
+ zn = (zn + lambda) / 4;
+ fn /= 4;
+ Q /= 4;
+ if(Q < An)
+ break;
+ }
+
+ policies::check_series_iterations<T, Policy>(function, k, pol);
+
+ T X = fn * (A0 - x) / An;
+ T Y = fn * (A0 - y) / An;
+ T Z = -(X + Y) / 3;
+ T E2 = X * Y - 6 * Z * Z;
+ T E3 = (3 * X * Y - 8 * Z * Z) * Z;
+ T E4 = 3 * (X * Y - Z * Z) * Z * Z;
+ T E5 = X * Y * Z * Z * Z;
+
+ T result = fn * pow(An, T(-3) / 2) *
+ (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16
+ + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68);
+ result += 3 * RD_sum;
+
+ return result;
}
} // namespace detail