cloudy
trunk
|
00001 /* This file contains routines (perhaps in modified form) written by third parties. 00002 * Use and distribution of these works are determined by their respective copyrights. */ 00003 #include "cddefines.h" 00004 #include "thirdparty.h" 00005 #include "physconst.h" 00006 00007 inline double polevl(double x, const double coef[], int N); 00008 inline double p1evl(double x, const double coef[], int N); 00009 inline double chbevl(double, const double[], int); 00010 00011 00012 /* the routine linfit was derived from the program slopes.f */ 00013 00014 /********************************************************************/ 00015 /* */ 00016 /* PROGRAM SLOPES */ 00017 /* */ 00018 /* PROGRAM TO COMPUTE THE THEORETICAL REGRESSION SLOPES */ 00019 /* AND UNCERTAINTIES (VIA DELTA METHOD), AND UNCERTAINTIES */ 00020 /* VIA BOOTSTRAP AND BIVARIATE NORMAL SIMULATION FOR A */ 00021 /* (X(I),Y(I)) DATA SET WITH UNKNOWN POPULATION DISTRIBUTION. */ 00022 /* */ 00023 /* WRITTEN BY ERIC FEIGELSON, PENN STATE. JUNE 1991 */ 00024 /* */ 00025 /* */ 00026 /* A full description of these methods can be found in: */ 00027 /* Isobe, T., Feigelson, E. D., Akritas, M. and Babu, G. J., */ 00028 /* Linear Regression in Astronomy I, Astrophys. J. 364, 104 */ 00029 /* (1990) */ 00030 /* Babu, G. J. and Feigelson, E. D., Analytical and Monte Carlo */ 00031 /* Comparisons of Six Different Linear Least Squares Fits, */ 00032 /* Communications in Statistics, Simulation & Computation, */ 00033 /* 21, 533 (1992) */ 00034 /* Feigelson, E. D. and Babu, G. J., Linear Regression in */ 00035 /* Astronomy II, Astrophys. J. 397, 55 (1992). */ 00036 /* */ 00037 /********************************************************************/ 00038 00039 /* this used to be sixlin, but only the first fit has been retained */ 00040 00041 /********************************************************************/ 00042 /************************* routine linfit ***************************/ 00043 /********************************************************************/ 00044 00045 bool linfit( 00046 long n, 00047 double x[], /* x[n] */ 00048 double y[], /* y[n] */ 00049 double &a, 00050 double &siga, 00051 double &b, 00052 double &sigb 00053 ) 00054 { 00055 00056 /* 00057 * linear regression 00058 * written by t. isobe, g. j. babu and e. d. feigelson 00059 * center for space research, m.i.t. 00060 * and 00061 * the pennsylvania state university 00062 * 00063 * rev. 1.0, september 1990 00064 * 00065 * this subroutine provides linear regression coefficients 00066 * computed by six different methods described in isobe, 00067 * feigelson, akritas, and babu 1990, astrophysical journal 00068 * and babu and feigelson 1990, subm. to technometrics. 00069 * the methods are ols(y/x), ols(x/y), ols bisector, orthogonal, 00070 * reduced major axis, and mean-ols regressions. 00071 * 00072 * [Modified and translated to C/C++ by Peter van Hoof, Royal 00073 * Observatory of Belgium; only the first method has been retained] 00074 * 00075 * 00076 * input 00077 * x[i] : independent variable 00078 * y[i] : dependent variable 00079 * n : number of data points 00080 * 00081 * output 00082 * a : intercept coefficients 00083 * b : slope coefficients 00084 * siga : standard deviations of intercepts 00085 * sigb : standard deviations of slopes 00086 * 00087 * error returns 00088 * returns true when division by zero will 00089 * occur (i.e. when sxy is zero) 00090 */ 00091 00092 DEBUG_ENTRY( "linfit()" ); 00093 00094 /* initializations */ 00095 a = 0.0; 00096 siga = 0.0; 00097 b = 0.0; 00098 sigb = 0.0; 00099 00100 /* compute averages and sums */ 00101 double s1 = 0.0; 00102 double s2 = 0.0; 00103 for( long i=0; i < n; i++ ) 00104 { 00105 s1 += x[i]; 00106 s2 += y[i]; 00107 } 00108 double rn = (double)n; 00109 double xavg = s1/rn; 00110 double yavg = s2/rn; 00111 double sxx = 0.0; 00112 double sxy = 0.0; 00113 for( long i=0; i < n; i++ ) 00114 { 00115 x[i] -= xavg; 00116 y[i] -= yavg; 00117 sxx += pow2(x[i]); 00118 sxy += x[i]*y[i]; 00119 } 00120 00121 if( sxy == 0.0 ) 00122 { 00123 return true; 00124 } 00125 00126 /* compute the slope coefficient */ 00127 b = sxy/sxx; 00128 00129 /* compute intercept coefficient */ 00130 a = yavg - b*xavg; 00131 00132 /* prepare for computation of variances */ 00133 double sum1 = 0.0; 00134 for( long i=0; i < n; i++ ) 00135 sum1 += pow2(x[i]*(y[i] - b*x[i])); 00136 00137 /* compute variance of the slope coefficient */ 00138 sigb = sum1/pow2(sxx); 00139 00140 /* compute variance of the intercept coefficient */ 00141 for( long i=0; i < n; i++ ) 00142 siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx)); 00143 00144 /* convert variances to standard deviations */ 00145 sigb = sqrt(sigb); 00146 siga = sqrt(siga)/rn; 00147 00148 /* return data arrays to their original form */ 00149 for( long i=0; i < n; i++ ) 00150 { 00151 x[i] += xavg; 00152 y[i] += yavg; 00153 } 00154 return false; 00155 } 00156 00157 00158 /* the routines factorial and lfactorial came originally from hydro_bauman.cpp 00159 * and were written by Robert Paul Bauman. lfactorial was modified by Peter van Hoof */ 00160 00161 /************************************************************************************************/ 00162 /* pre-calculated factorials */ 00163 /************************************************************************************************/ 00164 00165 static const double pre_factorial[NPRE_FACTORIAL] = 00166 { 00167 1.00000000000000000000e+00, 00168 1.00000000000000000000e+00, 00169 2.00000000000000000000e+00, 00170 6.00000000000000000000e+00, 00171 2.40000000000000000000e+01, 00172 1.20000000000000000000e+02, 00173 7.20000000000000000000e+02, 00174 5.04000000000000000000e+03, 00175 4.03200000000000000000e+04, 00176 3.62880000000000000000e+05, 00177 3.62880000000000000000e+06, 00178 3.99168000000000000000e+07, 00179 4.79001600000000000000e+08, 00180 6.22702080000000000000e+09, 00181 8.71782912000000000000e+10, 00182 1.30767436800000000000e+12, 00183 2.09227898880000000000e+13, 00184 3.55687428096000000000e+14, 00185 6.40237370572800000000e+15, 00186 1.21645100408832000000e+17, 00187 2.43290200817664000000e+18, 00188 5.10909421717094400000e+19, 00189 1.12400072777760768000e+21, 00190 2.58520167388849766400e+22, 00191 6.20448401733239439360e+23, 00192 1.55112100433309859840e+25, 00193 4.03291461126605635592e+26, 00194 1.08888694504183521614e+28, 00195 3.04888344611713860511e+29, 00196 8.84176199373970195470e+30, 00197 2.65252859812191058647e+32, 00198 8.22283865417792281807e+33, 00199 2.63130836933693530178e+35, 00200 8.68331761881188649615e+36, 00201 2.95232799039604140861e+38, 00202 1.03331479663861449300e+40, 00203 3.71993326789901217463e+41, 00204 1.37637530912263450457e+43, 00205 5.23022617466601111726e+44, 00206 2.03978820811974433568e+46, 00207 8.15915283247897734264e+47, 00208 3.34525266131638071044e+49, 00209 1.40500611775287989839e+51, 00210 6.04152630633738356321e+52, 00211 2.65827157478844876773e+54, 00212 1.19622220865480194551e+56, 00213 5.50262215981208894940e+57, 00214 2.58623241511168180614e+59, 00215 1.24139155925360726691e+61, 00216 6.08281864034267560801e+62, 00217 3.04140932017133780398e+64, 00218 1.55111875328738227999e+66, 00219 8.06581751709438785585e+67, 00220 4.27488328406002556374e+69, 00221 2.30843697339241380441e+71, 00222 1.26964033536582759243e+73, 00223 7.10998587804863451749e+74, 00224 4.05269195048772167487e+76, 00225 2.35056133128287857145e+78, 00226 1.38683118545689835713e+80, 00227 8.32098711274139014271e+81, 00228 5.07580213877224798711e+83, 00229 3.14699732603879375200e+85, 00230 1.98260831540444006372e+87, 00231 1.26886932185884164078e+89, 00232 8.24765059208247066512e+90, 00233 5.44344939077443063905e+92, 00234 3.64711109181886852801e+94, 00235 2.48003554243683059915e+96, 00236 1.71122452428141311337e+98, 00237 1.19785716699698917933e+100, 00238 8.50478588567862317347e+101, 00239 6.12344583768860868500e+103, 00240 4.47011546151268434004e+105, 00241 3.30788544151938641157e+107, 00242 2.48091408113953980872e+109, 00243 1.88549470166605025458e+111, 00244 1.45183092028285869606e+113, 00245 1.13242811782062978295e+115, 00246 8.94618213078297528506e+116, 00247 7.15694570462638022794e+118, 00248 5.79712602074736798470e+120, 00249 4.75364333701284174746e+122, 00250 3.94552396972065865030e+124, 00251 3.31424013456535326627e+126, 00252 2.81710411438055027626e+128, 00253 2.42270953836727323750e+130, 00254 2.10775729837952771662e+132, 00255 1.85482642257398439069e+134, 00256 1.65079551609084610774e+136, 00257 1.48571596448176149700e+138, 00258 1.35200152767840296226e+140, 00259 1.24384140546413072522e+142, 00260 1.15677250708164157442e+144, 00261 1.08736615665674307994e+146, 00262 1.03299784882390592592e+148, 00263 9.91677934870949688836e+149, 00264 9.61927596824821198159e+151, 00265 9.42689044888324774164e+153, 00266 9.33262154439441526381e+155, 00267 9.33262154439441526388e+157, 00268 9.42594775983835941673e+159, 00269 9.61446671503512660515e+161, 00270 9.90290071648618040340e+163, 00271 1.02990167451456276198e+166, 00272 1.08139675824029090008e+168, 00273 1.14628056373470835406e+170, 00274 1.22652020319613793888e+172, 00275 1.32464181945182897396e+174, 00276 1.44385958320249358163e+176, 00277 1.58824554152274293982e+178, 00278 1.76295255109024466316e+180, 00279 1.97450685722107402277e+182, 00280 2.23119274865981364576e+184, 00281 2.54355973347218755612e+186, 00282 2.92509369349301568964e+188, 00283 3.39310868445189820004e+190, 00284 3.96993716080872089396e+192, 00285 4.68452584975429065488e+194, 00286 5.57458576120760587943e+196, 00287 6.68950291344912705515e+198, 00288 8.09429852527344373681e+200, 00289 9.87504420083360135884e+202, 00290 1.21463043670253296712e+205, 00291 1.50614174151114087918e+207, 00292 1.88267717688892609901e+209, 00293 2.37217324288004688470e+211, 00294 3.01266001845765954361e+213, 00295 3.85620482362580421582e+215, 00296 4.97450422247728743840e+217, 00297 6.46685548922047366972e+219, 00298 8.47158069087882050755e+221, 00299 1.11824865119600430699e+224, 00300 1.48727070609068572828e+226, 00301 1.99294274616151887582e+228, 00302 2.69047270731805048244e+230, 00303 3.65904288195254865604e+232, 00304 5.01288874827499165889e+234, 00305 6.91778647261948848943e+236, 00306 9.61572319694108900019e+238, 00307 1.34620124757175246000e+241, 00308 1.89814375907617096864e+243, 00309 2.69536413788816277557e+245, 00310 3.85437071718007276916e+247, 00311 5.55029383273930478744e+249, 00312 8.04792605747199194159e+251, 00313 1.17499720439091082343e+254, 00314 1.72724589045463891049e+256, 00315 2.55632391787286558753e+258, 00316 3.80892263763056972532e+260, 00317 5.71338395644585458806e+262, 00318 8.62720977423324042775e+264, 00319 1.31133588568345254503e+267, 00320 2.00634390509568239384e+269, 00321 3.08976961384735088657e+271, 00322 4.78914290146339387432e+273, 00323 7.47106292628289444390e+275, 00324 1.17295687942641442768e+278, 00325 1.85327186949373479574e+280, 00326 2.94670227249503832518e+282, 00327 4.71472363599206132029e+284, 00328 7.59070505394721872577e+286, 00329 1.22969421873944943358e+289, 00330 2.00440157654530257674e+291, 00331 3.28721858553429622598e+293, 00332 5.42391066613158877297e+295, 00333 9.00369170577843736335e+297, 00334 1.50361651486499903974e+300, 00335 2.52607574497319838672e+302, 00336 4.26906800900470527345e+304, 00337 7.25741561530799896496e+306 00338 }; 00339 00340 double factorial( long n ) 00341 { 00342 DEBUG_ENTRY( "factorial()" ); 00343 00344 if( n < 0 || n >= NPRE_FACTORIAL ) 00345 { 00346 fprintf( ioQQQ, "factorial: domain error\n" ); 00347 cdEXIT(EXIT_FAILURE); 00348 } 00349 00350 return pre_factorial[n]; 00351 } 00352 00353 /* NB - this implementation is not thread-safe! */ 00354 00355 class t_lfact : public Singleton<t_lfact> 00356 { 00357 friend class Singleton<t_lfact>; 00358 protected: 00359 t_lfact() 00360 { 00361 p_lf.reserve( 512 ); 00362 p_lf.push_back( 0. ); /* log10( 0! ) */ 00363 p_lf.push_back( 0. ); /* log10( 1! ) */ 00364 } 00365 private: 00366 vector<double> p_lf; 00367 public: 00368 double get_lfact( unsigned long n ) 00369 { 00370 if( n < p_lf.size() ) 00371 { 00372 return p_lf[n]; 00373 } 00374 else 00375 { 00376 for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ ) 00377 p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) ); 00378 return p_lf[n]; 00379 } 00380 } 00381 }; 00382 00383 double lfactorial( long n ) 00384 { 00385 /******************************/ 00386 /* n */ 00387 /* --- */ 00388 /* log( n! ) = > log(j) */ 00389 /* --- */ 00390 /* j=1 */ 00391 /******************************/ 00392 00393 DEBUG_ENTRY( "lfactorial()" ); 00394 00395 if( n < 0 ) 00396 { 00397 fprintf( ioQQQ, "lfactorial: domain error\n" ); 00398 cdEXIT(EXIT_FAILURE); 00399 } 00400 00401 return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) ); 00402 } 00403 00404 /* complex Gamma function in double precision */ 00405 /* this routine is a slightly modified version of the one found 00406 * at http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html 00407 * The following copyright applies: 00408 * Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp). 00409 * You may use, copy, modify this code for any purpose and 00410 * without fee. You may distribute this ORIGINAL package. */ 00411 complex<double> cdgamma(complex<double> x) 00412 { 00413 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t; 00414 00415 DEBUG_ENTRY( "cdgamma()" ); 00416 00417 xr = x.real(); 00418 xi = x.imag(); 00419 if(xr < 0) 00420 { 00421 wr = 1. - xr; 00422 wi = -xi; 00423 } 00424 else 00425 { 00426 wr = xr; 00427 wi = xi; 00428 } 00429 ur = wr + 6.00009857740312429; 00430 vr = ur * (wr + 4.99999857982434025) - wi * wi; 00431 vi = wi * (wr + 4.99999857982434025) + ur * wi; 00432 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 + 00433 0.293729529320536228; 00434 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521; 00435 ur = vr * (wr + 4.00000003016801681) - vi * wi; 00436 ui = vi * (wr + 4.00000003016801681) + vr * wi; 00437 vr = ur * (wr + 2.99999999944915534) - ui * wi; 00438 vi = ui * (wr + 2.99999999944915534) + ur * wi; 00439 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063; 00440 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063; 00441 ur = vr * (wr + 2.00000000000603851) - vi * wi; 00442 ui = vi * (wr + 2.00000000000603851) + vr * wi; 00443 vr = ur * (wr + 0.999999999999975753) - ui * wi; 00444 vi = ui * (wr + 0.999999999999975753) + ur * wi; 00445 yr += ur * 10.5400280458730808 + vr; 00446 yi += ui * 10.5400280458730808 + vi; 00447 ur = vr * wr - vi * wi; 00448 ui = vi * wr + vr * wi; 00449 t = ur * ur + ui * ui; 00450 vr = yr * ur + yi * ui + t * 0.0327673720261526849; 00451 vi = yi * ur - yr * ui; 00452 yr = wr + 7.31790632447016203; 00453 ur = log(yr * yr + wi * wi) * 0.5 - 1; 00454 ui = atan2(wi, yr); 00455 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t; 00456 yi = ui * (wr - 0.5) + ur * wi; 00457 ur = yr * cos(yi); 00458 ui = yr * sin(yi); 00459 yr = ur * vr - ui * vi; 00460 yi = ui * vr + ur * vi; 00461 if(xr < 0) 00462 { 00463 wr = xr * 3.14159265358979324; 00464 wi = exp(xi * 3.14159265358979324); 00465 vi = 1 / wi; 00466 ur = (vi + wi) * sin(wr); 00467 ui = (vi - wi) * cos(wr); 00468 vr = ur * yr + ui * yi; 00469 vi = ui * yr - ur * yi; 00470 ur = 6.2831853071795862 / (vr * vr + vi * vi); 00471 yr = ur * vr; 00472 yi = ur * vi; 00473 } 00474 return complex<double>( yr, yi ); 00475 } 00476 00477 /*==================================================================== 00478 * 00479 * Below are routines from the Cephes library. 00480 * 00481 * This is the copyright statement included with the library: 00482 * 00483 * Some software in this archive may be from the book _Methods and 00484 * Programs for Mathematical Functions_ (Prentice-Hall, 1989) or 00485 * from the Cephes Mathematical Library, a commercial product. In 00486 * either event, it is copyrighted by the author. What you see here 00487 * may be used freely but it comes with no support or guarantee. 00488 * 00489 * The two known misprints in the book are repaired here in the 00490 * source listings for the gamma function and the incomplete beta 00491 * integral. 00492 * 00493 * Stephen L. Moshier 00494 * moshier@world.std.com 00495 * 00496 *====================================================================*/ 00497 00498 /* j0.c 00499 * 00500 * Bessel function of order zero 00501 * 00502 * 00503 * 00504 * SYNOPSIS: 00505 * 00506 * double x, y, j0(); 00507 * 00508 * y = j0( x ); 00509 * 00510 * 00511 * 00512 * DESCRIPTION: 00513 * 00514 * Returns Bessel function of order zero of the argument. 00515 * 00516 * The domain is divided into the intervals [0, 5] and 00517 * (5, infinity). In the first interval the following rational 00518 * approximation is used: 00519 * 00520 * 00521 * 2 2 00522 * (w - r ) (w - r ) P (w) / Q (w) 00523 * 1 2 3 8 00524 * 00525 * 2 00526 * where w = x and the two r's are zeros of the function. 00527 * 00528 * In the second interval, the Hankel asymptotic expansion 00529 * is employed with two rational functions of degree 6/6 00530 * and 7/7. 00531 * 00532 * 00533 * 00534 * ACCURACY: 00535 * 00536 * Absolute error: 00537 * arithmetic domain # trials peak rms 00538 * DEC 0, 30 10000 4.4e-17 6.3e-18 00539 * IEEE 0, 30 60000 4.2e-16 1.1e-16 00540 * 00541 */ 00542 /* y0.c 00543 * 00544 * Bessel function of the second kind, order zero 00545 * 00546 * 00547 * 00548 * SYNOPSIS: 00549 * 00550 * double x, y, y0(); 00551 * 00552 * y = y0( x ); 00553 * 00554 * 00555 * 00556 * DESCRIPTION: 00557 * 00558 * Returns Bessel function of the second kind, of order 00559 * zero, of the argument. 00560 * 00561 * The domain is divided into the intervals [0, 5] and 00562 * (5, infinity). In the first interval a rational approximation 00563 * R(x) is employed to compute 00564 * y0(x) = R(x) + 2 * log(x) * j0(x) / PI. 00565 * Thus a call to j0() is required. 00566 * 00567 * In the second interval, the Hankel asymptotic expansion 00568 * is employed with two rational functions of degree 6/6 00569 * and 7/7. 00570 * 00571 * 00572 * 00573 * ACCURACY: 00574 * 00575 * Absolute error, when y0(x) < 1; else relative error: 00576 * 00577 * arithmetic domain # trials peak rms 00578 * DEC 0, 30 9400 7.0e-17 7.9e-18 00579 * IEEE 0, 30 30000 1.3e-15 1.6e-16 00580 * 00581 */ 00582 00583 /* 00584 Cephes Math Library Release 2.8: June, 2000 00585 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 00586 */ 00587 00588 /* Note: all coefficients satisfy the relative error criterion 00589 * except YP, YQ which are designed for absolute error. */ 00590 00591 static const double b0_PP[7] = { 00592 7.96936729297347051624e-4, 00593 8.28352392107440799803e-2, 00594 1.23953371646414299388e0, 00595 5.44725003058768775090e0, 00596 8.74716500199817011941e0, 00597 5.30324038235394892183e0, 00598 9.99999999999999997821e-1, 00599 }; 00600 00601 static const double b0_PQ[7] = { 00602 9.24408810558863637013e-4, 00603 8.56288474354474431428e-2, 00604 1.25352743901058953537e0, 00605 5.47097740330417105182e0, 00606 8.76190883237069594232e0, 00607 5.30605288235394617618e0, 00608 1.00000000000000000218e0, 00609 }; 00610 00611 static const double b0_QP[8] = { 00612 -1.13663838898469149931e-2, 00613 -1.28252718670509318512e0, 00614 -1.95539544257735972385e1, 00615 -9.32060152123768231369e1, 00616 -1.77681167980488050595e2, 00617 -1.47077505154951170175e2, 00618 -5.14105326766599330220e1, 00619 -6.05014350600728481186e0, 00620 }; 00621 00622 static const double b0_QQ[7] = { 00623 /* 1.00000000000000000000e0,*/ 00624 6.43178256118178023184e1, 00625 8.56430025976980587198e2, 00626 3.88240183605401609683e3, 00627 7.24046774195652478189e3, 00628 5.93072701187316984827e3, 00629 2.06209331660327847417e3, 00630 2.42005740240291393179e2, 00631 }; 00632 00633 static const double b0_YP[8] = { 00634 1.55924367855235737965e4, 00635 -1.46639295903971606143e7, 00636 5.43526477051876500413e9, 00637 -9.82136065717911466409e11, 00638 8.75906394395366999549e13, 00639 -3.46628303384729719441e15, 00640 4.42733268572569800351e16, 00641 -1.84950800436986690637e16, 00642 }; 00643 00644 static const double b0_YQ[7] = { 00645 /* 1.00000000000000000000e0,*/ 00646 1.04128353664259848412e3, 00647 6.26107330137134956842e5, 00648 2.68919633393814121987e8, 00649 8.64002487103935000337e10, 00650 2.02979612750105546709e13, 00651 3.17157752842975028269e15, 00652 2.50596256172653059228e17, 00653 }; 00654 00655 /* 5.783185962946784521175995758455807035071 */ 00656 static const double DR1 = 5.78318596294678452118e0; 00657 /* 30.47126234366208639907816317502275584842 */ 00658 static const double DR2 = 3.04712623436620863991e1; 00659 00660 static double b0_RP[4] = { 00661 -4.79443220978201773821e9, 00662 1.95617491946556577543e12, 00663 -2.49248344360967716204e14, 00664 9.70862251047306323952e15, 00665 }; 00666 00667 static double b0_RQ[8] = { 00668 /* 1.00000000000000000000e0,*/ 00669 4.99563147152651017219e2, 00670 1.73785401676374683123e5, 00671 4.84409658339962045305e7, 00672 1.11855537045356834862e10, 00673 2.11277520115489217587e12, 00674 3.10518229857422583814e14, 00675 3.18121955943204943306e16, 00676 1.71086294081043136091e18, 00677 }; 00678 00679 static const double TWOOPI = 2./PI; 00680 static const double SQ2OPI = sqrt(2./PI); 00681 static const double PIO4 = PI/4.; 00682 00683 double bessel_j0(double x) 00684 { 00685 double w, z, p, q, xn; 00686 00687 DEBUG_ENTRY( "bessel_j0()" ); 00688 00689 if( x < 0 ) 00690 x = -x; 00691 00692 if( x <= 5.0 ) 00693 { 00694 z = x * x; 00695 if( x < 1.0e-5 ) 00696 return 1.0 - z/4.0; 00697 00698 p = (z - DR1) * (z - DR2); 00699 p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 ); 00700 return p; 00701 } 00702 00703 w = 5.0/x; 00704 q = 25.0/(x*x); 00705 p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 ); 00706 q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 ); 00707 xn = x - PIO4; 00708 p = p * cos(xn) - w * q * sin(xn); 00709 return p * SQ2OPI / sqrt(x); 00710 } 00711 00712 /* y0() 2 */ 00713 /* Bessel function of second kind, order zero */ 00714 00715 /* Rational approximation coefficients YP[], YQ[] are used here. 00716 * The function computed is y0(x) - 2 * log(x) * j0(x) / PI, 00717 * whose value at x = 0 is 2 * ( log(0.5) + EULER ) / PI 00718 * = 0.073804295108687225. 00719 */ 00720 00721 double bessel_y0(double x) 00722 { 00723 double w, z, p, q, xn; 00724 00725 DEBUG_ENTRY( "bessel_y0()" ); 00726 00727 if( x <= 5.0 ) 00728 { 00729 if( x <= 0.0 ) 00730 { 00731 fprintf( ioQQQ, "bessel_y0: domain error\n" ); 00732 cdEXIT(EXIT_FAILURE); 00733 } 00734 z = x * x; 00735 w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 ); 00736 w += TWOOPI * log(x) * bessel_j0(x); 00737 return w; 00738 } 00739 00740 w = 5.0/x; 00741 z = 25.0 / (x * x); 00742 p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 ); 00743 q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 ); 00744 xn = x - PIO4; 00745 p = p * sin(xn) + w * q * cos(xn); 00746 return p * SQ2OPI / sqrt(x); 00747 } 00748 00749 /* j1.c 00750 * 00751 * Bessel function of order one 00752 * 00753 * 00754 * 00755 * SYNOPSIS: 00756 * 00757 * double x, y, j1(); 00758 * 00759 * y = j1( x ); 00760 * 00761 * 00762 * 00763 * DESCRIPTION: 00764 * 00765 * Returns Bessel function of order one of the argument. 00766 * 00767 * The domain is divided into the intervals [0, 8] and 00768 * (8, infinity). In the first interval a 24 term Chebyshev 00769 * expansion is used. In the second, the asymptotic 00770 * trigonometric representation is employed using two 00771 * rational functions of degree 5/5. 00772 * 00773 * 00774 * 00775 * ACCURACY: 00776 * 00777 * Absolute error: 00778 * arithmetic domain # trials peak rms 00779 * DEC 0, 30 10000 4.0e-17 1.1e-17 00780 * IEEE 0, 30 30000 2.6e-16 1.1e-16 00781 * 00782 * 00783 */ 00784 /* y1.c 00785 * 00786 * Bessel function of second kind of order one 00787 * 00788 * 00789 * 00790 * SYNOPSIS: 00791 * 00792 * double x, y, y1(); 00793 * 00794 * y = y1( x ); 00795 * 00796 * 00797 * 00798 * DESCRIPTION: 00799 * 00800 * Returns Bessel function of the second kind of order one 00801 * of the argument. 00802 * 00803 * The domain is divided into the intervals [0, 8] and 00804 * (8, infinity). In the first interval a 25 term Chebyshev 00805 * expansion is used, and a call to j1() is required. 00806 * In the second, the asymptotic trigonometric representation 00807 * is employed using two rational functions of degree 5/5. 00808 * 00809 * 00810 * 00811 * ACCURACY: 00812 * 00813 * Absolute error: 00814 * arithmetic domain # trials peak rms 00815 * DEC 0, 30 10000 8.6e-17 1.3e-17 00816 * IEEE 0, 30 30000 1.0e-15 1.3e-16 00817 * 00818 * (error criterion relative when |y1| > 1). 00819 * 00820 */ 00821 00822 /* 00823 Cephes Math Library Release 2.8: June, 2000 00824 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 00825 */ 00826 00827 static const double b1_RP[4] = { 00828 -8.99971225705559398224e8, 00829 4.52228297998194034323e11, 00830 -7.27494245221818276015e13, 00831 3.68295732863852883286e15, 00832 }; 00833 00834 static const double b1_RQ[8] = { 00835 /* 1.00000000000000000000E0,*/ 00836 6.20836478118054335476e2, 00837 2.56987256757748830383e5, 00838 8.35146791431949253037e7, 00839 2.21511595479792499675e10, 00840 4.74914122079991414898e12, 00841 7.84369607876235854894e14, 00842 8.95222336184627338078e16, 00843 5.32278620332680085395e18, 00844 }; 00845 00846 static const double b1_PP[7] = { 00847 7.62125616208173112003e-4, 00848 7.31397056940917570436e-2, 00849 1.12719608129684925192e0, 00850 5.11207951146807644818e0, 00851 8.42404590141772420927e0, 00852 5.21451598682361504063e0, 00853 1.00000000000000000254e0, 00854 }; 00855 00856 static const double b1_PQ[7] = { 00857 5.71323128072548699714e-4, 00858 6.88455908754495404082e-2, 00859 1.10514232634061696926e0, 00860 5.07386386128601488557e0, 00861 8.39985554327604159757e0, 00862 5.20982848682361821619e0, 00863 9.99999999999999997461e-1, 00864 }; 00865 00866 static const double b1_QP[8] = { 00867 5.10862594750176621635e-2, 00868 4.98213872951233449420e0, 00869 7.58238284132545283818e1, 00870 3.66779609360150777800e2, 00871 7.10856304998926107277e2, 00872 5.97489612400613639965e2, 00873 2.11688757100572135698e2, 00874 2.52070205858023719784e1, 00875 }; 00876 00877 static const double b1_QQ[7] = { 00878 /* 1.00000000000000000000e0,*/ 00879 7.42373277035675149943e1, 00880 1.05644886038262816351e3, 00881 4.98641058337653607651e3, 00882 9.56231892404756170795e3, 00883 7.99704160447350683650e3, 00884 2.82619278517639096600e3, 00885 3.36093607810698293419e2, 00886 }; 00887 00888 static const double b1_YP[6] = { 00889 1.26320474790178026440e9, 00890 -6.47355876379160291031e11, 00891 1.14509511541823727583e14, 00892 -8.12770255501325109621e15, 00893 2.02439475713594898196e17, 00894 -7.78877196265950026825e17, 00895 }; 00896 00897 static const double b1_YQ[8] = { 00898 /* 1.00000000000000000000E0,*/ 00899 5.94301592346128195359E2, 00900 2.35564092943068577943E5, 00901 7.34811944459721705660E7, 00902 1.87601316108706159478E10, 00903 3.88231277496238566008E12, 00904 6.20557727146953693363E14, 00905 6.87141087355300489866E16, 00906 3.97270608116560655612E18, 00907 }; 00908 00909 static const double Z1 = 1.46819706421238932572E1; 00910 static const double Z2 = 4.92184563216946036703E1; 00911 00912 static const double THPIO4 = 3.*PI/4.; 00913 00914 double bessel_j1(double x) 00915 { 00916 double w, z, p, q, xn; 00917 00918 DEBUG_ENTRY( "bessel_j1()" ); 00919 00920 w = x; 00921 if( x < 0 ) 00922 w = -x; 00923 00924 if( w <= 5.0 ) 00925 { 00926 z = x * x; 00927 w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 ); 00928 w = w * x * (z - Z1) * (z - Z2); 00929 return w; 00930 } 00931 00932 w = 5.0/x; 00933 z = w * w; 00934 p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 ); 00935 q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 ); 00936 xn = x - THPIO4; 00937 p = p * cos(xn) - w * q * sin(xn); 00938 return p * SQ2OPI / sqrt(x); 00939 } 00940 00941 double bessel_y1(double x) 00942 { 00943 double w, z, p, q, xn; 00944 00945 DEBUG_ENTRY( "bessel_y1()" ); 00946 00947 if( x <= 5.0 ) 00948 { 00949 if( x <= 0.0 ) 00950 { 00951 fprintf( ioQQQ, "bessel_y1: domain error\n" ); 00952 cdEXIT(EXIT_FAILURE); 00953 } 00954 z = x * x; 00955 w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 )); 00956 w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x ); 00957 return w; 00958 } 00959 00960 w = 5.0/x; 00961 z = w * w; 00962 p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 ); 00963 q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 ); 00964 xn = x - THPIO4; 00965 p = p * sin(xn) + w * q * cos(xn); 00966 return p * SQ2OPI / sqrt(x); 00967 } 00968 00969 /* jn.c 00970 * 00971 * Bessel function of integer order 00972 * 00973 * 00974 * 00975 * SYNOPSIS: 00976 * 00977 * int n; 00978 * double x, y, jn(); 00979 * 00980 * y = jn( n, x ); 00981 * 00982 * 00983 * 00984 * DESCRIPTION: 00985 * 00986 * Returns Bessel function of order n, where n is a 00987 * (possibly negative) integer. 00988 * 00989 * The ratio of jn(x) to j0(x) is computed by backward 00990 * recurrence. First the ratio jn/jn-1 is found by a 00991 * continued fraction expansion. Then the recurrence 00992 * relating successive orders is applied until j0 or j1 is 00993 * reached. 00994 * 00995 * If n = 0 or 1 the routine for j0 or j1 is called 00996 * directly. 00997 * 00998 * 00999 * 01000 * ACCURACY: 01001 * 01002 * Absolute error: 01003 * arithmetic range # trials peak rms 01004 * DEC 0, 30 5500 6.9e-17 9.3e-18 01005 * IEEE 0, 30 5000 4.4e-16 7.9e-17 01006 * 01007 * 01008 * Not suitable for large n or x. Use jv() instead. 01009 * 01010 */ 01011 01012 /* jn.c 01013 Cephes Math Library Release 2.8: June, 2000 01014 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01015 */ 01016 01017 double bessel_jn(int n, double x) 01018 { 01019 double pkm2, pkm1, pk, xk, r, ans; 01020 int k, sign; 01021 01022 DEBUG_ENTRY( "bessel_jn()" ); 01023 01024 if( n < 0 ) 01025 { 01026 n = -n; 01027 if( (n & 1) == 0 ) /* -1**n */ 01028 sign = 1; 01029 else 01030 sign = -1; 01031 } 01032 else 01033 sign = 1; 01034 01035 if( x < 0.0 ) 01036 { 01037 if( n & 1 ) 01038 sign = -sign; 01039 x = -x; 01040 } 01041 01042 if( n == 0 ) 01043 { 01044 return sign * bessel_j0(x); 01045 } 01046 if( n == 1 ) 01047 { 01048 return sign * bessel_j1(x); 01049 } 01050 if( n == 2 ) 01051 { 01052 return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x)); 01053 } 01054 01055 if( x < DBL_EPSILON ) 01056 { 01057 return 0.0; 01058 } 01059 01060 /* continued fraction */ 01061 k = 52; 01062 01063 pk = 2 * (n + k); 01064 ans = pk; 01065 xk = x * x; 01066 01067 do 01068 { 01069 pk -= 2.0; 01070 ans = pk - (xk/ans); 01071 } 01072 while( --k > 0 ); 01073 ans = x/ans; 01074 01075 /* backward recurrence */ 01076 pk = 1.0; 01077 pkm1 = 1.0/ans; 01078 k = n-1; 01079 r = 2 * k; 01080 01081 do 01082 { 01083 pkm2 = (pkm1 * r - pk * x) / x; 01084 pk = pkm1; 01085 pkm1 = pkm2; 01086 r -= 2.0; 01087 } 01088 while( --k > 0 ); 01089 01090 if( fabs(pk) > fabs(pkm1) ) 01091 ans = bessel_j1(x)/pk; 01092 else 01093 ans = bessel_j0(x)/pkm1; 01094 return sign * ans; 01095 } 01096 01097 /* yn.c 01098 * 01099 * Bessel function of second kind of integer order 01100 * 01101 * 01102 * 01103 * SYNOPSIS: 01104 * 01105 * double x, y, yn(); 01106 * int n; 01107 * 01108 * y = yn( n, x ); 01109 * 01110 * 01111 * 01112 * DESCRIPTION: 01113 * 01114 * Returns Bessel function of order n, where n is a 01115 * (possibly negative) integer. 01116 * 01117 * The function is evaluated by forward recurrence on 01118 * n, starting with values computed by the routines 01119 * y0() and y1(). 01120 * 01121 * If n = 0 or 1 the routine for y0 or y1 is called 01122 * directly. 01123 * 01124 * 01125 * 01126 * ACCURACY: 01127 * 01128 * 01129 * Absolute error, except relative 01130 * when y > 1: 01131 * arithmetic domain # trials peak rms 01132 * DEC 0, 30 2200 2.9e-16 5.3e-17 01133 * IEEE 0, 30 30000 3.4e-15 4.3e-16 01134 * 01135 * 01136 * ERROR MESSAGES: 01137 * 01138 * message condition value returned 01139 * yn singularity x = 0 MAXNUM 01140 * yn overflow MAXNUM 01141 * 01142 * Spot checked against tables for x, n between 0 and 100. 01143 * 01144 */ 01145 01146 /* 01147 Cephes Math Library Release 2.8: June, 2000 01148 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01149 */ 01150 01151 double bessel_yn(int n, double x) 01152 { 01153 double an, anm1, anm2, r; 01154 int k, sign; 01155 01156 DEBUG_ENTRY( "bessel_yn()" ); 01157 01158 if( n < 0 ) 01159 { 01160 n = -n; 01161 if( (n & 1) == 0 ) /* -1**n */ 01162 sign = 1; 01163 else 01164 sign = -1; 01165 } 01166 else 01167 sign = 1; 01168 01169 if( n == 0 ) 01170 { 01171 return sign * bessel_y0(x); 01172 } 01173 if( n == 1 ) 01174 { 01175 return sign * bessel_y1(x); 01176 } 01177 01178 /* test for overflow */ 01179 if( x <= 0.0 ) 01180 { 01181 fprintf( ioQQQ, "bessel_yn: domain error\n" ); 01182 cdEXIT(EXIT_FAILURE); 01183 } 01184 01185 /* forward recurrence on n */ 01186 anm2 = bessel_y0(x); 01187 anm1 = bessel_y1(x); 01188 k = 1; 01189 r = 2 * k; 01190 do 01191 { 01192 an = r * anm1 / x - anm2; 01193 anm2 = anm1; 01194 anm1 = an; 01195 r += 2.0; 01196 ++k; 01197 } 01198 while( k < n ); 01199 return sign * an; 01200 } 01201 01202 /* k0.c 01203 * 01204 * Modified Bessel function, third kind, order zero 01205 * 01206 * 01207 * 01208 * SYNOPSIS: 01209 * 01210 * double x, y, k0(); 01211 * 01212 * y = k0( x ); 01213 * 01214 * 01215 * 01216 * DESCRIPTION: 01217 * 01218 * Returns modified Bessel function of the third kind 01219 * of order zero of the argument. 01220 * 01221 * The range is partitioned into the two intervals [0,8] and 01222 * (8, infinity). Chebyshev polynomial expansions are employed 01223 * in each interval. 01224 * 01225 * 01226 * 01227 * ACCURACY: 01228 * 01229 * Tested at 2000 random points between 0 and 8. Peak absolute 01230 * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. 01231 * Relative error: 01232 * arithmetic domain # trials peak rms 01233 * DEC 0, 30 3100 1.3e-16 2.1e-17 01234 * IEEE 0, 30 30000 1.2e-15 1.6e-16 01235 * 01236 * ERROR MESSAGES: 01237 * 01238 * message condition value returned 01239 * K0 domain x <= 0 MAXNUM 01240 * 01241 */ 01242 /* k0e() 01243 * 01244 * Modified Bessel function, third kind, order zero, 01245 * exponentially scaled 01246 * 01247 * 01248 * 01249 * SYNOPSIS: 01250 * 01251 * double x, y, k0e(); 01252 * 01253 * y = k0e( x ); 01254 * 01255 * 01256 * 01257 * DESCRIPTION: 01258 * 01259 * Returns exponentially scaled modified Bessel function 01260 * of the third kind of order zero of the argument. 01261 * 01262 * 01263 * 01264 * ACCURACY: 01265 * 01266 * Relative error: 01267 * arithmetic domain # trials peak rms 01268 * IEEE 0, 30 30000 1.4e-15 1.4e-16 01269 * See k0(). 01270 * 01271 */ 01272 01273 /* 01274 Cephes Math Library Release 2.8: June, 2000 01275 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01276 */ 01277 01278 /* Chebyshev coefficients for K0(x) + log(x/2) I0(x) 01279 * in the interval [0,2]. The odd order coefficients are all 01280 * zero; only the even order coefficients are listed. 01281 * 01282 * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EULER. 01283 */ 01284 01285 static const double k0_A[] = 01286 { 01287 1.37446543561352307156e-16, 01288 4.25981614279661018399e-14, 01289 1.03496952576338420167e-11, 01290 1.90451637722020886025e-9, 01291 2.53479107902614945675e-7, 01292 2.28621210311945178607e-5, 01293 1.26461541144692592338e-3, 01294 3.59799365153615016266e-2, 01295 3.44289899924628486886e-1, 01296 -5.35327393233902768720e-1 01297 }; 01298 01299 /* Chebyshev coefficients for exp(x) sqrt(x) K0(x) 01300 * in the inverted interval [2,infinity]. 01301 * 01302 * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2). 01303 */ 01304 01305 static const double k0_B[] = { 01306 5.30043377268626276149e-18, 01307 -1.64758043015242134646e-17, 01308 5.21039150503902756861e-17, 01309 -1.67823109680541210385e-16, 01310 5.51205597852431940784e-16, 01311 -1.84859337734377901440e-15, 01312 6.34007647740507060557e-15, 01313 -2.22751332699166985548e-14, 01314 8.03289077536357521100e-14, 01315 -2.98009692317273043925e-13, 01316 1.14034058820847496303e-12, 01317 -4.51459788337394416547e-12, 01318 1.85594911495471785253e-11, 01319 -7.95748924447710747776e-11, 01320 3.57739728140030116597e-10, 01321 -1.69753450938905987466e-9, 01322 8.57403401741422608519e-9, 01323 -4.66048989768794782956e-8, 01324 2.76681363944501510342e-7, 01325 -1.83175552271911948767e-6, 01326 1.39498137188764993662e-5, 01327 -1.28495495816278026384e-4, 01328 1.56988388573005337491e-3, 01329 -3.14481013119645005427e-2, 01330 2.44030308206595545468e0 01331 }; 01332 01333 double bessel_k0(double x) 01334 { 01335 double y, z; 01336 01337 DEBUG_ENTRY( "bessel_k0()" ); 01338 01339 if( x <= 0.0 ) 01340 { 01341 fprintf( ioQQQ, "bessel_k0: domain error\n" ); 01342 cdEXIT(EXIT_FAILURE); 01343 } 01344 01345 if( x <= 2.0 ) 01346 { 01347 y = x * x - 2.0; 01348 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x); 01349 return y; 01350 } 01351 z = 8.0/x - 2.0; 01352 y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x); 01353 return y; 01354 } 01355 01356 double bessel_k0_scaled(double x) 01357 { 01358 double y; 01359 01360 DEBUG_ENTRY( "bessel_k0_scaled()" ); 01361 01362 if( x <= 0.0 ) 01363 { 01364 fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" ); 01365 cdEXIT(EXIT_FAILURE); 01366 } 01367 01368 if( x <= 2.0 ) 01369 { 01370 y = x * x - 2.0; 01371 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x); 01372 return y * exp(x); 01373 } 01374 return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x); 01375 } 01376 01377 /* k1.c 01378 * 01379 * Modified Bessel function, third kind, order one 01380 * 01381 * 01382 * 01383 * SYNOPSIS: 01384 * 01385 * double x, y, k1(); 01386 * 01387 * y = k1( x ); 01388 * 01389 * 01390 * 01391 * DESCRIPTION: 01392 * 01393 * Computes the modified Bessel function of the third kind 01394 * of order one of the argument. 01395 * 01396 * The range is partitioned into the two intervals [0,2] and 01397 * (2, infinity). Chebyshev polynomial expansions are employed 01398 * in each interval. 01399 * 01400 * 01401 * 01402 * ACCURACY: 01403 * 01404 * Relative error: 01405 * arithmetic domain # trials peak rms 01406 * DEC 0, 30 3300 8.9e-17 2.2e-17 01407 * IEEE 0, 30 30000 1.2e-15 1.6e-16 01408 * 01409 * ERROR MESSAGES: 01410 * 01411 * message condition value returned 01412 * k1 domain x <= 0 MAXNUM 01413 * 01414 */ 01415 /* k1e.c 01416 * 01417 * Modified Bessel function, third kind, order one, 01418 * exponentially scaled 01419 * 01420 * 01421 * 01422 * SYNOPSIS: 01423 * 01424 * double x, y, k1e(); 01425 * 01426 * y = k1e( x ); 01427 * 01428 * 01429 * 01430 * DESCRIPTION: 01431 * 01432 * Returns exponentially scaled modified Bessel function 01433 * of the third kind of order one of the argument: 01434 * 01435 * k1e(x) = exp(x) * k1(x). 01436 * 01437 * 01438 * 01439 * ACCURACY: 01440 * 01441 * Relative error: 01442 * arithmetic domain # trials peak rms 01443 * IEEE 0, 30 30000 7.8e-16 1.2e-16 01444 * See k1(). 01445 * 01446 */ 01447 01448 /* 01449 Cephes Math Library Release 2.8: June, 2000 01450 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01451 */ 01452 01453 /* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x)) 01454 * in the interval [0,2]. 01455 * 01456 * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1. 01457 */ 01458 01459 static const double k1_A[] = 01460 { 01461 -7.02386347938628759343e-18, 01462 -2.42744985051936593393e-15, 01463 -6.66690169419932900609e-13, 01464 -1.41148839263352776110e-10, 01465 -2.21338763073472585583e-8, 01466 -2.43340614156596823496e-6, 01467 -1.73028895751305206302e-4, 01468 -6.97572385963986435018e-3, 01469 -1.22611180822657148235e-1, 01470 -3.53155960776544875667e-1, 01471 1.52530022733894777053e0 01472 }; 01473 01474 /* Chebyshev coefficients for exp(x) sqrt(x) K1(x) 01475 * in the interval [2,infinity]. 01476 * 01477 * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2). 01478 */ 01479 01480 static const double k1_B[] = 01481 { 01482 -5.75674448366501715755e-18, 01483 1.79405087314755922667e-17, 01484 -5.68946255844285935196e-17, 01485 1.83809354436663880070e-16, 01486 -6.05704724837331885336e-16, 01487 2.03870316562433424052e-15, 01488 -7.01983709041831346144e-15, 01489 2.47715442448130437068e-14, 01490 -8.97670518232499435011e-14, 01491 3.34841966607842919884e-13, 01492 -1.28917396095102890680e-12, 01493 5.13963967348173025100e-12, 01494 -2.12996783842756842877e-11, 01495 9.21831518760500529508e-11, 01496 -4.19035475934189648750e-10, 01497 2.01504975519703286596e-9, 01498 -1.03457624656780970260e-8, 01499 5.74108412545004946722e-8, 01500 -3.50196060308781257119e-7, 01501 2.40648494783721712015e-6, 01502 -1.93619797416608296024e-5, 01503 1.95215518471351631108e-4, 01504 -2.85781685962277938680e-3, 01505 1.03923736576817238437e-1, 01506 2.72062619048444266945e0 01507 }; 01508 01509 double bessel_k1(double x) 01510 { 01511 double y, z; 01512 01513 DEBUG_ENTRY( "bessel_k1()" ); 01514 01515 z = 0.5 * x; 01516 if( z <= 0.0 ) 01517 { 01518 fprintf( ioQQQ, "bessel_k1: domain error\n" ); 01519 cdEXIT(EXIT_FAILURE); 01520 } 01521 01522 if( x <= 2.0 ) 01523 { 01524 y = x * x - 2.0; 01525 y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x; 01526 return y; 01527 } 01528 return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x); 01529 } 01530 01531 double bessel_k1_scaled(double x) 01532 { 01533 double y; 01534 01535 DEBUG_ENTRY( "bessel_k1_scaled()" ); 01536 01537 if( x <= 0.0 ) 01538 { 01539 fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" ); 01540 cdEXIT(EXIT_FAILURE); 01541 } 01542 01543 if( x <= 2.0 ) 01544 { 01545 y = x * x - 2.0; 01546 y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x; 01547 return y * exp(x); 01548 } 01549 return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x); 01550 } 01551 01552 /* i0.c 01553 * 01554 * Modified Bessel function of order zero 01555 * 01556 * 01557 * 01558 * SYNOPSIS: 01559 * 01560 * double x, y, i0(); 01561 * 01562 * y = i0( x ); 01563 * 01564 * 01565 * 01566 * DESCRIPTION: 01567 * 01568 * Returns modified Bessel function of order zero of the 01569 * argument. 01570 * 01571 * The function is defined as i0(x) = j0( ix ). 01572 * 01573 * The range is partitioned into the two intervals [0,8] and 01574 * (8, infinity). Chebyshev polynomial expansions are employed 01575 * in each interval. 01576 * 01577 * 01578 * 01579 * ACCURACY: 01580 * 01581 * Relative error: 01582 * arithmetic domain # trials peak rms 01583 * DEC 0,30 6000 8.2e-17 1.9e-17 01584 * IEEE 0,30 30000 5.8e-16 1.4e-16 01585 * 01586 */ 01587 /* i0e.c 01588 * 01589 * Modified Bessel function of order zero, 01590 * exponentially scaled 01591 * 01592 * 01593 * 01594 * SYNOPSIS: 01595 * 01596 * double x, y, i0e(); 01597 * 01598 * y = i0e( x ); 01599 * 01600 * 01601 * 01602 * DESCRIPTION: 01603 * 01604 * Returns exponentially scaled modified Bessel function 01605 * of order zero of the argument. 01606 * 01607 * The function is defined as i0e(x) = exp(-|x|) j0( ix ). 01608 * 01609 * 01610 * 01611 * ACCURACY: 01612 * 01613 * Relative error: 01614 * arithmetic domain # trials peak rms 01615 * IEEE 0,30 30000 5.4e-16 1.2e-16 01616 * See i0(). 01617 * 01618 */ 01619 01620 /* 01621 Cephes Math Library Release 2.8: June, 2000 01622 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01623 */ 01624 01625 /* Chebyshev coefficients for exp(-x) I0(x) 01626 * in the interval [0,8]. 01627 * 01628 * lim(x->0){ exp(-x) I0(x) } = 1. 01629 */ 01630 01631 static const double i0_A[] = 01632 { 01633 -4.41534164647933937950e-18, 01634 3.33079451882223809783e-17, 01635 -2.43127984654795469359e-16, 01636 1.71539128555513303061e-15, 01637 -1.16853328779934516808e-14, 01638 7.67618549860493561688e-14, 01639 -4.85644678311192946090e-13, 01640 2.95505266312963983461e-12, 01641 -1.72682629144155570723e-11, 01642 9.67580903537323691224e-11, 01643 -5.18979560163526290666e-10, 01644 2.65982372468238665035e-9, 01645 -1.30002500998624804212e-8, 01646 6.04699502254191894932e-8, 01647 -2.67079385394061173391e-7, 01648 1.11738753912010371815e-6, 01649 -4.41673835845875056359e-6, 01650 1.64484480707288970893e-5, 01651 -5.75419501008210370398e-5, 01652 1.88502885095841655729e-4, 01653 -5.76375574538582365885e-4, 01654 1.63947561694133579842e-3, 01655 -4.32430999505057594430e-3, 01656 1.05464603945949983183e-2, 01657 -2.37374148058994688156e-2, 01658 4.93052842396707084878e-2, 01659 -9.49010970480476444210e-2, 01660 1.71620901522208775349e-1, 01661 -3.04682672343198398683e-1, 01662 6.76795274409476084995e-1 01663 }; 01664 01665 /* Chebyshev coefficients for exp(-x) sqrt(x) I0(x) 01666 * in the inverted interval [8,infinity]. 01667 * 01668 * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). 01669 */ 01670 01671 static const double i0_B[] = 01672 { 01673 -7.23318048787475395456e-18, 01674 -4.83050448594418207126e-18, 01675 4.46562142029675999901e-17, 01676 3.46122286769746109310e-17, 01677 -2.82762398051658348494e-16, 01678 -3.42548561967721913462e-16, 01679 1.77256013305652638360e-15, 01680 3.81168066935262242075e-15, 01681 -9.55484669882830764870e-15, 01682 -4.15056934728722208663e-14, 01683 1.54008621752140982691e-14, 01684 3.85277838274214270114e-13, 01685 7.18012445138366623367e-13, 01686 -1.79417853150680611778e-12, 01687 -1.32158118404477131188e-11, 01688 -3.14991652796324136454e-11, 01689 1.18891471078464383424e-11, 01690 4.94060238822496958910e-10, 01691 3.39623202570838634515e-9, 01692 2.26666899049817806459e-8, 01693 2.04891858946906374183e-7, 01694 2.89137052083475648297e-6, 01695 6.88975834691682398426e-5, 01696 3.36911647825569408990e-3, 01697 8.04490411014108831608e-1 01698 }; 01699 01700 double bessel_i0(double x) 01701 { 01702 double y; 01703 01704 DEBUG_ENTRY( "bessel_i0()" ); 01705 01706 if( x < 0 ) 01707 x = -x; 01708 01709 if( x <= 8.0 ) 01710 { 01711 y = (x/2.0) - 2.0; 01712 return exp(x) * chbevl( y, i0_A, 30 ); 01713 } 01714 return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x); 01715 } 01716 01717 double bessel_i0_scaled(double x) 01718 { 01719 double y; 01720 01721 DEBUG_ENTRY( "bessel_i0_scaled()" ); 01722 01723 if( x < 0 ) 01724 x = -x; 01725 01726 if( x <= 8.0 ) 01727 { 01728 y = (x/2.0) - 2.0; 01729 return chbevl( y, i0_A, 30 ); 01730 } 01731 return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x); 01732 } 01733 01734 /* i1.c 01735 * 01736 * Modified Bessel function of order one 01737 * 01738 * 01739 * 01740 * SYNOPSIS: 01741 * 01742 * double x, y, i1(); 01743 * 01744 * y = i1( x ); 01745 * 01746 * 01747 * 01748 * DESCRIPTION: 01749 * 01750 * Returns modified Bessel function of order one of the 01751 * argument. 01752 * 01753 * The function is defined as i1(x) = -i j1( ix ). 01754 * 01755 * The range is partitioned into the two intervals [0,8] and 01756 * (8, infinity). Chebyshev polynomial expansions are employed 01757 * in each interval. 01758 * 01759 * 01760 * 01761 * ACCURACY: 01762 * 01763 * Relative error: 01764 * arithmetic domain # trials peak rms 01765 * DEC 0, 30 3400 1.2e-16 2.3e-17 01766 * IEEE 0, 30 30000 1.9e-15 2.1e-16 01767 * 01768 * 01769 */ 01770 /* i1e.c 01771 * 01772 * Modified Bessel function of order one, 01773 * exponentially scaled 01774 * 01775 * 01776 * 01777 * SYNOPSIS: 01778 * 01779 * double x, y, i1e(); 01780 * 01781 * y = i1e( x ); 01782 * 01783 * 01784 * 01785 * DESCRIPTION: 01786 * 01787 * Returns exponentially scaled modified Bessel function 01788 * of order one of the argument. 01789 * 01790 * The function is defined as i1(x) = -i exp(-|x|) j1( ix ). 01791 * 01792 * 01793 * 01794 * ACCURACY: 01795 * 01796 * Relative error: 01797 * arithmetic domain # trials peak rms 01798 * IEEE 0, 30 30000 2.0e-15 2.0e-16 01799 * See i1(). 01800 * 01801 */ 01802 01803 /* 01804 Cephes Math Library Release 2.8: June, 2000 01805 Copyright 1985, 1987, 2000 by Stephen L. Moshier 01806 */ 01807 01808 /* Chebyshev coefficients for exp(-x) I1(x) / x 01809 * in the interval [0,8]. 01810 * 01811 * lim(x->0){ exp(-x) I1(x) / x } = 1/2. 01812 */ 01813 01814 static double i1_A[] = 01815 { 01816 2.77791411276104639959e-18, 01817 -2.11142121435816608115e-17, 01818 1.55363195773620046921e-16, 01819 -1.10559694773538630805e-15, 01820 7.60068429473540693410e-15, 01821 -5.04218550472791168711e-14, 01822 3.22379336594557470981e-13, 01823 -1.98397439776494371520e-12, 01824 1.17361862988909016308e-11, 01825 -6.66348972350202774223e-11, 01826 3.62559028155211703701e-10, 01827 -1.88724975172282928790e-9, 01828 9.38153738649577178388e-9, 01829 -4.44505912879632808065e-8, 01830 2.00329475355213526229e-7, 01831 -8.56872026469545474066e-7, 01832 3.47025130813767847674e-6, 01833 -1.32731636560394358279e-5, 01834 4.78156510755005422638e-5, 01835 -1.61760815825896745588e-4, 01836 5.12285956168575772895e-4, 01837 -1.51357245063125314899e-3, 01838 4.15642294431288815669e-3, 01839 -1.05640848946261981558e-2, 01840 2.47264490306265168283e-2, 01841 -5.29459812080949914269e-2, 01842 1.02643658689847095384e-1, 01843 -1.76416518357834055153e-1, 01844 2.52587186443633654823e-1 01845 }; 01846 01847 /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x) 01848 * in the inverted interval [8,infinity]. 01849 * 01850 * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi). 01851 */ 01852 01853 static double i1_B[] = 01854 { 01855 7.51729631084210481353e-18, 01856 4.41434832307170791151e-18, 01857 -4.65030536848935832153e-17, 01858 -3.20952592199342395980e-17, 01859 2.96262899764595013876e-16, 01860 3.30820231092092828324e-16, 01861 -1.88035477551078244854e-15, 01862 -3.81440307243700780478e-15, 01863 1.04202769841288027642e-14, 01864 4.27244001671195135429e-14, 01865 -2.10154184277266431302e-14, 01866 -4.08355111109219731823e-13, 01867 -7.19855177624590851209e-13, 01868 2.03562854414708950722e-12, 01869 1.41258074366137813316e-11, 01870 3.25260358301548823856e-11, 01871 -1.89749581235054123450e-11, 01872 -5.58974346219658380687e-10, 01873 -3.83538038596423702205e-9, 01874 -2.63146884688951950684e-8, 01875 -2.51223623787020892529e-7, 01876 -3.88256480887769039346e-6, 01877 -1.10588938762623716291e-4, 01878 -9.76109749136146840777e-3, 01879 7.78576235018280120474e-1 01880 }; 01881 01882 double bessel_i1(double x) 01883 { 01884 double y, z; 01885 01886 DEBUG_ENTRY( "bessel_i1()" ); 01887 01888 z = fabs(x); 01889 if( z <= 8.0 ) 01890 { 01891 y = (z/2.0) - 2.0; 01892 z = chbevl( y, i1_A, 29 ) * z * exp(z); 01893 } 01894 else 01895 { 01896 z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z); 01897 } 01898 if( x < 0.0 ) 01899 z = -z; 01900 return z; 01901 } 01902 01903 double bessel_i1_scaled(double x) 01904 { 01905 double y, z; 01906 01907 DEBUG_ENTRY( "bessel_i1_scaled()" ); 01908 01909 z = fabs(x); 01910 if( z <= 8.0 ) 01911 { 01912 y = (z/2.0) - 2.0; 01913 z = chbevl( y, i1_A, 29 ) * z; 01914 } 01915 else 01916 { 01917 z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z); 01918 } 01919 if( x < 0.0 ) 01920 z = -z; 01921 return z; 01922 } 01923 01924 /* ellpk.c 01925 * 01926 * Complete elliptic integral of the first kind 01927 * 01928 * 01929 * 01930 * SYNOPSIS: 01931 * 01932 * double m1, y, ellpk(); 01933 * 01934 * y = ellpk( m1 ); 01935 * 01936 * 01937 * 01938 * DESCRIPTION: 01939 * 01940 * Approximates the integral 01941 * 01942 * 01943 * 01944 * pi/2 01945 * - 01946 * | | 01947 * | dt 01948 * K(m) = | ------------------ 01949 * | 2 01950 * | | sqrt( 1 - m sin t ) 01951 * - 01952 * 0 01953 * 01954 * where m = 1 - m1, using the approximation 01955 * 01956 * P(x) - log x Q(x). 01957 * 01958 * The argument m1 is used rather than m so that the logarithmic 01959 * singularity at m = 1 will be shifted to the origin; this 01960 * preserves maximum accuracy. 01961 * 01962 * K(0) = pi/2. 01963 * 01964 * ACCURACY: 01965 * 01966 * Relative error: 01967 * arithmetic domain # trials peak rms 01968 * DEC 0,1 16000 3.5e-17 1.1e-17 01969 * IEEE 0,1 30000 2.5e-16 6.8e-17 01970 * 01971 * ERROR MESSAGES: 01972 * 01973 * message condition value returned 01974 * ellpk domain x<0, x>1 0.0 01975 * 01976 */ 01977 01978 /* 01979 Cephes Math Library, Release 2.8: June, 2000 01980 Copyright 1984, 1987, 2000 by Stephen L. Moshier 01981 */ 01982 01983 static const double elk_P[] = 01984 { 01985 1.37982864606273237150e-4, 01986 2.28025724005875567385e-3, 01987 7.97404013220415179367e-3, 01988 9.85821379021226008714e-3, 01989 6.87489687449949877925e-3, 01990 6.18901033637687613229e-3, 01991 8.79078273952743772254e-3, 01992 1.49380448916805252718e-2, 01993 3.08851465246711995998e-2, 01994 9.65735902811690126535e-2, 01995 1.38629436111989062502e0 01996 }; 01997 01998 static const double elk_Q[] = 01999 { 02000 2.94078955048598507511e-5, 02001 9.14184723865917226571e-4, 02002 5.94058303753167793257e-3, 02003 1.54850516649762399335e-2, 02004 2.39089602715924892727e-2, 02005 3.01204715227604046988e-2, 02006 3.73774314173823228969e-2, 02007 4.88280347570998239232e-2, 02008 7.03124996963957469739e-2, 02009 1.24999999999870820058e-1, 02010 4.99999999999999999821e-1 02011 }; 02012 02013 static const double C1 = 1.3862943611198906188e0; /* log(4) */ 02014 02015 double ellpk(double x) 02016 { 02017 DEBUG_ENTRY( "ellpk()" ); 02018 02019 if( x < 0.0 || x > 1.0 ) 02020 { 02021 fprintf( ioQQQ, "ellpk: domain error\n" ); 02022 cdEXIT(EXIT_FAILURE); 02023 } 02024 02025 if( x > DBL_EPSILON ) 02026 { 02027 return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10); 02028 } 02029 else 02030 { 02031 if( x == 0.0 ) 02032 { 02033 fprintf( ioQQQ, "ellpk: domain error\n" ); 02034 cdEXIT(EXIT_FAILURE); 02035 } 02036 else 02037 { 02038 return C1 - 0.5 * log(x); 02039 } 02040 } 02041 } 02042 02043 /* expn.c 02044 * 02045 * Exponential integral En 02046 * 02047 * 02048 * 02049 * SYNOPSIS: 02050 * 02051 * int n; 02052 * double x, y, expn(); 02053 * 02054 * y = expn( n, x ); 02055 * 02056 * 02057 * 02058 * DESCRIPTION: 02059 * 02060 * Evaluates the exponential integral 02061 * 02062 * inf. 02063 * - 02064 * | | -xt 02065 * | e 02066 * E (x) = | ---- dt. 02067 * n | n 02068 * | | t 02069 * - 02070 * 1 02071 * 02072 * 02073 * Both n and x must be nonnegative. 02074 * 02075 * The routine employs either a power series, a continued 02076 * fraction, or an asymptotic formula depending on the 02077 * relative values of n and x. 02078 * 02079 * ACCURACY: 02080 * 02081 * Relative error: 02082 * arithmetic domain # trials peak rms 02083 * DEC 0, 30 5000 2.0e-16 4.6e-17 02084 * IEEE 0, 30 10000 1.7e-15 3.6e-16 02085 * 02086 */ 02087 02088 /* Cephes Math Library Release 2.8: June, 2000 02089 Copyright 1985, 2000 by Stephen L. Moshier */ 02090 02091 static const double MAXLOG = log(DBL_MAX); 02092 static const double BIG = 1.44115188075855872E+17; /* 2^57 */ 02093 02094 double expn(int n, double x) 02095 { 02096 double ans, r, t, yk, xk; 02097 double pk, pkm1, pkm2, qk, qkm1, qkm2; 02098 double psi, z; 02099 int i, k; 02100 02101 DEBUG_ENTRY( "expn()" ); 02102 02103 if( n < 0 || x < 0. ) 02104 { 02105 fprintf( ioQQQ, "expn: domain error\n" ); 02106 cdEXIT(EXIT_FAILURE); 02107 } 02108 02109 if( x > MAXLOG ) 02110 { 02111 return 0.0; 02112 } 02113 02114 if( x == 0.0 ) 02115 { 02116 if( n < 2 ) 02117 { 02118 fprintf( ioQQQ, "expn: domain error\n" ); 02119 cdEXIT(EXIT_FAILURE); 02120 } 02121 else 02122 { 02123 return 1.0/((double)n-1.0); 02124 } 02125 } 02126 02127 if( n == 0 ) 02128 { 02129 return exp(-x)/x; 02130 } 02131 02132 /* Expansion for large n */ 02133 if( n > 5000 ) 02134 { 02135 xk = x + n; 02136 yk = 1.0 / (xk * xk); 02137 t = n; 02138 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t); 02139 ans = yk * (ans + t * (t - 2.0 * x)); 02140 ans = yk * (ans + t); 02141 ans = (ans + 1.0) * exp( -x ) / xk; 02142 return ans; 02143 } 02144 02145 if( x <= 1.0 ) 02146 { 02147 /* Power series expansion */ 02148 psi = -EULER - log(x); 02149 for( i=1; i < n; i++ ) 02150 psi = psi + 1.0/i; 02151 02152 z = -x; 02153 xk = 0.0; 02154 yk = 1.0; 02155 pk = 1.0 - n; 02156 if( n == 1 ) 02157 ans = 0.0; 02158 else 02159 ans = 1.0/pk; 02160 do 02161 { 02162 xk += 1.0; 02163 yk *= z/xk; 02164 pk += 1.0; 02165 if( pk != 0.0 ) 02166 { 02167 ans += yk/pk; 02168 } 02169 if( ans != 0.0 ) 02170 t = fabs(yk/ans); 02171 else 02172 t = 1.0; 02173 } 02174 while( t > DBL_EPSILON ); 02175 ans = powi(z,n-1)*psi/factorial(n-1) - ans; 02176 return ans; 02177 } 02178 else 02179 { 02180 /* continued fraction */ 02181 k = 1; 02182 pkm2 = 1.0; 02183 qkm2 = x; 02184 pkm1 = 1.0; 02185 qkm1 = x + n; 02186 ans = pkm1/qkm1; 02187 02188 do 02189 { 02190 k += 1; 02191 if( is_odd(k) ) 02192 { 02193 yk = 1.0; 02194 xk = static_cast<double>(n + (k-1)/2); 02195 } 02196 else 02197 { 02198 yk = x; 02199 xk = static_cast<double>(k/2); 02200 } 02201 pk = pkm1 * yk + pkm2 * xk; 02202 qk = qkm1 * yk + qkm2 * xk; 02203 if( qk != 0 ) 02204 { 02205 r = pk/qk; 02206 t = fabs( (ans - r)/r ); 02207 ans = r; 02208 } 02209 else 02210 t = 1.0; 02211 pkm2 = pkm1; 02212 pkm1 = pk; 02213 qkm2 = qkm1; 02214 qkm1 = qk; 02215 if( fabs(pk) > BIG ) 02216 { 02217 pkm2 /= BIG; 02218 pkm1 /= BIG; 02219 qkm2 /= BIG; 02220 qkm1 /= BIG; 02221 } 02222 } 02223 while( t > DBL_EPSILON ); 02224 02225 ans *= exp( -x ); 02226 return ans; 02227 } 02228 } 02229 02230 /* polevl.c 02231 * p1evl.c 02232 * 02233 * Evaluate polynomial 02234 * 02235 * 02236 * 02237 * SYNOPSIS: 02238 * 02239 * int N; 02240 * double x, y, coef[N+1], polevl[]; 02241 * 02242 * y = polevl( x, coef, N ); 02243 * 02244 * 02245 * 02246 * DESCRIPTION: 02247 * 02248 * Evaluates polynomial of degree N: 02249 * 02250 * 2 N 02251 * y = C + C x + C x +...+ C x 02252 * 0 1 2 N 02253 * 02254 * Coefficients are stored in reverse order: 02255 * 02256 * coef[0] = C , ..., coef[N] = C . 02257 * N 0 02258 * 02259 * The function p1evl() assumes that coef[N] = 1.0 and is 02260 * omitted from the array. Its calling arguments are 02261 * otherwise the same as polevl(). 02262 * 02263 * 02264 * SPEED: 02265 * 02266 * In the interest of speed, there are no checks for out 02267 * of bounds arithmetic. This routine is used by most of 02268 * the functions in the library. Depending on available 02269 * equipment features, the user may wish to rewrite the 02270 * program in microcode or assembly language. 02271 * 02272 */ 02273 02274 /* 02275 Cephes Math Library Release 2.1: December, 1988 02276 Copyright 1984, 1987, 1988 by Stephen L. Moshier 02277 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 02278 */ 02279 02280 inline double polevl(double x, const double coef[], int N) 02281 { 02282 double ans; 02283 int i; 02284 const double *p = coef; 02285 02286 ans = *p++; 02287 i = N; 02288 02289 do 02290 ans = ans * x + *p++; 02291 while( --i ); 02292 02293 return ans; 02294 } 02295 02296 /* p1evl() */ 02297 /* N 02298 * Evaluate polynomial when coefficient of x is 1.0. 02299 * Otherwise same as polevl. 02300 */ 02301 02302 inline double p1evl(double x, const double coef[], int N) 02303 { 02304 double ans; 02305 const double *p = coef; 02306 int i; 02307 02308 ans = x + *p++; 02309 i = N-1; 02310 02311 do 02312 ans = ans * x + *p++; 02313 while( --i ); 02314 02315 return ans; 02316 } 02317 02318 /* chbevl.c 02319 * 02320 * Evaluate Chebyshev series 02321 * 02322 * 02323 * 02324 * SYNOPSIS: 02325 * 02326 * int N; 02327 * double x, y, coef[N], chebevl(); 02328 * 02329 * y = chbevl( x, coef, N ); 02330 * 02331 * 02332 * 02333 * DESCRIPTION: 02334 * 02335 * Evaluates the series 02336 * 02337 * N-1 02338 * - ' 02339 * y = > coef[i] T (x/2) 02340 * - i 02341 * i=0 02342 * 02343 * of Chebyshev polynomials Ti at argument x/2. 02344 * 02345 * Coefficients are stored in reverse order, i.e. the zero 02346 * order term is last in the array. Note N is the number of 02347 * coefficients, not the order. 02348 * 02349 * If coefficients are for the interval a to b, x must 02350 * have been transformed to x -> 2(2x - b - a)/(b-a) before 02351 * entering the routine. This maps x from (a, b) to (-1, 1), 02352 * over which the Chebyshev polynomials are defined. 02353 * 02354 * If the coefficients are for the inverted interval, in 02355 * which (a, b) is mapped to (1/b, 1/a), the transformation 02356 * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity, 02357 * this becomes x -> 4a/x - 1. 02358 * 02359 * 02360 * 02361 * SPEED: 02362 * 02363 * Taking advantage of the recurrence properties of the 02364 * Chebyshev polynomials, the routine requires one more 02365 * addition per loop than evaluating a nested polynomial of 02366 * the same degree. 02367 * 02368 */ 02369 02370 /* 02371 Cephes Math Library Release 2.0: April, 1987 02372 Copyright 1985, 1987 by Stephen L. Moshier 02373 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 02374 */ 02375 02376 inline double chbevl(double x, const double array[], int n) 02377 { 02378 double b0, b1, b2; 02379 const double *p = array; 02380 int i; 02381 02382 b0 = *p++; 02383 b1 = 0.0; 02384 i = n - 1; 02385 02386 do 02387 { 02388 b2 = b1; 02389 b1 = b0; 02390 b0 = x * b1 - b2 + *p++; 02391 } 02392 while( --i ); 02393 02394 return 0.5*(b0-b2); 02395 } 02396 02397 02398 /* >>refer Mersenne Twister Matsumoto, M. & Nishimura, T. 1998, ACM Transactions on Modeling 02399 * >>refercon and Computer Simulation (TOMACS), 8, 1998 */ 02400 02401 /******************************************************************** 02402 * This copyright notice must accompany the following block of code * 02403 *******************************************************************/ 02404 02405 /* 02406 A C-program for MT19937, with initialization improved 2002/2/10. 02407 Coded by Takuji Nishimura and Makoto Matsumoto. 02408 This is a faster version by taking Shawn Cokus's optimization, 02409 Matthe Bellew's simplification, Isaku Wada's real version. 02410 02411 Before using, initialize the state by using init_genrand(seed) 02412 or init_by_array(init_key, key_length). 02413 02414 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 02415 All rights reserved. 02416 02417 Redistribution and use in source and binary forms, with or without 02418 modification, are permitted provided that the following conditions 02419 are met: 02420 02421 1. Redistributions of source code must retain the above copyright 02422 notice, this list of conditions and the following disclaimer. 02423 02424 2. Redistributions in binary form must reproduce the above copyright 02425 notice, this list of conditions and the following disclaimer in the 02426 documentation and/or other materials provided with the distribution. 02427 02428 3. The names of its contributors may not be used to endorse or promote 02429 products derived from this software without specific prior written 02430 permission. 02431 02432 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 02433 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 02434 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 02435 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 02436 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 02437 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 02438 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 02439 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 02440 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 02441 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 02442 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 02443 02444 02445 Any feedback is very welcome. 02446 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 02447 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) 02448 */ 02449 02450 /* Period parameters */ 02451 #define N 624 02452 #define M 397 02453 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ 02454 #define UMASK 0x80000000UL /* most significant w-r bits */ 02455 #define LMASK 0x7fffffffUL /* least significant r bits */ 02456 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) 02457 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) 02458 02459 static unsigned long state[N]; /* the array for the state vector */ 02460 static int nleft = 1; 02461 static int initf = 0; 02462 static unsigned long *next; 02463 02464 /* initializes state[N] with a seed */ 02465 void init_genrand(unsigned long s) 02466 { 02467 int j; 02468 state[0]= s & 0xffffffffUL; 02469 for (j=1; j<N; j++) { 02470 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j); 02471 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 02472 /* In the previous versions, MSBs of the seed affect */ 02473 /* only MSBs of the array state[]. */ 02474 /* 2002/01/09 modified by Makoto Matsumoto */ 02475 state[j] &= 0xffffffffUL; /* for >32 bit machines */ 02476 } 02477 nleft = 1; initf = 1; 02478 } 02479 02480 /* initialize by an array with array-length */ 02481 /* init_key is the array for initializing keys */ 02482 /* key_length is its length */ 02483 /* slight change for C++, 2004/2/26 */ 02484 void init_by_array(unsigned long init_key[], int key_length) 02485 { 02486 int i, j, k; 02487 init_genrand(19650218UL); 02488 i=1; j=0; 02489 k = (N>key_length ? N : key_length); 02490 for (; k; k--) { 02491 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) 02492 + init_key[j] + j; /* non linear */ 02493 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 02494 i++; j++; 02495 if (i>=N) { state[0] = state[N-1]; i=1; } 02496 if (j>=key_length) j=0; 02497 } 02498 for (k=N-1; k; k--) { 02499 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL)) 02500 - i; /* non linear */ 02501 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 02502 i++; 02503 if (i>=N) { state[0] = state[N-1]; i=1; } 02504 } 02505 02506 state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 02507 nleft = 1; initf = 1; 02508 } 02509 02510 static void next_state(void) 02511 { 02512 unsigned long *p=state; 02513 int j; 02514 02515 /* if init_genrand() has not been called, */ 02516 /* a default initial seed is used */ 02517 if (initf==0) init_genrand(5489UL); 02518 02519 nleft = N; 02520 next = state; 02521 02522 for (j=N-M+1; --j; p++) 02523 *p = p[M] ^ TWIST(p[0], p[1]); 02524 02525 for (j=M; --j; p++) 02526 *p = p[M-N] ^ TWIST(p[0], p[1]); 02527 02528 *p = p[M-N] ^ TWIST(p[0], state[0]); 02529 } 02530 02531 /* generates a random number on [0,0xffffffff]-interval */ 02532 unsigned long genrand_int32(void) 02533 { 02534 unsigned long y; 02535 02536 if (--nleft == 0) next_state(); 02537 y = *next++; 02538 02539 /* Tempering */ 02540 y ^= (y >> 11); 02541 y ^= (y << 7) & 0x9d2c5680UL; 02542 y ^= (y << 15) & 0xefc60000UL; 02543 y ^= (y >> 18); 02544 02545 return y; 02546 } 02547 02548 /* generates a random number on [0,0x7fffffff]-interval */ 02549 long genrand_int31(void) 02550 { 02551 unsigned long y; 02552 02553 if (--nleft == 0) next_state(); 02554 y = *next++; 02555 02556 /* Tempering */ 02557 y ^= (y >> 11); 02558 y ^= (y << 7) & 0x9d2c5680UL; 02559 y ^= (y << 15) & 0xefc60000UL; 02560 y ^= (y >> 18); 02561 02562 return (long)(y>>1); 02563 } 02564 02565 /* generates a random number on [0,1]-real-interval */ 02566 double genrand_real1(void) 02567 { 02568 unsigned long y; 02569 02570 if (--nleft == 0) next_state(); 02571 y = *next++; 02572 02573 /* Tempering */ 02574 y ^= (y >> 11); 02575 y ^= (y << 7) & 0x9d2c5680UL; 02576 y ^= (y << 15) & 0xefc60000UL; 02577 y ^= (y >> 18); 02578 02579 return (double)y * (1.0/4294967295.0); 02580 /* divided by 2^32-1 */ 02581 } 02582 02583 /* generates a random number on [0,1)-real-interval */ 02584 double genrand_real2(void) 02585 { 02586 unsigned long y; 02587 02588 if (--nleft == 0) next_state(); 02589 y = *next++; 02590 02591 /* Tempering */ 02592 y ^= (y >> 11); 02593 y ^= (y << 7) & 0x9d2c5680UL; 02594 y ^= (y << 15) & 0xefc60000UL; 02595 y ^= (y >> 18); 02596 02597 return (double)y * (1.0/4294967296.0); 02598 /* divided by 2^32 */ 02599 } 02600 02601 /* generates a random number on (0,1)-real-interval */ 02602 double genrand_real3(void) 02603 { 02604 unsigned long y; 02605 02606 if (--nleft == 0) next_state(); 02607 y = *next++; 02608 02609 /* Tempering */ 02610 y ^= (y >> 11); 02611 y ^= (y << 7) & 0x9d2c5680UL; 02612 y ^= (y << 15) & 0xefc60000UL; 02613 y ^= (y >> 18); 02614 02615 return ((double)y + 0.5) * (1.0/4294967296.0); 02616 /* divided by 2^32 */ 02617 } 02618 02619 /* generates a random number on [0,1) with 53-bit resolution*/ 02620 double genrand_res53(void) 02621 { 02622 unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; 02623 return(a*67108864.0+b)*(1.0/9007199254740992.0); 02624 } 02625 02626 /* These real versions are due to Isaku Wada, 2002/01/09 added */ 02627 02628 /************************************************************************ 02629 * This marks the end of the block of code from Matsumoto and Nishimura * 02630 ************************************************************************/