1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006-2013 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)
49 // [5.2] Special functions
51 // Implementation-space details.
54 _GLIBCXX_BEGIN_NAMESPACE_VERSION
57 * @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$
60 * The Carlson elliptic function of the first kind is defined by:
62 * R_F(x,y,z) = \frac{1}{2} \int_0^\infty
63 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}}
66 * @param __x The first of three symmetric arguments.
67 * @param __y The second of three symmetric arguments.
68 * @param __z The third of three symmetric arguments.
69 * @return The Carlson elliptic function of the first kind.
71 template<typename _Tp>
73 __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
75 const _Tp __min = std::numeric_limits<_Tp>::min();
76 const _Tp __max = std::numeric_limits<_Tp>::max();
77 const _Tp __lolim = _Tp(5) * __min;
78 const _Tp __uplim = __max / _Tp(5);
80 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
81 std::__throw_domain_error(__N("Argument less than zero "
83 else if (__x + __y < __lolim || __x + __z < __lolim
84 || __y + __z < __lolim)
85 std::__throw_domain_error(__N("Argument too small in __ellint_rf"));
88 const _Tp __c0 = _Tp(1) / _Tp(4);
89 const _Tp __c1 = _Tp(1) / _Tp(24);
90 const _Tp __c2 = _Tp(1) / _Tp(10);
91 const _Tp __c3 = _Tp(3) / _Tp(44);
92 const _Tp __c4 = _Tp(1) / _Tp(14);
98 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
99 const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
101 _Tp __xndev, __yndev, __zndev;
103 const unsigned int __max_iter = 100;
104 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
106 __mu = (__xn + __yn + __zn) / _Tp(3);
107 __xndev = 2 - (__mu + __xn) / __mu;
108 __yndev = 2 - (__mu + __yn) / __mu;
109 __zndev = 2 - (__mu + __zn) / __mu;
110 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
111 __epsilon = std::max(__epsilon, std::abs(__zndev));
112 if (__epsilon < __errtol)
114 const _Tp __xnroot = std::sqrt(__xn);
115 const _Tp __ynroot = std::sqrt(__yn);
116 const _Tp __znroot = std::sqrt(__zn);
117 const _Tp __lambda = __xnroot * (__ynroot + __znroot)
118 + __ynroot * __znroot;
119 __xn = __c0 * (__xn + __lambda);
120 __yn = __c0 * (__yn + __lambda);
121 __zn = __c0 * (__zn + __lambda);
124 const _Tp __e2 = __xndev * __yndev - __zndev * __zndev;
125 const _Tp __e3 = __xndev * __yndev * __zndev;
126 const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2
129 return __s / std::sqrt(__mu);
135 * @brief Return the complete elliptic integral of the first kind
136 * @f$ K(k) @f$ by series expansion.
138 * The complete elliptic integral of the first kind is defined as
140 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
141 * {\sqrt{1 - k^2sin^2\theta}}
144 * This routine is not bad as long as |k| is somewhat smaller than 1
145 * but is not is good as the Carlson elliptic integral formulation.
147 * @param __k The argument of the complete elliptic function.
148 * @return The complete elliptic function of the first kind.
150 template<typename _Tp>
152 __comp_ellint_1_series(const _Tp __k)
155 const _Tp __kk = __k * __k;
157 _Tp __term = __kk / _Tp(4);
158 _Tp __sum = _Tp(1) + __term;
160 const unsigned int __max_iter = 1000;
161 for (unsigned int __i = 2; __i < __max_iter; ++__i)
163 __term *= (2 * __i - 1) * __kk / (2 * __i);
164 if (__term < std::numeric_limits<_Tp>::epsilon())
169 return __numeric_constants<_Tp>::__pi_2() * __sum;
174 * @brief Return the complete elliptic integral of the first kind
175 * @f$ K(k) @f$ using the Carlson formulation.
177 * The complete elliptic integral of the first kind is defined as
179 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
180 * {\sqrt{1 - k^2 sin^2\theta}}
182 * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
185 * @param __k The argument of the complete elliptic function.
186 * @return The complete elliptic function of the first kind.
188 template<typename _Tp>
190 __comp_ellint_1(const _Tp __k)
194 return std::numeric_limits<_Tp>::quiet_NaN();
195 else if (std::abs(__k) >= _Tp(1))
196 return std::numeric_limits<_Tp>::quiet_NaN();
198 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
203 * @brief Return the incomplete elliptic integral of the first kind
204 * @f$ F(k,\phi) @f$ using the Carlson formulation.
206 * The incomplete elliptic integral of the first kind is defined as
208 * F(k,\phi) = \int_0^{\phi}\frac{d\theta}
209 * {\sqrt{1 - k^2 sin^2\theta}}
212 * @param __k The argument of the elliptic function.
213 * @param __phi The integral limit argument of the elliptic function.
214 * @return The elliptic function of the first kind.
216 template<typename _Tp>
218 __ellint_1(const _Tp __k, const _Tp __phi)
221 if (__isnan(__k) || __isnan(__phi))
222 return std::numeric_limits<_Tp>::quiet_NaN();
223 else if (std::abs(__k) > _Tp(1))
224 std::__throw_domain_error(__N("Bad argument in __ellint_1."));
227 // Reduce phi to -pi/2 < phi < +pi/2.
228 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
230 const _Tp __phi_red = __phi
231 - __n * __numeric_constants<_Tp>::__pi();
233 const _Tp __s = std::sin(__phi_red);
234 const _Tp __c = std::cos(__phi_red);
237 * __ellint_rf(__c * __c,
238 _Tp(1) - __k * __k * __s * __s, _Tp(1));
243 return __F + _Tp(2) * __n * __comp_ellint_1(__k);
249 * @brief Return the complete elliptic integral of the second kind
250 * @f$ E(k) @f$ by series expansion.
252 * The complete elliptic integral of the second kind is defined as
254 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
257 * This routine is not bad as long as |k| is somewhat smaller than 1
258 * but is not is good as the Carlson elliptic integral formulation.
260 * @param __k The argument of the complete elliptic function.
261 * @return The complete elliptic function of the second kind.
263 template<typename _Tp>
265 __comp_ellint_2_series(const _Tp __k)
268 const _Tp __kk = __k * __k;
273 const unsigned int __max_iter = 1000;
274 for (unsigned int __i = 2; __i < __max_iter; ++__i)
276 const _Tp __i2m = 2 * __i - 1;
277 const _Tp __i2 = 2 * __i;
278 __term *= __i2m * __i2m * __kk / (__i2 * __i2);
279 if (__term < std::numeric_limits<_Tp>::epsilon())
281 __sum += __term / __i2m;
284 return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
289 * @brief Return the Carlson elliptic function of the second kind
290 * @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where
291 * @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function
294 * The Carlson elliptic function of the second kind is defined by:
296 * R_D(x,y,z) = \frac{3}{2} \int_0^\infty
297 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}}
300 * Based on Carlson's algorithms:
301 * - B. C. Carlson Numer. Math. 33, 1 (1979)
302 * - B. C. Carlson, Special Functions of Applied Mathematics (1977)
303 * - Numerical Recipes in C, 2nd ed, pp. 261-269,
304 * by Press, Teukolsky, Vetterling, Flannery (1992)
306 * @param __x The first of two symmetric arguments.
307 * @param __y The second of two symmetric arguments.
308 * @param __z The third argument.
309 * @return The Carlson elliptic function of the second kind.
311 template<typename _Tp>
313 __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
315 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
316 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
317 const _Tp __min = std::numeric_limits<_Tp>::min();
318 const _Tp __max = std::numeric_limits<_Tp>::max();
319 const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3));
320 const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3));
322 if (__x < _Tp(0) || __y < _Tp(0))
323 std::__throw_domain_error(__N("Argument less than zero "
325 else if (__x + __y < __lolim || __z < __lolim)
326 std::__throw_domain_error(__N("Argument too small "
330 const _Tp __c0 = _Tp(1) / _Tp(4);
331 const _Tp __c1 = _Tp(3) / _Tp(14);
332 const _Tp __c2 = _Tp(1) / _Tp(6);
333 const _Tp __c3 = _Tp(9) / _Tp(22);
334 const _Tp __c4 = _Tp(3) / _Tp(26);
339 _Tp __sigma = _Tp(0);
340 _Tp __power4 = _Tp(1);
343 _Tp __xndev, __yndev, __zndev;
345 const unsigned int __max_iter = 100;
346 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
348 __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5);
349 __xndev = (__mu - __xn) / __mu;
350 __yndev = (__mu - __yn) / __mu;
351 __zndev = (__mu - __zn) / __mu;
352 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
353 __epsilon = std::max(__epsilon, std::abs(__zndev));
354 if (__epsilon < __errtol)
356 _Tp __xnroot = std::sqrt(__xn);
357 _Tp __ynroot = std::sqrt(__yn);
358 _Tp __znroot = std::sqrt(__zn);
359 _Tp __lambda = __xnroot * (__ynroot + __znroot)
360 + __ynroot * __znroot;
361 __sigma += __power4 / (__znroot * (__zn + __lambda));
363 __xn = __c0 * (__xn + __lambda);
364 __yn = __c0 * (__yn + __lambda);
365 __zn = __c0 * (__zn + __lambda);
368 // Note: __ea is an SPU badname.
369 _Tp __eaa = __xndev * __yndev;
370 _Tp __eb = __zndev * __zndev;
371 _Tp __ec = __eaa - __eb;
372 _Tp __ed = __eaa - _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 - __eaa));
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(const _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(const _Tp __k, const _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(const _Tp __x, const _Tp __y)
496 const _Tp __min = std::numeric_limits<_Tp>::min();
497 const _Tp __max = std::numeric_limits<_Tp>::max();
498 const _Tp __lolim = _Tp(5) * __min;
499 const _Tp __uplim = __max / _Tp(5);
501 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
502 std::__throw_domain_error(__N("Argument less than zero "
506 const _Tp __c0 = _Tp(1) / _Tp(4);
507 const _Tp __c1 = _Tp(1) / _Tp(7);
508 const _Tp __c2 = _Tp(9) / _Tp(22);
509 const _Tp __c3 = _Tp(3) / _Tp(10);
510 const _Tp __c4 = _Tp(3) / _Tp(8);
515 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
516 const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
520 const unsigned int __max_iter = 100;
521 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
523 __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
524 __sn = (__yn + __mu) / __mu - _Tp(2);
525 if (std::abs(__sn) < __errtol)
527 const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
529 __xn = __c0 * (__xn + __lambda);
530 __yn = __c0 * (__yn + __lambda);
533 _Tp __s = __sn * __sn
534 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
536 return (_Tp(1) + __s) / std::sqrt(__mu);
542 * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$
545 * The Carlson elliptic function of the third kind is defined by:
547 * R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty
548 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)}
551 * Based on Carlson's algorithms:
552 * - B. C. Carlson Numer. Math. 33, 1 (1979)
553 * - B. C. Carlson, Special Functions of Applied Mathematics (1977)
554 * - Numerical Recipes in C, 2nd ed, pp. 261-269,
555 * by Press, Teukolsky, Vetterling, Flannery (1992)
557 * @param __x The first of three symmetric arguments.
558 * @param __y The second of three symmetric arguments.
559 * @param __z The third of three symmetric arguments.
560 * @param __p The fourth argument.
561 * @return The Carlson elliptic function of the fourth kind.
563 template<typename _Tp>
565 __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
567 const _Tp __min = std::numeric_limits<_Tp>::min();
568 const _Tp __max = std::numeric_limits<_Tp>::max();
569 const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3));
570 const _Tp __uplim = _Tp(0.3L)
571 * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3));
573 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
574 std::__throw_domain_error(__N("Argument less than zero "
576 else if (__x + __y < __lolim || __x + __z < __lolim
577 || __y + __z < __lolim || __p < __lolim)
578 std::__throw_domain_error(__N("Argument too small "
582 const _Tp __c0 = _Tp(1) / _Tp(4);
583 const _Tp __c1 = _Tp(3) / _Tp(14);
584 const _Tp __c2 = _Tp(1) / _Tp(3);
585 const _Tp __c3 = _Tp(3) / _Tp(22);
586 const _Tp __c4 = _Tp(3) / _Tp(26);
592 _Tp __sigma = _Tp(0);
593 _Tp __power4 = _Tp(1);
595 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
596 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
599 _Tp __xndev, __yndev, __zndev, __pndev;
601 const unsigned int __max_iter = 100;
602 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
604 __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5);
605 __xndev = (__mu - __xn) / __mu;
606 __yndev = (__mu - __yn) / __mu;
607 __zndev = (__mu - __zn) / __mu;
608 __pndev = (__mu - __pn) / __mu;
609 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
610 __epsilon = std::max(__epsilon, std::abs(__zndev));
611 __epsilon = std::max(__epsilon, std::abs(__pndev));
612 if (__epsilon < __errtol)
614 const _Tp __xnroot = std::sqrt(__xn);
615 const _Tp __ynroot = std::sqrt(__yn);
616 const _Tp __znroot = std::sqrt(__zn);
617 const _Tp __lambda = __xnroot * (__ynroot + __znroot)
618 + __ynroot * __znroot;
619 const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot)
620 + __xnroot * __ynroot * __znroot;
621 const _Tp __alpha2 = __alpha1 * __alpha1;
622 const _Tp __beta = __pn * (__pn + __lambda)
624 __sigma += __power4 * __ellint_rc(__alpha2, __beta);
626 __xn = __c0 * (__xn + __lambda);
627 __yn = __c0 * (__yn + __lambda);
628 __zn = __c0 * (__zn + __lambda);
629 __pn = __c0 * (__pn + __lambda);
632 // Note: __ea is an SPU badname.
633 _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev;
634 _Tp __eb = __xndev * __yndev * __zndev;
635 _Tp __ec = __pndev * __pndev;
636 _Tp __e2 = __eaa - _Tp(3) * __ec;
637 _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec);
638 _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4)
639 - _Tp(3) * __c4 * __e3 / _Tp(2));
640 _Tp __s2 = __eb * (__c2 / _Tp(2)
641 + __pndev * (-__c3 - __c3 + __pndev * __c4));
642 _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3)
643 - __c2 * __pndev * __ec;
645 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
646 / (__mu * std::sqrt(__mu));
652 * @brief Return the complete elliptic integral of the third kind
653 * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the
654 * Carlson formulation.
656 * The complete elliptic integral of the third kind is defined as
658 * \Pi(k,\nu) = \int_0^{\pi/2}
660 * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
663 * @param __k The argument of the elliptic function.
664 * @param __nu The second argument of the elliptic function.
665 * @return The complete elliptic function of the third kind.
667 template<typename _Tp>
669 __comp_ellint_3(const _Tp __k, const _Tp __nu)
672 if (__isnan(__k) || __isnan(__nu))
673 return std::numeric_limits<_Tp>::quiet_NaN();
674 else if (__nu == _Tp(1))
675 return std::numeric_limits<_Tp>::infinity();
676 else if (std::abs(__k) > _Tp(1))
677 std::__throw_domain_error(__N("Bad argument in __comp_ellint_3."));
680 const _Tp __kk = __k * __k;
682 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
684 * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
691 * @brief Return the incomplete elliptic integral of the third kind
692 * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation.
694 * The incomplete elliptic integral of the third kind is defined as
696 * \Pi(k,\nu,\phi) = \int_0^{\phi}
698 * {(1 - \nu \sin^2\theta)
699 * \sqrt{1 - k^2 \sin^2\theta}}
702 * @param __k The argument of the elliptic function.
703 * @param __nu The second argument of the elliptic function.
704 * @param __phi The integral limit argument of the elliptic function.
705 * @return The elliptic function of the third kind.
707 template<typename _Tp>
709 __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
712 if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))
713 return std::numeric_limits<_Tp>::quiet_NaN();
714 else if (std::abs(__k) > _Tp(1))
715 std::__throw_domain_error(__N("Bad argument in __ellint_3."));
718 // Reduce phi to -pi/2 < phi < +pi/2.
719 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
721 const _Tp __phi_red = __phi
722 - __n * __numeric_constants<_Tp>::__pi();
724 const _Tp __kk = __k * __k;
725 const _Tp __s = std::sin(__phi_red);
726 const _Tp __ss = __s * __s;
727 const _Tp __sss = __ss * __s;
728 const _Tp __c = std::cos(__phi_red);
729 const _Tp __cc = __c * __c;
732 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
734 * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
735 _Tp(1) + __nu * __ss) / _Tp(3);
740 return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
744 _GLIBCXX_END_NAMESPACE_VERSION
745 } // namespace std::tr1::__detail
749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC