1 | // Special functions -*- C++ -*- |
2 | |
3 | // Copyright (C) 2006-2021 Free Software Foundation, Inc. |
4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free |
6 | // software; you can redistribute it and/or modify it under the |
7 | // terms of the GNU General Public License as published by the |
8 | // Free Software Foundation; either version 3, or (at your option) |
9 | // any later version. |
10 | // |
11 | // This library is distributed in the hope that it will be useful, |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | // GNU General Public License for more details. |
15 | // |
16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version |
18 | // 3.1, as published by the Free Software Foundation. |
19 | |
20 | // You should have received a copy of the GNU General Public License and |
21 | // a copy of the GCC Runtime Library Exception along with this program; |
22 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
23 | // <http://www.gnu.org/licenses/>. |
24 | |
25 | /** @file tr1/ell_integral.tcc |
26 | * This is an internal header file, included by other library headers. |
27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
28 | */ |
29 | |
30 | // |
31 | // ISO C++ 14882 TR1: 5.2 Special functions |
32 | // |
33 | |
34 | // Written by Edward Smith-Rowland based on: |
35 | // (1) B. C. Carlson Numer. Math. 33, 1 (1979) |
36 | // (2) B. C. Carlson, Special Functions of Applied Mathematics (1977) |
37 | // (3) The Gnu Scientific Library, http://www.gnu.org/software/gsl |
38 | // (4) Numerical Recipes in C, 2nd ed, by W. H. Press, S. A. Teukolsky, |
39 | // W. T. Vetterling, B. P. Flannery, Cambridge University Press |
40 | // (1992), pp. 261-269 |
41 | |
42 | #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
43 | #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1 |
44 | |
45 | namespace std _GLIBCXX_VISIBILITY(default) |
46 | { |
47 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
48 | |
49 | #if _GLIBCXX_USE_STD_SPEC_FUNCS |
50 | #elif defined(_GLIBCXX_TR1_CMATH) |
51 | namespace tr1 |
52 | { |
53 | #else |
54 | # error do not include this header directly, use <cmath> or <tr1/cmath> |
55 | #endif |
56 | // [5.2] Special functions |
57 | |
58 | // Implementation-space details. |
59 | namespace __detail |
60 | { |
61 | /** |
62 | * @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$ |
63 | * of the first kind. |
64 | * |
65 | * The Carlson elliptic function of the first kind is defined by: |
66 | * @f[ |
67 | * R_F(x,y,z) = \frac{1}{2} \int_0^\infty |
68 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}} |
69 | * @f] |
70 | * |
71 | * @param __x The first of three symmetric arguments. |
72 | * @param __y The second of three symmetric arguments. |
73 | * @param __z The third of three symmetric arguments. |
74 | * @return The Carlson elliptic function of the first kind. |
75 | */ |
76 | template<typename _Tp> |
77 | _Tp |
78 | __ellint_rf(_Tp __x, _Tp __y, _Tp __z) |
79 | { |
80 | const _Tp __min = std::numeric_limits<_Tp>::min(); |
81 | const _Tp __lolim = _Tp(5) * __min; |
82 | |
83 | if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) |
84 | std::__throw_domain_error(__N("Argument less than zero " |
85 | "in __ellint_rf." )); |
86 | else if (__x + __y < __lolim || __x + __z < __lolim |
87 | || __y + __z < __lolim) |
88 | std::__throw_domain_error(__N("Argument too small in __ellint_rf" )); |
89 | else |
90 | { |
91 | const _Tp __c0 = _Tp(1) / _Tp(4); |
92 | const _Tp __c1 = _Tp(1) / _Tp(24); |
93 | const _Tp __c2 = _Tp(1) / _Tp(10); |
94 | const _Tp __c3 = _Tp(3) / _Tp(44); |
95 | const _Tp __c4 = _Tp(1) / _Tp(14); |
96 | |
97 | _Tp __xn = __x; |
98 | _Tp __yn = __y; |
99 | _Tp __zn = __z; |
100 | |
101 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); |
102 | const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6)); |
103 | _Tp __mu; |
104 | _Tp __xndev, __yndev, __zndev; |
105 | |
106 | const unsigned int __max_iter = 100; |
107 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) |
108 | { |
109 | __mu = (__xn + __yn + __zn) / _Tp(3); |
110 | __xndev = 2 - (__mu + __xn) / __mu; |
111 | __yndev = 2 - (__mu + __yn) / __mu; |
112 | __zndev = 2 - (__mu + __zn) / __mu; |
113 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); |
114 | __epsilon = std::max(__epsilon, std::abs(__zndev)); |
115 | if (__epsilon < __errtol) |
116 | break; |
117 | const _Tp __xnroot = std::sqrt(__xn); |
118 | const _Tp __ynroot = std::sqrt(__yn); |
119 | const _Tp __znroot = std::sqrt(__zn); |
120 | const _Tp __lambda = __xnroot * (__ynroot + __znroot) |
121 | + __ynroot * __znroot; |
122 | __xn = __c0 * (__xn + __lambda); |
123 | __yn = __c0 * (__yn + __lambda); |
124 | __zn = __c0 * (__zn + __lambda); |
125 | } |
126 | |
127 | const _Tp __e2 = __xndev * __yndev - __zndev * __zndev; |
128 | const _Tp __e3 = __xndev * __yndev * __zndev; |
129 | const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2 |
130 | + __c4 * __e3; |
131 | |
132 | return __s / std::sqrt(__mu); |
133 | } |
134 | } |
135 | |
136 | |
137 | /** |
138 | * @brief Return the complete elliptic integral of the first kind |
139 | * @f$ K(k) @f$ by series expansion. |
140 | * |
141 | * The complete elliptic integral of the first kind is defined as |
142 | * @f[ |
143 | * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} |
144 | * {\sqrt{1 - k^2sin^2\theta}} |
145 | * @f] |
146 | * |
147 | * This routine is not bad as long as |k| is somewhat smaller than 1 |
148 | * but is not is good as the Carlson elliptic integral formulation. |
149 | * |
150 | * @param __k The argument of the complete elliptic function. |
151 | * @return The complete elliptic function of the first kind. |
152 | */ |
153 | template<typename _Tp> |
154 | _Tp |
155 | __comp_ellint_1_series(_Tp __k) |
156 | { |
157 | |
158 | const _Tp __kk = __k * __k; |
159 | |
160 | _Tp __term = __kk / _Tp(4); |
161 | _Tp __sum = _Tp(1) + __term; |
162 | |
163 | const unsigned int __max_iter = 1000; |
164 | for (unsigned int __i = 2; __i < __max_iter; ++__i) |
165 | { |
166 | __term *= (2 * __i - 1) * __kk / (2 * __i); |
167 | if (__term < std::numeric_limits<_Tp>::epsilon()) |
168 | break; |
169 | __sum += __term; |
170 | } |
171 | |
172 | return __numeric_constants<_Tp>::__pi_2() * __sum; |
173 | } |
174 | |
175 | |
176 | /** |
177 | * @brief Return the complete elliptic integral of the first kind |
178 | * @f$ K(k) @f$ using the Carlson formulation. |
179 | * |
180 | * The complete elliptic integral of the first kind is defined as |
181 | * @f[ |
182 | * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} |
183 | * {\sqrt{1 - k^2 sin^2\theta}} |
184 | * @f] |
185 | * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the |
186 | * first kind. |
187 | * |
188 | * @param __k The argument of the complete elliptic function. |
189 | * @return The complete elliptic function of the first kind. |
190 | */ |
191 | template<typename _Tp> |
192 | _Tp |
193 | __comp_ellint_1(_Tp __k) |
194 | { |
195 | |
196 | if (__isnan(__k)) |
197 | return std::numeric_limits<_Tp>::quiet_NaN(); |
198 | else if (std::abs(__k) >= _Tp(1)) |
199 | return std::numeric_limits<_Tp>::quiet_NaN(); |
200 | else |
201 | return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1)); |
202 | } |
203 | |
204 | |
205 | /** |
206 | * @brief Return the incomplete elliptic integral of the first kind |
207 | * @f$ F(k,\phi) @f$ using the Carlson formulation. |
208 | * |
209 | * The incomplete elliptic integral of the first kind is defined as |
210 | * @f[ |
211 | * F(k,\phi) = \int_0^{\phi}\frac{d\theta} |
212 | * {\sqrt{1 - k^2 sin^2\theta}} |
213 | * @f] |
214 | * |
215 | * @param __k The argument of the elliptic function. |
216 | * @param __phi The integral limit argument of the elliptic function. |
217 | * @return The elliptic function of the first kind. |
218 | */ |
219 | template<typename _Tp> |
220 | _Tp |
221 | __ellint_1(_Tp __k, _Tp __phi) |
222 | { |
223 | |
224 | if (__isnan(__k) || __isnan(__phi)) |
225 | return std::numeric_limits<_Tp>::quiet_NaN(); |
226 | else if (std::abs(__k) > _Tp(1)) |
227 | std::__throw_domain_error(__N("Bad argument in __ellint_1." )); |
228 | else |
229 | { |
230 | // Reduce phi to -pi/2 < phi < +pi/2. |
231 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() |
232 | + _Tp(0.5L)); |
233 | const _Tp __phi_red = __phi |
234 | - __n * __numeric_constants<_Tp>::__pi(); |
235 | |
236 | const _Tp __s = std::sin(__phi_red); |
237 | const _Tp __c = std::cos(__phi_red); |
238 | |
239 | const _Tp __F = __s |
240 | * __ellint_rf(__c * __c, |
241 | _Tp(1) - __k * __k * __s * __s, _Tp(1)); |
242 | |
243 | if (__n == 0) |
244 | return __F; |
245 | else |
246 | return __F + _Tp(2) * __n * __comp_ellint_1(__k); |
247 | } |
248 | } |
249 | |
250 | |
251 | /** |
252 | * @brief Return the complete elliptic integral of the second kind |
253 | * @f$ E(k) @f$ by series expansion. |
254 | * |
255 | * The complete elliptic integral of the second kind is defined as |
256 | * @f[ |
257 | * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} |
258 | * @f] |
259 | * |
260 | * This routine is not bad as long as |k| is somewhat smaller than 1 |
261 | * but is not is good as the Carlson elliptic integral formulation. |
262 | * |
263 | * @param __k The argument of the complete elliptic function. |
264 | * @return The complete elliptic function of the second kind. |
265 | */ |
266 | template<typename _Tp> |
267 | _Tp |
268 | __comp_ellint_2_series(_Tp __k) |
269 | { |
270 | |
271 | const _Tp __kk = __k * __k; |
272 | |
273 | _Tp __term = __kk; |
274 | _Tp __sum = __term; |
275 | |
276 | const unsigned int __max_iter = 1000; |
277 | for (unsigned int __i = 2; __i < __max_iter; ++__i) |
278 | { |
279 | const _Tp __i2m = 2 * __i - 1; |
280 | const _Tp __i2 = 2 * __i; |
281 | __term *= __i2m * __i2m * __kk / (__i2 * __i2); |
282 | if (__term < std::numeric_limits<_Tp>::epsilon()) |
283 | break; |
284 | __sum += __term / __i2m; |
285 | } |
286 | |
287 | return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum); |
288 | } |
289 | |
290 | |
291 | /** |
292 | * @brief Return the Carlson elliptic function of the second kind |
293 | * @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where |
294 | * @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function |
295 | * of the third kind. |
296 | * |
297 | * The Carlson elliptic function of the second kind is defined by: |
298 | * @f[ |
299 | * R_D(x,y,z) = \frac{3}{2} \int_0^\infty |
300 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}} |
301 | * @f] |
302 | * |
303 | * Based on Carlson's algorithms: |
304 | * - B. C. Carlson Numer. Math. 33, 1 (1979) |
305 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) |
306 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
307 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
308 | * |
309 | * @param __x The first of two symmetric arguments. |
310 | * @param __y The second of two symmetric arguments. |
311 | * @param __z The third argument. |
312 | * @return The Carlson elliptic function of the second kind. |
313 | */ |
314 | template<typename _Tp> |
315 | _Tp |
316 | __ellint_rd(_Tp __x, _Tp __y, _Tp __z) |
317 | { |
318 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); |
319 | const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); |
320 | const _Tp __max = std::numeric_limits<_Tp>::max(); |
321 | const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3)); |
322 | |
323 | if (__x < _Tp(0) || __y < _Tp(0)) |
324 | std::__throw_domain_error(__N("Argument less than zero " |
325 | "in __ellint_rd." )); |
326 | else if (__x + __y < __lolim || __z < __lolim) |
327 | std::__throw_domain_error(__N("Argument too small " |
328 | "in __ellint_rd." )); |
329 | else |
330 | { |
331 | const _Tp __c0 = _Tp(1) / _Tp(4); |
332 | const _Tp __c1 = _Tp(3) / _Tp(14); |
333 | const _Tp __c2 = _Tp(1) / _Tp(6); |
334 | const _Tp __c3 = _Tp(9) / _Tp(22); |
335 | const _Tp __c4 = _Tp(3) / _Tp(26); |
336 | |
337 | _Tp __xn = __x; |
338 | _Tp __yn = __y; |
339 | _Tp __zn = __z; |
340 | _Tp __sigma = _Tp(0); |
341 | _Tp __power4 = _Tp(1); |
342 | |
343 | _Tp __mu; |
344 | _Tp __xndev, __yndev, __zndev; |
345 | |
346 | const unsigned int __max_iter = 100; |
347 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) |
348 | { |
349 | __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5); |
350 | __xndev = (__mu - __xn) / __mu; |
351 | __yndev = (__mu - __yn) / __mu; |
352 | __zndev = (__mu - __zn) / __mu; |
353 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); |
354 | __epsilon = std::max(__epsilon, std::abs(__zndev)); |
355 | if (__epsilon < __errtol) |
356 | break; |
357 | _Tp __xnroot = std::sqrt(__xn); |
358 | _Tp __ynroot = std::sqrt(__yn); |
359 | _Tp __znroot = std::sqrt(__zn); |
360 | _Tp __lambda = __xnroot * (__ynroot + __znroot) |
361 | + __ynroot * __znroot; |
362 | __sigma += __power4 / (__znroot * (__zn + __lambda)); |
363 | __power4 *= __c0; |
364 | __xn = __c0 * (__xn + __lambda); |
365 | __yn = __c0 * (__yn + __lambda); |
366 | __zn = __c0 * (__zn + __lambda); |
367 | } |
368 | |
369 | _Tp __ea = __xndev * __yndev; |
370 | _Tp __eb = __zndev * __zndev; |
371 | _Tp __ec = __ea - __eb; |
372 | _Tp __ed = __ea - _Tp(6) * __eb; |
373 | _Tp __ef = __ed + __ec + __ec; |
374 | _Tp __s1 = __ed * (-__c1 + __c3 * __ed |
375 | / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef |
376 | / _Tp(2)); |
377 | _Tp __s2 = __zndev |
378 | * (__c2 * __ef |
379 | + __zndev * (-__c3 * __ec - __zndev * __c4 - __ea)); |
380 | |
381 | return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2) |
382 | / (__mu * std::sqrt(__mu)); |
383 | } |
384 | } |
385 | |
386 | |
387 | /** |
388 | * @brief Return the complete elliptic integral of the second kind |
389 | * @f$ E(k) @f$ using the Carlson formulation. |
390 | * |
391 | * The complete elliptic integral of the second kind is defined as |
392 | * @f[ |
393 | * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} |
394 | * @f] |
395 | * |
396 | * @param __k The argument of the complete elliptic function. |
397 | * @return The complete elliptic function of the second kind. |
398 | */ |
399 | template<typename _Tp> |
400 | _Tp |
401 | __comp_ellint_2(_Tp __k) |
402 | { |
403 | |
404 | if (__isnan(__k)) |
405 | return std::numeric_limits<_Tp>::quiet_NaN(); |
406 | else if (std::abs(__k) == 1) |
407 | return _Tp(1); |
408 | else if (std::abs(__k) > _Tp(1)) |
409 | std::__throw_domain_error(__N("Bad argument in __comp_ellint_2." )); |
410 | else |
411 | { |
412 | const _Tp __kk = __k * __k; |
413 | |
414 | return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) |
415 | - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3); |
416 | } |
417 | } |
418 | |
419 | |
420 | /** |
421 | * @brief Return the incomplete elliptic integral of the second kind |
422 | * @f$ E(k,\phi) @f$ using the Carlson formulation. |
423 | * |
424 | * The incomplete elliptic integral of the second kind is defined as |
425 | * @f[ |
426 | * E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta} |
427 | * @f] |
428 | * |
429 | * @param __k The argument of the elliptic function. |
430 | * @param __phi The integral limit argument of the elliptic function. |
431 | * @return The elliptic function of the second kind. |
432 | */ |
433 | template<typename _Tp> |
434 | _Tp |
435 | __ellint_2(_Tp __k, _Tp __phi) |
436 | { |
437 | |
438 | if (__isnan(__k) || __isnan(__phi)) |
439 | return std::numeric_limits<_Tp>::quiet_NaN(); |
440 | else if (std::abs(__k) > _Tp(1)) |
441 | std::__throw_domain_error(__N("Bad argument in __ellint_2." )); |
442 | else |
443 | { |
444 | // Reduce phi to -pi/2 < phi < +pi/2. |
445 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() |
446 | + _Tp(0.5L)); |
447 | const _Tp __phi_red = __phi |
448 | - __n * __numeric_constants<_Tp>::__pi(); |
449 | |
450 | const _Tp __kk = __k * __k; |
451 | const _Tp __s = std::sin(__phi_red); |
452 | const _Tp __ss = __s * __s; |
453 | const _Tp __sss = __ss * __s; |
454 | const _Tp __c = std::cos(__phi_red); |
455 | const _Tp __cc = __c * __c; |
456 | |
457 | const _Tp __E = __s |
458 | * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) |
459 | - __kk * __sss |
460 | * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1)) |
461 | / _Tp(3); |
462 | |
463 | if (__n == 0) |
464 | return __E; |
465 | else |
466 | return __E + _Tp(2) * __n * __comp_ellint_2(__k); |
467 | } |
468 | } |
469 | |
470 | |
471 | /** |
472 | * @brief Return the Carlson elliptic function |
473 | * @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$ |
474 | * is the Carlson elliptic function of the first kind. |
475 | * |
476 | * The Carlson elliptic function is defined by: |
477 | * @f[ |
478 | * R_C(x,y) = \frac{1}{2} \int_0^\infty |
479 | * \frac{dt}{(t + x)^{1/2}(t + y)} |
480 | * @f] |
481 | * |
482 | * Based on Carlson's algorithms: |
483 | * - B. C. Carlson Numer. Math. 33, 1 (1979) |
484 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) |
485 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
486 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
487 | * |
488 | * @param __x The first argument. |
489 | * @param __y The second argument. |
490 | * @return The Carlson elliptic function. |
491 | */ |
492 | template<typename _Tp> |
493 | _Tp |
494 | __ellint_rc(_Tp __x, _Tp __y) |
495 | { |
496 | const _Tp __min = std::numeric_limits<_Tp>::min(); |
497 | const _Tp __lolim = _Tp(5) * __min; |
498 | |
499 | if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim) |
500 | std::__throw_domain_error(__N("Argument less than zero " |
501 | "in __ellint_rc." )); |
502 | else |
503 | { |
504 | const _Tp __c0 = _Tp(1) / _Tp(4); |
505 | const _Tp __c1 = _Tp(1) / _Tp(7); |
506 | const _Tp __c2 = _Tp(9) / _Tp(22); |
507 | const _Tp __c3 = _Tp(3) / _Tp(10); |
508 | const _Tp __c4 = _Tp(3) / _Tp(8); |
509 | |
510 | _Tp __xn = __x; |
511 | _Tp __yn = __y; |
512 | |
513 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); |
514 | const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6)); |
515 | _Tp __mu; |
516 | _Tp __sn; |
517 | |
518 | const unsigned int __max_iter = 100; |
519 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) |
520 | { |
521 | __mu = (__xn + _Tp(2) * __yn) / _Tp(3); |
522 | __sn = (__yn + __mu) / __mu - _Tp(2); |
523 | if (std::abs(__sn) < __errtol) |
524 | break; |
525 | const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn) |
526 | + __yn; |
527 | __xn = __c0 * (__xn + __lambda); |
528 | __yn = __c0 * (__yn + __lambda); |
529 | } |
530 | |
531 | _Tp __s = __sn * __sn |
532 | * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2))); |
533 | |
534 | return (_Tp(1) + __s) / std::sqrt(__mu); |
535 | } |
536 | } |
537 | |
538 | |
539 | /** |
540 | * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$ |
541 | * of the third kind. |
542 | * |
543 | * The Carlson elliptic function of the third kind is defined by: |
544 | * @f[ |
545 | * R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty |
546 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)} |
547 | * @f] |
548 | * |
549 | * Based on Carlson's algorithms: |
550 | * - B. C. Carlson Numer. Math. 33, 1 (1979) |
551 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) |
552 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
553 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
554 | * |
555 | * @param __x The first of three symmetric arguments. |
556 | * @param __y The second of three symmetric arguments. |
557 | * @param __z The third of three symmetric arguments. |
558 | * @param __p The fourth argument. |
559 | * @return The Carlson elliptic function of the fourth kind. |
560 | */ |
561 | template<typename _Tp> |
562 | _Tp |
563 | __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p) |
564 | { |
565 | const _Tp __min = std::numeric_limits<_Tp>::min(); |
566 | const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3)); |
567 | |
568 | if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) |
569 | std::__throw_domain_error(__N("Argument less than zero " |
570 | "in __ellint_rj." )); |
571 | else if (__x + __y < __lolim || __x + __z < __lolim |
572 | || __y + __z < __lolim || __p < __lolim) |
573 | std::__throw_domain_error(__N("Argument too small " |
574 | "in __ellint_rj" )); |
575 | else |
576 | { |
577 | const _Tp __c0 = _Tp(1) / _Tp(4); |
578 | const _Tp __c1 = _Tp(3) / _Tp(14); |
579 | const _Tp __c2 = _Tp(1) / _Tp(3); |
580 | const _Tp __c3 = _Tp(3) / _Tp(22); |
581 | const _Tp __c4 = _Tp(3) / _Tp(26); |
582 | |
583 | _Tp __xn = __x; |
584 | _Tp __yn = __y; |
585 | _Tp __zn = __z; |
586 | _Tp __pn = __p; |
587 | _Tp __sigma = _Tp(0); |
588 | _Tp __power4 = _Tp(1); |
589 | |
590 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); |
591 | const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); |
592 | |
593 | _Tp __mu; |
594 | _Tp __xndev, __yndev, __zndev, __pndev; |
595 | |
596 | const unsigned int __max_iter = 100; |
597 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) |
598 | { |
599 | __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5); |
600 | __xndev = (__mu - __xn) / __mu; |
601 | __yndev = (__mu - __yn) / __mu; |
602 | __zndev = (__mu - __zn) / __mu; |
603 | __pndev = (__mu - __pn) / __mu; |
604 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); |
605 | __epsilon = std::max(__epsilon, std::abs(__zndev)); |
606 | __epsilon = std::max(__epsilon, std::abs(__pndev)); |
607 | if (__epsilon < __errtol) |
608 | break; |
609 | const _Tp __xnroot = std::sqrt(__xn); |
610 | const _Tp __ynroot = std::sqrt(__yn); |
611 | const _Tp __znroot = std::sqrt(__zn); |
612 | const _Tp __lambda = __xnroot * (__ynroot + __znroot) |
613 | + __ynroot * __znroot; |
614 | const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) |
615 | + __xnroot * __ynroot * __znroot; |
616 | const _Tp __alpha2 = __alpha1 * __alpha1; |
617 | const _Tp __beta = __pn * (__pn + __lambda) |
618 | * (__pn + __lambda); |
619 | __sigma += __power4 * __ellint_rc(__alpha2, __beta); |
620 | __power4 *= __c0; |
621 | __xn = __c0 * (__xn + __lambda); |
622 | __yn = __c0 * (__yn + __lambda); |
623 | __zn = __c0 * (__zn + __lambda); |
624 | __pn = __c0 * (__pn + __lambda); |
625 | } |
626 | |
627 | _Tp __ea = __xndev * (__yndev + __zndev) + __yndev * __zndev; |
628 | _Tp __eb = __xndev * __yndev * __zndev; |
629 | _Tp __ec = __pndev * __pndev; |
630 | _Tp __e2 = __ea - _Tp(3) * __ec; |
631 | _Tp __e3 = __eb + _Tp(2) * __pndev * (__ea - __ec); |
632 | _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4) |
633 | - _Tp(3) * __c4 * __e3 / _Tp(2)); |
634 | _Tp __s2 = __eb * (__c2 / _Tp(2) |
635 | + __pndev * (-__c3 - __c3 + __pndev * __c4)); |
636 | _Tp __s3 = __pndev * __ea * (__c2 - __pndev * __c3) |
637 | - __c2 * __pndev * __ec; |
638 | |
639 | return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3) |
640 | / (__mu * std::sqrt(__mu)); |
641 | } |
642 | } |
643 | |
644 | |
645 | /** |
646 | * @brief Return the complete elliptic integral of the third kind |
647 | * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the |
648 | * Carlson formulation. |
649 | * |
650 | * The complete elliptic integral of the third kind is defined as |
651 | * @f[ |
652 | * \Pi(k,\nu) = \int_0^{\pi/2} |
653 | * \frac{d\theta} |
654 | * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}} |
655 | * @f] |
656 | * |
657 | * @param __k The argument of the elliptic function. |
658 | * @param __nu The second argument of the elliptic function. |
659 | * @return The complete elliptic function of the third kind. |
660 | */ |
661 | template<typename _Tp> |
662 | _Tp |
663 | __comp_ellint_3(_Tp __k, _Tp __nu) |
664 | { |
665 | |
666 | if (__isnan(__k) || __isnan(__nu)) |
667 | return std::numeric_limits<_Tp>::quiet_NaN(); |
668 | else if (__nu == _Tp(1)) |
669 | return std::numeric_limits<_Tp>::infinity(); |
670 | else if (std::abs(__k) > _Tp(1)) |
671 | std::__throw_domain_error(__N("Bad argument in __comp_ellint_3." )); |
672 | else |
673 | { |
674 | const _Tp __kk = __k * __k; |
675 | |
676 | return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) |
677 | + __nu |
678 | * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) - __nu) |
679 | / _Tp(3); |
680 | } |
681 | } |
682 | |
683 | |
684 | /** |
685 | * @brief Return the incomplete elliptic integral of the third kind |
686 | * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation. |
687 | * |
688 | * The incomplete elliptic integral of the third kind is defined as |
689 | * @f[ |
690 | * \Pi(k,\nu,\phi) = \int_0^{\phi} |
691 | * \frac{d\theta} |
692 | * {(1 - \nu \sin^2\theta) |
693 | * \sqrt{1 - k^2 \sin^2\theta}} |
694 | * @f] |
695 | * |
696 | * @param __k The argument of the elliptic function. |
697 | * @param __nu The second argument of the elliptic function. |
698 | * @param __phi The integral limit argument of the elliptic function. |
699 | * @return The elliptic function of the third kind. |
700 | */ |
701 | template<typename _Tp> |
702 | _Tp |
703 | __ellint_3(_Tp __k, _Tp __nu, _Tp __phi) |
704 | { |
705 | |
706 | if (__isnan(__k) || __isnan(__nu) || __isnan(__phi)) |
707 | return std::numeric_limits<_Tp>::quiet_NaN(); |
708 | else if (std::abs(__k) > _Tp(1)) |
709 | std::__throw_domain_error(__N("Bad argument in __ellint_3." )); |
710 | else |
711 | { |
712 | // Reduce phi to -pi/2 < phi < +pi/2. |
713 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() |
714 | + _Tp(0.5L)); |
715 | const _Tp __phi_red = __phi |
716 | - __n * __numeric_constants<_Tp>::__pi(); |
717 | |
718 | const _Tp __kk = __k * __k; |
719 | const _Tp __s = std::sin(__phi_red); |
720 | const _Tp __ss = __s * __s; |
721 | const _Tp __sss = __ss * __s; |
722 | const _Tp __c = std::cos(__phi_red); |
723 | const _Tp __cc = __c * __c; |
724 | |
725 | const _Tp __Pi = __s |
726 | * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) |
727 | + __nu * __sss |
728 | * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1), |
729 | _Tp(1) - __nu * __ss) / _Tp(3); |
730 | |
731 | if (__n == 0) |
732 | return __Pi; |
733 | else |
734 | return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu); |
735 | } |
736 | } |
737 | } // namespace __detail |
738 | #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
739 | } // namespace tr1 |
740 | #endif |
741 | |
742 | _GLIBCXX_END_NAMESPACE_VERSION |
743 | } |
744 | |
745 | #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
746 | |
747 | |