diff --git a/include/boost/math/special_functions/ellint_1.hpp b/include/boost/math/special_functions/ellint_1.hpp
index 33d84c281..b4f550dd3 100644
--- a/include/boost/math/special_functions/ellint_1.hpp
+++ b/include/boost/math/special_functions/ellint_1.hpp
@@ -41,10 +41,12 @@ template <typename T, typename Policy>
 T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&);
 template <typename T, typename Policy>
 T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&);
+template <typename T, typename Policy>
+T ellint_k_imp(T k, const Policy& pol, T one_minus_k2);
 
 // Elliptic integral (Legendre form) of the first kind
 template <typename T, typename Policy>
-T ellint_f_imp(T phi, T k, const Policy& pol)
+T ellint_f_imp(T phi, T k, const Policy& pol, T one_minus_k2)
 {
     BOOST_MATH_STD_USING
     using namespace boost::math::tools;
@@ -80,7 +82,7 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
           std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
        > precision_tag_type;
 
-       result = 2 * phi * ellint_k_imp(k, pol, precision_tag_type()) / constants::pi<T>();
+       result = 2 * phi * ellint_k_imp(k, pol, one_minus_k2) / constants::pi<T>();
        BOOST_MATH_INSTRUMENT_VARIABLE(result);
     }
     else
@@ -121,31 +123,40 @@ T ellint_f_imp(T phi, T k, const Policy& pol)
           BOOST_MATH_ASSERT(rphi != 0); // precondition, can't be true if sin(rphi) != 0.
           //
           // Use http://dlmf.nist.gov/19.25#E5, note that
-          // c-1 simplifies to cot^2(rphi) which avoid cancellation:
+          // c-1 simplifies to cot^2(rphi) which avoids cancellation.
+          // Likewise c - k^2 is the same as (c - 1) + (1 - k^2).
           //
           T c = 1 / sinp;
-          result = static_cast<T>(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol));
+          T c_minus_one = cosp / sinp;
+          T cross = fabs(c / (k * k));
+          T arg2;
+          if ((cross > 0.9) && (cross < 1.1))
+             arg2 = c_minus_one + one_minus_k2;
+          else
+             arg2 = c - k * k;
+          result = static_cast<T>(s * ellint_rf_imp(c_minus_one, arg2, c, pol));
        }
        else
           result = s * sin(rphi);
        BOOST_MATH_INSTRUMENT_VARIABLE(result);
        if(m != 0)
        {
-          typedef std::integral_constant<int,
-             std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
-             std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
-          > precision_tag_type;
-
-          result += m * ellint_k_imp(k, pol, precision_tag_type());
+          result += m * ellint_k_imp(k, pol, one_minus_k2);
           BOOST_MATH_INSTRUMENT_VARIABLE(result);
        }
     }
     return invert ? T(-result) : result;
 }
 
+template <typename T, typename Policy>
+inline T ellint_f_imp(T phi, T k, const Policy& pol)
+{
+   return ellint_f_imp(phi, k, pol, T(1 - k * k));
+}
+
 // Complete elliptic integral (Legendre form) of the first kind
 template <typename T, typename Policy>
-T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
+T ellint_k_imp(T k, const Policy& pol, T one_minus_k2)
 {
     BOOST_MATH_STD_USING
     using namespace boost::math::tools;
@@ -162,12 +173,16 @@ T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
     }
 
     T x = 0;
-    T y = 1 - k * k;
     T z = 1;
-    T value = ellint_rf_imp(x, y, z, pol);
+    T value = ellint_rf_imp(x, one_minus_k2, z, pol);
 
     return value;
 }
+template <typename T, typename Policy>
+inline T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
+{
+   return ellint_k_imp(k, pol, T(1 - k * k));
+}
 
 //
 // Special versions for double and 80-bit long double precision,
diff --git a/include/boost/math/special_functions/heuman_lambda.hpp b/include/boost/math/special_functions/heuman_lambda.hpp
index c2d16d46e..0fbf4a980 100644
--- a/include/boost/math/special_functions/heuman_lambda.hpp
+++ b/include/boost/math/special_functions/heuman_lambda.hpp
@@ -63,8 +63,8 @@ T heuman_lambda_imp(T phi, T k, const Policy& pol)
           return policies::raise_domain_error<T>(function, "When 1-k^2 == 1 then phi must be < Pi/2, but got phi = %1%", phi, pol);
        }
        else
-          ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol, precision_tag_type());
-       result = ratio + ellint_k_imp(k, pol, precision_tag_type()) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi<T>();
+          ratio = ellint_f_imp(phi, rkp, pol, k2) / ellint_k_imp(rkp, pol, k2);
+       result = ratio + ellint_k_imp(k, pol, precision_tag_type()) * jacobi_zeta_imp(phi, rkp, pol, k2) / constants::half_pi<T>();
     }
     return result;
 }
diff --git a/include/boost/math/special_functions/jacobi_zeta.hpp b/include/boost/math/special_functions/jacobi_zeta.hpp
index aca9ff6d4..4586992af 100644
--- a/include/boost/math/special_functions/jacobi_zeta.hpp
+++ b/include/boost/math/special_functions/jacobi_zeta.hpp
@@ -27,7 +27,7 @@ namespace detail{
 
 // Elliptic integral - Jacobi Zeta
 template <typename T, typename Policy>
-T jacobi_zeta_imp(T phi, T k, const Policy& pol)
+T jacobi_zeta_imp(T phi, T k, const Policy& pol, T kp)
 {
     BOOST_MATH_STD_USING
     using namespace boost::math::tools;
@@ -44,8 +44,9 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol)
     T sinp = sin(phi);
     T cosp = cos(phi);
     T s2 = sinp * sinp;
+    T c2 = cosp * cosp;
+    T one_minus_ks2 = kp + c2 - kp * c2;
     T k2 = k * k;
-    T kp = 1 - k2;
     if(k == 1)
        result = sinp * (boost::math::sign)(cosp);  // We get here by simplifying JacobiZeta[w, 1] in Mathematica, and the fact that 0 <= phi.
     else
@@ -54,11 +55,15 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol)
           std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
           std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
        > precision_tag_type;
-       result = k2 * sinp * cosp * sqrt(1 - k2 * s2) * ellint_rj_imp(T(0), kp, T(1), T(1 - k2 * s2), pol) / (3 * ellint_k_imp(k, pol, precision_tag_type()));
+       result = k2 * sinp * cosp * sqrt(one_minus_ks2) * ellint_rj_imp(T(0), kp, T(1), one_minus_ks2, pol) / (3 * ellint_k_imp(k, pol, kp));
     }
     return invert ? T(-result) : result;
 }
-
+template <typename T, typename Policy>
+inline T jacobi_zeta_imp(T phi, T k, const Policy& pol)
+{
+   return jacobi_zeta_imp(phi, k, pol, T(1 - k * k));
+}
 } // detail
 
 template <class T1, class T2, class Policy>