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