STLdoc
STLdocumentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Macros | Functions
ell_integral.tcc File Reference

Macros

#define _GLIBCXX_TR1_ELL_INTEGRAL_TCC   1
 

Functions

namespace std _GLIBCXX_VISIBILITY (default)
 

Detailed Description

This is an internal header file, included by other library headers. Do not attempt to use it directly. {tr1/cmath}

Macro Definition Documentation

#define _GLIBCXX_TR1_ELL_INTEGRAL_TCC   1

Function Documentation

namespace std _GLIBCXX_VISIBILITY ( default  )

Return the Carlson elliptic function $ R_F(x,y,z) $ of the first kind.

The Carlson elliptic function of the first kind is defined by:

\[ R_F(x,y,z) = \frac{1}{2} \int_0^\infty \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}} \]

Parameters
__xThe first of three symmetric arguments.
__yThe second of three symmetric arguments.
__zThe third of three symmetric arguments.
Returns
The Carlson elliptic function of the first kind.

Return the complete elliptic integral of the first kind $ K(k) $ by series expansion.

The complete elliptic integral of the first kind is defined as

\[ K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} {\sqrt{1 - k^2sin^2\theta}} \]

This routine is not bad as long as |k| is somewhat smaller than 1 but is not is good as the Carlson elliptic integral formulation.

Parameters
__kThe argument of the complete elliptic function.
Returns
The complete elliptic function of the first kind.

Return the complete elliptic integral of the first kind $ K(k) $ using the Carlson formulation.

The complete elliptic integral of the first kind is defined as

\[ K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} {\sqrt{1 - k^2 sin^2\theta}} \]

where $ F(k,\phi) $ is the incomplete elliptic integral of the first kind.

Parameters
__kThe argument of the complete elliptic function.
Returns
The complete elliptic function of the first kind.

Return the incomplete elliptic integral of the first kind $ F(k,\phi) $ using the Carlson formulation.

The incomplete elliptic integral of the first kind is defined as

\[ F(k,\phi) = \int_0^{\phi}\frac{d\theta} {\sqrt{1 - k^2 sin^2\theta}} \]

Parameters
__kThe argument of the elliptic function.
__phiThe integral limit argument of the elliptic function.
Returns
The elliptic function of the first kind.

Return the complete elliptic integral of the second kind $ E(k) $ by series expansion.

The complete elliptic integral of the second kind is defined as

\[ E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} \]

This routine is not bad as long as |k| is somewhat smaller than 1 but is not is good as the Carlson elliptic integral formulation.

Parameters
__kThe argument of the complete elliptic function.
Returns
The complete elliptic function of the second kind.

Return the Carlson elliptic function of the second kind $ R_D(x,y,z) = R_J(x,y,z,z) $ where $ R_J(x,y,z,p) $ is the Carlson elliptic function of the third kind.

The Carlson elliptic function of the second kind is defined by:

\[ R_D(x,y,z) = \frac{3}{2} \int_0^\infty \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}} \]

Based on Carlson's algorithms:

  • B. C. Carlson Numer. Math. 33, 1 (1979)
  • B. C. Carlson, Special Functions of Applied Mathematics (1977)
  • Numerical Recipes in C, 2nd ed, pp. 261-269, by Press, Teukolsky, Vetterling, Flannery (1992)
Parameters
__xThe first of two symmetric arguments.
__yThe second of two symmetric arguments.
__zThe third argument.
Returns
The Carlson elliptic function of the second kind.

Return the complete elliptic integral of the second kind $ E(k) $ using the Carlson formulation.

The complete elliptic integral of the second kind is defined as

\[ E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} \]

Parameters
__kThe argument of the complete elliptic function.
Returns
The complete elliptic function of the second kind.

Return the incomplete elliptic integral of the second kind $ E(k,\phi) $ using the Carlson formulation.

The incomplete elliptic integral of the second kind is defined as

\[ E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta} \]

Parameters
__kThe argument of the elliptic function.
__phiThe integral limit argument of the elliptic function.
Returns
The elliptic function of the second kind.

Return the Carlson elliptic function $ R_C(x,y) = R_F(x,y,y) $ where $ R_F(x,y,z) $ is the Carlson elliptic function of the first kind.

The Carlson elliptic function is defined by:

\[ R_C(x,y) = \frac{1}{2} \int_0^\infty \frac{dt}{(t + x)^{1/2}(t + y)} \]

Based on Carlson's algorithms:

  • B. C. Carlson Numer. Math. 33, 1 (1979)
  • B. C. Carlson, Special Functions of Applied Mathematics (1977)
  • Numerical Recipes in C, 2nd ed, pp. 261-269, by Press, Teukolsky, Vetterling, Flannery (1992)
Parameters
__xThe first argument.
__yThe second argument.
Returns
The Carlson elliptic function.

Return the Carlson elliptic function $ R_J(x,y,z,p) $ of the third kind.

The Carlson elliptic function of the third kind is defined by:

\[ R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)} \]

