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