summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/ellint_rj.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_rj.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_rj.hpp')
-rw-r--r--boost/math/special_functions/ellint_rj.hpp356
1 files changed, 239 insertions, 117 deletions
diff --git a/boost/math/special_functions/ellint_rj.hpp b/boost/math/special_functions/ellint_rj.hpp
index 8a242f0..ac39bed 100644
--- a/boost/math/special_functions/ellint_rj.hpp
+++ b/boost/math/special_functions/ellint_rj.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 p < 0 case.
+// Updated 2015 to use Carlson's latest methods.
//
#ifndef BOOST_MATH_ELLINT_RJ_HPP
@@ -22,6 +23,7 @@
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/ellint_rc.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
+#include <boost/math/special_functions/ellint_rd.hpp>
// Carlson's elliptic integral of the third kind
// R_J(x, y, z, p) = 1.5 * \int_{0}^{\infty} (t+p)^{-1} [(t+x)(t+y)(t+z)]^{-1/2} dt
@@ -30,124 +32,244 @@
namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
+T ellint_rc1p_imp(T y, const Policy& pol)
+{
+ using namespace boost::math;
+ // Calculate RC(1, 1 + x)
+ BOOST_MATH_STD_USING
+
+ static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)";
+
+ if(y == -1)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument y must not be zero but got %1%", y, pol);
+ }
+
+ // for 1 + y < 0, the integral is singular, return Cauchy principal value
+ T result;
+ if(y < -1)
+ {
+ result = sqrt(1 / -y) * detail::ellint_rc_imp(T(-y), T(-1 - y), pol);
+ }
+ else if(y == 0)
+ {
+ result = 1;
+ }
+ else if(y > 0)
+ {
+ result = atan(sqrt(y)) / sqrt(y);
+ }
+ else
+ {
+ if(y > -0.5)
+ {
+ T arg = sqrt(-y);
+ result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(-y));
+ }
+ else
+ {
+ result = log((1 + sqrt(-y)) / sqrt(1 + y)) / sqrt(-y);
+ }
+ }
+ return result;
+}
+
+template <typename T, typename Policy>
T ellint_rj_imp(T x, T y, T z, T p, const Policy& pol)
{
- T value, u, lambda, alpha, beta, sigma, factor, tolerance;
- T X, Y, Z, P, EA, EB, EC, E2, E3, S1, S2, S3;
- unsigned long k;
-
- BOOST_MATH_STD_USING
- using namespace boost::math::tools;
-
- static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)";
-
- if (x < 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument x must be non-negative, but got x = %1%", x, pol);
- }
- if(y < 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument y must be non-negative, but got y = %1%", y, pol);
- }
- if(z < 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument z must be non-negative, but got z = %1%", z, pol);
- }
- if(p == 0)
- {
- return policies::raise_domain_error<T>(function,
- "Argument p must not be zero, but got p = %1%", p, pol);
- }
- if (x + y == 0 || y + z == 0 || z + x == 0)
- {
- return policies::raise_domain_error<T>(function,
- "At most one argument can be zero, "
- "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
- }
-
- // error scales as the 6th power of tolerance
- tolerance = pow(T(1) * tools::epsilon<T>() / 3, T(1) / 6);
-
- // for p < 0, the integral is singular, return Cauchy principal value
- if (p < 0)
- {
- //
- // We must ensure that (z - y) * (y - x) is positive.
- // Since the integral is symmetrical in x, y and z
- // we can just permute the values:
- //
- if(x > y)
- std::swap(x, y);
- if(y > z)
- std::swap(y, z);
- if(x > y)
- std::swap(x, y);
-
- T q = -p;
- T pmy = (z - y) * (y - x) / (y + q); // p - y
-
- BOOST_ASSERT(pmy >= 0);
-
- p = pmy + y;
- value = boost::math::ellint_rj(x, y, z, p, pol);
- value *= pmy;
- value -= 3 * boost::math::ellint_rf(x, y, z, pol);
- value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(x * z + p * q, p * q, pol);
- value /= (y + q);
- return value;
- }
-
- // duplication
- sigma = 0;
- factor = 1;
- k = 1;
- do
- {
- u = (x + y + z + p + p) / 5;
- X = (u - x) / u;
- Y = (u - y) / u;
- Z = (u - z) / u;
- P = (u - p) / u;
-
- if ((tools::max)(abs(X), abs(Y), abs(Z), abs(P)) < tolerance)
- break;
-
- T sx = sqrt(x);
- T sy = sqrt(y);
- T sz = sqrt(z);
-
- lambda = sy * (sx + sz) + sz * sx;
- alpha = p * (sx + sy + sz) + sx * sy * sz;
- alpha *= alpha;
- beta = p * (p + lambda) * (p + lambda);
- sigma += factor * boost::math::ellint_rc(alpha, beta, pol);
- factor /= 4;
- x = (x + lambda) / 4;
- y = (y + lambda) / 4;
- z = (z + lambda) / 4;
- p = (p + 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 + Y * Z + Z * X;
- EB = X * Y * Z;
- EC = P * P;
- E2 = EA - 3 * EC;
- E3 = EB + 2 * P * (EA - EC);
- S1 = 1 + E2 * (E2 * T(9) / 88 - E3 * T(9) / 52 - T(3) / 14);
- S2 = EB * (T(1) / 6 + P * (T(-6) / 22 + P * T(3) / 26));
- S3 = P * ((EA - EC) / 3 - P * EA * T(3) / 22);
- value = 3 * sigma + factor * (S1 + S2 + S3) / (u * sqrt(u));
-
- return value;
+ BOOST_MATH_STD_USING
+
+ static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)";
+
+ if(x < 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument x must be non-negative, but got x = %1%", x, pol);
+ }
+ if(y < 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument y must be non-negative, but got y = %1%", y, pol);
+ }
+ if(z < 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument z must be non-negative, but got z = %1%", z, pol);
+ }
+ if(p == 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "Argument p must not be zero, but got p = %1%", p, pol);
+ }
+ if(x + y == 0 || y + z == 0 || z + x == 0)
+ {
+ return policies::raise_domain_error<T>(function,
+ "At most one argument can be zero, "
+ "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
+ }
+
+ // for p < 0, the integral is singular, return Cauchy principal value
+ if(p < 0)
+ {
+ //
+ // We must ensure that x < y < z.
+ // Since the integral is symmetrical in x, y and z
+ // we can just permute the values:
+ //
+ if(x > y)
+ std::swap(x, y);
+ if(y > z)
+ std::swap(y, z);
+ if(x > y)
+ std::swap(x, y);
+
+ BOOST_ASSERT(x <= y);
+ BOOST_ASSERT(y <= z);
+
+ T q = -p;
+ p = (z * (x + y + q) - x * y) / (z + q);
+
+ BOOST_ASSERT(p >= 0);
+
+ T value = (p - z) * ellint_rj_imp(x, y, z, p, pol);
+ value -= 3 * ellint_rf_imp(x, y, z, pol);
+ value += 3 * sqrt((x * y * z) / (x * y + p * q)) * ellint_rc_imp(T(x * y + p * q), T(p * q), pol);
+ value /= (z + q);
+ return value;
+ }
+
+ //
+ // Special cases from http://dlmf.nist.gov/19.20#iii
+ //
+ if(x == y)
+ {
+ if(x == z)
+ {
+ if(x == p)
+ {
+ // All values equal:
+ return 1 / (x * sqrt(x));
+ }
+ else
+ {
+ // x = y = z:
+ return 3 * (ellint_rc_imp(x, p, pol) - 1 / sqrt(x)) / (x - p);
+ }
+ }
+ else
+ {
+ // x = y only, permute so y = z:
+ using std::swap;
+ swap(x, z);
+ if(y == p)
+ {
+ return ellint_rd_imp(x, y, y, pol);
+ }
+ else if((std::max)(y, p) / (std::min)(y, p) > 1.2)
+ {
+ return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y);
+ }
+ // Otherwise fall through to normal method, special case above will suffer too much cancellation...
+ }
+ }
+ if(y == z)
+ {
+ if(y == p)
+ {
+ // y = z = p:
+ return ellint_rd_imp(x, y, y, pol);
+ }
+ else if((std::max)(y, p) / (std::min)(y, p) > 1.2)
+ {
+ // y = z:
+ return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y);
+ }
+ // Otherwise fall through to normal method, special case above will suffer too much cancellation...
+ }
+ if(z == p)
+ {
+ return ellint_rd_imp(x, y, z, pol);
+ }
+
+ T xn = x;
+ T yn = y;
+ T zn = z;
+ T pn = p;
+ T An = (x + y + z + 2 * p) / 5;
+ T A0 = An;
+ T delta = (p - x) * (p - y) * (p - z);
+ T Q = pow(tools::epsilon<T>() / 5, -T(1) / 8) * (std::max)((std::max)(fabs(An - x), fabs(An - y)), (std::max)(fabs(An - z), fabs(An - p)));
+
+ unsigned n;
+ T lambda;
+ T Dn;
+ T En;
+ T rx, ry, rz, rp;
+ T fmn = 1; // 4^-n
+ T RC_sum = 0;
+
+ for(n = 0; n < policies::get_max_series_iterations<Policy>(); ++n)
+ {
+ rx = sqrt(xn);
+ ry = sqrt(yn);
+ rz = sqrt(zn);
+ rp = sqrt(pn);
+ Dn = (rp + rx) * (rp + ry) * (rp + rz);
+ En = delta / Dn;
+ En /= Dn;
+ if((En < -0.5) && (En > -1.5))
+ {
+ //
+ // Occationally En ~ -1, we then have no means of calculating
+ // RC(1, 1+En) without terrible cancellation error, so we
+ // need to get to 1+En directly. By substitution we have
+ //
+ // 1+E_0 = 1 + (p-x)*(p-y)*(p-z)/((sqrt(p) + sqrt(x))*(sqrt(p)+sqrt(y))*(sqrt(p)+sqrt(z)))^2
+ // = 2*sqrt(p)*(p+sqrt(x) * (sqrt(y)+sqrt(z)) + sqrt(y)*sqrt(z)) / ((sqrt(p) + sqrt(x))*(sqrt(p) + sqrt(y)*(sqrt(p)+sqrt(z))))
+ //
+ // And since this is just an application of the duplication formula for RJ, the same
+ // expression works for 1+En if we use x,y,z,p_n etc.
+ // This branch is taken only once or twice at the start of iteration,
+ // after than En reverts to it's usual very small values.
+ //
+ T b = 2 * rp * (pn + rx * (ry + rz) + ry * rz) / Dn;
+ RC_sum += fmn / Dn * detail::ellint_rc_imp(T(1), b, pol);
+ }
+ else
+ {
+ RC_sum += fmn / Dn * ellint_rc1p_imp(En, pol);
+ }
+ lambda = rx * ry + rx * rz + ry * rz;
+
+ // From here on we move to n+1:
+ An = (An + lambda) / 4;
+ fmn /= 4;
+
+ if(fmn * Q < An)
+ break;
+
+ xn = (xn + lambda) / 4;
+ yn = (yn + lambda) / 4;
+ zn = (zn + lambda) / 4;
+ pn = (pn + lambda) / 4;
+ delta /= 64;
+ }
+
+ T X = fmn * (A0 - x) / An;
+ T Y = fmn * (A0 - y) / An;
+ T Z = fmn * (A0 - z) / An;
+ T P = (-X - Y - Z) / 2;
+ T E2 = X * Y + X * Z + Y * Z - 3 * P * P;
+ T E3 = X * Y * Z + 2 * E2 * P + 4 * P * P * P;
+ T E4 = (2 * X * Y * Z + E2 * P + 3 * P * P * P) * P;
+ T E5 = X * Y * Z * P * P;
+ T result = fmn * 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 += 6 * RC_sum;
+ return result;
}
} // namespace detail