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