Based on Carlson's algorithms:

  • B. C. Carlson Numer. Math. 33, 1 (1979)
  • B. C. Carlson, Special Functions of Applied Mathematics (1977)
  • Numerical Recipes in C, 2nd ed, pp. 261-269, by Press, Teukolsky, Vetterling, Flannery (1992)
Parameters
__xThe first of three symmetric arguments.
__yThe second of three symmetric arguments.
__zThe third of three symmetric arguments.
__pThe fourth argument.
Returns
The Carlson elliptic function of the fourth kind.

Return the complete elliptic integral of the third kind $ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) $ using the Carlson formulation.

The complete elliptic integral of the third kind is defined as

\[ \Pi(k,\nu) = \int_0^{\pi/2} \frac{d\theta} {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}} \]

Parameters
__kThe argument of the elliptic function.
__nuThe second argument of the elliptic function.
Returns
The complete elliptic function of the third kind.

Return the incomplete elliptic integral of the third kind $ \Pi(k,\nu,\phi) $ using the Carlson formulation.

The incomplete elliptic integral of the third kind is defined as

\[ \Pi(k,\nu,\phi) = \int_0^{\phi} \frac{d\theta} {(1 - \nu \sin^2\theta) \sqrt{1 - k^2 \sin^2\theta}} \]

