1 // Copyright (c) 2006 Xiaogang Zhang
2 // Copyright (c) 2006 John Maddock
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0. (See accompanying file
5 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
8 // XZ wrote the original of this file as part of the Google
9 // Summer of Code 2006. JM modified it to fit into the
10 // Boost.Math conceptual framework better, and to ensure
11 // that the code continues to work no matter how many digits
14 #ifndef BOOST_MATH_ELLINT_D_HPP
15 #define BOOST_MATH_ELLINT_D_HPP
21 #include <boost/math/special_functions/math_fwd.hpp>
22 #include <boost/math/special_functions/ellint_rf.hpp>
23 #include <boost/math/special_functions/ellint_rd.hpp>
24 #include <boost/math/special_functions/ellint_rg.hpp>
25 #include <boost/math/constants/constants.hpp>
26 #include <boost/math/policies/error_handling.hpp>
27 #include <boost/math/tools/workaround.hpp>
28 #include <boost/math/special_functions/round.hpp>
30 // Elliptic integrals (complete and incomplete) of the second kind
31 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
33 namespace boost { namespace math {
35 template <class T1, class T2, class Policy>
36 typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol);
40 template <typename T, typename Policy>
41 T ellint_d_imp(T k, const Policy& pol);
43 // Elliptic integral (Legendre form) of the second kind
44 template <typename T, typename Policy>
45 T ellint_d_imp(T phi, T k, const Policy& pol)
48 using namespace boost::math::tools;
49 using namespace boost::math::constants;
60 if(phi >= tools::max_value<T>())
62 // Need to handle infinity as a special case:
63 result = policies::raise_overflow_error<T>("boost::math::ellint_d<%1%>(%1%,%1%)", 0, pol);
65 else if(phi > 1 / tools::epsilon<T>())
67 // Phi is so large that phi%pi is necessarily zero (or garbage),
68 // just return the second part of the duplication formula:
69 result = 2 * phi * ellint_d_imp(k, pol) / constants::pi<T>();
73 // Carlson's algorithm works only for |phi| <= pi/2,
74 // use the integrand's periodicity to normalize phi
76 T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
77 T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
79 if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
83 rphi = constants::half_pi<T>() - rphi;
87 T c = 1 / (sinp * sinp);
88 T cm1 = cosp * cosp / (sinp * sinp); // c - 1
92 return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
100 // http://dlmf.nist.gov/19.25#E10
101 result = s * ellint_rd_imp(cm1, T(c - k2), c, pol) / 3;
104 result += m * ellint_d_imp(k, pol);
106 return invert ? T(-result) : result;
109 // Complete elliptic integral (Legendre form) of the second kind
110 template <typename T, typename Policy>
111 T ellint_d_imp(T k, const Policy& pol)
114 using namespace boost::math::tools;
118 return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)",
119 "Got k = %1%, function requires |k| <= 1", k, pol);
121 if(fabs(k) <= tools::root_epsilon<T>())
122 return constants::pi<T>() / 4;
128 T value = ellint_rd_imp(x, y, z, pol) / 3;
133 template <typename T, typename Policy>
134 inline typename tools::promote_args<T>::type ellint_d(T k, const Policy& pol, const mpl::true_&)
136 typedef typename tools::promote_args<T>::type result_type;
137 typedef typename policies::evaluation<result_type, Policy>::type value_type;
138 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(k), pol), "boost::math::ellint_d<%1%>(%1%)");
141 // Elliptic integral (Legendre form) of the second kind
142 template <class T1, class T2>
143 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const mpl::false_&)
145 return boost::math::ellint_d(k, phi, policies::policy<>());
150 // Complete elliptic integral (Legendre form) of the second kind
151 template <typename T>
152 inline typename tools::promote_args<T>::type ellint_d(T k)
154 return ellint_d(k, policies::policy<>());
157 // Elliptic integral (Legendre form) of the second kind
158 template <class T1, class T2>
159 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi)
161 typedef typename policies::is_policy<T2>::type tag_type;
162 return detail::ellint_d(k, phi, tag_type());
165 template <class T1, class T2, class Policy>
166 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol)
168 typedef typename tools::promote_args<T1, T2>::type result_type;
169 typedef typename policies::evaluation<result_type, Policy>::type value_type;
170 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");
175 #endif // BOOST_MATH_ELLINT_D_HPP