From 5eb3bbbfcccc8b0bd53f54e38887a28a77e2bf6c Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Mon, 16 Feb 2015 17:49:56 +0000
Subject: [PATCH] Apply Carlson's latest relations for the Legendre Elliptic
 integrals.

---
 .../boost/math/special_functions/ellint_1.hpp |  9 +++++-
 .../boost/math/special_functions/ellint_2.hpp | 32 +++++++++++++++----
 2 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/include/boost/math/special_functions/ellint_1.hpp b/include/boost/math/special_functions/ellint_1.hpp
index da16bc6f2..3590da6d1 100644
--- a/include/boost/math/special_functions/ellint_1.hpp
+++ b/include/boost/math/special_functions/ellint_1.hpp
@@ -103,10 +103,17 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
           BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
        }
        T sinp = sin(rphi);
+       sinp *= sinp;
        T cosp = cos(rphi);
+       cosp *= cosp;
+       T c = 1 / sinp;
        BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
        BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
-       result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol);
+       //
+       // Use http://dlmf.nist.gov/19.25#E5, note that
+       // c-1 simplifies to cot^2(rphi) which avoid cancellation:
+       //
+       result = rphi == 0 ? 0 : s * ellint_rf_imp(T(cosp / sinp), c - k * k, c, pol);
        BOOST_MATH_INSTRUMENT_VARIABLE(result);
        if(m != 0)
        {
diff --git a/include/boost/math/special_functions/ellint_2.hpp b/include/boost/math/special_functions/ellint_2.hpp
index 72caf3eb1..80ca0c931 100644
--- a/include/boost/math/special_functions/ellint_2.hpp
+++ b/include/boost/math/special_functions/ellint_2.hpp
@@ -21,6 +21,7 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/ellint_rf.hpp>
 #include <boost/math/special_functions/ellint_rd.hpp>
+#include <boost/math/special_functions/ellint_rg.hpp>
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/tools/workaround.hpp>
@@ -87,11 +88,30 @@ T ellint_e_imp(T phi, T k, const Policy& pol)
        }
        T sinp = sin(rphi);
        T cosp = cos(rphi);
-       T x = cosp * cosp;
-       T t = k * k * sinp * sinp;
-       T y = 1 - t;
-       T z = 1;
-       result = s * sinp * (ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3);
+       T c = 1 / (sinp * sinp);
+       T cm1 = cosp * cosp / (sinp * sinp);  // c - 1
+       T k2 = k * k;
+       if(k2 > 1)
+       {
+          return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
+       }
+       else if(rphi == 0)
+       {
+          result = 0;
+       }
+       else if(k == 0)
+       {
+          result = rphi;
+       }
+       else if(k2 == 1)
+       {
+          result = sinp;
+       }
+       else
+       {
+          // http://dlmf.nist.gov/19.25#E10
+          result = s * ((1 - k2) * ellint_rf_imp(cm1, c - k2, c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, c - k2, pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2))));
+       }
        if(m != 0)
           result += m * ellint_e_imp(k, pol);
     }
@@ -119,7 +139,7 @@ T ellint_e_imp(T k, const Policy& pol)
     T t = k * k;
     T y = 1 - t;
     T z = 1;
-    T value = ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3;
+    T value = 2 * ellint_rg_imp(x, y, z, pol);
 
     return value;
 }