00001
00002
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 inline double dawson(double x, int order);
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 bool linfit(
00047 long n,
00048 const double xorg[],
00049 const double yorg[],
00050 double &a,
00051 double &siga,
00052 double &b,
00053 double &sigb
00054 )
00055 {
00056
00057
00058
00059
00060
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
00089
00090
00091
00092
00093 DEBUG_ENTRY( "linfit()" );
00094
00095 ASSERT( n >= 2 );
00096
00097 valarray<double> x(n);
00098 valarray<double> y(n);
00099
00100 for( long i=0; i < n; i++ )
00101 {
00102 x[i] = xorg[i];
00103 y[i] = yorg[i];
00104 }
00105
00106
00107 a = 0.0;
00108 siga = 0.0;
00109 b = 0.0;
00110 sigb = 0.0;
00111
00112
00113 double s1 = 0.0;
00114 double s2 = 0.0;
00115 for( long i=0; i < n; i++ )
00116 {
00117 s1 += x[i];
00118 s2 += y[i];
00119 }
00120 double rn = (double)n;
00121 double xavg = s1/rn;
00122 double yavg = s2/rn;
00123 double sxx = 0.0;
00124 double sxy = 0.0;
00125 for( long i=0; i < n; i++ )
00126 {
00127 x[i] -= xavg;
00128 y[i] -= yavg;
00129 sxx += pow2(x[i]);
00130 sxy += x[i]*y[i];
00131 }
00132
00133 if( pow2(sxx) == 0.0 )
00134 {
00135 return true;
00136 }
00137
00138
00139 b = sxy/sxx;
00140
00141
00142 a = yavg - b*xavg;
00143
00144
00145 double sum1 = 0.0;
00146 for( long i=0; i < n; i++ )
00147 sum1 += pow2(x[i]*(y[i] - b*x[i]));
00148
00149
00150 sigb = sum1/pow2(sxx);
00151
00152
00153 for( long i=0; i < n; i++ )
00154 siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
00155
00156
00157 sigb = sqrt(sigb);
00158 siga = sqrt(siga)/rn;
00159
00160
00161 for( long i=0; i < n; i++ )
00162 {
00163 x[i] += xavg;
00164 y[i] += yavg;
00165 }
00166 return false;
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181 static const double pre_factorial[NPRE_FACTORIAL] =
00182 {
00183 1.00000000000000000000e+00,
00184 1.00000000000000000000e+00,
00185 2.00000000000000000000e+00,
00186 6.00000000000000000000e+00,
00187 2.40000000000000000000e+01,
00188 1.20000000000000000000e+02,
00189 7.20000000000000000000e+02,
00190 5.04000000000000000000e+03,
00191 4.03200000000000000000e+04,
00192 3.62880000000000000000e+05,
00193 3.62880000000000000000e+06,
00194 3.99168000000000000000e+07,
00195 4.79001600000000000000e+08,
00196 6.22702080000000000000e+09,
00197 8.71782912000000000000e+10,
00198 1.30767436800000000000e+12,
00199 2.09227898880000000000e+13,
00200 3.55687428096000000000e+14,
00201 6.40237370572800000000e+15,
00202 1.21645100408832000000e+17,
00203 2.43290200817664000000e+18,
00204 5.10909421717094400000e+19,
00205 1.12400072777760768000e+21,
00206 2.58520167388849766400e+22,
00207 6.20448401733239439360e+23,
00208 1.55112100433309859840e+25,
00209 4.03291461126605635592e+26,
00210 1.08888694504183521614e+28,
00211 3.04888344611713860511e+29,
00212 8.84176199373970195470e+30,
00213 2.65252859812191058647e+32,
00214 8.22283865417792281807e+33,
00215 2.63130836933693530178e+35,
00216 8.68331761881188649615e+36,
00217 2.95232799039604140861e+38,
00218 1.03331479663861449300e+40,
00219 3.71993326789901217463e+41,
00220 1.37637530912263450457e+43,
00221 5.23022617466601111726e+44,
00222 2.03978820811974433568e+46,
00223 8.15915283247897734264e+47,
00224 3.34525266131638071044e+49,
00225 1.40500611775287989839e+51,
00226 6.04152630633738356321e+52,
00227 2.65827157478844876773e+54,
00228 1.19622220865480194551e+56,
00229 5.50262215981208894940e+57,
00230 2.58623241511168180614e+59,
00231 1.24139155925360726691e+61,
00232 6.08281864034267560801e+62,
00233 3.04140932017133780398e+64,
00234 1.55111875328738227999e+66,
00235 8.06581751709438785585e+67,
00236 4.27488328406002556374e+69,
00237 2.30843697339241380441e+71,
00238 1.26964033536582759243e+73,
00239 7.10998587804863451749e+74,
00240 4.05269195048772167487e+76,
00241 2.35056133128287857145e+78,
00242 1.38683118545689835713e+80,
00243 8.32098711274139014271e+81,
00244 5.07580213877224798711e+83,
00245 3.14699732603879375200e+85,
00246 1.98260831540444006372e+87,
00247 1.26886932185884164078e+89,
00248 8.24765059208247066512e+90,
00249 5.44344939077443063905e+92,
00250 3.64711109181886852801e+94,
00251 2.48003554243683059915e+96,
00252 1.71122452428141311337e+98,
00253 1.19785716699698917933e+100,
00254 8.50478588567862317347e+101,
00255 6.12344583768860868500e+103,
00256 4.47011546151268434004e+105,
00257 3.30788544151938641157e+107,
00258 2.48091408113953980872e+109,
00259 1.88549470166605025458e+111,
00260 1.45183092028285869606e+113,
00261 1.13242811782062978295e+115,
00262 8.94618213078297528506e+116,
00263 7.15694570462638022794e+118,
00264 5.79712602074736798470e+120,
00265 4.75364333701284174746e+122,
00266 3.94552396972065865030e+124,
00267 3.31424013456535326627e+126,
00268 2.81710411438055027626e+128,
00269 2.42270953836727323750e+130,
00270 2.10775729837952771662e+132,
00271 1.85482642257398439069e+134,
00272 1.65079551609084610774e+136,
00273 1.48571596448176149700e+138,
00274 1.35200152767840296226e+140,
00275 1.24384140546413072522e+142,
00276 1.15677250708164157442e+144,
00277 1.08736615665674307994e+146,
00278 1.03299784882390592592e+148,
00279 9.91677934870949688836e+149,
00280 9.61927596824821198159e+151,
00281 9.42689044888324774164e+153,
00282 9.33262154439441526381e+155,
00283 9.33262154439441526388e+157,
00284 9.42594775983835941673e+159,
00285 9.61446671503512660515e+161,
00286 9.90290071648618040340e+163,
00287 1.02990167451456276198e+166,
00288 1.08139675824029090008e+168,
00289 1.14628056373470835406e+170,
00290 1.22652020319613793888e+172,
00291 1.32464181945182897396e+174,
00292 1.44385958320249358163e+176,
00293 1.58824554152274293982e+178,
00294 1.76295255109024466316e+180,
00295 1.97450685722107402277e+182,
00296 2.23119274865981364576e+184,
00297 2.54355973347218755612e+186,
00298 2.92509369349301568964e+188,
00299 3.39310868445189820004e+190,
00300 3.96993716080872089396e+192,
00301 4.68452584975429065488e+194,
00302 5.57458576120760587943e+196,
00303 6.68950291344912705515e+198,
00304 8.09429852527344373681e+200,
00305 9.87504420083360135884e+202,
00306 1.21463043670253296712e+205,
00307 1.50614174151114087918e+207,
00308 1.88267717688892609901e+209,
00309 2.37217324288004688470e+211,
00310 3.01266001845765954361e+213,
00311 3.85620482362580421582e+215,
00312 4.97450422247728743840e+217,
00313 6.46685548922047366972e+219,
00314 8.47158069087882050755e+221,
00315 1.11824865119600430699e+224,
00316 1.48727070609068572828e+226,
00317 1.99294274616151887582e+228,
00318 2.69047270731805048244e+230,
00319 3.65904288195254865604e+232,
00320 5.01288874827499165889e+234,
00321 6.91778647261948848943e+236,
00322 9.61572319694108900019e+238,
00323 1.34620124757175246000e+241,
00324 1.89814375907617096864e+243,
00325 2.69536413788816277557e+245,
00326 3.85437071718007276916e+247,
00327 5.55029383273930478744e+249,
00328 8.04792605747199194159e+251,
00329 1.17499720439091082343e+254,
00330 1.72724589045463891049e+256,
00331 2.55632391787286558753e+258,
00332 3.80892263763056972532e+260,
00333 5.71338395644585458806e+262,
00334 8.62720977423324042775e+264,
00335 1.31133588568345254503e+267,
00336 2.00634390509568239384e+269,
00337 3.08976961384735088657e+271,
00338 4.78914290146339387432e+273,
00339 7.47106292628289444390e+275,
00340 1.17295687942641442768e+278,
00341 1.85327186949373479574e+280,
00342 2.94670227249503832518e+282,
00343 4.71472363599206132029e+284,
00344 7.59070505394721872577e+286,
00345 1.22969421873944943358e+289,
00346 2.00440157654530257674e+291,
00347 3.28721858553429622598e+293,
00348 5.42391066613158877297e+295,
00349 9.00369170577843736335e+297,
00350 1.50361651486499903974e+300,
00351 2.52607574497319838672e+302,
00352 4.26906800900470527345e+304,
00353 7.25741561530799896496e+306
00354 };
00355
00356 double factorial( long n )
00357 {
00358 DEBUG_ENTRY( "factorial()" );
00359
00360 if( n < 0 || n >= NPRE_FACTORIAL )
00361 {
00362 fprintf( ioQQQ, "factorial: domain error\n" );
00363 cdEXIT(EXIT_FAILURE);
00364 }
00365
00366 return pre_factorial[n];
00367 }
00368
00369
00370
00371 class t_lfact : public Singleton<t_lfact>
00372 {
00373 friend class Singleton<t_lfact>;
00374 protected:
00375 t_lfact()
00376 {
00377 p_lf.reserve( 512 );
00378 p_lf.push_back( 0. );
00379 p_lf.push_back( 0. );
00380 }
00381 private:
00382 vector<double> p_lf;
00383 public:
00384 double get_lfact( unsigned long n )
00385 {
00386 if( n < p_lf.size() )
00387 {
00388 return p_lf[n];
00389 }
00390 else
00391 {
00392 for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ )
00393 p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) );
00394 return p_lf[n];
00395 }
00396 }
00397 };
00398
00399 double lfactorial( long n )
00400 {
00401
00402
00403
00404
00405
00406
00407
00408
00409 DEBUG_ENTRY( "lfactorial()" );
00410
00411 if( n < 0 )
00412 {
00413 fprintf( ioQQQ, "lfactorial: domain error\n" );
00414 cdEXIT(EXIT_FAILURE);
00415 }
00416
00417 return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 complex<double> cdgamma(complex<double> x)
00433 {
00434 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
00435
00436 DEBUG_ENTRY( "cdgamma()" );
00437
00438 xr = x.real();
00439 xi = x.imag();
00440 if(xr < 0)
00441 {
00442 wr = 1. - xr;
00443 wi = -xi;
00444 }
00445 else
00446 {
00447 wr = xr;
00448 wi = xi;
00449 }
00450 ur = wr + 6.00009857740312429;
00451 vr = ur * (wr + 4.99999857982434025) - wi * wi;
00452 vi = wi * (wr + 4.99999857982434025) + ur * wi;
00453 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
00454 0.293729529320536228;
00455 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
00456 ur = vr * (wr + 4.00000003016801681) - vi * wi;
00457 ui = vi * (wr + 4.00000003016801681) + vr * wi;
00458 vr = ur * (wr + 2.99999999944915534) - ui * wi;
00459 vi = ui * (wr + 2.99999999944915534) + ur * wi;
00460 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
00461 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
00462 ur = vr * (wr + 2.00000000000603851) - vi * wi;
00463 ui = vi * (wr + 2.00000000000603851) + vr * wi;
00464 vr = ur * (wr + 0.999999999999975753) - ui * wi;
00465 vi = ui * (wr + 0.999999999999975753) + ur * wi;
00466 yr += ur * 10.5400280458730808 + vr;
00467 yi += ui * 10.5400280458730808 + vi;
00468 ur = vr * wr - vi * wi;
00469 ui = vi * wr + vr * wi;
00470 t = ur * ur + ui * ui;
00471 vr = yr * ur + yi * ui + t * 0.0327673720261526849;
00472 vi = yi * ur - yr * ui;
00473 yr = wr + 7.31790632447016203;
00474 ur = log(yr * yr + wi * wi) * 0.5 - 1;
00475 ui = atan2(wi, yr);
00476 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
00477 yi = ui * (wr - 0.5) + ur * wi;
00478 ur = yr * cos(yi);
00479 ui = yr * sin(yi);
00480 yr = ur * vr - ui * vi;
00481 yi = ui * vr + ur * vi;
00482 if(xr < 0)
00483 {
00484 wr = xr * 3.14159265358979324;
00485 wi = exp(xi * 3.14159265358979324);
00486 vi = 1 / wi;
00487 ur = (vi + wi) * sin(wr);
00488 ui = (vi - wi) * cos(wr);
00489 vr = ur * yr + ui * yi;
00490 vi = ui * yr - ur * yi;
00491 ur = 6.2831853071795862 / (vr * vr + vi * vi);
00492 yr = ur * vr;
00493 yi = ur * vi;
00494 }
00495 return complex<double>( yr, yi );
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 static const double b0_PP[7] = {
00617 7.96936729297347051624e-4,
00618 8.28352392107440799803e-2,
00619 1.23953371646414299388e0,
00620 5.44725003058768775090e0,
00621 8.74716500199817011941e0,
00622 5.30324038235394892183e0,
00623 9.99999999999999997821e-1,
00624 };
00625
00626 static const double b0_PQ[7] = {
00627 9.24408810558863637013e-4,
00628 8.56288474354474431428e-2,
00629 1.25352743901058953537e0,
00630 5.47097740330417105182e0,
00631 8.76190883237069594232e0,
00632 5.30605288235394617618e0,
00633 1.00000000000000000218e0,
00634 };
00635
00636 static const double b0_QP[8] = {
00637 -1.13663838898469149931e-2,
00638 -1.28252718670509318512e0,
00639 -1.95539544257735972385e1,
00640 -9.32060152123768231369e1,
00641 -1.77681167980488050595e2,
00642 -1.47077505154951170175e2,
00643 -5.14105326766599330220e1,
00644 -6.05014350600728481186e0,
00645 };
00646
00647 static const double b0_QQ[7] = {
00648
00649 6.43178256118178023184e1,
00650 8.56430025976980587198e2,
00651 3.88240183605401609683e3,
00652 7.24046774195652478189e3,
00653 5.93072701187316984827e3,
00654 2.06209331660327847417e3,
00655 2.42005740240291393179e2,
00656 };
00657
00658 static const double b0_YP[8] = {
00659 1.55924367855235737965e4,
00660 -1.46639295903971606143e7,
00661 5.43526477051876500413e9,
00662 -9.82136065717911466409e11,
00663 8.75906394395366999549e13,
00664 -3.46628303384729719441e15,
00665 4.42733268572569800351e16,
00666 -1.84950800436986690637e16,
00667 };
00668
00669 static const double b0_YQ[7] = {
00670
00671 1.04128353664259848412e3,
00672 6.26107330137134956842e5,
00673 2.68919633393814121987e8,
00674 8.64002487103935000337e10,
00675 2.02979612750105546709e13,
00676 3.17157752842975028269e15,
00677 2.50596256172653059228e17,
00678 };
00679
00680
00681 static const double DR1 = 5.78318596294678452118e0;
00682
00683 static const double DR2 = 3.04712623436620863991e1;
00684
00685 static double b0_RP[4] = {
00686 -4.79443220978201773821e9,
00687 1.95617491946556577543e12,
00688 -2.49248344360967716204e14,
00689 9.70862251047306323952e15,
00690 };
00691
00692 static double b0_RQ[8] = {
00693
00694 4.99563147152651017219e2,
00695 1.73785401676374683123e5,
00696 4.84409658339962045305e7,
00697 1.11855537045356834862e10,
00698 2.11277520115489217587e12,
00699 3.10518229857422583814e14,
00700 3.18121955943204943306e16,
00701 1.71086294081043136091e18,
00702 };
00703
00704 static const double TWOOPI = 2./PI;
00705 static const double SQ2OPI = sqrt(2./PI);
00706 static const double PIO4 = PI/4.;
00707
00708 double bessel_j0(double x)
00709 {
00710 double w, z, p, q, xn;
00711
00712 DEBUG_ENTRY( "bessel_j0()" );
00713
00714 if( x < 0 )
00715 x = -x;
00716
00717 if( x <= 5.0 )
00718 {
00719 z = x * x;
00720 if( x < 1.0e-5 )
00721 return 1.0 - z/4.0;
00722
00723 p = (z - DR1) * (z - DR2);
00724 p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
00725 return p;
00726 }
00727
00728 w = 5.0/x;
00729 q = 25.0/(x*x);
00730 p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
00731 q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
00732 xn = x - PIO4;
00733 p = p * cos(xn) - w * q * sin(xn);
00734 return p * SQ2OPI / sqrt(x);
00735 }
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746 double bessel_y0(double x)
00747 {
00748 double w, z, p, q, xn;
00749
00750 DEBUG_ENTRY( "bessel_y0()" );
00751
00752 if( x <= 5.0 )
00753 {
00754 if( x <= 0.0 )
00755 {
00756 fprintf( ioQQQ, "bessel_y0: domain error\n" );
00757 cdEXIT(EXIT_FAILURE);
00758 }
00759 z = x * x;
00760 w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
00761 w += TWOOPI * log(x) * bessel_j0(x);
00762 return w;
00763 }
00764
00765 w = 5.0/x;
00766 z = 25.0 / (x * x);
00767 p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
00768 q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
00769 xn = x - PIO4;
00770 p = p * sin(xn) + w * q * cos(xn);
00771 return p * SQ2OPI / sqrt(x);
00772 }
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852 static const double b1_RP[4] = {
00853 -8.99971225705559398224e8,
00854 4.52228297998194034323e11,
00855 -7.27494245221818276015e13,
00856 3.68295732863852883286e15,
00857 };
00858
00859 static const double b1_RQ[8] = {
00860
00861 6.20836478118054335476e2,
00862 2.56987256757748830383e5,
00863 8.35146791431949253037e7,
00864 2.21511595479792499675e10,
00865 4.74914122079991414898e12,
00866 7.84369607876235854894e14,
00867 8.95222336184627338078e16,
00868 5.32278620332680085395e18,
00869 };
00870
00871 static const double b1_PP[7] = {
00872 7.62125616208173112003e-4,
00873 7.31397056940917570436e-2,
00874 1.12719608129684925192e0,
00875 5.11207951146807644818e0,
00876 8.42404590141772420927e0,
00877 5.21451598682361504063e0,
00878 1.00000000000000000254e0,
00879 };
00880
00881 static const double b1_PQ[7] = {
00882 5.71323128072548699714e-4,
00883 6.88455908754495404082e-2,
00884 1.10514232634061696926e0,
00885 5.07386386128601488557e0,
00886 8.39985554327604159757e0,
00887 5.20982848682361821619e0,
00888 9.99999999999999997461e-1,
00889 };
00890
00891 static const double b1_QP[8] = {
00892 5.10862594750176621635e-2,
00893 4.98213872951233449420e0,
00894 7.58238284132545283818e1,
00895 3.66779609360150777800e2,
00896 7.10856304998926107277e2,
00897 5.97489612400613639965e2,
00898 2.11688757100572135698e2,
00899 2.52070205858023719784e1,
00900 };
00901
00902 static const double b1_QQ[7] = {
00903
00904 7.42373277035675149943e1,
00905 1.05644886038262816351e3,
00906 4.98641058337653607651e3,
00907 9.56231892404756170795e3,
00908 7.99704160447350683650e3,
00909 2.82619278517639096600e3,
00910 3.36093607810698293419e2,
00911 };
00912
00913 static const double b1_YP[6] = {
00914 1.26320474790178026440e9,
00915 -6.47355876379160291031e11,
00916 1.14509511541823727583e14,
00917 -8.12770255501325109621e15,
00918 2.02439475713594898196e17,
00919 -7.78877196265950026825e17,
00920 };
00921
00922 static const double b1_YQ[8] = {
00923
00924 5.94301592346128195359E2,
00925 2.35564092943068577943E5,
00926 7.34811944459721705660E7,
00927 1.87601316108706159478E10,
00928 3.88231277496238566008E12,
00929 6.20557727146953693363E14,
00930 6.87141087355300489866E16,
00931 3.97270608116560655612E18,
00932 };
00933
00934 static const double Z1 = 1.46819706421238932572E1;
00935 static const double Z2 = 4.92184563216946036703E1;
00936
00937 static const double THPIO4 = 3.*PI/4.;
00938
00939 double bessel_j1(double x)
00940 {
00941 double w, z, p, q, xn;
00942
00943 DEBUG_ENTRY( "bessel_j1()" );
00944
00945 w = x;
00946 if( x < 0 )
00947 w = -x;
00948
00949 if( w <= 5.0 )
00950 {
00951 z = x * x;
00952 w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
00953 w = w * x * (z - Z1) * (z - Z2);
00954 return w;
00955 }
00956
00957 w = 5.0/x;
00958 z = w * w;
00959 p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
00960 q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
00961 xn = x - THPIO4;
00962 p = p * cos(xn) - w * q * sin(xn);
00963 return p * SQ2OPI / sqrt(x);
00964 }
00965
00966 double bessel_y1(double x)
00967 {
00968 double w, z, p, q, xn;
00969
00970 DEBUG_ENTRY( "bessel_y1()" );
00971
00972 if( x <= 5.0 )
00973 {
00974 if( x <= 0.0 )
00975 {
00976 fprintf( ioQQQ, "bessel_y1: domain error\n" );
00977 cdEXIT(EXIT_FAILURE);
00978 }
00979 z = x * x;
00980 w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
00981 w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
00982 return w;
00983 }
00984
00985 w = 5.0/x;
00986 z = w * w;
00987 p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
00988 q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
00989 xn = x - THPIO4;
00990 p = p * sin(xn) + w * q * cos(xn);
00991 return p * SQ2OPI / sqrt(x);
00992 }
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 double bessel_jn(int n, double x)
01043 {
01044 double pkm2, pkm1, pk, xk, r, ans;
01045 int k, sign;
01046
01047 DEBUG_ENTRY( "bessel_jn()" );
01048
01049 if( n < 0 )
01050 {
01051 n = -n;
01052 if( (n & 1) == 0 )
01053 sign = 1;
01054 else
01055 sign = -1;
01056 }
01057 else
01058 sign = 1;
01059
01060 if( x < 0.0 )
01061 {
01062 if( n & 1 )
01063 sign = -sign;
01064 x = -x;
01065 }
01066
01067 if( x < DBL_EPSILON )
01068 {
01069 return sign * powi(x/2.,n)/factorial(n);
01070 }
01071
01072 if( n == 0 )
01073 {
01074 return sign * bessel_j0(x);
01075 }
01076 if( n == 1 )
01077 {
01078 return sign * bessel_j1(x);
01079 }
01080
01081 if( n == 2 && x > 0.1 )
01082 {
01083 return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
01084 }
01085
01086
01087 k = 52;
01088
01089 pk = 2 * (n + k);
01090 ans = pk;
01091 xk = x * x;
01092
01093 do
01094 {
01095 pk -= 2.0;
01096 ans = pk - (xk/ans);
01097 }
01098 while( --k > 0 );
01099 ans = x/ans;
01100
01101
01102 pk = 1.0;
01103 pkm1 = 1.0/ans;
01104 k = n-1;
01105 r = 2 * k;
01106
01107 do
01108 {
01109 pkm2 = (pkm1 * r - pk * x) / x;
01110 pk = pkm1;
01111 pkm1 = pkm2;
01112 r -= 2.0;
01113 }
01114 while( --k > 0 );
01115
01116 if( fabs(pk) > fabs(pkm1) )
01117 ans = bessel_j1(x)/pk;
01118 else
01119 ans = bessel_j0(x)/pkm1;
01120 return sign * ans;
01121 }
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177 double bessel_yn(int n, double x)
01178 {
01179 double an, anm1, anm2, r;
01180 int k, sign;
01181
01182 DEBUG_ENTRY( "bessel_yn()" );
01183
01184 if( n < 0 )
01185 {
01186 n = -n;
01187 if( (n & 1) == 0 )
01188 sign = 1;
01189 else
01190 sign = -1;
01191 }
01192 else
01193 sign = 1;
01194
01195 if( n == 0 )
01196 {
01197 return sign * bessel_y0(x);
01198 }
01199 if( n == 1 )
01200 {
01201 return sign * bessel_y1(x);
01202 }
01203
01204
01205 if( x <= 0.0 )
01206 {
01207 fprintf( ioQQQ, "bessel_yn: domain error\n" );
01208 cdEXIT(EXIT_FAILURE);
01209 }
01210
01211
01212 anm2 = bessel_y0(x);
01213 anm1 = bessel_y1(x);
01214 k = 1;
01215 r = 2 * k;
01216 do
01217 {
01218 an = r * anm1 / x - anm2;
01219 anm2 = anm1;
01220 anm1 = an;
01221 r += 2.0;
01222 ++k;
01223 }
01224 while( k < n );
01225 return sign * an;
01226 }
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311 static const double k0_A[] =
01312 {
01313 1.37446543561352307156e-16,
01314 4.25981614279661018399e-14,
01315 1.03496952576338420167e-11,
01316 1.90451637722020886025e-9,
01317 2.53479107902614945675e-7,
01318 2.28621210311945178607e-5,
01319 1.26461541144692592338e-3,
01320 3.59799365153615016266e-2,
01321 3.44289899924628486886e-1,
01322 -5.35327393233902768720e-1
01323 };
01324
01325
01326
01327
01328
01329
01330
01331 static const double k0_B[] = {
01332 5.30043377268626276149e-18,
01333 -1.64758043015242134646e-17,
01334 5.21039150503902756861e-17,
01335 -1.67823109680541210385e-16,
01336 5.51205597852431940784e-16,
01337 -1.84859337734377901440e-15,
01338 6.34007647740507060557e-15,
01339 -2.22751332699166985548e-14,
01340 8.03289077536357521100e-14,
01341 -2.98009692317273043925e-13,
01342 1.14034058820847496303e-12,
01343 -4.51459788337394416547e-12,
01344 1.85594911495471785253e-11,
01345 -7.95748924447710747776e-11,
01346 3.57739728140030116597e-10,
01347 -1.69753450938905987466e-9,
01348 8.57403401741422608519e-9,
01349 -4.66048989768794782956e-8,
01350 2.76681363944501510342e-7,
01351 -1.83175552271911948767e-6,
01352 1.39498137188764993662e-5,
01353 -1.28495495816278026384e-4,
01354 1.56988388573005337491e-3,
01355 -3.14481013119645005427e-2,
01356 2.44030308206595545468e0
01357 };
01358
01359 double bessel_k0(double x)
01360 {
01361 double y, z;
01362
01363 DEBUG_ENTRY( "bessel_k0()" );
01364
01365 if( x <= 0.0 )
01366 {
01367 fprintf( ioQQQ, "bessel_k0: domain error\n" );
01368 cdEXIT(EXIT_FAILURE);
01369 }
01370
01371 if( x <= 2.0 )
01372 {
01373 y = x * x - 2.0;
01374 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01375 return y;
01376 }
01377 z = 8.0/x - 2.0;
01378 y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x);
01379 return y;
01380 }
01381
01382 double bessel_k0_scaled(double x)
01383 {
01384 double y;
01385
01386 DEBUG_ENTRY( "bessel_k0_scaled()" );
01387
01388 if( x <= 0.0 )
01389 {
01390 fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" );
01391 cdEXIT(EXIT_FAILURE);
01392 }
01393
01394 if( x <= 2.0 )
01395 {
01396 y = x * x - 2.0;
01397 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01398 return y * exp(x);
01399 }
01400 return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x);
01401 }
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485 static const double k1_A[] =
01486 {
01487 -7.02386347938628759343e-18,
01488 -2.42744985051936593393e-15,
01489 -6.66690169419932900609e-13,
01490 -1.41148839263352776110e-10,
01491 -2.21338763073472585583e-8,
01492 -2.43340614156596823496e-6,
01493 -1.73028895751305206302e-4,
01494 -6.97572385963986435018e-3,
01495 -1.22611180822657148235e-1,
01496 -3.53155960776544875667e-1,
01497 1.52530022733894777053e0
01498 };
01499
01500
01501
01502
01503
01504
01505
01506 static const double k1_B[] =
01507 {
01508 -5.75674448366501715755e-18,
01509 1.79405087314755922667e-17,
01510 -5.68946255844285935196e-17,
01511 1.83809354436663880070e-16,
01512 -6.05704724837331885336e-16,
01513 2.03870316562433424052e-15,
01514 -7.01983709041831346144e-15,
01515 2.47715442448130437068e-14,
01516 -8.97670518232499435011e-14,
01517 3.34841966607842919884e-13,
01518 -1.28917396095102890680e-12,
01519 5.13963967348173025100e-12,
01520 -2.12996783842756842877e-11,
01521 9.21831518760500529508e-11,
01522 -4.19035475934189648750e-10,
01523 2.01504975519703286596e-9,
01524 -1.03457624656780970260e-8,
01525 5.74108412545004946722e-8,
01526 -3.50196060308781257119e-7,
01527 2.40648494783721712015e-6,
01528 -1.93619797416608296024e-5,
01529 1.95215518471351631108e-4,
01530 -2.85781685962277938680e-3,
01531 1.03923736576817238437e-1,
01532 2.72062619048444266945e0
01533 };
01534
01535 double bessel_k1(double x)
01536 {
01537 double y, z;
01538
01539 DEBUG_ENTRY( "bessel_k1()" );
01540
01541 z = 0.5 * x;
01542 if( z <= 0.0 )
01543 {
01544 fprintf( ioQQQ, "bessel_k1: domain error\n" );
01545 cdEXIT(EXIT_FAILURE);
01546 }
01547
01548 if( x <= 2.0 )
01549 {
01550 y = x * x - 2.0;
01551 y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01552 return y;
01553 }
01554 return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01555 }
01556
01557 double bessel_k1_scaled(double x)
01558 {
01559 double y;
01560
01561 DEBUG_ENTRY( "bessel_k1_scaled()" );
01562
01563 if( x <= 0.0 )
01564 {
01565 fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" );
01566 cdEXIT(EXIT_FAILURE);
01567 }
01568
01569 if( x <= 2.0 )
01570 {
01571 y = x * x - 2.0;
01572 y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01573 return y * exp(x);
01574 }
01575 return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01576 }
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657 static const double i0_A[] =
01658 {
01659 -4.41534164647933937950e-18,
01660 3.33079451882223809783e-17,
01661 -2.43127984654795469359e-16,
01662 1.71539128555513303061e-15,
01663 -1.16853328779934516808e-14,
01664 7.67618549860493561688e-14,
01665 -4.85644678311192946090e-13,
01666 2.95505266312963983461e-12,
01667 -1.72682629144155570723e-11,
01668 9.67580903537323691224e-11,
01669 -5.18979560163526290666e-10,
01670 2.65982372468238665035e-9,
01671 -1.30002500998624804212e-8,
01672 6.04699502254191894932e-8,
01673 -2.67079385394061173391e-7,
01674 1.11738753912010371815e-6,
01675 -4.41673835845875056359e-6,
01676 1.64484480707288970893e-5,
01677 -5.75419501008210370398e-5,
01678 1.88502885095841655729e-4,
01679 -5.76375574538582365885e-4,
01680 1.63947561694133579842e-3,
01681 -4.32430999505057594430e-3,
01682 1.05464603945949983183e-2,
01683 -2.37374148058994688156e-2,
01684 4.93052842396707084878e-2,
01685 -9.49010970480476444210e-2,
01686 1.71620901522208775349e-1,
01687 -3.04682672343198398683e-1,
01688 6.76795274409476084995e-1
01689 };
01690
01691
01692
01693
01694
01695
01696
01697 static const double i0_B[] =
01698 {
01699 -7.23318048787475395456e-18,
01700 -4.83050448594418207126e-18,
01701 4.46562142029675999901e-17,
01702 3.46122286769746109310e-17,
01703 -2.82762398051658348494e-16,
01704 -3.42548561967721913462e-16,
01705 1.77256013305652638360e-15,
01706 3.81168066935262242075e-15,
01707 -9.55484669882830764870e-15,
01708 -4.15056934728722208663e-14,
01709 1.54008621752140982691e-14,
01710 3.85277838274214270114e-13,
01711 7.18012445138366623367e-13,
01712 -1.79417853150680611778e-12,
01713 -1.32158118404477131188e-11,
01714 -3.14991652796324136454e-11,
01715 1.18891471078464383424e-11,
01716 4.94060238822496958910e-10,
01717 3.39623202570838634515e-9,
01718 2.26666899049817806459e-8,
01719 2.04891858946906374183e-7,
01720 2.89137052083475648297e-6,
01721 6.88975834691682398426e-5,
01722 3.36911647825569408990e-3,
01723 8.04490411014108831608e-1
01724 };
01725
01726 double bessel_i0(double x)
01727 {
01728 double y;
01729
01730 DEBUG_ENTRY( "bessel_i0()" );
01731
01732 if( x < 0 )
01733 x = -x;
01734
01735 if( x <= 8.0 )
01736 {
01737 y = (x/2.0) - 2.0;
01738 return exp(x) * chbevl( y, i0_A, 30 );
01739 }
01740 return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01741 }
01742
01743 double bessel_i0_scaled(double x)
01744 {
01745 double y;
01746
01747 DEBUG_ENTRY( "bessel_i0_scaled()" );
01748
01749 if( x < 0 )
01750 x = -x;
01751
01752 if( x <= 8.0 )
01753 {
01754 y = (x/2.0) - 2.0;
01755 return chbevl( y, i0_A, 30 );
01756 }
01757 return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01758 }
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840 static double i1_A[] =
01841 {
01842 2.77791411276104639959e-18,
01843 -2.11142121435816608115e-17,
01844 1.55363195773620046921e-16,
01845 -1.10559694773538630805e-15,
01846 7.60068429473540693410e-15,
01847 -5.04218550472791168711e-14,
01848 3.22379336594557470981e-13,
01849 -1.98397439776494371520e-12,
01850 1.17361862988909016308e-11,
01851 -6.66348972350202774223e-11,
01852 3.62559028155211703701e-10,
01853 -1.88724975172282928790e-9,
01854 9.38153738649577178388e-9,
01855 -4.44505912879632808065e-8,
01856 2.00329475355213526229e-7,
01857 -8.56872026469545474066e-7,
01858 3.47025130813767847674e-6,
01859 -1.32731636560394358279e-5,
01860 4.78156510755005422638e-5,
01861 -1.61760815825896745588e-4,
01862 5.12285956168575772895e-4,
01863 -1.51357245063125314899e-3,
01864 4.15642294431288815669e-3,
01865 -1.05640848946261981558e-2,
01866 2.47264490306265168283e-2,
01867 -5.29459812080949914269e-2,
01868 1.02643658689847095384e-1,
01869 -1.76416518357834055153e-1,
01870 2.52587186443633654823e-1
01871 };
01872
01873
01874
01875
01876
01877
01878
01879 static double i1_B[] =
01880 {
01881 7.51729631084210481353e-18,
01882 4.41434832307170791151e-18,
01883 -4.65030536848935832153e-17,
01884 -3.20952592199342395980e-17,
01885 2.96262899764595013876e-16,
01886 3.30820231092092828324e-16,
01887 -1.88035477551078244854e-15,
01888 -3.81440307243700780478e-15,
01889 1.04202769841288027642e-14,
01890 4.27244001671195135429e-14,
01891 -2.10154184277266431302e-14,
01892 -4.08355111109219731823e-13,
01893 -7.19855177624590851209e-13,
01894 2.03562854414708950722e-12,
01895 1.41258074366137813316e-11,
01896 3.25260358301548823856e-11,
01897 -1.89749581235054123450e-11,
01898 -5.58974346219658380687e-10,
01899 -3.83538038596423702205e-9,
01900 -2.63146884688951950684e-8,
01901 -2.51223623787020892529e-7,
01902 -3.88256480887769039346e-6,
01903 -1.10588938762623716291e-4,
01904 -9.76109749136146840777e-3,
01905 7.78576235018280120474e-1
01906 };
01907
01908 double bessel_i1(double x)
01909 {
01910 double y, z;
01911
01912 DEBUG_ENTRY( "bessel_i1()" );
01913
01914 z = fabs(x);
01915 if( z <= 8.0 )
01916 {
01917 y = (z/2.0) - 2.0;
01918 z = chbevl( y, i1_A, 29 ) * z * exp(z);
01919 }
01920 else
01921 {
01922 z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01923 }
01924 if( x < 0.0 )
01925 z = -z;
01926 return z;
01927 }
01928
01929 double bessel_i1_scaled(double x)
01930 {
01931 double y, z;
01932
01933 DEBUG_ENTRY( "bessel_i1_scaled()" );
01934
01935 z = fabs(x);
01936 if( z <= 8.0 )
01937 {
01938 y = (z/2.0) - 2.0;
01939 z = chbevl( y, i1_A, 29 ) * z;
01940 }
01941 else
01942 {
01943 z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01944 }
01945 if( x < 0.0 )
01946 z = -z;
01947 return z;
01948 }
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009 static const double elk_P[] =
02010 {
02011 1.37982864606273237150e-4,
02012 2.28025724005875567385e-3,
02013 7.97404013220415179367e-3,
02014 9.85821379021226008714e-3,
02015 6.87489687449949877925e-3,
02016 6.18901033637687613229e-3,
02017 8.79078273952743772254e-3,
02018 1.49380448916805252718e-2,
02019 3.08851465246711995998e-2,
02020 9.65735902811690126535e-2,
02021 1.38629436111989062502e0
02022 };
02023
02024 static const double elk_Q[] =
02025 {
02026 2.94078955048598507511e-5,
02027 9.14184723865917226571e-4,
02028 5.94058303753167793257e-3,
02029 1.54850516649762399335e-2,
02030 2.39089602715924892727e-2,
02031 3.01204715227604046988e-2,
02032 3.73774314173823228969e-2,
02033 4.88280347570998239232e-2,
02034 7.03124996963957469739e-2,
02035 1.24999999999870820058e-1,
02036 4.99999999999999999821e-1
02037 };
02038
02039 static const double C1 = 1.3862943611198906188e0;
02040
02041 double ellpk(double x)
02042 {
02043 DEBUG_ENTRY( "ellpk()" );
02044
02045 if( x < 0.0 || x > 1.0 )
02046 {
02047 fprintf( ioQQQ, "ellpk: domain error\n" );
02048 cdEXIT(EXIT_FAILURE);
02049 }
02050
02051 if( x > DBL_EPSILON )
02052 {
02053 return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
02054 }
02055 else
02056 {
02057 if( x == 0.0 )
02058 {
02059 fprintf( ioQQQ, "ellpk: domain error\n" );
02060 cdEXIT(EXIT_FAILURE);
02061 }
02062 else
02063 {
02064 return C1 - 0.5 * log(x);
02065 }
02066 }
02067 }
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117 static const double MAXLOG = log(DBL_MAX);
02118 static const double BIG = 1.44115188075855872E+17;
02119
02120
02121 double expn(int n, double x)
02122 {
02123 double ans, r, t, yk, xk;
02124 double pk, pkm1, pkm2, qk, qkm1, qkm2;
02125 double psi, z;
02126 int i, k;
02127
02128 DEBUG_ENTRY( "expn()" );
02129
02130 if( n < 0 || x < 0. )
02131 {
02132 fprintf( ioQQQ, "expn: domain error\n" );
02133 cdEXIT(EXIT_FAILURE);
02134 }
02135
02136 if( x > MAXLOG )
02137 {
02138 return 0.0;
02139 }
02140
02141 if( x == 0.0 )
02142 {
02143 if( n < 2 )
02144 {
02145 fprintf( ioQQQ, "expn: domain error\n" );
02146 cdEXIT(EXIT_FAILURE);
02147 }
02148 else
02149 {
02150 return 1.0/((double)n-1.0);
02151 }
02152 }
02153
02154 if( n == 0 )
02155 {
02156 return exp(-x)/x;
02157 }
02158
02159
02160 if( n > 5000 )
02161 {
02162 xk = x + n;
02163 yk = 1.0 / (xk * xk);
02164 t = n;
02165 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
02166 ans = yk * (ans + t * (t - 2.0 * x));
02167 ans = yk * (ans + t);
02168 ans = (ans + 1.0) * exp( -x ) / xk;
02169 return ans;
02170 }
02171
02172 if( x <= 1.0 )
02173 {
02174
02175 psi = -EULER - log(x);
02176 for( i=1; i < n; i++ )
02177 psi = psi + 1.0/i;
02178
02179 z = -x;
02180 xk = 0.0;
02181 yk = 1.0;
02182 pk = 1.0 - n;
02183 if( n == 1 )
02184 ans = 0.0;
02185 else
02186 ans = 1.0/pk;
02187 do
02188 {
02189 xk += 1.0;
02190 yk *= z/xk;
02191 pk += 1.0;
02192 if( pk != 0.0 )
02193 {
02194 ans += yk/pk;
02195 }
02196 if( ans != 0.0 )
02197 t = fabs(yk/ans);
02198 else
02199 t = 1.0;
02200 }
02201 while( t > DBL_EPSILON );
02202 ans = powi(z,n-1)*psi/factorial(n-1) - ans;
02203 return ans;
02204 }
02205 else
02206 {
02207
02208 k = 1;
02209 pkm2 = 1.0;
02210 qkm2 = x;
02211 pkm1 = 1.0;
02212 qkm1 = x + n;
02213 ans = pkm1/qkm1;
02214
02215 do
02216 {
02217 k += 1;
02218 if( is_odd(k) )
02219 {
02220 yk = 1.0;
02221 xk = static_cast<double>(n + (k-1)/2);
02222 }
02223 else
02224 {
02225 yk = x;
02226 xk = static_cast<double>(k/2);
02227 }
02228 pk = pkm1 * yk + pkm2 * xk;
02229 qk = qkm1 * yk + qkm2 * xk;
02230 if( qk != 0 )
02231 {
02232 r = pk/qk;
02233 t = fabs( (ans - r)/r );
02234 ans = r;
02235 }
02236 else
02237 t = 1.0;
02238 pkm2 = pkm1;
02239 pkm1 = pk;
02240 qkm2 = qkm1;
02241 qkm1 = qk;
02242 if( fabs(pk) > BIG )
02243 {
02244 pkm2 /= BIG;
02245 pkm1 /= BIG;
02246 qkm2 /= BIG;
02247 qkm1 /= BIG;
02248 }
02249 }
02250 while( t > DBL_EPSILON );
02251
02252 ans *= exp( -x );
02253 return ans;
02254 }
02255 }
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354 static double erf_P[] = {
02355 2.46196981473530512524e-10,
02356 5.64189564831068821977e-1,
02357 7.46321056442269912687e0,
02358 4.86371970985681366614e1,
02359 1.96520832956077098242e2,
02360 5.26445194995477358631e2,
02361 9.34528527171957607540e2,
02362 1.02755188689515710272e3,
02363 5.57535335369399327526e2
02364 };
02365 static double erf_Q[] = {
02366
02367 1.32281951154744992508e1,
02368 8.67072140885989742329e1,
02369 3.54937778887819891062e2,
02370 9.75708501743205489753e2,
02371 1.82390916687909736289e3,
02372 2.24633760818710981792e3,
02373 1.65666309194161350182e3,
02374 5.57535340817727675546e2
02375 };
02376 static double erf_R[] = {
02377 5.64189583547755073984e-1,
02378 1.27536670759978104416e0,
02379 5.01905042251180477414e0,
02380 6.16021097993053585195e0,
02381 7.40974269950448939160e0,
02382 2.97886665372100240670e0
02383 };
02384 static double erf_S[] = {
02385
02386 2.26052863220117276590e0,
02387 9.39603524938001434673e0,
02388 1.20489539808096656605e1,
02389 1.70814450747565897222e1,
02390 9.60896809063285878198e0,
02391 3.36907645100081516050e0
02392 };
02393
02394
02395 #ifndef HAVE_ERFC
02396
02397 STATIC double expx2(double x, int sign);
02398
02399
02400
02401
02402 const bool lgUSE_EXPXSQ = true;
02403
02404 double erfc(double a)
02405 {
02406 double p,q,x,y,z;
02407
02408 DEBUG_ENTRY( "erfc()" );
02409
02410 x = abs(a);
02411
02412 if( x < 1.0 )
02413 return 1.0 - erf(a);
02414
02415 z = -a * a;
02416
02417 if( z < -MAXLOG )
02418 return ( a < 0.0 ) ? 2.0 : 0.0;
02419
02420 if( lgUSE_EXPXSQ )
02421
02422 z = expx2(a, -1);
02423 else
02424 z = exp(z);
02425
02426 if( x < 8.0 )
02427 {
02428 p = polevl( x, erf_P, 8 );
02429 q = p1evl( x, erf_Q, 8 );
02430 }
02431 else
02432 {
02433 p = polevl( x, erf_R, 5 );
02434 q = p1evl( x, erf_S, 6 );
02435 }
02436 y = (z * p)/q;
02437
02438 if( a < 0 )
02439 y = 2.0 - y;
02440
02441 if( y == 0.0 )
02442 return ( a < 0. ) ? 2.0 : 0.0;
02443
02444 return y;
02445 }
02446
02447 #endif
02448
02449
02450
02451
02452
02453
02454 double erfce(double x)
02455 {
02456 double p,q;
02457
02458 DEBUG_ENTRY( "erfce()" );
02459
02460 if( x < 8.0 )
02461 {
02462 p = polevl( x, erf_P, 8 );
02463 q = p1evl( x, erf_Q, 8 );
02464 }
02465 else
02466 {
02467 p = polevl( x, erf_R, 5 );
02468 q = p1evl( x, erf_S, 6 );
02469 }
02470 return p/q;
02471 }
02472
02473
02474 #ifndef HAVE_ERF
02475
02476 static double erf_T[] = {
02477 9.60497373987051638749e0,
02478 9.00260197203842689217e1,
02479 2.23200534594684319226e3,
02480 7.00332514112805075473e3,
02481 5.55923013010394962768e4
02482 };
02483 static double erf_U[] = {
02484
02485 3.35617141647503099647e1,
02486 5.21357949780152679795e2,
02487 4.59432382970980127987e3,
02488 2.26290000613890934246e4,
02489 4.92673942608635921086e4
02490 };
02491
02492 double erf(double x)
02493 {
02494 double y, z;
02495
02496 DEBUG_ENTRY( "erf()" );
02497
02498 if( abs(x) > 1.0 )
02499 return 1.0 - erfc(x);
02500 z = x * x;
02501 y = x * polevl( z, erf_T, 4 ) / p1evl( z, erf_U, 5 );
02502 return y;
02503 }
02504
02505 #endif
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546 #ifndef HAVE_ERFC
02547
02548 const double exp_M = 128.0;
02549 const double exp_MINV = .0078125;
02550
02551 STATIC double expx2(double x, int sign)
02552 {
02553 double u, u1, m, f;
02554
02555 DEBUG_ENTRY( "expx2()" );
02556
02557 x = abs(x);
02558 if( sign < 0 )
02559 x = -x;
02560
02561
02562
02563
02564 m = exp_MINV * floor(exp_M * x + 0.5);
02565 f = x - m;
02566
02567
02568 u = m * m;
02569 u1 = 2 * m * f + f * f;
02570
02571 if( sign < 0 )
02572 {
02573 u = -u;
02574 u1 = -u1;
02575 }
02576
02577 if( (u+u1) > MAXLOG )
02578 return DBL_MAX;
02579
02580
02581 u = exp(u) * exp(u1);
02582 return u;
02583 }
02584
02585 #endif
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638 inline double polevl(double x, const double coef[], int N)
02639 {
02640 double ans;
02641 int i;
02642 const double *p = coef;
02643
02644 ans = *p++;
02645 i = N;
02646
02647 do
02648 ans = ans * x + *p++;
02649 while( --i );
02650
02651 return ans;
02652 }
02653
02654
02655
02656
02657
02658
02659
02660 inline double p1evl(double x, const double coef[], int N)
02661 {
02662 double ans;
02663 const double *p = coef;
02664 int i;
02665
02666 ans = x + *p++;
02667 i = N-1;
02668
02669 do
02670 ans = ans * x + *p++;
02671 while( --i );
02672
02673 return ans;
02674 }
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734 inline double chbevl(double x, const double array[], int n)
02735 {
02736 double b0, b1, b2;
02737 const double *p = array;
02738 int i;
02739
02740 b0 = *p++;
02741 b1 = 0.0;
02742 i = n - 1;
02743
02744 do
02745 {
02746 b2 = b1;
02747 b1 = b0;
02748 b0 = x * b1 - b2 + *p++;
02749 }
02750 while( --i );
02751
02752 return 0.5*(b0-b2);
02753 }
02754
02755
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807
02808
02809
02810
02811
02812
02813
02814 static const int N = 624;
02815 static const int M = 397;
02816 static const unsigned long MATRIX_A = 0x9908b0dfUL;
02817 static const unsigned long UMASK = 0x80000000UL;
02818 static const unsigned long LMASK = 0x7fffffffUL;
02819 inline unsigned long MIXBITS(unsigned long u, unsigned long v)
02820 {
02821 return (u & UMASK) | (v & LMASK);
02822 }
02823 inline unsigned long TWIST(unsigned long u, unsigned long v)
02824 {
02825 return (MIXBITS(u,v) >> 1) ^ (v&1UL ? MATRIX_A : 0UL);
02826 }
02827
02828 static unsigned long state[N];
02829 static int nleft = 1;
02830 static int initf = 0;
02831 static unsigned long *nexxt;
02832
02833
02834 void init_genrand(unsigned long s)
02835 {
02836 int j;
02837 state[0]= s & 0xffffffffUL;
02838 for (j=1; j<N; j++) {
02839 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
02840
02841
02842
02843
02844 state[j] &= 0xffffffffUL;
02845 }
02846 nleft = 1; initf = 1;
02847 }
02848
02849
02850
02851
02852
02853 void init_by_array(unsigned long init_key[], int key_length)
02854 {
02855 int i, j, k;
02856 init_genrand(19650218UL);
02857 i=1; j=0;
02858 k = (N>key_length ? N : key_length);
02859 for (; k; k--) {
02860 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
02861 + init_key[j] + j;
02862 state[i] &= 0xffffffffUL;
02863 i++; j++;
02864 if (i>=N) { state[0] = state[N-1]; i=1; }
02865 if (j>=key_length) j=0;
02866 }
02867 for (k=N-1; k; k--) {
02868 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
02869 - i;
02870 state[i] &= 0xffffffffUL;
02871 i++;
02872 if (i>=N) { state[0] = state[N-1]; i=1; }
02873 }
02874
02875 state[0] = 0x80000000UL;
02876 nleft = 1; initf = 1;
02877 }
02878
02879 static void next_state()
02880 {
02881 unsigned long *p=state;
02882 int j;
02883
02884
02885
02886 if (initf==0) init_genrand(5489UL);
02887
02888 nleft = N;
02889 nexxt = state;
02890
02891 for (j=N-M+1; --j; p++)
02892 *p = p[M] ^ TWIST(p[0], p[1]);
02893
02894 for (j=M; --j; p++)
02895 *p = p[M-N] ^ TWIST(p[0], p[1]);
02896
02897 *p = p[M-N] ^ TWIST(p[0], state[0]);
02898 }
02899
02900
02901 unsigned long genrand_int32()
02902 {
02903 unsigned long y;
02904
02905 if (--nleft == 0) next_state();
02906 y = *nexxt++;
02907
02908
02909 y ^= (y >> 11);
02910 y ^= (y << 7) & 0x9d2c5680UL;
02911 y ^= (y << 15) & 0xefc60000UL;
02912 y ^= (y >> 18);
02913
02914 return y;
02915 }
02916
02917
02918 long genrand_int31()
02919 {
02920 unsigned long y;
02921
02922 if (--nleft == 0) next_state();
02923 y = *nexxt++;
02924
02925
02926 y ^= (y >> 11);
02927 y ^= (y << 7) & 0x9d2c5680UL;
02928 y ^= (y << 15) & 0xefc60000UL;
02929 y ^= (y >> 18);
02930
02931 return (long)(y>>1);
02932 }
02933
02934
02935 double genrand_real1()
02936 {
02937 unsigned long y;
02938
02939 if (--nleft == 0) next_state();
02940 y = *nexxt++;
02941
02942
02943 y ^= (y >> 11);
02944 y ^= (y << 7) & 0x9d2c5680UL;
02945 y ^= (y << 15) & 0xefc60000UL;
02946 y ^= (y >> 18);
02947
02948 return (double)y * (1.0/4294967295.0);
02949
02950 }
02951
02952
02953 double genrand_real2()
02954 {
02955 unsigned long y;
02956
02957 if (--nleft == 0) next_state();
02958 y = *nexxt++;
02959
02960
02961 y ^= (y >> 11);
02962 y ^= (y << 7) & 0x9d2c5680UL;
02963 y ^= (y << 15) & 0xefc60000UL;
02964 y ^= (y >> 18);
02965
02966 return (double)y * (1.0/4294967296.0);
02967
02968 }
02969
02970
02971 double genrand_real3()
02972 {
02973 unsigned long y;
02974
02975 if (--nleft == 0) next_state();
02976 y = *nexxt++;
02977
02978
02979 y ^= (y >> 11);
02980 y ^= (y << 7) & 0x9d2c5680UL;
02981 y ^= (y << 15) & 0xefc60000UL;
02982 y ^= (y >> 18);
02983
02984 return ((double)y + 0.5) * (1.0/4294967296.0);
02985
02986 }
02987
02988
02989 double genrand_res53()
02990 {
02991 unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
02992 return (a*67108864.0+b)*(1.0/9007199254740992.0);
02993 }
02994
02995
02996
02997
02998
02999
03000
03001
03002 const int N_DAWSON = 100;
03003
03004
03005
03006 static const double tbl_dawson[N_DAWSON+1] = {
03007 0.000000000000000000e+0,
03008 9.933599239785286114e-2,
03009 1.947510333680280496e-1,
03010 2.826316650213119286e-1,
03011 3.599434819348881042e-1,
03012 4.244363835020222959e-1,
03013 4.747632036629779306e-1,
03014 5.105040575592317787e-1,
03015 5.321017070563654290e-1,
03016 5.407243187262986750e-1,
03017 5.380794985696822201e-1,
03018 5.262066799705525356e-1,
03019 5.072734964077396141e-1,
03020 4.833975173848241360e-1,
03021 4.565072375268972572e-1,
03022 4.282490710853986254e-1,
03023 3.999398943230814126e-1,
03024 3.725593489740788250e-1,
03025 3.467727691148722451e-1,
03026 3.229743193228178382e-1,
03027 3.013403889237919660e-1,
03028 2.818849389255277938e-1,
03029 2.645107599508319576e-1,
03030 2.490529568377666955e-1,
03031 2.353130556638425762e-1,
03032 2.230837221674354811e-1,
03033 2.121651242424990041e-1,
03034 2.023745109105139857e-1,
03035 1.935507238593667923e-1,
03036 1.855552345354997718e-1,
03037 1.782710306105582873e-1,
03038 1.716003557161446928e-1,
03039 1.654619998786752031e-1,
03040 1.597885804744950500e-1,
03041 1.545240577369634525e-1,
03042 1.496215930807564847e-1,
03043 1.450417730540888593e-1,
03044 1.407511741154154125e-1,
03045 1.367212216746364963e-1,
03046 1.329272910810892667e-1,
03047 1.293480012360051155e-1,
03048 1.259646584343461329e-1,
03049 1.227608160065229225e-1,
03050 1.197219228088390280e-1,
03051 1.168350399532972540e-1,
03052 1.140886102268249801e-1,
03053 1.114722685321253072e-1,
03054 1.089766845891902214e-1,
03055 1.065934312832810744e-1,
03056 1.043148736220704706e-1,
03057 1.021340744242768354e-1,
03058 1.000447137201476355e-1,
03059 9.804101948507806734e-2,
03060 9.611770781195023073e-2,
03061 9.426993099823683440e-2,
03062 9.249323231075475996e-2,
03063 9.078350641561113352e-2,
03064 8.913696463869524546e-2,
03065 8.755010436413927499e-2,
03066 8.601968199264808016e-2,
03067 8.454268897454385223e-2,
03068 8.311633050835148859e-2,
03069 8.173800655824702378e-2,
03070 8.040529489538834118e-2,
03071 7.911593591113373223e-2,
03072 7.786781898606987138e-2,
03073 7.665897022891428948e-2,
03074 7.548754142476270211e-2,
03075 7.435180005364808165e-2,
03076 7.325012025863538754e-2,
03077 7.218097465823629202e-2,
03078 7.114292691123472774e-2,
03079 7.013462495342931107e-2,
03080 6.915479483562112754e-2,
03081 6.820223510065167384e-2,
03082 6.727581164463061598e-2,
03083 6.637445301385693290e-2,
03084 6.549714609447248675e-2,
03085 6.464293215671364666e-2,
03086 6.381090321984490301e-2,
03087 6.300019870755338791e-2,
03088 6.221000236682679049e-2,
03089 6.143953942619040819e-2,
03090 6.068807397169402555e-2,
03091 5.995490652126037542e-2,
03092 5.923937177997213955e-2,
03093 5.854083656061641403e-2,
03094 5.785869785535237086e-2,
03095 5.719238104574369009e-2,
03096 5.654133823962313062e-2,
03097 5.590504672435046070e-2,
03098 5.528300752700259690e-2,
03099 5.467474407290986634e-2,
03100 5.407980093473671650e-2,
03101 5.349774266500934015e-2,
03102 5.292815270562564644e-2,
03103 5.237063236845275277e-2,
03104 5.182479988163068060e-2,
03105 5.129028949666435102e-2,
03106 5.076675065180469937e-2,
03107 5.025384718759852810e-2
03108 };
03109
03110
03111
03112
03113
03114
03115
03116
03117
03118
03119
03120
03121
03122
03123
03124
03125
03126
03127
03128
03129
03130
03131
03132
03133
03134
03135
03136
03137 inline double dawson(double x, int order)
03138 {
03139 double x10 = x*10.;
03140 if( order == 1 )
03141 {
03142 int ind = min(max(int(x10),0),N_DAWSON-1);
03143 double p = x10 - double(ind);
03144 return tbl_dawson[ind] + p*(tbl_dawson[ind+1]-tbl_dawson[ind]);
03145 }
03146 else if( order == 3 )
03147 {
03148 int ind = min(max(int(x10-1.),0),N_DAWSON-3);
03149 double p = x10 - double(ind+1);
03150 double pm1 = p-1.;
03151 double pm2 = p-2.;
03152 double pp1 = p+1.;
03153
03154 return p*pm1*(pp1*tbl_dawson[ind+3]-pm2*tbl_dawson[ind])/6. +
03155 pm2*pp1*(pm1*tbl_dawson[ind+1]-p*tbl_dawson[ind+2])/2.;
03156 }
03157 else
03158 {
03159 TotalInsanity();
03160 }
03161 }
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171 realnum FastVoigtH(realnum a, realnum v)
03172 {
03173 DEBUG_ENTRY( "FastVoigtH()" );
03174
03175 ASSERT( a <= 0.101f );
03176
03177
03178
03179
03180
03181
03182
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198
03199 v = abs(v);
03200
03201 if( v > 9.f )
03202 {
03203
03204
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215 realnum vm2 = 1.f/pow2(v);
03216 return a*vm2/realnum(SQRTPI)*(((105.f/8.f*vm2 + 15.f/4.f)*vm2 + 3.f/2.f)*vm2 + 1.f);
03217 }
03218 else
03219 {
03220 realnum v2 = pow2(v);
03221
03222
03223
03224
03225 realnum emv2 = realnum(dsexp(v2));
03226 int order = ( a > 0.003f || v > 1.5f ) ? 3 : 1;
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240 return emv2*(1.f-pow2(a)*(2.f*v2-1.f)) +
03241 2.f*a/realnum(SQRTPI)*(-1.f + 2.f*v*realnum(dawson(v,order)));
03242 }
03243 }
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257 void humlik(int n, const realnum x[], realnum y, realnum k[])
03258 {
03259 DEBUG_ENTRY( "humlik()" );
03260
03261
03262
03263
03264
03265
03266
03267
03268
03269
03270 static const double c[6] = { 1.0117281, -.75197147, .012557727, .010022008, -2.4206814e-4, 5.0084806e-7 };
03271 static const double s[6] = { 1.393237, .23115241, -.15535147, .0062183662, 9.1908299e-5, -6.2752596e-7 };
03272 static const double t[6] = { .31424038, .94778839, 1.5976826, 2.2795071, 3.020637, 3.8897249 };
03273
03274 static const double R0 = 146.7, R1 = 14.67;
03275
03276
03277 double d, ki;
03278 int i, j;
03279 double a0, d0, e0, d2, e2, h0, e4, h2, h4, h6, p0,
03280 p2, p4, p6, p8, z0, z2, z4, z6, z8, mf[6], pf[6],
03281 mq[6], yf, pq[6], xm[6], ym[6], xp[6], xq, yq, yp[6];
03282 bool rg1, rg2, rg3;
03283 double abx, ypy0, xlim0, xlim1, xlim2, xlim3, xlim4, ypy0q, yrrtpi;
03284
03285 rg1 = rg2 = rg3 = true;
03286
03287 z0 = z2 = z4 = z6 = z8 = 0.;
03288 p0 = p2 = p4 = p6 = p8 = 0.;
03289 h0 = h2 = h4 = h6 = 0.;
03290 a0 = d0 = d2 = e0 = e2 = e4 = 0.;
03291
03292 yq = y * y;
03293 yrrtpi = y * .56418958;
03294
03295 xlim0 = R0 - y;
03296 xlim1 = R1 - y;
03297 xlim3 = y * 3.097 - .45;
03298 xlim2 = 6.8 - y;
03299 xlim4 = y * 18.1 + 1.65;
03300 if (y <= 1e-6)
03301 {
03302 xlim1 = xlim0;
03303 xlim2 = xlim0;
03304 }
03305
03306 for (i = 0; i < n; ++i)
03307 {
03308 abx = fabs(x[i]);
03309 xq = abx * abx;
03310 if (abx > xlim0)
03311 {
03312 k[i] = realnum(yrrtpi / (xq + yq));
03313 }
03314 else if (abx > xlim1)
03315 {
03316 if (rg1)
03317 {
03318
03319 rg1 = false;
03320 a0 = yq + .5;
03321 d0 = a0 * a0;
03322 d2 = yq + yq - 1.;
03323 }
03324 d = .56418958 / (d0 + xq * (d2 + xq));
03325 k[i] = realnum(d * y * (a0 + xq));
03326 }
03327 else if (abx > xlim2)
03328 {
03329 if (rg2)
03330 {
03331 rg2 = false;
03332 h0 = yq * (yq * (yq * (yq + 6.) + 10.5) + 4.5) + .5625;
03333 h2 = yq * (yq * (yq * 4. + 6.) + 9.) - 4.5;
03334 h4 = 10.5 - yq * (6. - yq * 6.);
03335 h6 = yq * 4. - 6.;
03336 e0 = yq * (yq * (yq + 5.5) + 8.25) + 1.875;
03337 e2 = yq * (yq * 3. + 1.) + 5.25;
03338 e4 = h6 * .75;
03339 }
03340 d = .56418958 / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
03341 k[i] = realnum(d * y * (e0 + xq * (e2 + xq * (e4 + xq))));
03342 }
03343 else if (abx < xlim3)
03344 {
03345 if (rg3)
03346 {
03347
03348 rg3 = false;
03349 z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 13.3988) +
03350 88.26741) + 369.1989) + 1074.409) + 2256.981) + 3447.629) +
03351 3764.966) + 2802.87) + 1280.829) + 272.1014;
03352 z2 = y * (y * (y * (y * (y * (y * (y * (y * 5. + 53.59518) +
03353 266.2987) + 793.4273) + 1549.675) + 2037.31) + 1758.336) +
03354 902.3066) + 211.678;
03355 z4 = y * (y * (y * (y * (y * (y * 10. + 80.39278) + 269.2916) +
03356 479.2576) + 497.3014) + 308.1852) + 78.86585;
03357 z6 = y * (y * (y * (y * 10. + 53.59518) + 92.75679) + 55.02933) +
03358 22.03523;
03359 z8 = y * (y * 5. + 13.3988) + 1.49646;
03360 p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * .3183291 +
03361 4.264678) + 27.93941) + 115.3772) + 328.2151) + 662.8097) +
03362 946.897) + 919.4955) + 549.3954) + 153.5168;
03363 p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163 + 12.79458) +
03364 56.81652) + 139.4665) + 189.773) + 124.5975) - 1.322256) -
03365 34.16955;
03366 p4 = y * (y * (y * (y * (y * 1.9099744 + 12.79568) + 29.81482) +
03367 24.01655) + 10.46332) + 2.584042;
03368 p6 = y * (y * (y * 1.273316 + 4.266322) + .9377051) - .07272979;
03369 p8 = y * .3183291 + 5.480304e-4;
03370 }
03371 d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
03372 k[i] = realnum(d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8)))));
03373 }
03374 else
03375 {
03376 ypy0 = y + 1.5;
03377 ypy0q = ypy0 * ypy0;
03378 ki = 0.;
03379 for (j = 0; j <= 5; ++j)
03380 {
03381 d = x[i] - t[j];
03382 mq[j] = d * d;
03383 mf[j] = 1. / (mq[j] + ypy0q);
03384 xm[j] = mf[j] * d;
03385 ym[j] = mf[j] * ypy0;
03386 d = x[i] + t[j];
03387 pq[j] = d * d;
03388 pf[j] = 1. / (pq[j] + ypy0q);
03389 xp[j] = pf[j] * d;
03390 yp[j] = pf[j] * ypy0;
03391 }
03392 if (abx <= xlim4)
03393 {
03394 for (j = 0; j <= 5; ++j)
03395 {
03396 ki = ki + c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
03397 }
03398 }
03399 else
03400 {
03401 yf = y + 3.;
03402 for (j = 0; j <= 5; ++j)
03403 {
03404 ki = ki + (c[j] * (mq[j] * mf[j] - ym[j] * 1.5) + s[j] * yf * xm[j]) / (mq[j] + 2.25) +
03405 (c[j] * (pq[j] * pf[j] - yp[j] * 1.5) - s[j] * yf * xp[j]) / (pq[j] + 2.25);
03406 }
03407 ki = y * ki + exp(-xq);
03408 }
03409 k[i] = realnum(ki);
03410 }
03411 }
03412 }
03413
03414
03415
03416
03417
03418 STATIC uint32 MD5swap( uint32 word );
03419 STATIC void MD5_Transform (uint32 *digest, const uint32 *in);
03420
03421
03422
03423
03424
03425 string MD5file(const char* fnam, access_scheme scheme)
03426 {
03427 DEBUG_ENTRY( "MD5file()" );
03428
03429 fstream ioFile;
03430 open_data( ioFile, fnam, mode_r, scheme );
03431
03432 char c;
03433 string content;
03434 while( ioFile.get(c) )
03435 content += c;
03436
03437 return MD5string( content );
03438 }
03439
03440
03441
03442
03443
03444 string MD5datafile(const char* fnam, access_scheme scheme)
03445 {
03446 DEBUG_ENTRY( "MD5datafile()" );
03447
03448 fstream ioFile;
03449 open_data( ioFile, fnam, mode_r, scheme );
03450
03451 string line, content;
03452 while( getline( ioFile, line ) )
03453 if( line[0] != '#' )
03454 content += line;
03455
03456 return MD5string( content );
03457 }
03458
03459
03460
03461 string MD5string(const string& str)
03462 {
03463 DEBUG_ENTRY( "MD5string()" );
03464
03465 uint32 state[4];
03466
03467 state[0] = 0x67452301L;
03468 state[1] = 0xefcdab89L;
03469 state[2] = 0x98badcfeL;
03470 state[3] = 0x10325476L;
03471
03472 string lstr = str;
03473
03474
03475 uint32 bytes = str.length()%64;
03476 uint32 padlen = ( bytes >= 56 ) ? 64 + 56 - bytes : 56 - bytes;
03477 lstr += '\x80';
03478 for( uint32 i=1; i < padlen; ++i )
03479 lstr += '\0';
03480
03481 ASSERT( lstr.length()%64 == 56 );
03482
03483 uint32 i;
03484 union {
03485 uint32 in[16];
03486 unsigned char chr[64];
03487 } u;
03488
03489 for( i=0; i < lstr.length()/64; ++i )
03490 {
03491 for( uint32 j=0; j < 64; ++j )
03492 {
03493 if( cpu.i().little_endian() )
03494 u.chr[j] = lstr[i*64+j];
03495 else
03496 {
03497 uint32 jr = j%4;
03498 uint32 j4 = j-jr;
03499 u.chr[j4+3-jr] = lstr[i*64+j];
03500 }
03501 }
03502 MD5_Transform( state, u.in );
03503 }
03504 for( uint32 j=0; j < 56; ++j )
03505 {
03506 if( cpu.i().little_endian() )
03507 u.chr[j] = lstr[i*64+j];
03508 else
03509 {
03510 uint32 jr = j%4;
03511 uint32 j4 = j-jr;
03512 u.chr[j4+3-jr] = lstr[i*64+j];
03513 }
03514 }
03515
03516 u.in[14] = (str.length()<<3)&0xffffffff;
03517 u.in[15] = (str.length()>>29)&0xffffffff;
03518
03519 MD5_Transform( state, u.in );
03520
03521 ostringstream hash;
03522 for( uint32 i=0; i < 4; ++i )
03523 hash << hex << setfill('0') << setw(8) << MD5swap(state[i]);
03524
03525 return hash.str();
03526 }
03527
03528 STATIC uint32 MD5swap( uint32 word )
03529 {
03530 DEBUG_ENTRY( "MD5swap()" );
03531
03532 union {
03533 uint32 word;
03534 unsigned char byte[4];
03535 } ui, uo;
03536
03537 ui.word = word;
03538 uo.byte[0] = ui.byte[3];
03539 uo.byte[1] = ui.byte[2];
03540 uo.byte[2] = ui.byte[1];
03541 uo.byte[3] = ui.byte[0];
03542
03543 return uo.word;
03544 }
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564
03565
03566
03567
03568
03569
03570
03571
03572
03573
03574
03575
03576
03577
03578
03579
03580
03581
03582
03583
03584
03585
03586
03587
03588
03589
03590
03591
03592
03593
03594
03595
03596
03597
03598
03599
03600
03601
03602
03603
03604
03605
03606
03607
03608
03609
03610
03611
03612
03613
03614
03615
03616
03617
03618
03619 inline uint32 rotlFixed(uint32 x, unsigned int y)
03620 {
03621 return uint32((x<<y) | (x>>(32-y)));
03622 }
03623
03624 STATIC void MD5_Transform (uint32 *digest, const uint32 *in)
03625 {
03626 DEBUG_ENTRY( "MD5_Transform()" );
03627
03628 #define F1(x, y, z) (z ^ (x & (y ^ z)))
03629 #define F2(x, y, z) F1(z, x, y)
03630 #define F3(x, y, z) (x ^ y ^ z)
03631 #define F4(x, y, z) (y ^ (x | ~z))
03632
03633 #define MD5STEP(f, w, x, y, z, data, s) \
03634 w = rotlFixed(w + f(x, y, z) + data, s) + x
03635
03636 uint32 a, b, c, d;
03637
03638 a=digest[0];
03639 b=digest[1];
03640 c=digest[2];
03641 d=digest[3];
03642
03643 MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
03644 MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
03645 MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
03646 MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
03647 MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
03648 MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
03649 MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
03650 MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
03651 MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
03652 MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
03653 MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
03654 MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
03655 MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
03656 MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
03657 MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
03658 MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
03659
03660 MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
03661 MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
03662 MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
03663 MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
03664 MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
03665 MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
03666 MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
03667 MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
03668 MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
03669 MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
03670 MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
03671 MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
03672 MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
03673 MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
03674 MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
03675 MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
03676
03677 MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
03678 MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
03679 MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
03680 MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
03681 MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
03682 MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
03683 MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
03684 MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
03685 MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
03686 MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
03687 MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
03688 MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
03689 MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
03690 MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
03691 MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
03692 MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
03693
03694 MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
03695 MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
03696 MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
03697 MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
03698 MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
03699 MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
03700 MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
03701 MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
03702 MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
03703 MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
03704 MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
03705 MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
03706 MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
03707 MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
03708 MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
03709 MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
03710
03711 digest[0] += a;
03712 digest[1] += b;
03713 digest[2] += c;
03714 digest[3] += d;
03715 }
03716
03717
03718
03719