diff options
Diffstat (limited to 'boost/math/special_functions/ellint_2.hpp')
-rw-r--r-- | boost/math/special_functions/ellint_2.hpp | 40 |
1 files changed, 34 insertions, 6 deletions
diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hpp index 72caf3eb11..f4f65cc11d 100644 --- a/boost/math/special_functions/ellint_2.hpp +++ b/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> @@ -67,6 +68,14 @@ T ellint_e_imp(T phi, T k, const Policy& pol) // just return the second part of the duplication formula: result = 2 * phi * ellint_e_imp(k, pol) / constants::pi<T>(); } + else if(k == 0) + { + return invert ? T(-phi) : phi; + } + else if(fabs(k) == 1) + { + return invert ? T(-sin(phi)) : sin(phi); + } else { // Carlson's algorithm works only for |phi| <= pi/2, @@ -87,11 +96,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(sinp * sinp < tools::min_value<T>()) + { + 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); + } + else + { + // http://dlmf.nist.gov/19.25#E10 + result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); + } if(m != 0) result += m * ellint_e_imp(k, pol); } @@ -119,7 +147,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; } |