00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC
00044 #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1
00045
00046 namespace std
00047 {
00048 namespace tr1
00049 {
00050
00051
00052
00053
00054 namespace __detail
00055 {
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 template<typename _Tp>
00073 _Tp
00074 __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
00075 {
00076 const _Tp __min = std::numeric_limits<_Tp>::min();
00077 const _Tp __max = std::numeric_limits<_Tp>::max();
00078 const _Tp __lolim = _Tp(5) * __min;
00079 const _Tp __uplim = __max / _Tp(5);
00080
00081 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
00082 std::__throw_domain_error(__N("Argument less than zero "
00083 "in __ellint_rf."));
00084 else if (__x + __y < __lolim || __x + __z < __lolim
00085 || __y + __z < __lolim)
00086 std::__throw_domain_error(__N("Argument too small in __ellint_rf"));
00087 else
00088 {
00089 const _Tp __c0 = _Tp(1) / _Tp(4);
00090 const _Tp __c1 = _Tp(1) / _Tp(24);
00091 const _Tp __c2 = _Tp(1) / _Tp(10);
00092 const _Tp __c3 = _Tp(3) / _Tp(44);
00093 const _Tp __c4 = _Tp(1) / _Tp(14);
00094
00095 _Tp __xn = __x;
00096 _Tp __yn = __y;
00097 _Tp __zn = __z;
00098
00099 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00100 const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
00101 _Tp __mu;
00102 _Tp __xndev, __yndev, __zndev;
00103
00104 const unsigned int __max_iter = 100;
00105 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00106 {
00107 __mu = (__xn + __yn + __zn) / _Tp(3);
00108 __xndev = 2 - (__mu + __xn) / __mu;
00109 __yndev = 2 - (__mu + __yn) / __mu;
00110 __zndev = 2 - (__mu + __zn) / __mu;
00111 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00112 __epsilon = std::max(__epsilon, std::abs(__zndev));
00113 if (__epsilon < __errtol)
00114 break;
00115 const _Tp __xnroot = std::sqrt(__xn);
00116 const _Tp __ynroot = std::sqrt(__yn);
00117 const _Tp __znroot = std::sqrt(__zn);
00118 const _Tp __lambda = __xnroot * (__ynroot + __znroot)
00119 + __ynroot * __znroot;
00120 __xn = __c0 * (__xn + __lambda);
00121 __yn = __c0 * (__yn + __lambda);
00122 __zn = __c0 * (__zn + __lambda);
00123 }
00124
00125 const _Tp __e2 = __xndev * __yndev - __zndev * __zndev;
00126 const _Tp __e3 = __xndev * __yndev * __zndev;
00127 const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2
00128 + __c4 * __e3;
00129
00130 return __s / std::sqrt(__mu);
00131 }
00132 }
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 template<typename _Tp>
00152 _Tp
00153 __comp_ellint_1_series(const _Tp __k)
00154 {
00155
00156 const _Tp __kk = __k * __k;
00157
00158 _Tp __term = __kk / _Tp(4);
00159 _Tp __sum = _Tp(1) + __term;
00160
00161 const unsigned int __max_iter = 1000;
00162 for (unsigned int __i = 2; __i < __max_iter; ++__i)
00163 {
00164 __term *= (2 * __i - 1) * __kk / (2 * __i);
00165 if (__term < std::numeric_limits<_Tp>::epsilon())
00166 break;
00167 __sum += __term;
00168 }
00169
00170 return __numeric_constants<_Tp>::__pi_2() * __sum;
00171 }
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 template<typename _Tp>
00190 _Tp
00191 __comp_ellint_1(const _Tp __k)
00192 {
00193
00194 if (__isnan(__k))
00195 return std::numeric_limits<_Tp>::quiet_NaN();
00196 else if (std::abs(__k) >= _Tp(1))
00197 return std::numeric_limits<_Tp>::quiet_NaN();
00198 else
00199 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
00200 }
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 template<typename _Tp>
00218 _Tp
00219 __ellint_1(const _Tp __k, const _Tp __phi)
00220 {
00221
00222 if (__isnan(__k) || __isnan(__phi))
00223 return std::numeric_limits<_Tp>::quiet_NaN();
00224 else if (std::abs(__k) > _Tp(1))
00225 std::__throw_domain_error(__N("Bad argument in __ellint_1."));
00226 else
00227 {
00228
00229 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00230 + _Tp(0.5L));
00231 const _Tp __phi_red = __phi
00232 - __n * __numeric_constants<_Tp>::__pi();
00233
00234 const _Tp __s = std::sin(__phi_red);
00235 const _Tp __c = std::cos(__phi_red);
00236
00237 const _Tp __F = __s
00238 * __ellint_rf(__c * __c,
00239 _Tp(1) - __k * __k * __s * __s, _Tp(1));
00240
00241 if (__n == 0)
00242 return __F;
00243 else
00244 return __F + _Tp(2) * __n * __comp_ellint_1(__k);
00245 }
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264 template<typename _Tp>
00265 _Tp
00266 __comp_ellint_2_series(const _Tp __k)
00267 {
00268
00269 const _Tp __kk = __k * __k;
00270
00271 _Tp __term = __kk;
00272 _Tp __sum = __term;
00273
00274 const unsigned int __max_iter = 1000;
00275 for (unsigned int __i = 2; __i < __max_iter; ++__i)
00276 {
00277 const _Tp __i2m = 2 * __i - 1;
00278 const _Tp __i2 = 2 * __i;
00279 __term *= __i2m * __i2m * __kk / (__i2 * __i2);
00280 if (__term < std::numeric_limits<_Tp>::epsilon())
00281 break;
00282 __sum += __term / __i2m;
00283 }
00284
00285 return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312 template<typename _Tp>
00313 _Tp
00314 __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
00315 {
00316 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00317 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
00318 const _Tp __min = std::numeric_limits<_Tp>::min();
00319 const _Tp __max = std::numeric_limits<_Tp>::max();
00320 const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3));
00321 const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3));
00322
00323 if (__x < _Tp(0) || __y < _Tp(0))
00324 std::__throw_domain_error(__N("Argument less than zero "
00325 "in __ellint_rd."));
00326 else if (__x + __y < __lolim || __z < __lolim)
00327 std::__throw_domain_error(__N("Argument too small "
00328 "in __ellint_rd."));
00329 else
00330 {
00331 const _Tp __c0 = _Tp(1) / _Tp(4);
00332 const _Tp __c1 = _Tp(3) / _Tp(14);
00333 const _Tp __c2 = _Tp(1) / _Tp(6);
00334 const _Tp __c3 = _Tp(9) / _Tp(22);
00335 const _Tp __c4 = _Tp(3) / _Tp(26);
00336
00337 _Tp __xn = __x;
00338 _Tp __yn = __y;
00339 _Tp __zn = __z;
00340 _Tp __sigma = _Tp(0);
00341 _Tp __power4 = _Tp(1);
00342
00343 _Tp __mu;
00344 _Tp __xndev, __yndev, __zndev;
00345
00346 const unsigned int __max_iter = 100;
00347 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00348 {
00349 __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5);
00350 __xndev = (__mu - __xn) / __mu;
00351 __yndev = (__mu - __yn) / __mu;
00352 __zndev = (__mu - __zn) / __mu;
00353 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00354 __epsilon = std::max(__epsilon, std::abs(__zndev));
00355 if (__epsilon < __errtol)
00356 break;
00357 _Tp __xnroot = std::sqrt(__xn);
00358 _Tp __ynroot = std::sqrt(__yn);
00359 _Tp __znroot = std::sqrt(__zn);
00360 _Tp __lambda = __xnroot * (__ynroot + __znroot)
00361 + __ynroot * __znroot;
00362 __sigma += __power4 / (__znroot * (__zn + __lambda));
00363 __power4 *= __c0;
00364 __xn = __c0 * (__xn + __lambda);
00365 __yn = __c0 * (__yn + __lambda);
00366 __zn = __c0 * (__zn + __lambda);
00367 }
00368
00369
00370 _Tp __eaa = __xndev * __yndev;
00371 _Tp __eb = __zndev * __zndev;
00372 _Tp __ec = __eaa - __eb;
00373 _Tp __ed = __eaa - _Tp(6) * __eb;
00374 _Tp __ef = __ed + __ec + __ec;
00375 _Tp __s1 = __ed * (-__c1 + __c3 * __ed
00376 / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef
00377 / _Tp(2));
00378 _Tp __s2 = __zndev
00379 * (__c2 * __ef
00380 + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
00381
00382 return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
00383 / (__mu * std::sqrt(__mu));
00384 }
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 template<typename _Tp>
00401 _Tp
00402 __comp_ellint_2(const _Tp __k)
00403 {
00404
00405 if (__isnan(__k))
00406 return std::numeric_limits<_Tp>::quiet_NaN();
00407 else if (std::abs(__k) == 1)
00408 return _Tp(1);
00409 else if (std::abs(__k) > _Tp(1))
00410 std::__throw_domain_error(__N("Bad argument in __comp_ellint_2."));
00411 else
00412 {
00413 const _Tp __kk = __k * __k;
00414
00415 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
00416 - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
00417 }
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434 template<typename _Tp>
00435 _Tp
00436 __ellint_2(const _Tp __k, const _Tp __phi)
00437 {
00438
00439 if (__isnan(__k) || __isnan(__phi))
00440 return std::numeric_limits<_Tp>::quiet_NaN();
00441 else if (std::abs(__k) > _Tp(1))
00442 std::__throw_domain_error(__N("Bad argument in __ellint_2."));
00443 else
00444 {
00445
00446 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00447 + _Tp(0.5L));
00448 const _Tp __phi_red = __phi
00449 - __n * __numeric_constants<_Tp>::__pi();
00450
00451 const _Tp __kk = __k * __k;
00452 const _Tp __s = std::sin(__phi_red);
00453 const _Tp __ss = __s * __s;
00454 const _Tp __sss = __ss * __s;
00455 const _Tp __c = std::cos(__phi_red);
00456 const _Tp __cc = __c * __c;
00457
00458 const _Tp __E = __s
00459 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00460 - __kk * __sss
00461 * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00462 / _Tp(3);
00463
00464 if (__n == 0)
00465 return __E;
00466 else
00467 return __E + _Tp(2) * __n * __comp_ellint_2(__k);
00468 }
00469 }
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 template<typename _Tp>
00494 _Tp
00495 __ellint_rc(const _Tp __x, const _Tp __y)
00496 {
00497 const _Tp __min = std::numeric_limits<_Tp>::min();
00498 const _Tp __max = std::numeric_limits<_Tp>::max();
00499 const _Tp __lolim = _Tp(5) * __min;
00500 const _Tp __uplim = __max / _Tp(5);
00501
00502 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
00503 std::__throw_domain_error(__N("Argument less than zero "
00504 "in __ellint_rc."));
00505 else
00506 {
00507 const _Tp __c0 = _Tp(1) / _Tp(4);
00508 const _Tp __c1 = _Tp(1) / _Tp(7);
00509 const _Tp __c2 = _Tp(9) / _Tp(22);
00510 const _Tp __c3 = _Tp(3) / _Tp(10);
00511 const _Tp __c4 = _Tp(3) / _Tp(8);
00512
00513 _Tp __xn = __x;
00514 _Tp __yn = __y;
00515
00516 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00517 const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
00518 _Tp __mu;
00519 _Tp __sn;
00520
00521 const unsigned int __max_iter = 100;
00522 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00523 {
00524 __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
00525 __sn = (__yn + __mu) / __mu - _Tp(2);
00526 if (std::abs(__sn) < __errtol)
00527 break;
00528 const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
00529 + __yn;
00530 __xn = __c0 * (__xn + __lambda);
00531 __yn = __c0 * (__yn + __lambda);
00532 }
00533
00534 _Tp __s = __sn * __sn
00535 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
00536
00537 return (_Tp(1) + __s) / std::sqrt(__mu);
00538 }
00539 }
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564 template<typename _Tp>
00565 _Tp
00566 __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
00567 {
00568 const _Tp __min = std::numeric_limits<_Tp>::min();
00569 const _Tp __max = std::numeric_limits<_Tp>::max();
00570 const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3));
00571 const _Tp __uplim = _Tp(0.3L)
00572 * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3));
00573
00574 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
00575 std::__throw_domain_error(__N("Argument less than zero "
00576 "in __ellint_rj."));
00577 else if (__x + __y < __lolim || __x + __z < __lolim
00578 || __y + __z < __lolim || __p < __lolim)
00579 std::__throw_domain_error(__N("Argument too small "
00580 "in __ellint_rj"));
00581 else
00582 {
00583 const _Tp __c0 = _Tp(1) / _Tp(4);
00584 const _Tp __c1 = _Tp(3) / _Tp(14);
00585 const _Tp __c2 = _Tp(1) / _Tp(3);
00586 const _Tp __c3 = _Tp(3) / _Tp(22);
00587 const _Tp __c4 = _Tp(3) / _Tp(26);
00588
00589 _Tp __xn = __x;
00590 _Tp __yn = __y;
00591 _Tp __zn = __z;
00592 _Tp __pn = __p;
00593 _Tp __sigma = _Tp(0);
00594 _Tp __power4 = _Tp(1);
00595
00596 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00597 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
00598
00599 _Tp __lambda, __mu;
00600 _Tp __xndev, __yndev, __zndev, __pndev;
00601
00602 const unsigned int __max_iter = 100;
00603 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00604 {
00605 __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5);
00606 __xndev = (__mu - __xn) / __mu;
00607 __yndev = (__mu - __yn) / __mu;
00608 __zndev = (__mu - __zn) / __mu;
00609 __pndev = (__mu - __pn) / __mu;
00610 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00611 __epsilon = std::max(__epsilon, std::abs(__zndev));
00612 __epsilon = std::max(__epsilon, std::abs(__pndev));
00613 if (__epsilon < __errtol)
00614 break;
00615 const _Tp __xnroot = std::sqrt(__xn);
00616 const _Tp __ynroot = std::sqrt(__yn);
00617 const _Tp __znroot = std::sqrt(__zn);
00618 const _Tp __lambda = __xnroot * (__ynroot + __znroot)
00619 + __ynroot * __znroot;
00620 const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot)
00621 + __xnroot * __ynroot * __znroot;
00622 const _Tp __alpha2 = __alpha1 * __alpha1;
00623 const _Tp __beta = __pn * (__pn + __lambda)
00624 * (__pn + __lambda);
00625 __sigma += __power4 * __ellint_rc(__alpha2, __beta);
00626 __power4 *= __c0;
00627 __xn = __c0 * (__xn + __lambda);
00628 __yn = __c0 * (__yn + __lambda);
00629 __zn = __c0 * (__zn + __lambda);
00630 __pn = __c0 * (__pn + __lambda);
00631 }
00632
00633
00634 _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev;
00635 _Tp __eb = __xndev * __yndev * __zndev;
00636 _Tp __ec = __pndev * __pndev;
00637 _Tp __e2 = __eaa - _Tp(3) * __ec;
00638 _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec);
00639 _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4)
00640 - _Tp(3) * __c4 * __e3 / _Tp(2));
00641 _Tp __s2 = __eb * (__c2 / _Tp(2)
00642 + __pndev * (-__c3 - __c3 + __pndev * __c4));
00643 _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3)
00644 - __c2 * __pndev * __ec;
00645
00646 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
00647 / (__mu * std::sqrt(__mu));
00648 }
00649 }
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668 template<typename _Tp>
00669 _Tp
00670 __comp_ellint_3(const _Tp __k, const _Tp __nu)
00671 {
00672
00673 if (__isnan(__k) || __isnan(__nu))
00674 return std::numeric_limits<_Tp>::quiet_NaN();
00675 else if (__nu == _Tp(1))
00676 return std::numeric_limits<_Tp>::infinity();
00677 else if (std::abs(__k) > _Tp(1))
00678 std::__throw_domain_error(__N("Bad argument in __comp_ellint_3."));
00679 else
00680 {
00681 const _Tp __kk = __k * __k;
00682
00683 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
00684 - __nu
00685 * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
00686 / _Tp(3);
00687 }
00688 }
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708 template<typename _Tp>
00709 _Tp
00710 __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
00711 {
00712
00713 if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))
00714 return std::numeric_limits<_Tp>::quiet_NaN();
00715 else if (std::abs(__k) > _Tp(1))
00716 std::__throw_domain_error(__N("Bad argument in __ellint_3."));
00717 else
00718 {
00719
00720 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00721 + _Tp(0.5L));
00722 const _Tp __phi_red = __phi
00723 - __n * __numeric_constants<_Tp>::__pi();
00724
00725 const _Tp __kk = __k * __k;
00726 const _Tp __s = std::sin(__phi_red);
00727 const _Tp __ss = __s * __s;
00728 const _Tp __sss = __ss * __s;
00729 const _Tp __c = std::cos(__phi_red);
00730 const _Tp __cc = __c * __c;
00731
00732 const _Tp __Pi = __s
00733 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00734 - __nu * __sss
00735 * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
00736 _Tp(1) + __nu * __ss) / _Tp(3);
00737
00738 if (__n == 0)
00739 return __Pi;
00740 else
00741 return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
00742 }
00743 }
00744
00745 }
00746 }
00747 }
00748
00749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC
00750