diff --git a/include/boost/math/special_functions/ellint_3.hpp b/include/boost/math/special_functions/ellint_3.hpp
index 26b429c87..d90659110 100644
--- a/include/boost/math/special_functions/ellint_3.hpp
+++ b/include/boost/math/special_functions/ellint_3.hpp
@@ -259,39 +259,14 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol)
 
     if(v < 0)
     {
+       // Apply A&S 17.7.17:
        T k2 = k * k;
        T N = (k2 - v) / (1 - v);
        T Nm1 = (1 - k2) / (1 - v);
-       T p2 = sqrt(-v * (k2 - v) / (1 - v));
        T result = 0;
-       if(k2 < N)
-       {
-          result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
-          result *= sqrt(Nm1 * (1 - k2 / N));
-       }
-       /* else result = 0;
-          Result so far must be zero if k^2/N == 1, note that
-          since k^2 <= 1 it is necessarily true that k^2/N <= 1
-          even though slight roundoff error may actually yield
-          values > 1
-          */
-
-       if(k != 0)
-       {
-          if((k2 <= tools::min_value<T>()) || (p2 <= tools::min_value<T>()))
-          {
-             // We need to use logs to calculate k2 / p2:
-             T lk2 = 2 * log(fabs(k)) - (log(-v) + log(k2 - v) - log(1 - v)) / 2;
-             lk2 = exp(lk2);
-             if(lk2 != 0)
-               result += ellint_k_imp(k, pol) * lk2;
-          }
-          else
-          {
-             result += ellint_k_imp(k, pol) * k2 / p2;
-          }
-       }
-       result /= sqrt((1 - v) * (1 - k2 / v));
+       result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
+       result *= -v * (1 - k2) / ((1 - v) * (k2 - v));
+       result += ellint_k_imp(k, pol) * k2 / (k2 - v);
        return result;
     }
 
diff --git a/test/test_ellint_3.hpp b/test/test_ellint_3.hpp
index 813ef0876..79f5a8376 100644
--- a/test/test_ellint_3.hpp
+++ b/test/test_ellint_3.hpp
@@ -125,7 +125,7 @@ void test_spots(T, const char* type_name)
     do_test_ellint_pi3<T>(ellint_pi3_large_data, type_name, "Elliptic Integral PI: Large Random Data");
 
     // function values calculated on http://functions.wolfram.com/
-    static const boost::array<boost::array<T, 3>, 14> data2 = {{
+    static const boost::array<boost::array<T, 3>, 16> data2 = {{
         {{ SC_(0.0), SC_(0.2), SC_(1.586867847454166237308008033828114192951) }},
         {{ SC_(0.0), SC_(0.4), SC_(1.639999865864511206865258329748601457626) }},
         {{ SC_(0.0), SC_(0.0), SC_(1.57079632679489661923132169163975144209858469968755291048747) }},
@@ -141,6 +141,8 @@ void test_spots(T, const char* type_name)
         { { SC_(1e-170), SC_(-1E-164), SC_(1.57079632679489661923132169163975144209858469968755291048747) } },
         { { -1.5f * ldexp(T(1), -52), SC_(-0.9375), SC_(2.48840049140103464299631535211815755485846563527849342319632) } },
         { { -1.5f * ldexp(T(1), -52), SC_(0.9375), SC_(2.48840049140103464299631535211815755485846563527849342319632) } },
+        { { ldexp(T(1), -560), ldexp(T(1), -165), SC_(1.57079632679489661923132169163975144209858469968756130722545) } },
+        { { ldexp(T(1), -560), -ldexp(T(1), -165), SC_(1.57079632679489661923132169163975144209858469968754451374949) } },
     } };
 
     do_test_ellint_pi2<T>(data2, type_name, "Complete Elliptic Integral PI: Mathworld Data");