summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'boost/math/special_functions/ellint_2.hpp')
-rw-r--r--boost/math/special_functions/ellint_2.hpp40
1 files changed, 34 insertions, 6 deletions
 diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hppindex 72caf3e..f4f65cc 100644--- a/boost/math/special_functions/ellint_2.hpp+++ b/boost/math/special_functions/ellint_2.hpp@@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -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(); }+ 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("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 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; }