]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | // | |
6 | // History: | |
7 | // XZ wrote the original of this file as part of the Google | |
8 | // Summer of Code 2006. JM modified it to fit into the | |
9 | // Boost.Math conceptual framework better, and to correctly | |
10 | // handle the p < 0 case. | |
11 | // Updated 2015 to use Carlson's latest methods. | |
12 | // | |
13 | ||
14 | #ifndef BOOST_MATH_ELLINT_RJ_HPP | |
15 | #define BOOST_MATH_ELLINT_RJ_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/tools/config.hpp> | |
23 | #include <boost/math/policies/error_handling.hpp> | |
24 | #include <boost/math/special_functions/ellint_rc.hpp> | |
25 | #include <boost/math/special_functions/ellint_rf.hpp> | |
26 | #include <boost/math/special_functions/ellint_rd.hpp> | |
27 | ||
28 | // Carlson's elliptic integral of the third kind | |
29 | // R_J(x, y, z, p) = 1.5 * \int_{0}^{\infty} (t+p)^{-1} [(t+x)(t+y)(t+z)]^{-1/2} dt | |
30 | // Carlson, Numerische Mathematik, vol 33, 1 (1979) | |
31 | ||
32 | namespace boost { namespace math { namespace detail{ | |
33 | ||
34 | template <typename T, typename Policy> | |
35 | T ellint_rc1p_imp(T y, const Policy& pol) | |
36 | { | |
37 | using namespace boost::math; | |
38 | // Calculate RC(1, 1 + x) | |
39 | BOOST_MATH_STD_USING | |
40 | ||
41 | static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)"; | |
42 | ||
43 | if(y == -1) | |
44 | { | |
45 | return policies::raise_domain_error<T>(function, | |
46 | "Argument y must not be zero but got %1%", y, pol); | |
47 | } | |
48 | ||
49 | // for 1 + y < 0, the integral is singular, return Cauchy principal value | |
50 | T result; | |
51 | if(y < -1) | |
52 | { | |
53 | result = sqrt(1 / -y) * detail::ellint_rc_imp(T(-y), T(-1 - y), pol); | |
54 | } | |
55 | else if(y == 0) | |
56 | { | |
57 | result = 1; | |
58 | } | |
59 | else if(y > 0) | |
60 | { | |
61 | result = atan(sqrt(y)) / sqrt(y); | |
62 | } | |
63 | else | |
64 | { | |
65 | if(y > -0.5) | |
66 | { | |
67 | T arg = sqrt(-y); | |
68 | result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(-y)); | |
69 | } | |
70 | else | |
71 | { | |
72 | result = log((1 + sqrt(-y)) / sqrt(1 + y)) / sqrt(-y); | |
73 | } | |
74 | } | |
75 | return result; | |
76 | } | |
77 | ||
78 | template <typename T, typename Policy> | |
79 | T ellint_rj_imp(T x, T y, T z, T p, const Policy& pol) | |
80 | { | |
81 | BOOST_MATH_STD_USING | |
82 | ||
83 | static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)"; | |
84 | ||
85 | if(x < 0) | |
86 | { | |
87 | return policies::raise_domain_error<T>(function, | |
88 | "Argument x must be non-negative, but got x = %1%", x, pol); | |
89 | } | |
90 | if(y < 0) | |
91 | { | |
92 | return policies::raise_domain_error<T>(function, | |
93 | "Argument y must be non-negative, but got y = %1%", y, pol); | |
94 | } | |
95 | if(z < 0) | |
96 | { | |
97 | return policies::raise_domain_error<T>(function, | |
98 | "Argument z must be non-negative, but got z = %1%", z, pol); | |
99 | } | |
100 | if(p == 0) | |
101 | { | |
102 | return policies::raise_domain_error<T>(function, | |
103 | "Argument p must not be zero, but got p = %1%", p, pol); | |
104 | } | |
105 | if(x + y == 0 || y + z == 0 || z + x == 0) | |
106 | { | |
107 | return policies::raise_domain_error<T>(function, | |
108 | "At most one argument can be zero, " | |
109 | "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); | |
110 | } | |
111 | ||
112 | // for p < 0, the integral is singular, return Cauchy principal value | |
113 | if(p < 0) | |
114 | { | |
115 | // | |
116 | // We must ensure that x < y < z. | |
117 | // Since the integral is symmetrical in x, y and z | |
118 | // we can just permute the values: | |
119 | // | |
120 | if(x > y) | |
121 | std::swap(x, y); | |
122 | if(y > z) | |
123 | std::swap(y, z); | |
124 | if(x > y) | |
125 | std::swap(x, y); | |
126 | ||
127 | BOOST_ASSERT(x <= y); | |
128 | BOOST_ASSERT(y <= z); | |
129 | ||
130 | T q = -p; | |
131 | p = (z * (x + y + q) - x * y) / (z + q); | |
132 | ||
133 | BOOST_ASSERT(p >= 0); | |
134 | ||
135 | T value = (p - z) * ellint_rj_imp(x, y, z, p, pol); | |
136 | value -= 3 * ellint_rf_imp(x, y, z, pol); | |
137 | value += 3 * sqrt((x * y * z) / (x * y + p * q)) * ellint_rc_imp(T(x * y + p * q), T(p * q), pol); | |
138 | value /= (z + q); | |
139 | return value; | |
140 | } | |
141 | ||
142 | // | |
143 | // Special cases from http://dlmf.nist.gov/19.20#iii | |
144 | // | |
145 | if(x == y) | |
146 | { | |
147 | if(x == z) | |
148 | { | |
149 | if(x == p) | |
150 | { | |
151 | // All values equal: | |
152 | return 1 / (x * sqrt(x)); | |
153 | } | |
154 | else | |
155 | { | |
156 | // x = y = z: | |
157 | return 3 * (ellint_rc_imp(x, p, pol) - 1 / sqrt(x)) / (x - p); | |
158 | } | |
159 | } | |
160 | else | |
161 | { | |
162 | // x = y only, permute so y = z: | |
163 | using std::swap; | |
164 | swap(x, z); | |
165 | if(y == p) | |
166 | { | |
167 | return ellint_rd_imp(x, y, y, pol); | |
168 | } | |
169 | else if((std::max)(y, p) / (std::min)(y, p) > 1.2) | |
170 | { | |
171 | return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y); | |
172 | } | |
173 | // Otherwise fall through to normal method, special case above will suffer too much cancellation... | |
174 | } | |
175 | } | |
176 | if(y == z) | |
177 | { | |
178 | if(y == p) | |
179 | { | |
180 | // y = z = p: | |
181 | return ellint_rd_imp(x, y, y, pol); | |
182 | } | |
183 | else if((std::max)(y, p) / (std::min)(y, p) > 1.2) | |
184 | { | |
185 | // y = z: | |
186 | return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y); | |
187 | } | |
188 | // Otherwise fall through to normal method, special case above will suffer too much cancellation... | |
189 | } | |
190 | if(z == p) | |
191 | { | |
192 | return ellint_rd_imp(x, y, z, pol); | |
193 | } | |
194 | ||
195 | T xn = x; | |
196 | T yn = y; | |
197 | T zn = z; | |
198 | T pn = p; | |
199 | T An = (x + y + z + 2 * p) / 5; | |
200 | T A0 = An; | |
201 | T delta = (p - x) * (p - y) * (p - z); | |
202 | T Q = pow(tools::epsilon<T>() / 5, -T(1) / 8) * (std::max)((std::max)(fabs(An - x), fabs(An - y)), (std::max)(fabs(An - z), fabs(An - p))); | |
203 | ||
204 | unsigned n; | |
205 | T lambda; | |
206 | T Dn; | |
207 | T En; | |
208 | T rx, ry, rz, rp; | |
209 | T fmn = 1; // 4^-n | |
210 | T RC_sum = 0; | |
211 | ||
212 | for(n = 0; n < policies::get_max_series_iterations<Policy>(); ++n) | |
213 | { | |
214 | rx = sqrt(xn); | |
215 | ry = sqrt(yn); | |
216 | rz = sqrt(zn); | |
217 | rp = sqrt(pn); | |
218 | Dn = (rp + rx) * (rp + ry) * (rp + rz); | |
219 | En = delta / Dn; | |
220 | En /= Dn; | |
221 | if((En < -0.5) && (En > -1.5)) | |
222 | { | |
223 | // | |
224 | // Occationally En ~ -1, we then have no means of calculating | |
225 | // RC(1, 1+En) without terrible cancellation error, so we | |
226 | // need to get to 1+En directly. By substitution we have | |
227 | // | |
228 | // 1+E_0 = 1 + (p-x)*(p-y)*(p-z)/((sqrt(p) + sqrt(x))*(sqrt(p)+sqrt(y))*(sqrt(p)+sqrt(z)))^2 | |
229 | // = 2*sqrt(p)*(p+sqrt(x) * (sqrt(y)+sqrt(z)) + sqrt(y)*sqrt(z)) / ((sqrt(p) + sqrt(x))*(sqrt(p) + sqrt(y)*(sqrt(p)+sqrt(z)))) | |
230 | // | |
231 | // And since this is just an application of the duplication formula for RJ, the same | |
232 | // expression works for 1+En if we use x,y,z,p_n etc. | |
233 | // This branch is taken only once or twice at the start of iteration, | |
234 | // after than En reverts to it's usual very small values. | |
235 | // | |
236 | T b = 2 * rp * (pn + rx * (ry + rz) + ry * rz) / Dn; | |
237 | RC_sum += fmn / Dn * detail::ellint_rc_imp(T(1), b, pol); | |
238 | } | |
239 | else | |
240 | { | |
241 | RC_sum += fmn / Dn * ellint_rc1p_imp(En, pol); | |
242 | } | |
243 | lambda = rx * ry + rx * rz + ry * rz; | |
244 | ||
245 | // From here on we move to n+1: | |
246 | An = (An + lambda) / 4; | |
247 | fmn /= 4; | |
248 | ||
249 | if(fmn * Q < An) | |
250 | break; | |
251 | ||
252 | xn = (xn + lambda) / 4; | |
253 | yn = (yn + lambda) / 4; | |
254 | zn = (zn + lambda) / 4; | |
255 | pn = (pn + lambda) / 4; | |
256 | delta /= 64; | |
257 | } | |
258 | ||
259 | T X = fmn * (A0 - x) / An; | |
260 | T Y = fmn * (A0 - y) / An; | |
261 | T Z = fmn * (A0 - z) / An; | |
262 | T P = (-X - Y - Z) / 2; | |
263 | T E2 = X * Y + X * Z + Y * Z - 3 * P * P; | |
264 | T E3 = X * Y * Z + 2 * E2 * P + 4 * P * P * P; | |
265 | T E4 = (2 * X * Y * Z + E2 * P + 3 * P * P * P) * P; | |
266 | T E5 = X * Y * Z * P * P; | |
267 | T result = fmn * pow(An, T(-3) / 2) * | |
268 | (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16 | |
269 | + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68); | |
270 | ||
271 | result += 6 * RC_sum; | |
272 | return result; | |
273 | } | |
274 | ||
275 | } // namespace detail | |
276 | ||
277 | template <class T1, class T2, class T3, class T4, class Policy> | |
278 | inline typename tools::promote_args<T1, T2, T3, T4>::type | |
279 | ellint_rj(T1 x, T2 y, T3 z, T4 p, const Policy& pol) | |
280 | { | |
281 | typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type; | |
282 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
283 | return policies::checked_narrowing_cast<result_type, Policy>( | |
284 | detail::ellint_rj_imp( | |
285 | static_cast<value_type>(x), | |
286 | static_cast<value_type>(y), | |
287 | static_cast<value_type>(z), | |
288 | static_cast<value_type>(p), | |
289 | pol), "boost::math::ellint_rj<%1%>(%1%,%1%,%1%,%1%)"); | |
290 | } | |
291 | ||
292 | template <class T1, class T2, class T3, class T4> | |
293 | inline typename tools::promote_args<T1, T2, T3, T4>::type | |
294 | ellint_rj(T1 x, T2 y, T3 z, T4 p) | |
295 | { | |
296 | return ellint_rj(x, y, z, p, policies::policy<>()); | |
297 | } | |
298 | ||
299 | }} // namespaces | |
300 | ||
301 | #endif // BOOST_MATH_ELLINT_RJ_HPP | |
302 |