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
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 __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