summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/ellint_2.hpp
diff options
context:
space:
mode:
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.hpp
index 72caf3e..f4f65cc 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;
}