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