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)
6 //
7 // History:
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
12 // type T has.
13
14 #ifndef BOOST_MATH_ELLINT_D_HPP
15 #define BOOST_MATH_ELLINT_D_HPP
16
17 #ifdef _MSC_VER
18 #pragma once
19 #endif
20
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>
29
30 // Elliptic integrals (complete and incomplete) of the second kind
31 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
32
33 namespace boost { namespace math {
34
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);
37
38 namespace detail{
39
40 template <typename T, typename Policy>
41 T ellint_d_imp(T k, const Policy& pol);
42
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)
46 {
47 BOOST_MATH_STD_USING
48 using namespace boost::math::tools;
49 using namespace boost::math::constants;
50
51 bool invert = false;
52 if(phi < 0)
53 {
54 phi = fabs(phi);
55 invert = true;
56 }
57
58 T result;
59
60 if(phi >= tools::max_value<T>())
61 {
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);
64 }
65 else if(phi > 1 / tools::epsilon<T>())
66 {
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>();
70 }
71 else
72 {
73 // Carlson's algorithm works only for |phi| <= pi/2,
74 // use the integrand's periodicity to normalize phi
75 //
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>());
78 int s = 1;
79 if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
80 {
81 m += 1;
82 s = -1;
83 rphi = constants::half_pi<T>() - rphi;
84 }
85 T sinp = sin(rphi);
86 T cosp = cos(rphi);
87 T c = 1 / (sinp * sinp);
88 T cm1 = cosp * cosp / (sinp * sinp); // c - 1
89 T k2 = k * k;
90 if(k2 > 1)
91 {
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);
93 }
94 else if(rphi == 0)
95 {
96 result = 0;
97 }
98 else
99 {
100 // http://dlmf.nist.gov/19.25#E10
101 result = s * ellint_rd_imp(cm1, T(c - k2), c, pol) / 3;
102 }
103 if(m != 0)
104 result += m * ellint_d_imp(k, pol);
105 }
106 return invert ? T(-result) : result;
107 }
108
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)
112 {
113 BOOST_MATH_STD_USING
114 using namespace boost::math::tools;
115
116 if (abs(k) >= 1)
117 {
118 return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)",
119 "Got k = %1%, function requires |k| <= 1", k, pol);
120 }
121 if(fabs(k) <= tools::root_epsilon<T>())
122 return constants::pi<T>() / 4;
123
124 T x = 0;
125 T t = k * k;
126 T y = 1 - t;
127 T z = 1;
128 T value = ellint_rd_imp(x, y, z, pol) / 3;
129
130 return value;
131 }
132
133 template <typename T, typename Policy>
ellint_d(T k,const Policy & pol,const mpl::true_ &)134 inline typename tools::promote_args<T>::type ellint_d(T k, const Policy& pol, const mpl::true_&)
135 {
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%)");
139 }
140
141 // Elliptic integral (Legendre form) of the second kind
142 template <class T1, class T2>
ellint_d(T1 k,T2 phi,const mpl::false_ &)143 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const mpl::false_&)
144 {
145 return boost::math::ellint_d(k, phi, policies::policy<>());
146 }
147
148 } // detail
149
150 // Complete elliptic integral (Legendre form) of the second kind
151 template <typename T>
ellint_d(T k)152 inline typename tools::promote_args<T>::type ellint_d(T k)
153 {
154 return ellint_d(k, policies::policy<>());
155 }
156
157 // Elliptic integral (Legendre form) of the second kind
158 template <class T1, class T2>
ellint_d(T1 k,T2 phi)159 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi)
160 {
161 typedef typename policies::is_policy<T2>::type tag_type;
162 return detail::ellint_d(k, phi, tag_type());
163 }
164
165 template <class T1, class T2, class Policy>
ellint_d(T1 k,T2 phi,const Policy & pol)166 inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol)
167 {
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%)");
171 }
172
173 }} // namespaces
174
175 #endif // BOOST_MATH_ELLINT_D_HPP
176
177