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_2_HPP
15 #define BOOST_MATH_ELLINT_2_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_2(T1 k, T2 phi, const Policy& pol);
37
38 namespace detail{
39
40 template <typename T, typename Policy>
41 T ellint_e_imp(T k, const Policy& pol);
42
43 // Elliptic integral (Legendre form) of the second kind
44 template <typename T, typename Policy>
ellint_e_imp(T phi,T k,const Policy & pol)45 T ellint_e_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 return 0;
54
55 if(phi < 0)
56 {
57 phi = fabs(phi);
58 invert = true;
59 }
60
61 T result;
62
63 if(phi >= tools::max_value<T>())
64 {
65 // Need to handle infinity as a special case:
66 result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", 0, pol);
67 }
68 else if(phi > 1 / tools::epsilon<T>())
69 {
70 // Phi is so large that phi%pi is necessarily zero (or garbage),
71 // just return the second part of the duplication formula:
72 result = 2 * phi * ellint_e_imp(k, pol) / constants::pi<T>();
73 }
74 else if(k == 0)
75 {
76 return invert ? T(-phi) : phi;
77 }
78 else if(fabs(k) == 1)
79 {
80 return invert ? T(-sin(phi)) : T(sin(phi));
81 }
82 else
83 {
84 // Carlson's algorithm works only for |phi| <= pi/2,
85 // use the integrand's periodicity to normalize phi
86 //
87 // Xiaogang's original code used a cast to long long here
88 // but that fails if T has more digits than a long long,
89 // so rewritten to use fmod instead:
90 //
91 T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
92 T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
93 int s = 1;
94 if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
95 {
96 m += 1;
97 s = -1;
98 rphi = constants::half_pi<T>() - rphi;
99 }
100 T k2 = k * k;
101 if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon<T>() * fabs(rphi))
102 {
103 // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/
104 result = s * rphi;
105 }
106 else
107 {
108 // http://dlmf.nist.gov/19.25#E10
109 T sinp = sin(rphi);
110 if (k2 * sinp * sinp >= 1)
111 {
112 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);
113 }
114 T cosp = cos(rphi);
115 T c = 1 / (sinp * sinp);
116 T cm1 = cosp * cosp / (sinp * sinp); // c - 1
117 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))));
118 }
119 if(m != 0)
120 result += m * ellint_e_imp(k, pol);
121 }
122 return invert ? T(-result) : result;
123 }
124
125 // Complete elliptic integral (Legendre form) of the second kind
126 template <typename T, typename Policy>
ellint_e_imp(T k,const Policy & pol)127 T ellint_e_imp(T k, const Policy& pol)
128 {
129 BOOST_MATH_STD_USING
130 using namespace boost::math::tools;
131
132 if (abs(k) > 1)
133 {
134 return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
135 "Got k = %1%, function requires |k| <= 1", k, pol);
136 }
137 if (abs(k) == 1)
138 {
139 return static_cast<T>(1);
140 }
141
142 T x = 0;
143 T t = k * k;
144 T y = 1 - t;
145 T z = 1;
146 T value = 2 * ellint_rg_imp(x, y, z, pol);
147
148 return value;
149 }
150
151 template <typename T, typename Policy>
ellint_2(T k,const Policy & pol,const boost::true_type &)152 inline typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const boost::true_type&)
153 {
154 typedef typename tools::promote_args<T>::type result_type;
155 typedef typename policies::evaluation<result_type, Policy>::type value_type;
156 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%)");
157 }
158
159 // Elliptic integral (Legendre form) of the second kind
160 template <class T1, class T2>
ellint_2(T1 k,T2 phi,const boost::false_type &)161 inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const boost::false_type&)
162 {
163 return boost::math::ellint_2(k, phi, policies::policy<>());
164 }
165
166 } // detail
167
168 // Complete elliptic integral (Legendre form) of the second kind
169 template <typename T>
ellint_2(T k)170 inline typename tools::promote_args<T>::type ellint_2(T k)
171 {
172 return ellint_2(k, policies::policy<>());
173 }
174
175 // Elliptic integral (Legendre form) of the second kind
176 template <class T1, class T2>
ellint_2(T1 k,T2 phi)177 inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
178 {
179 typedef typename policies::is_policy<T2>::type tag_type;
180 return detail::ellint_2(k, phi, tag_type());
181 }
182
183 template <class T1, class T2, class Policy>
ellint_2(T1 k,T2 phi,const Policy & pol)184 inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)
185 {
186 typedef typename tools::promote_args<T1, T2>::type result_type;
187 typedef typename policies::evaluation<result_type, Policy>::type value_type;
188 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");
189 }
190
191 }} // namespaces
192
193 #endif // BOOST_MATH_ELLINT_2_HPP
194
195