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 61014d3866..03b73b159f 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