Parameters
__kThe argument of the elliptic function.
__nuThe second argument of the elliptic function.
__phiThe integral limit argument of the elliptic function.
Returns
The elliptic function of the third kind.
46 {
47 namespace tr1
48 {
49  // [5.2] Special functions
50 
51  // Implementation-space details.
52  namespace __detail
53  {
54  _GLIBCXX_BEGIN_NAMESPACE_VERSION
55 
71  template<typename _Tp>
72  _Tp
73  __ellint_rf(_Tp __x, _Tp __y, _Tp __z)
74  {
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);
79 
80  if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
81  std::__throw_domain_error(__N("Argument less than zero "
82  "in __ellint_rf."));
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"));
86  else
87  {
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);
93 
94  _Tp __xn = __x;
95  _Tp __yn = __y;
96  _Tp __zn = __z;
97 
98  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
99  const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
100  _Tp __mu;
101  _Tp __xndev, __yndev, __zndev;
102 
103  const unsigned int __max_iter = 100;
104  for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
105  {
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)
113  break;
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);
122  }
123 
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
127  + __c4 * __e3;
128 
129  return __s / std::sqrt(__mu);
130  }
131  }
132 
133 
150  template<typename _Tp>
151  _Tp
152  __comp_ellint_1_series(_Tp __k)
153  {
154 
155  const _Tp __kk = __k * __k;
156 
157  _Tp __term = __kk / _Tp(4);
158  _Tp __sum = _Tp(1) + __term;
159 
160  const unsigned int __max_iter = 1000;
161  for (unsigned int __i = 2; __i < __max_iter; ++__i)
162  {
163  __term *= (2 * __i - 1) * __kk / (2 * __i);
164  if (__term < std::numeric_limits<_Tp>::epsilon())
165  break;
166  __sum += __term;
167  }
168 
169  return __numeric_constants<_Tp>::__pi_2() * __sum;
170  }
171 
172 
188  template<typename _Tp>
189  _Tp
190  __comp_ellint_1(_Tp __k)
191  {
192 
193  if (__isnan(__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();
197  else
198  return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
199  }
200 
201 
216  template<typename _Tp>
217  _Tp
218  __ellint_1(_Tp __k, _Tp __phi)
219  {
220 
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."));
225  else
226  {
227  // Reduce phi to -pi/2 < phi < +pi/2.
228  const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
229  + _Tp(0.5L));
230  const _Tp __phi_red = __phi
231  - __n * __numeric_constants<_Tp>::__pi();
232 
233  const _Tp __s = std::sin(__phi_red);
234  const _Tp __c = std::cos(__phi_red);
235 
236  const _Tp __F = __s
237  * __ellint_rf(__c * __c,
238  _Tp(1) - __k * __k * __s * __s, _Tp(1));
239 
240  if (__n == 0)
241  return __F;
242  else
243  return __F + _Tp(2) * __n * __comp_ellint_1(__k);
244  }
245  }
246 
247 
263  template<typename _Tp>
264  _Tp
265  __comp_ellint_2_series(_Tp __k)
266  {
267 
268  const _Tp __kk = __k * __k;
269 
270  _Tp __term = __kk;
271  _Tp __sum = __term;
272 
273  const unsigned int __max_iter = 1000;
274  for (unsigned int __i = 2; __i < __max_iter; ++__i)
275  {
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())
280  break;
281  __sum += __term / __i2m;
282  }
283 
284  return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
285  }
286 
287 
311  template<typename _Tp>
312  _Tp
313  __ellint_rd(_Tp __x, _Tp __y, _Tp __z)
314  {
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));
321 
322  if (__x < _Tp(0) || __y < _Tp(0))
323  std::__throw_domain_error(__N("Argument less than zero "
324  "in __ellint_rd."));
325  else if (__x + __y < __lolim || __z < __lolim)
326  std::__throw_domain_error(__N("Argument too small "
327  "in __ellint_rd."));
328  else
329  {
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);
335 
336  _Tp __xn = __x;
337  _Tp __yn = __y;
338  _Tp __zn = __z;
339  _Tp __sigma = _Tp(0);
340  _Tp __power4 = _Tp(1);
341 
342  _Tp __mu;
343  _Tp __xndev, __yndev, __zndev;
344 
345  const unsigned int __max_iter = 100;
346  for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
347  {
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)
355  break;
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));
362  __power4 *= __c0;
363  __xn = __c0 * (__xn + __lambda);
364  __yn = __c0 * (__yn + __lambda);
365  __zn = __c0 * (__zn + __lambda);
366  }
367 
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
376  / _Tp(2));
377  _Tp __s2 = __zndev
378  * (__c2 * __ef
379  + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
380 
381  return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
382  / (__mu * std::sqrt(__mu));
383  }
384  }
385 
386 
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 
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 
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 __max = std::numeric_limits<_Tp>::max();
498  const _Tp __lolim = _Tp(5) * __min;
499  const _Tp __uplim = __max / _Tp(5);
500 
501  if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
502  std::__throw_domain_error(__N("Argument less than zero "
503  "in __ellint_rc."));
504  else
505  {
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);
511 
512  _Tp __xn = __x;
513  _Tp __yn = __y;
514 
515  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
516  const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
517  _Tp __mu;
518  _Tp __sn;
519 
520  const unsigned int __max_iter = 100;
521  for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
522  {
523  __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
524  __sn = (__yn + __mu) / __mu - _Tp(2);
525  if (std::abs(__sn) < __errtol)
526  break;
527  const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
528  + __yn;
529  __xn = __c0 * (__xn + __lambda);
530  __yn = __c0 * (__yn + __lambda);
531  }
532 
533  _Tp __s = __sn * __sn
534  * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
535 
536  return (_Tp(1) + __s) / std::sqrt(__mu);
537  }
538  }
539 
540 
563  template<typename _Tp>
564  _Tp
565  __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p)
566  {
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));
572 
573  if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
574  std::__throw_domain_error(__N("Argument less than zero "
575  "in __ellint_rj."));
576  else if (__x + __y < __lolim || __x + __z < __lolim
577  || __y + __z < __lolim || __p < __lolim)
578  std::__throw_domain_error(__N("Argument too small "
579  "in __ellint_rj"));
580  else
581  {
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);
587 
588  _Tp __xn = __x;
589  _Tp __yn = __y;
590  _Tp __zn = __z;
591  _Tp __pn = __p;
592  _Tp __sigma = _Tp(0);
593  _Tp __power4 = _Tp(1);
594 
595  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
596  const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
597 
598  _Tp __lambda, __mu;
599  _Tp __xndev, __yndev, __zndev, __pndev;
600 
601  const unsigned int __max_iter = 100;
602  for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
603  {
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)
613  break;
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)
623  * (__pn + __lambda);
624  __sigma += __power4 * __ellint_rc(__alpha2, __beta);
625  __power4 *= __c0;
626  __xn = __c0 * (__xn + __lambda);
627  __yn = __c0 * (__yn + __lambda);
628  __zn = __c0 * (__zn + __lambda);
629  __pn = __c0 * (__pn + __lambda);
630  }
631 
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;
644 
645  return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
646  / (__mu * std::sqrt(__mu));
647  }
648  }
649 
650 
667  template<typename _Tp>
668  _Tp
669  __comp_ellint_3(_Tp __k, _Tp __nu)
670  {
671 
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."));
678  else
679  {
680  const _Tp __kk = __k * __k;
681 
682  return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
683  - __nu
684  * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
685  / _Tp(3);
686  }
687  }
688 
689 
707  template<typename _Tp>
708  _Tp
709  __ellint_3(_Tp __k, _Tp __nu, _Tp __phi)
710  {
711 
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."));
716  else
717  {
718  // Reduce phi to -pi/2 < phi < +pi/2.
719  const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
720  + _Tp(0.5L));
721  const _Tp __phi_red = __phi
722  - __n * __numeric_constants<_Tp>::__pi();
723 
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;
730 
731  const _Tp __Pi = __s
732  * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
733  - __nu * __sss
734  * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
735  _Tp(1) + __nu * __ss) / _Tp(3);
736 
737  if (__n == 0)
738  return __Pi;
739  else
740  return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
741  }
742  }
743 
744  _GLIBCXX_END_NAMESPACE_VERSION
745  } // namespace std::tr1::__detail
746 }
747 }
__gnu_parallel::min
const _Tp & min(const _Tp &__a, const _Tp &__b)
Equivalent to std::min.
Definition: base.h:144
__F
__inline __m256 float float float float float __F
Definition: avxintrin.h:1189
__gnu_parallel::max
const _Tp & max(const _Tp &__a, const _Tp &__b)
Equivalent to std::max.
Definition: base.h:150
__E
__inline __m256 float float float float __E
Definition: avxintrin.h:1189
Generated on Sat Jul 22 2017 14:42:36 for STLdoc by   doxygen 1.8.6