1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006-2024 Free Software Foundation, Inc.
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)
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.
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.
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/>.
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}
31 // ISO C++ 14882 TR1: 5.2 Special functions
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
42 #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC
43 #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1
45 namespace std _GLIBCXX_VISIBILITY(default)
47 _GLIBCXX_BEGIN_NAMESPACE_VERSION
49 #if _GLIBCXX_USE_STD_SPEC_FUNCS
50 #elif defined(_GLIBCXX_TR1_CMATH)
54 # error do not include this header directly, use <cmath> or <tr1/cmath>
56 // [5.2] Special functions
58 // Implementation-space details.
62 * @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$
65 * The Carlson elliptic function of the first kind is defined by:
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}}
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.
76 template<typename _Tp>
78 __ellint_rf(_Tp __x, _Tp __y, _Tp __z)
80 const _Tp __min = std::numeric_limits<_Tp>::min();
81 const _Tp __lolim = _Tp(5) * __min;
83 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
84 std::__throw_domain_error(__N("Argument less than zero "
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"));
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);
101 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
102 const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
104 _Tp __xndev, __yndev, __zndev;
106 const unsigned int __max_iter = 100;
107 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
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)
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);
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
132 return __s / std::sqrt(__mu);
138 * @brief Return the complete elliptic integral of the first kind
139 * @f$ K(k) @f$ by series expansion.
141 * The complete elliptic integral of the first kind is defined as
143 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
144 * {\sqrt{1 - k^2sin^2\theta}}
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.
150 * @param __k The argument of the complete elliptic function.
151 * @return The complete elliptic function of the first kind.
153 template<typename _Tp>
155 __comp_ellint_1_series(_Tp __k)
158 const _Tp __kk = __k * __k;
160 _Tp __term = __kk / _Tp(4);
161 _Tp __sum = _Tp(1) + __term;
163 const unsigned int __max_iter = 1000;
164 for (unsigned int __i = 2; __i < __max_iter; ++__i)
166 __term *= (2 * __i - 1) * __kk / (2 * __i);
167 if (__term < std::numeric_limits<_Tp>::epsilon())
172 return __numeric_constants<_Tp>::__pi_2() * __sum;
177 * @brief Return the complete elliptic integral of the first kind
178 * @f$ K(k) @f$ using the Carlson formulation.
180 * The complete elliptic integral of the first kind is defined as
182 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
183 * {\sqrt{1 - k^2 sin^2\theta}}
185 * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
188 * @param __k The argument of the complete elliptic function.
189 * @return The complete elliptic function of the first kind.
191 template<typename _Tp>
193 __comp_ellint_1(_Tp __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();
201 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
206 * @brief Return the incomplete elliptic integral of the first kind
207 * @f$ F(k,\phi) @f$ using the Carlson formulation.
209 * The incomplete elliptic integral of the first kind is defined as
211 * F(k,\phi) = \int_0^{\phi}\frac{d\theta}
212 * {\sqrt{1 - k^2 sin^2\theta}}
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.
219 template<typename _Tp>
221 __ellint_1(_Tp __k, _Tp __phi)
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."));
230 // Reduce phi to -pi/2 < phi < +pi/2.
231 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
233 const _Tp __phi_red = __phi
234 - __n * __numeric_constants<_Tp>::__pi();
236 const _Tp __s = std::sin(__phi_red);
237 const _Tp __c = std::cos(__phi_red);
240 * __ellint_rf(__c * __c,
241 _Tp(1) - __k * __k * __s * __s, _Tp(1));
246 return __F + _Tp(2) * __n * __comp_ellint_1(__k);
252 * @brief Return the complete elliptic integral of the second kind
253 * @f$ E(k) @f$ by series expansion.
255 * The complete elliptic integral of the second kind is defined as
257 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
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.
263 * @param __k The argument of the complete elliptic function.
264 * @return The complete elliptic function of the second kind.
266 template<typename _Tp>
268 __comp_ellint_2_series(_Tp __k)
271 const _Tp __kk = __k * __k;
276 const unsigned int __max_iter = 1000;
277 for (unsigned int __i = 2; __i < __max_iter; ++__i)
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())
284 __sum += __term / __i2m;
287 return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
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
297 * The Carlson elliptic function of the second kind is defined by:
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}}
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)
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.
314 template<typename _Tp>
316 __ellint_rd(_Tp __x, _Tp __y, _Tp __z)
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));
323 if (__x < _Tp(0) || __y < _Tp(0))
324 std::__throw_domain_error(__N("Argument less than zero "
326 else if (__x + __y < __lolim || __z < __lolim)
327 std::__throw_domain_error(__N("Argument too small "
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);
340 _Tp __sigma = _Tp(0);
341 _Tp __power4 = _Tp(1);
344 _Tp __xndev, __yndev, __zndev;
346 const unsigned int __max_iter = 100;
347 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
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)
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));
364 __xn = __c0 * (__xn + __lambda);
365 __yn = __c0 * (__yn + __lambda);
366 __zn = __c0 * (__zn + __lambda);
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
379 + __zndev * (-__c3 * __ec - __zndev * __c4 - __ea));
381 return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
382 / (__mu * std::sqrt(__mu));
388 * @brief Return the complete elliptic integral of the second kind
389 * @f$ E(k) @f$ using the Carlson formulation.
391 * The complete elliptic integral of the second kind is defined as
393 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
396 * @param __k The argument of the complete elliptic function.
397 * @return The complete elliptic function of the second kind.
399 template<typename _Tp>
401 __comp_ellint_2(_Tp __k)
405 return std::numeric_limits<_Tp>::quiet_NaN();
406 else if (std::abs(__k) == 1)
408 else if (std::abs(__k) > _Tp(1))
409 std::__throw_domain_error(__N("Bad argument in __comp_ellint_2."));
412 const _Tp __kk = __k * __k;
414 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
415 - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
421 * @brief Return the incomplete elliptic integral of the second kind
422 * @f$ E(k,\phi) @f$ using the Carlson formulation.
424 * The incomplete elliptic integral of the second kind is defined as
426 * E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
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.
433 template<typename _Tp>
435 __ellint_2(_Tp __k, _Tp __phi)
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."));
444 // Reduce phi to -pi/2 < phi < +pi/2.
445 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
447 const _Tp __phi_red = __phi
448 - __n * __numeric_constants<_Tp>::__pi();
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;
458 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
460 * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1))
466 return __E + _Tp(2) * __n * __comp_ellint_2(__k);
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.
476 * The Carlson elliptic function is defined by:
478 * R_C(x,y) = \frac{1}{2} \int_0^\infty
479 * \frac{dt}{(t + x)^{1/2}(t + y)}
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)
488 * @param __x The first argument.
489 * @param __y The second argument.
490 * @return The Carlson elliptic function.
492 template<typename _Tp>
494 __ellint_rc(_Tp __x, _Tp __y)
496 const _Tp __min = std::numeric_limits<_Tp>::min();
497 const _Tp __lolim = _Tp(5) * __min;
499 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
500 std::__throw_domain_error(__N("Argument less than zero "
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);
513 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
514 const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
518 const unsigned int __max_iter = 100;
519 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
521 __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
522 __sn = (__yn + __mu) / __mu - _Tp(2);
523 if (std::abs(__sn) < __errtol)
525 const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
527 __xn = __c0 * (__xn + __lambda);
528 __yn = __c0 * (__yn + __lambda);
531 _Tp __s = __sn * __sn
532 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
534 return (_Tp(1) + __s) / std::sqrt(__mu);
540 * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$
543 * The Carlson elliptic function of the third kind is defined by:
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)}
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)
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.
561 template<typename _Tp>
563 __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p)
565 const _Tp __min = std::numeric_limits<_Tp>::min();
566 const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3));
568 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
569 std::__throw_domain_error(__N("Argument less than zero "
571 else if (__x + __y < __lolim || __x + __z < __lolim
572 || __y + __z < __lolim || __p < __lolim)
573 std::__throw_domain_error(__N("Argument too small "
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);
587 _Tp __sigma = _Tp(0);
588 _Tp __power4 = _Tp(1);
590 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
591 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
594 _Tp __xndev, __yndev, __zndev, __pndev;
596 const unsigned int __max_iter = 100;
597 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
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)
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)
619 __sigma += __power4 * __ellint_rc(__alpha2, __beta);
621 __xn = __c0 * (__xn + __lambda);
622 __yn = __c0 * (__yn + __lambda);
623 __zn = __c0 * (__zn + __lambda);
624 __pn = __c0 * (__pn + __lambda);
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;
639 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
640 / (__mu * std::sqrt(__mu));
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.
650 * The complete elliptic integral of the third kind is defined as
652 * \Pi(k,\nu) = \int_0^{\pi/2}
654 * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
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.
661 template<typename _Tp>
663 __comp_ellint_3(_Tp __k, _Tp __nu)
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."));
674 const _Tp __kk = __k * __k;
676 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
678 * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) - __nu)
685 * @brief Return the incomplete elliptic integral of the third kind
686 * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation.
688 * The incomplete elliptic integral of the third kind is defined as
690 * \Pi(k,\nu,\phi) = \int_0^{\phi}
692 * {(1 - \nu \sin^2\theta)
693 * \sqrt{1 - k^2 \sin^2\theta}}
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.
701 template<typename _Tp>
703 __ellint_3(_Tp __k, _Tp __nu, _Tp __phi)
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."));
712 // Reduce phi to -pi/2 < phi < +pi/2.
713 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
715 const _Tp __phi_red = __phi
716 - __n * __numeric_constants<_Tp>::__pi();
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;
726 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
728 * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
729 _Tp(1) - __nu * __ss) / _Tp(3);
734 return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
737 } // namespace __detail
738 #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)
742 _GLIBCXX_END_NAMESPACE_VERSION
745 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC