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
00044
00045
00046
00047
00048
00049
00050
00051
00052 #ifndef _GLIBCXX_TR1_BESSEL_FUNCTION_TCC
00053 #define _GLIBCXX_TR1_BESSEL_FUNCTION_TCC 1
00054
00055 #include "special_function_util.h"
00056
00057 namespace std
00058 {
00059 namespace tr1
00060 {
00061
00062
00063
00064
00065 namespace __detail
00066 {
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 template <typename _Tp>
00094 void
00095 __gamma_temme(const _Tp __mu,
00096 _Tp & __gam1, _Tp & __gam2, _Tp & __gampl, _Tp & __gammi)
00097 {
00098 #if _GLIBCXX_USE_C99_MATH_TR1
00099 __gampl = _Tp(1) / std::tr1::tgamma(_Tp(1) + __mu);
00100 __gammi = _Tp(1) / std::tr1::tgamma(_Tp(1) - __mu);
00101 #else
00102 __gampl = _Tp(1) / __gamma(_Tp(1) + __mu);
00103 __gammi = _Tp(1) / __gamma(_Tp(1) - __mu);
00104 #endif
00105
00106 if (std::abs(__mu) < std::numeric_limits<_Tp>::epsilon())
00107 __gam1 = -_Tp(__numeric_constants<_Tp>::__gamma_e());
00108 else
00109 __gam1 = (__gammi - __gampl) / (_Tp(2) * __mu);
00110
00111 __gam2 = (__gammi + __gampl) / (_Tp(2));
00112
00113 return;
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 template <typename _Tp>
00132 void
00133 __bessel_jn(const _Tp __nu, const _Tp __x,
00134 _Tp & __Jnu, _Tp & __Nnu, _Tp & __Jpnu, _Tp & __Npnu)
00135 {
00136 if (__x == _Tp(0))
00137 {
00138 if (__nu == _Tp(0))
00139 {
00140 __Jnu = _Tp(1);
00141 __Jpnu = _Tp(0);
00142 }
00143 else if (__nu == _Tp(1))
00144 {
00145 __Jnu = _Tp(0);
00146 __Jpnu = _Tp(0.5L);
00147 }
00148 else
00149 {
00150 __Jnu = _Tp(0);
00151 __Jpnu = _Tp(0);
00152 }
00153 __Nnu = -std::numeric_limits<_Tp>::infinity();
00154 __Npnu = std::numeric_limits<_Tp>::infinity();
00155 return;
00156 }
00157
00158 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00159
00160
00161
00162
00163 const _Tp __fp_min = std::sqrt(std::numeric_limits<_Tp>::min());
00164 const int __max_iter = 15000;
00165 const _Tp __x_min = _Tp(2);
00166
00167 const int __nl = (__x < __x_min
00168 ? static_cast<int>(__nu + _Tp(0.5L))
00169 : std::max(0, static_cast<int>(__nu - __x + _Tp(1.5L))));
00170
00171 const _Tp __mu = __nu - __nl;
00172 const _Tp __mu2 = __mu * __mu;
00173 const _Tp __xi = _Tp(1) / __x;
00174 const _Tp __xi2 = _Tp(2) * __xi;
00175 _Tp __w = __xi2 / __numeric_constants<_Tp>::__pi();
00176 int __isign = 1;
00177 _Tp __h = __nu * __xi;
00178 if (__h < __fp_min)
00179 __h = __fp_min;
00180 _Tp __b = __xi2 * __nu;
00181 _Tp __d = _Tp(0);
00182 _Tp __c = __h;
00183 int __i;
00184 for (__i = 1; __i <= __max_iter; ++__i)
00185 {
00186 __b += __xi2;
00187 __d = __b - __d;
00188 if (std::abs(__d) < __fp_min)
00189 __d = __fp_min;
00190 __c = __b - _Tp(1) / __c;
00191 if (std::abs(__c) < __fp_min)
00192 __c = __fp_min;
00193 __d = _Tp(1) / __d;
00194 const _Tp __del = __c * __d;
00195 __h *= __del;
00196 if (__d < _Tp(0))
00197 __isign = -__isign;
00198 if (std::abs(__del - _Tp(1)) < __eps)
00199 break;
00200 }
00201 if (__i > __max_iter)
00202 std::__throw_runtime_error(__N("Argument x too large in __bessel_jn; "
00203 "try asymptotic expansion."));
00204 _Tp __Jnul = __isign * __fp_min;
00205 _Tp __Jpnul = __h * __Jnul;
00206 _Tp __Jnul1 = __Jnul;
00207 _Tp __Jpnu1 = __Jpnul;
00208 _Tp __fact = __nu * __xi;
00209 for ( int __l = __nl; __l >= 1; --__l )
00210 {
00211 const _Tp __Jnutemp = __fact * __Jnul + __Jpnul;
00212 __fact -= __xi;
00213 __Jpnul = __fact * __Jnutemp - __Jnul;
00214 __Jnul = __Jnutemp;
00215 }
00216 if (__Jnul == _Tp(0))
00217 __Jnul = __eps;
00218 _Tp __f= __Jpnul / __Jnul;
00219 _Tp __Nmu, __Nnu1, __Npmu, __Jmu;
00220 if (__x < __x_min)
00221 {
00222 const _Tp __x2 = __x / _Tp(2);
00223 const _Tp __pimu = __numeric_constants<_Tp>::__pi() * __mu;
00224 _Tp __fact = (std::abs(__pimu) < __eps
00225 ? _Tp(1) : __pimu / std::sin(__pimu));
00226 _Tp __d = -std::log(__x2);
00227 _Tp __e = __mu * __d;
00228 _Tp __fact2 = (std::abs(__e) < __eps
00229 ? _Tp(1) : std::sinh(__e) / __e);
00230 _Tp __gam1, __gam2, __gampl, __gammi;
00231 __gamma_temme(__mu, __gam1, __gam2, __gampl, __gammi);
00232 _Tp __ff = (_Tp(2) / __numeric_constants<_Tp>::__pi())
00233 * __fact * (__gam1 * std::cosh(__e) + __gam2 * __fact2 * __d);
00234 __e = std::exp(__e);
00235 _Tp __p = __e / (__numeric_constants<_Tp>::__pi() * __gampl);
00236 _Tp __q = _Tp(1) / (__e * __numeric_constants<_Tp>::__pi() * __gammi);
00237 const _Tp __pimu2 = __pimu / _Tp(2);
00238 _Tp __fact3 = (std::abs(__pimu2) < __eps
00239 ? _Tp(1) : std::sin(__pimu2) / __pimu2 );
00240 _Tp __r = __numeric_constants<_Tp>::__pi() * __pimu2 * __fact3 * __fact3;
00241 _Tp __c = _Tp(1);
00242 __d = -__x2 * __x2;
00243 _Tp __sum = __ff + __r * __q;
00244 _Tp __sum1 = __p;
00245 for (__i = 1; __i <= __max_iter; ++__i)
00246 {
00247 __ff = (__i * __ff + __p + __q) / (__i * __i - __mu2);
00248 __c *= __d / _Tp(__i);
00249 __p /= _Tp(__i) - __mu;
00250 __q /= _Tp(__i) + __mu;
00251 const _Tp __del = __c * (__ff + __r * __q);
00252 __sum += __del;
00253 const _Tp __del1 = __c * __p - __i * __del;
00254 __sum1 += __del1;
00255 if ( std::abs(__del) < __eps * (_Tp(1) + std::abs(__sum)) )
00256 break;
00257 }
00258 if ( __i > __max_iter )
00259 std::__throw_runtime_error(__N("Bessel y series failed to converge "
00260 "in __bessel_jn."));
00261 __Nmu = -__sum;
00262 __Nnu1 = -__sum1 * __xi2;
00263 __Npmu = __mu * __xi * __Nmu - __Nnu1;
00264 __Jmu = __w / (__Npmu - __f * __Nmu);
00265 }
00266 else
00267 {
00268 _Tp __a = _Tp(0.25L) - __mu2;
00269 _Tp __q = _Tp(1);
00270 _Tp __p = -__xi / _Tp(2);
00271 _Tp __br = _Tp(2) * __x;
00272 _Tp __bi = _Tp(2);
00273 _Tp __fact = __a * __xi / (__p * __p + __q * __q);
00274 _Tp __cr = __br + __q * __fact;
00275 _Tp __ci = __bi + __p * __fact;
00276 _Tp __den = __br * __br + __bi * __bi;
00277 _Tp __dr = __br / __den;
00278 _Tp __di = -__bi / __den;
00279 _Tp __dlr = __cr * __dr - __ci * __di;
00280 _Tp __dli = __cr * __di + __ci * __dr;
00281 _Tp __temp = __p * __dlr - __q * __dli;
00282 __q = __p * __dli + __q * __dlr;
00283 __p = __temp;
00284 int __i;
00285 for (__i = 2; __i <= __max_iter; ++__i)
00286 {
00287 __a += _Tp(2 * (__i - 1));
00288 __bi += _Tp(2);
00289 __dr = __a * __dr + __br;
00290 __di = __a * __di + __bi;
00291 if (std::abs(__dr) + std::abs(__di) < __fp_min)
00292 __dr = __fp_min;
00293 __fact = __a / (__cr * __cr + __ci * __ci);
00294 __cr = __br + __cr * __fact;
00295 __ci = __bi - __ci * __fact;
00296 if (std::abs(__cr) + std::abs(__ci) < __fp_min)
00297 __cr = __fp_min;
00298 __den = __dr * __dr + __di * __di;
00299 __dr /= __den;
00300 __di /= -__den;
00301 __dlr = __cr * __dr - __ci * __di;
00302 __dli = __cr * __di + __ci * __dr;
00303 __temp = __p * __dlr - __q * __dli;
00304 __q = __p * __dli + __q * __dlr;
00305 __p = __temp;
00306 if (std::abs(__dlr - _Tp(1)) + std::abs(__dli) < __eps)
00307 break;
00308 }
00309 if (__i > __max_iter)
00310 std::__throw_runtime_error(__N("Lentz's method failed "
00311 "in __bessel_jn."));
00312 const _Tp __gam = (__p - __f) / __q;
00313 __Jmu = std::sqrt(__w / ((__p - __f) * __gam + __q));
00314 #if _GLIBCXX_USE_C99_MATH_TR1
00315 __Jmu = std::tr1::copysign(__Jmu, __Jnul);
00316 #else
00317 if (__Jmu * __Jnul < _Tp(0))
00318 __Jmu = -__Jmu;
00319 #endif
00320 __Nmu = __gam * __Jmu;
00321 __Npmu = (__p + __q / __gam) * __Nmu;
00322 __Nnu1 = __mu * __xi * __Nmu - __Npmu;
00323 }
00324 __fact = __Jmu / __Jnul;
00325 __Jnu = __fact * __Jnul1;
00326 __Jpnu = __fact * __Jpnu1;
00327 for (__i = 1; __i <= __nl; ++__i)
00328 {
00329 const _Tp __Nnutemp = (__mu + __i) * __xi2 * __Nnu1 - __Nmu;
00330 __Nmu = __Nnu1;
00331 __Nnu1 = __Nnutemp;
00332 }
00333 __Nnu = __Nmu;
00334 __Npnu = __nu * __xi * __Nmu - __Nnu1;
00335
00336 return;
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 template <typename _Tp>
00357 void
00358 __cyl_bessel_jn_asymp(const _Tp __nu, const _Tp __x,
00359 _Tp & __Jnu, _Tp & __Nnu)
00360 {
00361 const _Tp __coef = std::sqrt(_Tp(2)
00362 / (__numeric_constants<_Tp>::__pi() * __x));
00363 const _Tp __mu = _Tp(4) * __nu * __nu;
00364 const _Tp __mum1 = __mu - _Tp(1);
00365 const _Tp __mum9 = __mu - _Tp(9);
00366 const _Tp __mum25 = __mu - _Tp(25);
00367 const _Tp __mum49 = __mu - _Tp(49);
00368 const _Tp __xx = _Tp(64) * __x * __x;
00369 const _Tp __P = _Tp(1) - __mum1 * __mum9 / (_Tp(2) * __xx)
00370 * (_Tp(1) - __mum25 * __mum49 / (_Tp(12) * __xx));
00371 const _Tp __Q = __mum1 / (_Tp(8) * __x)
00372 * (_Tp(1) - __mum9 * __mum25 / (_Tp(6) * __xx));
00373
00374 const _Tp __chi = __x - (__nu + _Tp(0.5L))
00375 * __numeric_constants<_Tp>::__pi_2();
00376 const _Tp __c = std::cos(__chi);
00377 const _Tp __s = std::sin(__chi);
00378
00379 __Jnu = __coef * (__c * __P - __s * __Q);
00380 __Nnu = __coef * (__s * __P + __c * __Q);
00381
00382 return;
00383 }
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 template <typename _Tp>
00414 _Tp
00415 __cyl_bessel_ij_series(const _Tp __nu, const _Tp __x, const _Tp __sgn,
00416 const unsigned int __max_iter)
00417 {
00418
00419 const _Tp __x2 = __x / _Tp(2);
00420 _Tp __fact = __nu * std::log(__x2);
00421 #if _GLIBCXX_USE_C99_MATH_TR1
00422 __fact -= std::tr1::lgamma(__nu + _Tp(1));
00423 #else
00424 __fact -= __log_gamma(__nu + _Tp(1));
00425 #endif
00426 __fact = std::exp(__fact);
00427 const _Tp __xx4 = __sgn * __x2 * __x2;
00428 _Tp __Jn = _Tp(1);
00429 _Tp __term = _Tp(1);
00430
00431 for (unsigned int __i = 1; __i < __max_iter; ++__i)
00432 {
00433 __term *= __xx4 / (_Tp(__i) * (__nu + _Tp(__i)));
00434 __Jn += __term;
00435 if (std::abs(__term / __Jn) < std::numeric_limits<_Tp>::epsilon())
00436 break;
00437 }
00438
00439 return __fact * __Jn;
00440 }
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 template<typename _Tp>
00458 _Tp
00459 __cyl_bessel_j(const _Tp __nu, const _Tp __x)
00460 {
00461 if (__nu < _Tp(0) || __x < _Tp(0))
00462 std::__throw_domain_error(__N("Bad argument "
00463 "in __cyl_bessel_j."));
00464 else if (__isnan(__nu) || __isnan(__x))
00465 return std::numeric_limits<_Tp>::quiet_NaN();
00466 else if (__x * __x < _Tp(10) * (__nu + _Tp(1)))
00467 return __cyl_bessel_ij_series(__nu, __x, -_Tp(1), 200);
00468 else if (__x > _Tp(1000))
00469 {
00470 _Tp __J_nu, __N_nu;
00471 __cyl_bessel_jn_asymp(__nu, __x, __J_nu, __N_nu);
00472 return __J_nu;
00473 }
00474 else
00475 {
00476 _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;
00477 __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);
00478 return __J_nu;
00479 }
00480 }
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499 template<typename _Tp>
00500 _Tp
00501 __cyl_neumann_n(const _Tp __nu, const _Tp __x)
00502 {
00503 if (__nu < _Tp(0) || __x < _Tp(0))
00504 std::__throw_domain_error(__N("Bad argument "
00505 "in __cyl_neumann_n."));
00506 else if (__isnan(__nu) || __isnan(__x))
00507 return std::numeric_limits<_Tp>::quiet_NaN();
00508 else if (__x > _Tp(1000))
00509 {
00510 _Tp __J_nu, __N_nu;
00511 __cyl_bessel_jn_asymp(__nu, __x, __J_nu, __N_nu);
00512 return __N_nu;
00513 }
00514 else
00515 {
00516 _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;
00517 __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);
00518 return __N_nu;
00519 }
00520 }
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 template <typename _Tp>
00537 void
00538 __sph_bessel_jn(const unsigned int __n, const _Tp __x,
00539 _Tp & __j_n, _Tp & __n_n, _Tp & __jp_n, _Tp & __np_n)
00540 {
00541 const _Tp __nu = _Tp(__n) + _Tp(0.5L);
00542
00543 _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;
00544 __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);
00545
00546 const _Tp __factor = __numeric_constants<_Tp>::__sqrtpio2()
00547 / std::sqrt(__x);
00548
00549 __j_n = __factor * __J_nu;
00550 __n_n = __factor * __N_nu;
00551 __jp_n = __factor * __Jp_nu - __j_n / (_Tp(2) * __x);
00552 __np_n = __factor * __Np_nu - __n_n / (_Tp(2) * __x);
00553
00554 return;
00555 }
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 template <typename _Tp>
00572 _Tp
00573 __sph_bessel(const unsigned int __n, const _Tp __x)
00574 {
00575 if (__x < _Tp(0))
00576 std::__throw_domain_error(__N("Bad argument "
00577 "in __sph_bessel."));
00578 else if (__isnan(__x))
00579 return std::numeric_limits<_Tp>::quiet_NaN();
00580 else if (__x == _Tp(0))
00581 {
00582 if (__n == 0)
00583 return _Tp(1);
00584 else
00585 return _Tp(0);
00586 }
00587 else
00588 {
00589 _Tp __j_n, __n_n, __jp_n, __np_n;
00590 __sph_bessel_jn(__n, __x, __j_n, __n_n, __jp_n, __np_n);
00591 return __j_n;
00592 }
00593 }
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 template <typename _Tp>
00610 _Tp
00611 __sph_neumann(const unsigned int __n, const _Tp __x)
00612 {
00613 if (__x < _Tp(0))
00614 std::__throw_domain_error(__N("Bad argument "
00615 "in __sph_neumann."));
00616 else if (__isnan(__x))
00617 return std::numeric_limits<_Tp>::quiet_NaN();
00618 else if (__x == _Tp(0))
00619 return -std::numeric_limits<_Tp>::infinity();
00620 else
00621 {
00622 _Tp __j_n, __n_n, __jp_n, __np_n;
00623 __sph_bessel_jn(__n, __x, __j_n, __n_n, __jp_n, __np_n);
00624 return __n_n;
00625 }
00626 }
00627
00628 }
00629 }
00630 }
00631
00632 #endif // _GLIBCXX_TR1_BESSEL_FUNCTION_TCC