43 #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC
44 #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1
72 template<
typename _Tp>
78 const _Tp __lolim = _Tp(5) * __min;
79 const _Tp __uplim = __max / _Tp(5);
81 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
82 std::__throw_domain_error(__N(
"Argument less than zero "
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"));
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);
100 const _Tp __errtol =
std::pow(__eps, _Tp(1) / _Tp(6));
102 _Tp __xndev, __yndev, __zndev;
104 const unsigned int __max_iter = 100;
105 for (
unsigned int __iter = 0; __iter < __max_iter; ++__iter)
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;
113 if (__epsilon < __errtol)
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);
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
151 template<
typename _Tp>
156 const _Tp __kk = __k * __k;
158 _Tp __term = __kk / _Tp(4);
159 _Tp __sum = _Tp(1) + __term;
161 const unsigned int __max_iter = 1000;
162 for (
unsigned int __i = 2; __i < __max_iter; ++__i)
164 __term *= (2 * __i - 1) * __kk / (2 * __i);
189 template<
typename _Tp>
199 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
217 template<
typename _Tp>
222 if (__isnan(__k) || __isnan(__phi))
225 std::__throw_domain_error(__N(
"Bad argument in __ellint_1."));
231 const _Tp __phi_red = __phi
234 const _Tp __s =
std::sin(__phi_red);
235 const _Tp __c =
std::cos(__phi_red);
239 _Tp(1) - __k * __k * __s * __s, _Tp(1));
264 template<
typename _Tp>
269 const _Tp __kk = __k * __k;
274 const unsigned int __max_iter = 1000;
275 for (
unsigned int __i = 2; __i < __max_iter; ++__i)
277 const _Tp __i2m = 2 * __i - 1;
278 const _Tp __i2 = 2 * __i;
279 __term *= __i2m * __i2m * __kk / (__i2 * __i2);
282 __sum += __term / __i2m;
312 template<
typename _Tp>
317 const _Tp __errtol =
std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
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));
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;
355 if (__epsilon < __errtol)
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);
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
380 + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
382 return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
400 template<
typename _Tp>
410 std::__throw_domain_error(__N(
"Bad argument in __comp_ellint_2."));
413 const _Tp __kk = __k * __k;
416 - __kk *
__ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
434 template<
typename _Tp>
439 if (__isnan(__k) || __isnan(__phi))
442 std::__throw_domain_error(__N(
"Bad argument in __ellint_2."));
448 const _Tp __phi_red = __phi
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;
493 template<
typename _Tp>
499 const _Tp __lolim = _Tp(5) * __min;
500 const _Tp __uplim = __max / _Tp(5);
502 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
503 std::__throw_domain_error(__N(
"Argument less than zero "
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);
517 const _Tp __errtol =
std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
521 const unsigned int __max_iter = 100;
522 for (
unsigned int __iter = 0; __iter < __max_iter; ++__iter)
524 __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
525 __sn = (__yn + __mu) / __mu - _Tp(2);
530 __xn = __c0 * (__xn + __lambda);
531 __yn = __c0 * (__yn + __lambda);
534 _Tp __s = __sn * __sn
535 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
564 template<
typename _Tp>
566 __ellint_rj(
const _Tp __x,
const _Tp __y,
const _Tp __z,
const _Tp __p)
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));
574 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
575 std::__throw_domain_error(__N(
"Argument less than zero "
577 else if (__x + __y < __lolim || __x + __z < __lolim
578 || __y + __z < __lolim || __p < __lolim)
579 std::__throw_domain_error(__N(
"Argument too small "
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);
593 _Tp __sigma = _Tp(0);
594 _Tp __power4 = _Tp(1);
597 const _Tp __errtol =
std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
600 _Tp __xndev, __yndev, __zndev, __pndev;
602 const unsigned int __max_iter = 100;
603 for (
unsigned int __iter = 0; __iter < __max_iter; ++__iter)
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;
613 if (__epsilon < __errtol)
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)
625 __sigma += __power4 *
__ellint_rc(__alpha2, __beta);
627 __xn = __c0 * (__xn + __lambda);
628 __yn = __c0 * (__yn + __lambda);
629 __zn = __c0 * (__zn + __lambda);
630 __pn = __c0 * (__pn + __lambda);
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;
646 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
668 template<
typename _Tp>
673 if (__isnan(__k) || __isnan(__nu))
675 else if (__nu == _Tp(1))
678 std::__throw_domain_error(__N(
"Bad argument in __comp_ellint_3."));
681 const _Tp __kk = __k * __k;
685 *
__ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
708 template<
typename _Tp>
713 if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))
716 std::__throw_domain_error(__N(
"Bad argument in __ellint_3."));
722 const _Tp __phi_red = __phi
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;
736 _Tp(1) + __nu * __ss) / _Tp(3);
749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC
_Tp __ellint_2(const _Tp __k, const _Tp __phi)
Return the incomplete elliptic integral of the second kind using the Carlson formulation.
A structure for numeric constants.
static _Tp __pi_2()
Constant .
complex< _Tp > cos(const complex< _Tp > &)
Return complex cosine of z.
complex< _Tp > sqrt(const complex< _Tp > &)
Return complex square root of z.
_Tp __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
Return the incomplete elliptic integral of the third kind using the Carlson formulation.
_Tp __comp_ellint_3(const _Tp __k, const _Tp __nu)
Return the complete elliptic integral of the third kind using the Carlson formulation.
const _Tp & max(const _Tp &, const _Tp &)
This does what you think it does.
_Tp __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
Return the Carlson elliptic function of the first kind.
Properties of fundamental types.
_Tp __comp_ellint_1(const _Tp __k)
Return the complete elliptic integral of the first kind using the Carlson formulation.
_Tp __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
Return the Carlson elliptic function of the third kind.
_Tp __comp_ellint_1_series(const _Tp __k)
Return the complete elliptic integral of the first kind by series expansion.
_Tp __ellint_rc(const _Tp __x, const _Tp __y)
Return the Carlson elliptic function where is the Carlson elliptic function of the first kind...
_Tp __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
Return the Carlson elliptic function of the second kind where is the Carlson elliptic function of t...
_Tp __ellint_1(const _Tp __k, const _Tp __phi)
Return the incomplete elliptic integral of the first kind using the Carlson formulation.
complex< _Tp > sin(const complex< _Tp > &)
Return complex sine of z.
complex< _Tp > pow(const complex< _Tp > &, const _Tp &)
Return x to the y'th power.
_Tp __comp_ellint_2(const _Tp __k)
Return the complete elliptic integral of the second kind using the Carlson formulation.
_Tp __comp_ellint_2_series(const _Tp __k)
Return the complete elliptic integral of the second kind by series expansion.
static _Tp __pi()
Constant .
_Tp __beta(_Tp __x, _Tp __y)
Return the beta function .
_Tp abs(const complex< _Tp > &)
Return magnitude of z.