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_MODIFIED_BESSEL_FUNC_TCC
00053 #define _GLIBCXX_TR1_MODIFIED_BESSEL_FUNC_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 template <typename _Tp>
00085 void
00086 __bessel_ik(const _Tp __nu, const _Tp __x,
00087 _Tp & __Inu, _Tp & __Knu, _Tp & __Ipnu, _Tp & __Kpnu)
00088 {
00089 if (__x == _Tp(0))
00090 {
00091 if (__nu == _Tp(0))
00092 {
00093 __Inu = _Tp(1);
00094 __Ipnu = _Tp(0);
00095 }
00096 else if (__nu == _Tp(1))
00097 {
00098 __Inu = _Tp(0);
00099 __Ipnu = _Tp(0.5L);
00100 }
00101 else
00102 {
00103 __Inu = _Tp(0);
00104 __Ipnu = _Tp(0);
00105 }
00106 __Knu = std::numeric_limits<_Tp>::infinity();
00107 __Kpnu = -std::numeric_limits<_Tp>::infinity();
00108 return;
00109 }
00110
00111 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00112 const _Tp __fp_min = _Tp(10) * std::numeric_limits<_Tp>::epsilon();
00113 const int __max_iter = 15000;
00114 const _Tp __x_min = _Tp(2);
00115
00116 const int __nl = static_cast<int>(__nu + _Tp(0.5L));
00117
00118 const _Tp __mu = __nu - __nl;
00119 const _Tp __mu2 = __mu * __mu;
00120 const _Tp __xi = _Tp(1) / __x;
00121 const _Tp __xi2 = _Tp(2) * __xi;
00122 _Tp __h = __nu * __xi;
00123 if ( __h < __fp_min )
00124 __h = __fp_min;
00125 _Tp __b = __xi2 * __nu;
00126 _Tp __d = _Tp(0);
00127 _Tp __c = __h;
00128 int __i;
00129 for ( __i = 1; __i <= __max_iter; ++__i )
00130 {
00131 __b += __xi2;
00132 __d = _Tp(1) / (__b + __d);
00133 __c = __b + _Tp(1) / __c;
00134 const _Tp __del = __c * __d;
00135 __h *= __del;
00136 if (std::abs(__del - _Tp(1)) < __eps)
00137 break;
00138 }
00139 if (__i > __max_iter)
00140 std::__throw_runtime_error(__N("Argument x too large "
00141 "in __bessel_jn; "
00142 "try asymptotic expansion."));
00143 _Tp __Inul = __fp_min;
00144 _Tp __Ipnul = __h * __Inul;
00145 _Tp __Inul1 = __Inul;
00146 _Tp __Ipnu1 = __Ipnul;
00147 _Tp __fact = __nu * __xi;
00148 for (int __l = __nl; __l >= 1; --__l)
00149 {
00150 const _Tp __Inutemp = __fact * __Inul + __Ipnul;
00151 __fact -= __xi;
00152 __Ipnul = __fact * __Inutemp + __Inul;
00153 __Inul = __Inutemp;
00154 }
00155 _Tp __f = __Ipnul / __Inul;
00156 _Tp __Kmu, __Knu1;
00157 if (__x < __x_min)
00158 {
00159 const _Tp __x2 = __x / _Tp(2);
00160 const _Tp __pimu = __numeric_constants<_Tp>::__pi() * __mu;
00161 const _Tp __fact = (std::abs(__pimu) < __eps
00162 ? _Tp(1) : __pimu / std::sin(__pimu));
00163 _Tp __d = -std::log(__x2);
00164 _Tp __e = __mu * __d;
00165 const _Tp __fact2 = (std::abs(__e) < __eps
00166 ? _Tp(1) : std::sinh(__e) / __e);
00167 _Tp __gam1, __gam2, __gampl, __gammi;
00168 __gamma_temme(__mu, __gam1, __gam2, __gampl, __gammi);
00169 _Tp __ff = __fact
00170 * (__gam1 * std::cosh(__e) + __gam2 * __fact2 * __d);
00171 _Tp __sum = __ff;
00172 __e = std::exp(__e);
00173 _Tp __p = __e / (_Tp(2) * __gampl);
00174 _Tp __q = _Tp(1) / (_Tp(2) * __e * __gammi);
00175 _Tp __c = _Tp(1);
00176 __d = __x2 * __x2;
00177 _Tp __sum1 = __p;
00178 int __i;
00179 for (__i = 1; __i <= __max_iter; ++__i)
00180 {
00181 __ff = (__i * __ff + __p + __q) / (__i * __i - __mu2);
00182 __c *= __d / __i;
00183 __p /= __i - __mu;
00184 __q /= __i + __mu;
00185 const _Tp __del = __c * __ff;
00186 __sum += __del;
00187 const _Tp __del1 = __c * (__p - __i * __ff);
00188 __sum1 += __del1;
00189 if (std::abs(__del) < __eps * std::abs(__sum))
00190 break;
00191 }
00192 if (__i > __max_iter)
00193 std::__throw_runtime_error(__N("Bessel k series failed to converge "
00194 "in __bessel_jn."));
00195 __Kmu = __sum;
00196 __Knu1 = __sum1 * __xi2;
00197 }
00198 else
00199 {
00200 _Tp __b = _Tp(2) * (_Tp(1) + __x);
00201 _Tp __d = _Tp(1) / __b;
00202 _Tp __delh = __d;
00203 _Tp __h = __delh;
00204 _Tp __q1 = _Tp(0);
00205 _Tp __q2 = _Tp(1);
00206 _Tp __a1 = _Tp(0.25L) - __mu2;
00207 _Tp __q = __c = __a1;
00208 _Tp __a = -__a1;
00209 _Tp __s = _Tp(1) + __q * __delh;
00210 int __i;
00211 for (__i = 2; __i <= __max_iter; ++__i)
00212 {
00213 __a -= 2 * (__i - 1);
00214 __c = -__a * __c / __i;
00215 const _Tp __qnew = (__q1 - __b * __q2) / __a;
00216 __q1 = __q2;
00217 __q2 = __qnew;
00218 __q += __c * __qnew;
00219 __b += _Tp(2);
00220 __d = _Tp(1) / (__b + __a * __d);
00221 __delh = (__b * __d - _Tp(1)) * __delh;
00222 __h += __delh;
00223 const _Tp __dels = __q * __delh;
00224 __s += __dels;
00225 if ( std::abs(__dels / __s) < __eps )
00226 break;
00227 }
00228 if (__i > __max_iter)
00229 std::__throw_runtime_error(__N("Steed's method failed "
00230 "in __bessel_jn."));
00231 __h = __a1 * __h;
00232 __Kmu = std::sqrt(__numeric_constants<_Tp>::__pi() / (_Tp(2) * __x))
00233 * std::exp(-__x) / __s;
00234 __Knu1 = __Kmu * (__mu + __x + _Tp(0.5L) - __h) * __xi;
00235 }
00236
00237 _Tp __Kpmu = __mu * __xi * __Kmu - __Knu1;
00238 _Tp __Inumu = __xi / (__f * __Kmu - __Kpmu);
00239 __Inu = __Inumu * __Inul1 / __Inul;
00240 __Ipnu = __Inumu * __Ipnu1 / __Inul;
00241 for ( __i = 1; __i <= __nl; ++__i )
00242 {
00243 const _Tp __Knutemp = (__mu + __i) * __xi2 * __Knu1 + __Kmu;
00244 __Kmu = __Knu1;
00245 __Knu1 = __Knutemp;
00246 }
00247 __Knu = __Kmu;
00248 __Kpnu = __nu * __xi * __Kmu - __Knu1;
00249
00250 return;
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 template<typename _Tp>
00269 _Tp
00270 __cyl_bessel_i(const _Tp __nu, const _Tp __x)
00271 {
00272 if (__nu < _Tp(0) || __x < _Tp(0))
00273 std::__throw_domain_error(__N("Bad argument "
00274 "in __cyl_bessel_i."));
00275 else if (__isnan(__nu) || __isnan(__x))
00276 return std::numeric_limits<_Tp>::quiet_NaN();
00277 else if (__x * __x < _Tp(10) * (__nu + _Tp(1)))
00278 return __cyl_bessel_ij_series(__nu, __x, +_Tp(1), 200);
00279 else
00280 {
00281 _Tp __I_nu, __K_nu, __Ip_nu, __Kp_nu;
00282 __bessel_ik(__nu, __x, __I_nu, __K_nu, __Ip_nu, __Kp_nu);
00283 return __I_nu;
00284 }
00285 }
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 template<typename _Tp>
00305 _Tp
00306 __cyl_bessel_k(const _Tp __nu, const _Tp __x)
00307 {
00308 if (__nu < _Tp(0) || __x < _Tp(0))
00309 std::__throw_domain_error(__N("Bad argument "
00310 "in __cyl_bessel_k."));
00311 else if (__isnan(__nu) || __isnan(__x))
00312 return std::numeric_limits<_Tp>::quiet_NaN();
00313 else
00314 {
00315 _Tp __I_nu, __K_nu, __Ip_nu, __Kp_nu;
00316 __bessel_ik(__nu, __x, __I_nu, __K_nu, __Ip_nu, __Kp_nu);
00317 return __K_nu;
00318 }
00319 }
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 template <typename _Tp>
00339 void
00340 __sph_bessel_ik(const unsigned int __n, const _Tp __x,
00341 _Tp & __i_n, _Tp & __k_n, _Tp & __ip_n, _Tp & __kp_n)
00342 {
00343 const _Tp __nu = _Tp(__n) + _Tp(0.5L);
00344
00345 _Tp __I_nu, __Ip_nu, __K_nu, __Kp_nu;
00346 __bessel_ik(__nu, __x, __I_nu, __K_nu, __Ip_nu, __Kp_nu);
00347
00348 const _Tp __factor = __numeric_constants<_Tp>::__sqrtpio2()
00349 / std::sqrt(__x);
00350
00351 __i_n = __factor * __I_nu;
00352 __k_n = __factor * __K_nu;
00353 __ip_n = __factor * __Ip_nu - __i_n / (_Tp(2) * __x);
00354 __kp_n = __factor * __Kp_nu - __k_n / (_Tp(2) * __x);
00355
00356 return;
00357 }
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 template <typename _Tp>
00374 void
00375 __airy(const _Tp __x,
00376 _Tp & __Ai, _Tp & __Bi, _Tp & __Aip, _Tp & __Bip)
00377 {
00378 const _Tp __absx = std::abs(__x);
00379 const _Tp __rootx = std::sqrt(__absx);
00380 const _Tp __z = _Tp(2) * __absx * __rootx / _Tp(3);
00381
00382 if (__isnan(__x))
00383 return std::numeric_limits<_Tp>::quiet_NaN();
00384 else if (__x > _Tp(0))
00385 {
00386 _Tp __I_nu, __Ip_nu, __K_nu, __Kp_nu;
00387
00388 __bessel_ik(_Tp(1) / _Tp(3), __z, __I_nu, __K_nu, __Ip_nu, __Kp_nu);
00389 __Ai = __rootx * __K_nu
00390 / (__numeric_constants<_Tp>::__sqrt3()
00391 * __numeric_constants<_Tp>::__pi());
00392 __Bi = __rootx * (__K_nu / __numeric_constants<_Tp>::__pi()
00393 + _Tp(2) * __I_nu / __numeric_constants<_Tp>::__sqrt3());
00394
00395 __bessel_ik(_Tp(2) / _Tp(3), __z, __I_nu, __K_nu, __Ip_nu, __Kp_nu);
00396 __Aip = -__x * __K_nu
00397 / (__numeric_constants<_Tp>::__sqrt3()
00398 * __numeric_constants<_Tp>::__pi());
00399 __Bip = __x * (__K_nu / __numeric_constants<_Tp>::__pi()
00400 + _Tp(2) * __I_nu
00401 / __numeric_constants<_Tp>::__sqrt3());
00402 }
00403 else if (__x < _Tp(0))
00404 {
00405 _Tp __J_nu, __Jp_nu, __N_nu, __Np_nu;
00406
00407 __bessel_jn(_Tp(1) / _Tp(3), __z, __J_nu, __N_nu, __Jp_nu, __Np_nu);
00408 __Ai = __rootx * (__J_nu
00409 - __N_nu / __numeric_constants<_Tp>::__sqrt3()) / _Tp(2);
00410 __Bi = -__rootx * (__N_nu
00411 + __J_nu / __numeric_constants<_Tp>::__sqrt3()) / _Tp(2);
00412
00413 __bessel_jn(_Tp(2) / _Tp(3), __z, __J_nu, __N_nu, __Jp_nu, __Np_nu);
00414 __Aip = __absx * (__N_nu / __numeric_constants<_Tp>::__sqrt3()
00415 + __J_nu) / _Tp(2);
00416 __Bip = __absx * (__J_nu / __numeric_constants<_Tp>::__sqrt3()
00417 - __N_nu) / _Tp(2);
00418 }
00419 else
00420 {
00421
00422
00423
00424 __Ai = _Tp(0.35502805388781723926L);
00425 __Bi = __Ai * __numeric_constants<_Tp>::__sqrt3();
00426
00427
00428
00429
00430 __Aip = -_Tp(0.25881940379280679840L);
00431 __Bip = -__Aip * __numeric_constants<_Tp>::__sqrt3();
00432 }
00433
00434 return;
00435 }
00436
00437 }
00438 }
00439 }
00440
00441 #endif // _GLIBCXX_TR1_MODIFIED_BESSEL_FUNC_TCC