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 8a242f06a4..ac39bed678 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