summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
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.hppindex 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 #include #include +#include // 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 +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(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 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(function,- "Argument x must be non-negative, but got x = %1%", x, pol);- }- if(y < 0)- {- return policies::raise_domain_error(function,- "Argument y must be non-negative, but got y = %1%", y, pol);- }- if(z < 0)- {- return policies::raise_domain_error(function,- "Argument z must be non-negative, but got z = %1%", z, pol);- }- if(p == 0)- {- return policies::raise_domain_error(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(function,- "At most one argument can be zero, "- "only possible result is %1%.", std::numeric_limits::quiet_NaN(), pol);- }-- // error scales as the 6th power of tolerance- tolerance = pow(T(1) * tools::epsilon() / 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());-- // Check to see if we gave up too soon:- policies::check_series_iterations(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(function,+ "Argument x must be non-negative, but got x = %1%", x, pol);+ }+ if(y < 0)+ {+ return policies::raise_domain_error(function,+ "Argument y must be non-negative, but got y = %1%", y, pol);+ }+ if(z < 0)+ {+ return policies::raise_domain_error(function,+ "Argument z must be non-negative, but got z = %1%", z, pol);+ }+ if(p == 0)+ {+ return policies::raise_domain_error(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(function,+ "At most one argument can be zero, "+ "only possible result is %1%.", std::numeric_limits::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() / 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(); ++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