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
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 bool linfit(
00046 long n,
00047 const double xorg[],
00048 const double yorg[],
00049 double &a,
00050 double &siga,
00051 double &b,
00052 double &sigb
00053 )
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 DEBUG_ENTRY( "linfit()" );
00093
00094 ASSERT( n >= 2 );
00095
00096 valarray<double> x(n);
00097 valarray<double> y(n);
00098
00099 for( long i=0; i < n; i++ )
00100 {
00101 x[i] = xorg[i];
00102 y[i] = yorg[i];
00103 }
00104
00105
00106 a = 0.0;
00107 siga = 0.0;
00108 b = 0.0;
00109 sigb = 0.0;
00110
00111
00112 double s1 = 0.0;
00113 double s2 = 0.0;
00114 for( long i=0; i < n; i++ )
00115 {
00116 s1 += x[i];
00117 s2 += y[i];
00118 }
00119 double rn = (double)n;
00120 double xavg = s1/rn;
00121 double yavg = s2/rn;
00122 double sxx = 0.0;
00123 double sxy = 0.0;
00124 for( long i=0; i < n; i++ )
00125 {
00126 x[i] -= xavg;
00127 y[i] -= yavg;
00128 sxx += pow2(x[i]);
00129 sxy += x[i]*y[i];
00130 }
00131
00132 if( pow2(sxx) == 0.0 )
00133 {
00134 return true;
00135 }
00136
00137
00138 b = sxy/sxx;
00139
00140
00141 a = yavg - b*xavg;
00142
00143
00144 double sum1 = 0.0;
00145 for( long i=0; i < n; i++ )
00146 sum1 += pow2(x[i]*(y[i] - b*x[i]));
00147
00148
00149 sigb = sum1/pow2(sxx);
00150
00151
00152 for( long i=0; i < n; i++ )
00153 siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
00154
00155
00156 sigb = sqrt(sigb);
00157 siga = sqrt(siga)/rn;
00158
00159
00160 for( long i=0; i < n; i++ )
00161 {
00162 x[i] += xavg;
00163 y[i] += yavg;
00164 }
00165 return false;
00166 }
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 static const double pre_factorial[NPRE_FACTORIAL] =
00181 {
00182 1.00000000000000000000e+00,
00183 1.00000000000000000000e+00,
00184 2.00000000000000000000e+00,
00185 6.00000000000000000000e+00,
00186 2.40000000000000000000e+01,
00187 1.20000000000000000000e+02,
00188 7.20000000000000000000e+02,
00189 5.04000000000000000000e+03,
00190 4.03200000000000000000e+04,
00191 3.62880000000000000000e+05,
00192 3.62880000000000000000e+06,
00193 3.99168000000000000000e+07,
00194 4.79001600000000000000e+08,
00195 6.22702080000000000000e+09,
00196 8.71782912000000000000e+10,
00197 1.30767436800000000000e+12,
00198 2.09227898880000000000e+13,
00199 3.55687428096000000000e+14,
00200 6.40237370572800000000e+15,
00201 1.21645100408832000000e+17,
00202 2.43290200817664000000e+18,
00203 5.10909421717094400000e+19,
00204 1.12400072777760768000e+21,
00205 2.58520167388849766400e+22,
00206 6.20448401733239439360e+23,
00207 1.55112100433309859840e+25,
00208 4.03291461126605635592e+26,
00209 1.08888694504183521614e+28,
00210 3.04888344611713860511e+29,
00211 8.84176199373970195470e+30,
00212 2.65252859812191058647e+32,
00213 8.22283865417792281807e+33,
00214 2.63130836933693530178e+35,
00215 8.68331761881188649615e+36,
00216 2.95232799039604140861e+38,
00217 1.03331479663861449300e+40,
00218 3.71993326789901217463e+41,
00219 1.37637530912263450457e+43,
00220 5.23022617466601111726e+44,
00221 2.03978820811974433568e+46,
00222 8.15915283247897734264e+47,
00223 3.34525266131638071044e+49,
00224 1.40500611775287989839e+51,
00225 6.04152630633738356321e+52,
00226 2.65827157478844876773e+54,
00227 1.19622220865480194551e+56,
00228 5.50262215981208894940e+57,
00229 2.58623241511168180614e+59,
00230 1.24139155925360726691e+61,
00231 6.08281864034267560801e+62,
00232 3.04140932017133780398e+64,
00233 1.55111875328738227999e+66,
00234 8.06581751709438785585e+67,
00235 4.27488328406002556374e+69,
00236 2.30843697339241380441e+71,
00237 1.26964033536582759243e+73,
00238 7.10998587804863451749e+74,
00239 4.05269195048772167487e+76,
00240 2.35056133128287857145e+78,
00241 1.38683118545689835713e+80,
00242 8.32098711274139014271e+81,
00243 5.07580213877224798711e+83,
00244 3.14699732603879375200e+85,
00245 1.98260831540444006372e+87,
00246 1.26886932185884164078e+89,
00247 8.24765059208247066512e+90,
00248 5.44344939077443063905e+92,
00249 3.64711109181886852801e+94,
00250 2.48003554243683059915e+96,
00251 1.71122452428141311337e+98,
00252 1.19785716699698917933e+100,
00253 8.50478588567862317347e+101,
00254 6.12344583768860868500e+103,
00255 4.47011546151268434004e+105,
00256 3.30788544151938641157e+107,
00257 2.48091408113953980872e+109,
00258 1.88549470166605025458e+111,
00259 1.45183092028285869606e+113,
00260 1.13242811782062978295e+115,
00261 8.94618213078297528506e+116,
00262 7.15694570462638022794e+118,
00263 5.79712602074736798470e+120,
00264 4.75364333701284174746e+122,
00265 3.94552396972065865030e+124,
00266 3.31424013456535326627e+126,
00267 2.81710411438055027626e+128,
00268 2.42270953836727323750e+130,
00269 2.10775729837952771662e+132,
00270 1.85482642257398439069e+134,
00271 1.65079551609084610774e+136,
00272 1.48571596448176149700e+138,
00273 1.35200152767840296226e+140,
00274 1.24384140546413072522e+142,
00275 1.15677250708164157442e+144,
00276 1.08736615665674307994e+146,
00277 1.03299784882390592592e+148,
00278 9.91677934870949688836e+149,
00279 9.61927596824821198159e+151,
00280 9.42689044888324774164e+153,
00281 9.33262154439441526381e+155,
00282 9.33262154439441526388e+157,
00283 9.42594775983835941673e+159,
00284 9.61446671503512660515e+161,
00285 9.90290071648618040340e+163,
00286 1.02990167451456276198e+166,
00287 1.08139675824029090008e+168,
00288 1.14628056373470835406e+170,
00289 1.22652020319613793888e+172,
00290 1.32464181945182897396e+174,
00291 1.44385958320249358163e+176,
00292 1.58824554152274293982e+178,
00293 1.76295255109024466316e+180,
00294 1.97450685722107402277e+182,
00295 2.23119274865981364576e+184,
00296 2.54355973347218755612e+186,
00297 2.92509369349301568964e+188,
00298 3.39310868445189820004e+190,
00299 3.96993716080872089396e+192,
00300 4.68452584975429065488e+194,
00301 5.57458576120760587943e+196,
00302 6.68950291344912705515e+198,
00303 8.09429852527344373681e+200,
00304 9.87504420083360135884e+202,
00305 1.21463043670253296712e+205,
00306 1.50614174151114087918e+207,
00307 1.88267717688892609901e+209,
00308 2.37217324288004688470e+211,
00309 3.01266001845765954361e+213,
00310 3.85620482362580421582e+215,
00311 4.97450422247728743840e+217,
00312 6.46685548922047366972e+219,
00313 8.47158069087882050755e+221,
00314 1.11824865119600430699e+224,
00315 1.48727070609068572828e+226,
00316 1.99294274616151887582e+228,
00317 2.69047270731805048244e+230,
00318 3.65904288195254865604e+232,
00319 5.01288874827499165889e+234,
00320 6.91778647261948848943e+236,
00321 9.61572319694108900019e+238,
00322 1.34620124757175246000e+241,
00323 1.89814375907617096864e+243,
00324 2.69536413788816277557e+245,
00325 3.85437071718007276916e+247,
00326 5.55029383273930478744e+249,
00327 8.04792605747199194159e+251,
00328 1.17499720439091082343e+254,
00329 1.72724589045463891049e+256,
00330 2.55632391787286558753e+258,
00331 3.80892263763056972532e+260,
00332 5.71338395644585458806e+262,
00333 8.62720977423324042775e+264,
00334 1.31133588568345254503e+267,
00335 2.00634390509568239384e+269,
00336 3.08976961384735088657e+271,
00337 4.78914290146339387432e+273,
00338 7.47106292628289444390e+275,
00339 1.17295687942641442768e+278,
00340 1.85327186949373479574e+280,
00341 2.94670227249503832518e+282,
00342 4.71472363599206132029e+284,
00343 7.59070505394721872577e+286,
00344 1.22969421873944943358e+289,
00345 2.00440157654530257674e+291,
00346 3.28721858553429622598e+293,
00347 5.42391066613158877297e+295,
00348 9.00369170577843736335e+297,
00349 1.50361651486499903974e+300,
00350 2.52607574497319838672e+302,
00351 4.26906800900470527345e+304,
00352 7.25741561530799896496e+306
00353 };
00354
00355 double factorial( long n )
00356 {
00357 DEBUG_ENTRY( "factorial()" );
00358
00359 if( n < 0 || n >= NPRE_FACTORIAL )
00360 {
00361 fprintf( ioQQQ, "factorial: domain error\n" );
00362 cdEXIT(EXIT_FAILURE);
00363 }
00364
00365 return pre_factorial[n];
00366 }
00367
00368
00369
00370 class t_lfact : public Singleton<t_lfact>
00371 {
00372 friend class Singleton<t_lfact>;
00373 protected:
00374 t_lfact()
00375 {
00376 p_lf.reserve( 512 );
00377 p_lf.push_back( 0. );
00378 p_lf.push_back( 0. );
00379 }
00380 private:
00381 vector<double> p_lf;
00382 public:
00383 double get_lfact( unsigned long n )
00384 {
00385 if( n < p_lf.size() )
00386 {
00387 return p_lf[n];
00388 }
00389 else
00390 {
00391 for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ )
00392 p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) );
00393 return p_lf[n];
00394 }
00395 }
00396 };
00397
00398 double lfactorial( long n )
00399 {
00400
00401
00402
00403
00404
00405
00406
00407
00408 DEBUG_ENTRY( "lfactorial()" );
00409
00410 if( n < 0 )
00411 {
00412 fprintf( ioQQQ, "lfactorial: domain error\n" );
00413 cdEXIT(EXIT_FAILURE);
00414 }
00415
00416 return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 complex<double> cdgamma(complex<double> x)
00432 {
00433 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
00434
00435 DEBUG_ENTRY( "cdgamma()" );
00436
00437 xr = x.real();
00438 xi = x.imag();
00439 if(xr < 0)
00440 {
00441 wr = 1. - xr;
00442 wi = -xi;
00443 }
00444 else
00445 {
00446 wr = xr;
00447 wi = xi;
00448 }
00449 ur = wr + 6.00009857740312429;
00450 vr = ur * (wr + 4.99999857982434025) - wi * wi;
00451 vi = wi * (wr + 4.99999857982434025) + ur * wi;
00452 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
00453 0.293729529320536228;
00454 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
00455 ur = vr * (wr + 4.00000003016801681) - vi * wi;
00456 ui = vi * (wr + 4.00000003016801681) + vr * wi;
00457 vr = ur * (wr + 2.99999999944915534) - ui * wi;
00458 vi = ui * (wr + 2.99999999944915534) + ur * wi;
00459 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
00460 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
00461 ur = vr * (wr + 2.00000000000603851) - vi * wi;
00462 ui = vi * (wr + 2.00000000000603851) + vr * wi;
00463 vr = ur * (wr + 0.999999999999975753) - ui * wi;
00464 vi = ui * (wr + 0.999999999999975753) + ur * wi;
00465 yr += ur * 10.5400280458730808 + vr;
00466 yi += ui * 10.5400280458730808 + vi;
00467 ur = vr * wr - vi * wi;
00468 ui = vi * wr + vr * wi;
00469 t = ur * ur + ui * ui;
00470 vr = yr * ur + yi * ui + t * 0.0327673720261526849;
00471 vi = yi * ur - yr * ui;
00472 yr = wr + 7.31790632447016203;
00473 ur = log(yr * yr + wi * wi) * 0.5 - 1;
00474 ui = atan2(wi, yr);
00475 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
00476 yi = ui * (wr - 0.5) + ur * wi;
00477 ur = yr * cos(yi);
00478 ui = yr * sin(yi);
00479 yr = ur * vr - ui * vi;
00480 yi = ui * vr + ur * vi;
00481 if(xr < 0)
00482 {
00483 wr = xr * 3.14159265358979324;
00484 wi = exp(xi * 3.14159265358979324);
00485 vi = 1 / wi;
00486 ur = (vi + wi) * sin(wr);
00487 ui = (vi - wi) * cos(wr);
00488 vr = ur * yr + ui * yi;
00489 vi = ui * yr - ur * yi;
00490 ur = 6.2831853071795862 / (vr * vr + vi * vi);
00491 yr = ur * vr;
00492 yi = ur * vi;
00493 }
00494 return complex<double>( yr, yi );
00495 }
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 static const double b0_PP[7] = {
00616 7.96936729297347051624e-4,
00617 8.28352392107440799803e-2,
00618 1.23953371646414299388e0,
00619 5.44725003058768775090e0,
00620 8.74716500199817011941e0,
00621 5.30324038235394892183e0,
00622 9.99999999999999997821e-1,
00623 };
00624
00625 static const double b0_PQ[7] = {
00626 9.24408810558863637013e-4,
00627 8.56288474354474431428e-2,
00628 1.25352743901058953537e0,
00629 5.47097740330417105182e0,
00630 8.76190883237069594232e0,
00631 5.30605288235394617618e0,
00632 1.00000000000000000218e0,
00633 };
00634
00635 static const double b0_QP[8] = {
00636 -1.13663838898469149931e-2,
00637 -1.28252718670509318512e0,
00638 -1.95539544257735972385e1,
00639 -9.32060152123768231369e1,
00640 -1.77681167980488050595e2,
00641 -1.47077505154951170175e2,
00642 -5.14105326766599330220e1,
00643 -6.05014350600728481186e0,
00644 };
00645
00646 static const double b0_QQ[7] = {
00647
00648 6.43178256118178023184e1,
00649 8.56430025976980587198e2,
00650 3.88240183605401609683e3,
00651 7.24046774195652478189e3,
00652 5.93072701187316984827e3,
00653 2.06209331660327847417e3,
00654 2.42005740240291393179e2,
00655 };
00656
00657 static const double b0_YP[8] = {
00658 1.55924367855235737965e4,
00659 -1.46639295903971606143e7,
00660 5.43526477051876500413e9,
00661 -9.82136065717911466409e11,
00662 8.75906394395366999549e13,
00663 -3.46628303384729719441e15,
00664 4.42733268572569800351e16,
00665 -1.84950800436986690637e16,
00666 };
00667
00668 static const double b0_YQ[7] = {
00669
00670 1.04128353664259848412e3,
00671 6.26107330137134956842e5,
00672 2.68919633393814121987e8,
00673 8.64002487103935000337e10,
00674 2.02979612750105546709e13,
00675 3.17157752842975028269e15,
00676 2.50596256172653059228e17,
00677 };
00678
00679
00680 static const double DR1 = 5.78318596294678452118e0;
00681
00682 static const double DR2 = 3.04712623436620863991e1;
00683
00684 static double b0_RP[4] = {
00685 -4.79443220978201773821e9,
00686 1.95617491946556577543e12,
00687 -2.49248344360967716204e14,
00688 9.70862251047306323952e15,
00689 };
00690
00691 static double b0_RQ[8] = {
00692
00693 4.99563147152651017219e2,
00694 1.73785401676374683123e5,
00695 4.84409658339962045305e7,
00696 1.11855537045356834862e10,
00697 2.11277520115489217587e12,
00698 3.10518229857422583814e14,
00699 3.18121955943204943306e16,
00700 1.71086294081043136091e18,
00701 };
00702
00703 static const double TWOOPI = 2./PI;
00704 static const double SQ2OPI = sqrt(2./PI);
00705 static const double PIO4 = PI/4.;
00706
00707 double bessel_j0(double x)
00708 {
00709 double w, z, p, q, xn;
00710
00711 DEBUG_ENTRY( "bessel_j0()" );
00712
00713 if( x < 0 )
00714 x = -x;
00715
00716 if( x <= 5.0 )
00717 {
00718 z = x * x;
00719 if( x < 1.0e-5 )
00720 return 1.0 - z/4.0;
00721
00722 p = (z - DR1) * (z - DR2);
00723 p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
00724 return p;
00725 }
00726
00727 w = 5.0/x;
00728 q = 25.0/(x*x);
00729 p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
00730 q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
00731 xn = x - PIO4;
00732 p = p * cos(xn) - w * q * sin(xn);
00733 return p * SQ2OPI / sqrt(x);
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745 double bessel_y0(double x)
00746 {
00747 double w, z, p, q, xn;
00748
00749 DEBUG_ENTRY( "bessel_y0()" );
00750
00751 if( x <= 5.0 )
00752 {
00753 if( x <= 0.0 )
00754 {
00755 fprintf( ioQQQ, "bessel_y0: domain error\n" );
00756 cdEXIT(EXIT_FAILURE);
00757 }
00758 z = x * x;
00759 w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
00760 w += TWOOPI * log(x) * bessel_j0(x);
00761 return w;
00762 }
00763
00764 w = 5.0/x;
00765 z = 25.0 / (x * x);
00766 p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
00767 q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
00768 xn = x - PIO4;
00769 p = p * sin(xn) + w * q * cos(xn);
00770 return p * SQ2OPI / sqrt(x);
00771 }
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 static const double b1_RP[4] = {
00852 -8.99971225705559398224e8,
00853 4.52228297998194034323e11,
00854 -7.27494245221818276015e13,
00855 3.68295732863852883286e15,
00856 };
00857
00858 static const double b1_RQ[8] = {
00859
00860 6.20836478118054335476e2,
00861 2.56987256757748830383e5,
00862 8.35146791431949253037e7,
00863 2.21511595479792499675e10,
00864 4.74914122079991414898e12,
00865 7.84369607876235854894e14,
00866 8.95222336184627338078e16,
00867 5.32278620332680085395e18,
00868 };
00869
00870 static const double b1_PP[7] = {
00871 7.62125616208173112003e-4,
00872 7.31397056940917570436e-2,
00873 1.12719608129684925192e0,
00874 5.11207951146807644818e0,
00875 8.42404590141772420927e0,
00876 5.21451598682361504063e0,
00877 1.00000000000000000254e0,
00878 };
00879
00880 static const double b1_PQ[7] = {
00881 5.71323128072548699714e-4,
00882 6.88455908754495404082e-2,
00883 1.10514232634061696926e0,
00884 5.07386386128601488557e0,
00885 8.39985554327604159757e0,
00886 5.20982848682361821619e0,
00887 9.99999999999999997461e-1,
00888 };
00889
00890 static const double b1_QP[8] = {
00891 5.10862594750176621635e-2,
00892 4.98213872951233449420e0,
00893 7.58238284132545283818e1,
00894 3.66779609360150777800e2,
00895 7.10856304998926107277e2,
00896 5.97489612400613639965e2,
00897 2.11688757100572135698e2,
00898 2.52070205858023719784e1,
00899 };
00900
00901 static const double b1_QQ[7] = {
00902
00903 7.42373277035675149943e1,
00904 1.05644886038262816351e3,
00905 4.98641058337653607651e3,
00906 9.56231892404756170795e3,
00907 7.99704160447350683650e3,
00908 2.82619278517639096600e3,
00909 3.36093607810698293419e2,
00910 };
00911
00912 static const double b1_YP[6] = {
00913 1.26320474790178026440e9,
00914 -6.47355876379160291031e11,
00915 1.14509511541823727583e14,
00916 -8.12770255501325109621e15,
00917 2.02439475713594898196e17,
00918 -7.78877196265950026825e17,
00919 };
00920
00921 static const double b1_YQ[8] = {
00922
00923 5.94301592346128195359E2,
00924 2.35564092943068577943E5,
00925 7.34811944459721705660E7,
00926 1.87601316108706159478E10,
00927 3.88231277496238566008E12,
00928 6.20557727146953693363E14,
00929 6.87141087355300489866E16,
00930 3.97270608116560655612E18,
00931 };
00932
00933 static const double Z1 = 1.46819706421238932572E1;
00934 static const double Z2 = 4.92184563216946036703E1;
00935
00936 static const double THPIO4 = 3.*PI/4.;
00937
00938 double bessel_j1(double x)
00939 {
00940 double w, z, p, q, xn;
00941
00942 DEBUG_ENTRY( "bessel_j1()" );
00943
00944 w = x;
00945 if( x < 0 )
00946 w = -x;
00947
00948 if( w <= 5.0 )
00949 {
00950 z = x * x;
00951 w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
00952 w = w * x * (z - Z1) * (z - Z2);
00953 return w;
00954 }
00955
00956 w = 5.0/x;
00957 z = w * w;
00958 p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
00959 q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
00960 xn = x - THPIO4;
00961 p = p * cos(xn) - w * q * sin(xn);
00962 return p * SQ2OPI / sqrt(x);
00963 }
00964
00965 double bessel_y1(double x)
00966 {
00967 double w, z, p, q, xn;
00968
00969 DEBUG_ENTRY( "bessel_y1()" );
00970
00971 if( x <= 5.0 )
00972 {
00973 if( x <= 0.0 )
00974 {
00975 fprintf( ioQQQ, "bessel_y1: domain error\n" );
00976 cdEXIT(EXIT_FAILURE);
00977 }
00978 z = x * x;
00979 w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
00980 w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
00981 return w;
00982 }
00983
00984 w = 5.0/x;
00985 z = w * w;
00986 p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
00987 q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
00988 xn = x - THPIO4;
00989 p = p * sin(xn) + w * q * cos(xn);
00990 return p * SQ2OPI / sqrt(x);
00991 }
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 double bessel_jn(int n, double x)
01042 {
01043 double pkm2, pkm1, pk, xk, r, ans;
01044 int k, sign;
01045
01046 DEBUG_ENTRY( "bessel_jn()" );
01047
01048 if( n < 0 )
01049 {
01050 n = -n;
01051 if( (n & 1) == 0 )
01052 sign = 1;
01053 else
01054 sign = -1;
01055 }
01056 else
01057 sign = 1;
01058
01059 if( x < 0.0 )
01060 {
01061 if( n & 1 )
01062 sign = -sign;
01063 x = -x;
01064 }
01065
01066 if( x < DBL_EPSILON )
01067 {
01068 return sign * powi(x/2.,n)/factorial(n);
01069 }
01070
01071 if( n == 0 )
01072 {
01073 return sign * bessel_j0(x);
01074 }
01075 if( n == 1 )
01076 {
01077 return sign * bessel_j1(x);
01078 }
01079
01080 if( n == 2 && x > 0.1 )
01081 {
01082 return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
01083 }
01084
01085
01086 k = 52;
01087
01088 pk = 2 * (n + k);
01089 ans = pk;
01090 xk = x * x;
01091
01092 do
01093 {
01094 pk -= 2.0;
01095 ans = pk - (xk/ans);
01096 }
01097 while( --k > 0 );
01098 ans = x/ans;
01099
01100
01101 pk = 1.0;
01102 pkm1 = 1.0/ans;
01103 k = n-1;
01104 r = 2 * k;
01105
01106 do
01107 {
01108 pkm2 = (pkm1 * r - pk * x) / x;
01109 pk = pkm1;
01110 pkm1 = pkm2;
01111 r -= 2.0;
01112 }
01113 while( --k > 0 );
01114
01115 if( fabs(pk) > fabs(pkm1) )
01116 ans = bessel_j1(x)/pk;
01117 else
01118 ans = bessel_j0(x)/pkm1;
01119 return sign * ans;
01120 }
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 double bessel_yn(int n, double x)
01177 {
01178 double an, anm1, anm2, r;
01179 int k, sign;
01180
01181 DEBUG_ENTRY( "bessel_yn()" );
01182
01183 if( n < 0 )
01184 {
01185 n = -n;
01186 if( (n & 1) == 0 )
01187 sign = 1;
01188 else
01189 sign = -1;
01190 }
01191 else
01192 sign = 1;
01193
01194 if( n == 0 )
01195 {
01196 return sign * bessel_y0(x);
01197 }
01198 if( n == 1 )
01199 {
01200 return sign * bessel_y1(x);
01201 }
01202
01203
01204 if( x <= 0.0 )
01205 {
01206 fprintf( ioQQQ, "bessel_yn: domain error\n" );
01207 cdEXIT(EXIT_FAILURE);
01208 }
01209
01210
01211 anm2 = bessel_y0(x);
01212 anm1 = bessel_y1(x);
01213 k = 1;
01214 r = 2 * k;
01215 do
01216 {
01217 an = r * anm1 / x - anm2;
01218 anm2 = anm1;
01219 anm1 = an;
01220 r += 2.0;
01221 ++k;
01222 }
01223 while( k < n );
01224 return sign * an;
01225 }
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 static const double k0_A[] =
01311 {
01312 1.37446543561352307156e-16,
01313 4.25981614279661018399e-14,
01314 1.03496952576338420167e-11,
01315 1.90451637722020886025e-9,
01316 2.53479107902614945675e-7,
01317 2.28621210311945178607e-5,
01318 1.26461541144692592338e-3,
01319 3.59799365153615016266e-2,
01320 3.44289899924628486886e-1,
01321 -5.35327393233902768720e-1
01322 };
01323
01324
01325
01326
01327
01328
01329
01330 static const double k0_B[] = {
01331 5.30043377268626276149e-18,
01332 -1.64758043015242134646e-17,
01333 5.21039150503902756861e-17,
01334 -1.67823109680541210385e-16,
01335 5.51205597852431940784e-16,
01336 -1.84859337734377901440e-15,
01337 6.34007647740507060557e-15,
01338 -2.22751332699166985548e-14,
01339 8.03289077536357521100e-14,
01340 -2.98009692317273043925e-13,
01341 1.14034058820847496303e-12,
01342 -4.51459788337394416547e-12,
01343 1.85594911495471785253e-11,
01344 -7.95748924447710747776e-11,
01345 3.57739728140030116597e-10,
01346 -1.69753450938905987466e-9,
01347 8.57403401741422608519e-9,
01348 -4.66048989768794782956e-8,
01349 2.76681363944501510342e-7,
01350 -1.83175552271911948767e-6,
01351 1.39498137188764993662e-5,
01352 -1.28495495816278026384e-4,
01353 1.56988388573005337491e-3,
01354 -3.14481013119645005427e-2,
01355 2.44030308206595545468e0
01356 };
01357
01358 double bessel_k0(double x)
01359 {
01360 double y, z;
01361
01362 DEBUG_ENTRY( "bessel_k0()" );
01363
01364 if( x <= 0.0 )
01365 {
01366 fprintf( ioQQQ, "bessel_k0: domain error\n" );
01367 cdEXIT(EXIT_FAILURE);
01368 }
01369
01370 if( x <= 2.0 )
01371 {
01372 y = x * x - 2.0;
01373 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01374 return y;
01375 }
01376 z = 8.0/x - 2.0;
01377 y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x);
01378 return y;
01379 }
01380
01381 double bessel_k0_scaled(double x)
01382 {
01383 double y;
01384
01385 DEBUG_ENTRY( "bessel_k0_scaled()" );
01386
01387 if( x <= 0.0 )
01388 {
01389 fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" );
01390 cdEXIT(EXIT_FAILURE);
01391 }
01392
01393 if( x <= 2.0 )
01394 {
01395 y = x * x - 2.0;
01396 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
01397 return y * exp(x);
01398 }
01399 return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x);
01400 }
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 static const double k1_A[] =
01485 {
01486 -7.02386347938628759343e-18,
01487 -2.42744985051936593393e-15,
01488 -6.66690169419932900609e-13,
01489 -1.41148839263352776110e-10,
01490 -2.21338763073472585583e-8,
01491 -2.43340614156596823496e-6,
01492 -1.73028895751305206302e-4,
01493 -6.97572385963986435018e-3,
01494 -1.22611180822657148235e-1,
01495 -3.53155960776544875667e-1,
01496 1.52530022733894777053e0
01497 };
01498
01499
01500
01501
01502
01503
01504
01505 static const double k1_B[] =
01506 {
01507 -5.75674448366501715755e-18,
01508 1.79405087314755922667e-17,
01509 -5.68946255844285935196e-17,
01510 1.83809354436663880070e-16,
01511 -6.05704724837331885336e-16,
01512 2.03870316562433424052e-15,
01513 -7.01983709041831346144e-15,
01514 2.47715442448130437068e-14,
01515 -8.97670518232499435011e-14,
01516 3.34841966607842919884e-13,
01517 -1.28917396095102890680e-12,
01518 5.13963967348173025100e-12,
01519 -2.12996783842756842877e-11,
01520 9.21831518760500529508e-11,
01521 -4.19035475934189648750e-10,
01522 2.01504975519703286596e-9,
01523 -1.03457624656780970260e-8,
01524 5.74108412545004946722e-8,
01525 -3.50196060308781257119e-7,
01526 2.40648494783721712015e-6,
01527 -1.93619797416608296024e-5,
01528 1.95215518471351631108e-4,
01529 -2.85781685962277938680e-3,
01530 1.03923736576817238437e-1,
01531 2.72062619048444266945e0
01532 };
01533
01534 double bessel_k1(double x)
01535 {
01536 double y, z;
01537
01538 DEBUG_ENTRY( "bessel_k1()" );
01539
01540 z = 0.5 * x;
01541 if( z <= 0.0 )
01542 {
01543 fprintf( ioQQQ, "bessel_k1: domain error\n" );
01544 cdEXIT(EXIT_FAILURE);
01545 }
01546
01547 if( x <= 2.0 )
01548 {
01549 y = x * x - 2.0;
01550 y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01551 return y;
01552 }
01553 return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01554 }
01555
01556 double bessel_k1_scaled(double x)
01557 {
01558 double y;
01559
01560 DEBUG_ENTRY( "bessel_k1_scaled()" );
01561
01562 if( x <= 0.0 )
01563 {
01564 fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" );
01565 cdEXIT(EXIT_FAILURE);
01566 }
01567
01568 if( x <= 2.0 )
01569 {
01570 y = x * x - 2.0;
01571 y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
01572 return y * exp(x);
01573 }
01574 return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
01575 }
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 static const double i0_A[] =
01657 {
01658 -4.41534164647933937950e-18,
01659 3.33079451882223809783e-17,
01660 -2.43127984654795469359e-16,
01661 1.71539128555513303061e-15,
01662 -1.16853328779934516808e-14,
01663 7.67618549860493561688e-14,
01664 -4.85644678311192946090e-13,
01665 2.95505266312963983461e-12,
01666 -1.72682629144155570723e-11,
01667 9.67580903537323691224e-11,
01668 -5.18979560163526290666e-10,
01669 2.65982372468238665035e-9,
01670 -1.30002500998624804212e-8,
01671 6.04699502254191894932e-8,
01672 -2.67079385394061173391e-7,
01673 1.11738753912010371815e-6,
01674 -4.41673835845875056359e-6,
01675 1.64484480707288970893e-5,
01676 -5.75419501008210370398e-5,
01677 1.88502885095841655729e-4,
01678 -5.76375574538582365885e-4,
01679 1.63947561694133579842e-3,
01680 -4.32430999505057594430e-3,
01681 1.05464603945949983183e-2,
01682 -2.37374148058994688156e-2,
01683 4.93052842396707084878e-2,
01684 -9.49010970480476444210e-2,
01685 1.71620901522208775349e-1,
01686 -3.04682672343198398683e-1,
01687 6.76795274409476084995e-1
01688 };
01689
01690
01691
01692
01693
01694
01695
01696 static const double i0_B[] =
01697 {
01698 -7.23318048787475395456e-18,
01699 -4.83050448594418207126e-18,
01700 4.46562142029675999901e-17,
01701 3.46122286769746109310e-17,
01702 -2.82762398051658348494e-16,
01703 -3.42548561967721913462e-16,
01704 1.77256013305652638360e-15,
01705 3.81168066935262242075e-15,
01706 -9.55484669882830764870e-15,
01707 -4.15056934728722208663e-14,
01708 1.54008621752140982691e-14,
01709 3.85277838274214270114e-13,
01710 7.18012445138366623367e-13,
01711 -1.79417853150680611778e-12,
01712 -1.32158118404477131188e-11,
01713 -3.14991652796324136454e-11,
01714 1.18891471078464383424e-11,
01715 4.94060238822496958910e-10,
01716 3.39623202570838634515e-9,
01717 2.26666899049817806459e-8,
01718 2.04891858946906374183e-7,
01719 2.89137052083475648297e-6,
01720 6.88975834691682398426e-5,
01721 3.36911647825569408990e-3,
01722 8.04490411014108831608e-1
01723 };
01724
01725 double bessel_i0(double x)
01726 {
01727 double y;
01728
01729 DEBUG_ENTRY( "bessel_i0()" );
01730
01731 if( x < 0 )
01732 x = -x;
01733
01734 if( x <= 8.0 )
01735 {
01736 y = (x/2.0) - 2.0;
01737 return exp(x) * chbevl( y, i0_A, 30 );
01738 }
01739 return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01740 }
01741
01742 double bessel_i0_scaled(double x)
01743 {
01744 double y;
01745
01746 DEBUG_ENTRY( "bessel_i0_scaled()" );
01747
01748 if( x < 0 )
01749 x = -x;
01750
01751 if( x <= 8.0 )
01752 {
01753 y = (x/2.0) - 2.0;
01754 return chbevl( y, i0_A, 30 );
01755 }
01756 return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
01757 }
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 static double i1_A[] =
01840 {
01841 2.77791411276104639959e-18,
01842 -2.11142121435816608115e-17,
01843 1.55363195773620046921e-16,
01844 -1.10559694773538630805e-15,
01845 7.60068429473540693410e-15,
01846 -5.04218550472791168711e-14,
01847 3.22379336594557470981e-13,
01848 -1.98397439776494371520e-12,
01849 1.17361862988909016308e-11,
01850 -6.66348972350202774223e-11,
01851 3.62559028155211703701e-10,
01852 -1.88724975172282928790e-9,
01853 9.38153738649577178388e-9,
01854 -4.44505912879632808065e-8,
01855 2.00329475355213526229e-7,
01856 -8.56872026469545474066e-7,
01857 3.47025130813767847674e-6,
01858 -1.32731636560394358279e-5,
01859 4.78156510755005422638e-5,
01860 -1.61760815825896745588e-4,
01861 5.12285956168575772895e-4,
01862 -1.51357245063125314899e-3,
01863 4.15642294431288815669e-3,
01864 -1.05640848946261981558e-2,
01865 2.47264490306265168283e-2,
01866 -5.29459812080949914269e-2,
01867 1.02643658689847095384e-1,
01868 -1.76416518357834055153e-1,
01869 2.52587186443633654823e-1
01870 };
01871
01872
01873
01874
01875
01876
01877
01878 static double i1_B[] =
01879 {
01880 7.51729631084210481353e-18,
01881 4.41434832307170791151e-18,
01882 -4.65030536848935832153e-17,
01883 -3.20952592199342395980e-17,
01884 2.96262899764595013876e-16,
01885 3.30820231092092828324e-16,
01886 -1.88035477551078244854e-15,
01887 -3.81440307243700780478e-15,
01888 1.04202769841288027642e-14,
01889 4.27244001671195135429e-14,
01890 -2.10154184277266431302e-14,
01891 -4.08355111109219731823e-13,
01892 -7.19855177624590851209e-13,
01893 2.03562854414708950722e-12,
01894 1.41258074366137813316e-11,
01895 3.25260358301548823856e-11,
01896 -1.89749581235054123450e-11,
01897 -5.58974346219658380687e-10,
01898 -3.83538038596423702205e-9,
01899 -2.63146884688951950684e-8,
01900 -2.51223623787020892529e-7,
01901 -3.88256480887769039346e-6,
01902 -1.10588938762623716291e-4,
01903 -9.76109749136146840777e-3,
01904 7.78576235018280120474e-1
01905 };
01906
01907 double bessel_i1(double x)
01908 {
01909 double y, z;
01910
01911 DEBUG_ENTRY( "bessel_i1()" );
01912
01913 z = fabs(x);
01914 if( z <= 8.0 )
01915 {
01916 y = (z/2.0) - 2.0;
01917 z = chbevl( y, i1_A, 29 ) * z * exp(z);
01918 }
01919 else
01920 {
01921 z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01922 }
01923 if( x < 0.0 )
01924 z = -z;
01925 return z;
01926 }
01927
01928 double bessel_i1_scaled(double x)
01929 {
01930 double y, z;
01931
01932 DEBUG_ENTRY( "bessel_i1_scaled()" );
01933
01934 z = fabs(x);
01935 if( z <= 8.0 )
01936 {
01937 y = (z/2.0) - 2.0;
01938 z = chbevl( y, i1_A, 29 ) * z;
01939 }
01940 else
01941 {
01942 z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
01943 }
01944 if( x < 0.0 )
01945 z = -z;
01946 return z;
01947 }
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 static const double elk_P[] =
02009 {
02010 1.37982864606273237150e-4,
02011 2.28025724005875567385e-3,
02012 7.97404013220415179367e-3,
02013 9.85821379021226008714e-3,
02014 6.87489687449949877925e-3,
02015 6.18901033637687613229e-3,
02016 8.79078273952743772254e-3,
02017 1.49380448916805252718e-2,
02018 3.08851465246711995998e-2,
02019 9.65735902811690126535e-2,
02020 1.38629436111989062502e0
02021 };
02022
02023 static const double elk_Q[] =
02024 {
02025 2.94078955048598507511e-5,
02026 9.14184723865917226571e-4,
02027 5.94058303753167793257e-3,
02028 1.54850516649762399335e-2,
02029 2.39089602715924892727e-2,
02030 3.01204715227604046988e-2,
02031 3.73774314173823228969e-2,
02032 4.88280347570998239232e-2,
02033 7.03124996963957469739e-2,
02034 1.24999999999870820058e-1,
02035 4.99999999999999999821e-1
02036 };
02037
02038 static const double C1 = 1.3862943611198906188e0;
02039
02040 double ellpk(double x)
02041 {
02042 DEBUG_ENTRY( "ellpk()" );
02043
02044 if( x < 0.0 || x > 1.0 )
02045 {
02046 fprintf( ioQQQ, "ellpk: domain error\n" );
02047 cdEXIT(EXIT_FAILURE);
02048 }
02049
02050 if( x > DBL_EPSILON )
02051 {
02052 return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
02053 }
02054 else
02055 {
02056 if( x == 0.0 )
02057 {
02058 fprintf( ioQQQ, "ellpk: domain error\n" );
02059 cdEXIT(EXIT_FAILURE);
02060 }
02061 else
02062 {
02063 return C1 - 0.5 * log(x);
02064 }
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 static const double MAXLOG = log(DBL_MAX);
02117 static const double BIG = 1.44115188075855872E+17;
02118
02119
02120 double expn(int n, double x)
02121 {
02122 double ans, r, t, yk, xk;
02123 double pk, pkm1, pkm2, qk, qkm1, qkm2;
02124 double psi, z;
02125 int i, k;
02126
02127 DEBUG_ENTRY( "expn()" );
02128
02129 if( n < 0 || x < 0. )
02130 {
02131 fprintf( ioQQQ, "expn: domain error\n" );
02132 cdEXIT(EXIT_FAILURE);
02133 }
02134
02135 if( x > MAXLOG )
02136 {
02137 return 0.0;
02138 }
02139
02140 if( x == 0.0 )
02141 {
02142 if( n < 2 )
02143 {
02144 fprintf( ioQQQ, "expn: domain error\n" );
02145 cdEXIT(EXIT_FAILURE);
02146 }
02147 else
02148 {
02149 return 1.0/((double)n-1.0);
02150 }
02151 }
02152
02153 if( n == 0 )
02154 {
02155 return exp(-x)/x;
02156 }
02157
02158
02159 if( n > 5000 )
02160 {
02161 xk = x + n;
02162 yk = 1.0 / (xk * xk);
02163 t = n;
02164 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
02165 ans = yk * (ans + t * (t - 2.0 * x));
02166 ans = yk * (ans + t);
02167 ans = (ans + 1.0) * exp( -x ) / xk;
02168 return ans;
02169 }
02170
02171 if( x <= 1.0 )
02172 {
02173
02174 psi = -EULER - log(x);
02175 for( i=1; i < n; i++ )
02176 psi = psi + 1.0/i;
02177
02178 z = -x;
02179 xk = 0.0;
02180 yk = 1.0;
02181 pk = 1.0 - n;
02182 if( n == 1 )
02183 ans = 0.0;
02184 else
02185 ans = 1.0/pk;
02186 do
02187 {
02188 xk += 1.0;
02189 yk *= z/xk;
02190 pk += 1.0;
02191 if( pk != 0.0 )
02192 {
02193 ans += yk/pk;
02194 }
02195 if( ans != 0.0 )
02196 t = fabs(yk/ans);
02197 else
02198 t = 1.0;
02199 }
02200 while( t > DBL_EPSILON );
02201 ans = powi(z,n-1)*psi/factorial(n-1) - ans;
02202 return ans;
02203 }
02204 else
02205 {
02206
02207 k = 1;
02208 pkm2 = 1.0;
02209 qkm2 = x;
02210 pkm1 = 1.0;
02211 qkm1 = x + n;
02212 ans = pkm1/qkm1;
02213
02214 do
02215 {
02216 k += 1;
02217 if( is_odd(k) )
02218 {
02219 yk = 1.0;
02220 xk = static_cast<double>(n + (k-1)/2);
02221 }
02222 else
02223 {
02224 yk = x;
02225 xk = static_cast<double>(k/2);
02226 }
02227 pk = pkm1 * yk + pkm2 * xk;
02228 qk = qkm1 * yk + qkm2 * xk;
02229 if( qk != 0 )
02230 {
02231 r = pk/qk;
02232 t = fabs( (ans - r)/r );
02233 ans = r;
02234 }
02235 else
02236 t = 1.0;
02237 pkm2 = pkm1;
02238 pkm1 = pk;
02239 qkm2 = qkm1;
02240 qkm1 = qk;
02241 if( fabs(pk) > BIG )
02242 {
02243 pkm2 /= BIG;
02244 pkm1 /= BIG;
02245 qkm2 /= BIG;
02246 qkm1 /= BIG;
02247 }
02248 }
02249 while( t > DBL_EPSILON );
02250
02251 ans *= exp( -x );
02252 return ans;
02253 }
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 inline double polevl(double x, const double coef[], int N)
02307 {
02308 double ans;
02309 int i;
02310 const double *p = coef;
02311
02312 ans = *p++;
02313 i = N;
02314
02315 do
02316 ans = ans * x + *p++;
02317 while( --i );
02318
02319 return ans;
02320 }
02321
02322
02323
02324
02325
02326
02327
02328 inline double p1evl(double x, const double coef[], int N)
02329 {
02330 double ans;
02331 const double *p = coef;
02332 int i;
02333
02334 ans = x + *p++;
02335 i = N-1;
02336
02337 do
02338 ans = ans * x + *p++;
02339 while( --i );
02340
02341 return ans;
02342 }
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402 inline double chbevl(double x, const double array[], int n)
02403 {
02404 double b0, b1, b2;
02405 const double *p = array;
02406 int i;
02407
02408 b0 = *p++;
02409 b1 = 0.0;
02410 i = n - 1;
02411
02412 do
02413 {
02414 b2 = b1;
02415 b1 = b0;
02416 b0 = x * b1 - b2 + *p++;
02417 }
02418 while( --i );
02419
02420 return 0.5*(b0-b2);
02421 }
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482 static const int N = 624;
02483 static const int M = 397;
02484 static const unsigned long MATRIX_A = 0x9908b0dfUL;
02485 static const unsigned long UMASK = 0x80000000UL;
02486 static const unsigned long LMASK = 0x7fffffffUL;
02487 inline unsigned long MIXBITS(unsigned long u, unsigned long v)
02488 {
02489 return (u & UMASK) | (v & LMASK);
02490 }
02491 inline unsigned long TWIST(unsigned long u, unsigned long v)
02492 {
02493 return (MIXBITS(u,v) >> 1) ^ (v&1UL ? MATRIX_A : 0UL);
02494 }
02495
02496 static unsigned long state[N];
02497 static int nleft = 1;
02498 static int initf = 0;
02499 static unsigned long *nexxt;
02500
02501
02502 void init_genrand(unsigned long s)
02503 {
02504 int j;
02505 state[0]= s & 0xffffffffUL;
02506 for (j=1; j<N; j++) {
02507 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
02508
02509
02510
02511
02512 state[j] &= 0xffffffffUL;
02513 }
02514 nleft = 1; initf = 1;
02515 }
02516
02517
02518
02519
02520
02521 void init_by_array(unsigned long init_key[], int key_length)
02522 {
02523 int i, j, k;
02524 init_genrand(19650218UL);
02525 i=1; j=0;
02526 k = (N>key_length ? N : key_length);
02527 for (; k; k--) {
02528 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
02529 + init_key[j] + j;
02530 state[i] &= 0xffffffffUL;
02531 i++; j++;
02532 if (i>=N) { state[0] = state[N-1]; i=1; }
02533 if (j>=key_length) j=0;
02534 }
02535 for (k=N-1; k; k--) {
02536 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
02537 - i;
02538 state[i] &= 0xffffffffUL;
02539 i++;
02540 if (i>=N) { state[0] = state[N-1]; i=1; }
02541 }
02542
02543 state[0] = 0x80000000UL;
02544 nleft = 1; initf = 1;
02545 }
02546
02547 static void next_state()
02548 {
02549 unsigned long *p=state;
02550 int j;
02551
02552
02553
02554 if (initf==0) init_genrand(5489UL);
02555
02556 nleft = N;
02557 nexxt = state;
02558
02559 for (j=N-M+1; --j; p++)
02560 *p = p[M] ^ TWIST(p[0], p[1]);
02561
02562 for (j=M; --j; p++)
02563 *p = p[M-N] ^ TWIST(p[0], p[1]);
02564
02565 *p = p[M-N] ^ TWIST(p[0], state[0]);
02566 }
02567
02568
02569 unsigned long genrand_int32()
02570 {
02571 unsigned long y;
02572
02573 if (--nleft == 0) next_state();
02574 y = *nexxt++;
02575
02576
02577 y ^= (y >> 11);
02578 y ^= (y << 7) & 0x9d2c5680UL;
02579 y ^= (y << 15) & 0xefc60000UL;
02580 y ^= (y >> 18);
02581
02582 return y;
02583 }
02584
02585
02586 long genrand_int31()
02587 {
02588 unsigned long y;
02589
02590 if (--nleft == 0) next_state();
02591 y = *nexxt++;
02592
02593
02594 y ^= (y >> 11);
02595 y ^= (y << 7) & 0x9d2c5680UL;
02596 y ^= (y << 15) & 0xefc60000UL;
02597 y ^= (y >> 18);
02598
02599 return (long)(y>>1);
02600 }
02601
02602
02603 double genrand_real1()
02604 {
02605 unsigned long y;
02606
02607 if (--nleft == 0) next_state();
02608 y = *nexxt++;
02609
02610
02611 y ^= (y >> 11);
02612 y ^= (y << 7) & 0x9d2c5680UL;
02613 y ^= (y << 15) & 0xefc60000UL;
02614 y ^= (y >> 18);
02615
02616 return (double)y * (1.0/4294967295.0);
02617
02618 }
02619
02620
02621 double genrand_real2()
02622 {
02623 unsigned long y;
02624
02625 if (--nleft == 0) next_state();
02626 y = *nexxt++;
02627
02628
02629 y ^= (y >> 11);
02630 y ^= (y << 7) & 0x9d2c5680UL;
02631 y ^= (y << 15) & 0xefc60000UL;
02632 y ^= (y >> 18);
02633
02634 return (double)y * (1.0/4294967296.0);
02635
02636 }
02637
02638
02639 double genrand_real3()
02640 {
02641 unsigned long y;
02642
02643 if (--nleft == 0) next_state();
02644 y = *nexxt++;
02645
02646
02647 y ^= (y >> 11);
02648 y ^= (y << 7) & 0x9d2c5680UL;
02649 y ^= (y << 15) & 0xefc60000UL;
02650 y ^= (y >> 18);
02651
02652 return ((double)y + 0.5) * (1.0/4294967296.0);
02653
02654 }
02655
02656
02657 double genrand_res53()
02658 {
02659 unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
02660 return (a*67108864.0+b)*(1.0/9007199254740992.0);
02661 }
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679 void humlik(int n, realnum x[], realnum y, realnum k[])
02680 {
02681
02682
02683
02684
02685
02686
02687 const realnum c[6] = { 1.0117281f,-.75197147f,.012557727f,.010022008f,
02688 -2.4206814e-4f,5.0084806e-7f };
02689 const realnum s[6] = { 1.393237,.23115241f,-.15535147f,.0062183662f,
02690 9.1908299e-5f,-6.2752596e-7f };
02691 const realnum t[6] = { .31424038f,.94778839f,1.5976826f,2.2795071f,
02692 3.020637f,3.8897249f };
02693
02694 const realnum R0 = 146.7f, R1 = 14.67f;
02695
02696
02697 realnum d;
02698 int i, j;
02699 realnum a0, d0, e0, d2, e2, h0, e4, h2, h4, h6, p0,
02700 p2, p4, p6, p8, z0, z2, z4, z6, z8, mf[6], pf[6],
02701 mq[6], yf, pq[6], xm[6], ym[6], xp[6], xq, yq, yp[6];
02702 bool rg1, rg2, rg3;
02703 realnum abx, ypy0, xlim0, xlim1, xlim2, xlim3, xlim4, ypy0q, yrrtpi;
02704
02705 rg1 = rg2 = rg3 = true;
02706
02707 z0 = z2 = z4 = z6 = z8 = 0.;
02708 p0 = p2 = p4 = p6 = p8 = 0.;
02709 h0 = h2 = h4 = h6 = 0.;
02710 a0 = d0 = d2 = e0 = e2 = e4 = 0.;
02711
02712 yq = y * y;
02713 yrrtpi = y * .56418958;
02714
02715 xlim0 = R0 - y;
02716 xlim1 = R1 - y;
02717 xlim3 = y * 3.097f - .45f;
02718 xlim2 = 6.8f - y;
02719 xlim4 = y * 18.1f + 1.65f;
02720 if (y <= 1e-6f)
02721 {
02722 xlim1 = xlim0;
02723 xlim2 = xlim0;
02724 }
02725
02726 for (i = 0; i < n; ++i)
02727 {
02728 abx = fabs(x[i]);
02729 xq = abx * abx;
02730 if (abx > xlim0)
02731 {
02732 k[i] = yrrtpi / (xq + yq);
02733 }
02734 else if (abx > xlim1)
02735 {
02736 if (rg1)
02737 {
02738
02739 rg1 = false;
02740 a0 = yq + .5f;
02741 d0 = a0 * a0;
02742 d2 = yq + yq - 1.f;
02743 }
02744 d = .56418958f / (d0 + xq * (d2 + xq));
02745 k[i] = d * y * (a0 + xq);
02746 }
02747 else if (abx > xlim2)
02748 {
02749 if (rg2)
02750 {
02751 rg2 = false;
02752 h0 = yq * (yq * (yq * (yq + 6.f) + 10.5f) + 4.5f) + .5625f;
02753 h2 = yq * (yq * (yq * 4.f + 6.f) + 9.f) - 4.5f;
02754 h4 = 10.5f - yq * (6.f - yq * 6.f);
02755 h6 = yq * 4. - 6.;
02756 e0 = yq * (yq * (yq + 5.5f) + 8.25f) + 1.875f;
02757 e2 = yq * (yq * 3.f + 1.f) + 5.25f;
02758 e4 = h6 * .75f;
02759 }
02760 d = .56418958f / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
02761 k[i] = d * y * (e0 + xq * (e2 + xq * (e4 + xq)));
02762 }
02763 else if (abx < xlim3)
02764 {
02765 if (rg3)
02766 {
02767
02768 rg3 = false;
02769 z0 = y * (y * (y * (y * (y * (y *
02770 (y * (y * (y * (y + 13.3988f)
02771 + 88.26741f)
02772 + 369.1989f) + 1074.409f) +
02773 2256.981f)
02774 + 3447.629f) + 3764.966f) + 2802.87f) +
02775 1280.829f) + 272.1014f;
02776 z2 = y * (y * (y * (y * (y *
02777 (y * (y * (y * 5.f + 53.59518f)
02778 + 266.2987f) + 793.4273f)
02779 + 1549.675f) +
02780 2037.31f) + 1758.336f) + 902.3066f) + 211.678f;
02781 z4 = y * (y * (y * (y * (y * (y * 10.f + 80.39278f) +
02782 269.2916f) + 479.2576f) + 497.3014f)
02783 + 308.1852f) + 78.86585f;
02784 z6 = y * (y * (y * (y * 10.f + 53.59518f) + 92.75679f) +
02785 55.02933f) + 22.03523f;
02786 z8 = y * (y * 5.f + 13.3988f) + 1.49646f;
02787 p0 = y * (y * (y * (y * (y * (y *
02788 (y * (y * (y *
02789 .3183291f + 4.264678f)
02790 + 27.93941f) + 115.3772f) +
02791 328.2151f) + 662.8097f)
02792 + 946.897f) + 919.4955f) +
02793 549.3954f) + 153.5168f;
02794 p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163f +
02795 12.79458f) + 56.81652f)
02796 + 139.4665f) + 189.773f) +
02797 124.5975f) - 1.322256f) - 34.16955f;
02798 p4 = y * (y * (y * (y * (y * 1.9099744f + 12.79568f) +
02799 29.81482f) + 24.01655f)
02800 + 10.46332f) + 2.584042f;
02801 p6 = y * (y * (y * 1.273316f + 4.266322f) + .9377051f) -
02802 .07272979f;
02803 p8 = y * .3183291f + 5.480304e-4f;
02804 }
02805 d = 1.7724538f / (z0 + xq * (z2 + xq * (z4 + xq *
02806 (z6 + xq * ( z8 + xq)))));
02807 k[i] = d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
02808 }
02809 else
02810 {
02811 ypy0 = y + 1.5f;
02812 ypy0q = ypy0 * ypy0;
02813 k[i] = 0.f;
02814 for (j = 0; j <= 5; ++j)
02815 {
02816 d = x[i] - t[j];
02817 mq[j] = d * d;
02818 mf[j] = 1.f / (mq[j] + ypy0q);
02819 xm[j] = mf[j] * d;
02820 ym[j] = mf[j] * ypy0;
02821 d = x[i] + t[j];
02822 pq[j] = d * d;
02823 pf[j] = 1. / (pq[j] + ypy0q);
02824 xp[j] = pf[j] * d;
02825 yp[j] = pf[j] * ypy0;
02826 }
02827 if (abx <= xlim4)
02828 {
02829 for (j = 0; j <= 5; ++j)
02830 {
02831 k[i] = k[i] + c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
02832 }
02833 }
02834 else
02835 {
02836 yf = y + 3.f;
02837 for (j = 0; j <= 5; ++j)
02838 {
02839 k[i] = k[i] + (c[j] * (mq[j] * mf[j] - ym[j] * 1.5f)
02840 + s[j] * yf * xm[j]) / (mq[j] + 2.25f) +
02841 (c[j] * (pq[j] * pf[j] - yp[j] * 1.5f) - s[j] * yf * xp[j])
02842 / (pq[j] + 2.25f);
02843 }
02844 k[i] = y * k[i] + exp(-xq);
02845 }
02846 }
02847 }
02848 }
02849
02850
02851
02852
02853
02854 STATIC uint32 MD5swap( uint32 word );
02855 STATIC void MD5_Transform (uint32 *digest, const uint32 *in);
02856
02857
02858
02859
02860
02861 string MD5file(const char* fnam, access_scheme scheme)
02862 {
02863 DEBUG_ENTRY( "MD5file()" );
02864
02865 fstream ioFile;
02866 open_data( ioFile, fnam, mode_r, scheme );
02867
02868 char c;
02869 string content;
02870 while( ioFile.get(c) )
02871 content += c;
02872
02873 return MD5string( content );
02874 }
02875
02876
02877
02878
02879
02880 string MD5datafile(const char* fnam, access_scheme scheme)
02881 {
02882 DEBUG_ENTRY( "MD5datafile()" );
02883
02884 fstream ioFile;
02885 open_data( ioFile, fnam, mode_r, scheme );
02886
02887 string line, content;
02888 while( getline( ioFile, line ) )
02889 if( line[0] != '#' )
02890 content += line;
02891
02892 return MD5string( content );
02893 }
02894
02895
02896
02897 string MD5string(const string& str)
02898 {
02899 DEBUG_ENTRY( "MD5string()" );
02900
02901 uint32 state[4];
02902
02903 state[0] = 0x67452301L;
02904 state[1] = 0xefcdab89L;
02905 state[2] = 0x98badcfeL;
02906 state[3] = 0x10325476L;
02907
02908 string lstr = str;
02909
02910
02911 uint32 bytes = str.length()%64;
02912 uint32 padlen = ( bytes >= 56 ) ? 64 + 56 - bytes : 56 - bytes;
02913 lstr += '\x80';
02914 for( uint32 i=1; i < padlen; ++i )
02915 lstr += '\0';
02916
02917 ASSERT( lstr.length()%64 == 56 );
02918
02919 uint32 i;
02920 union {
02921 uint32 in[16];
02922 unsigned char chr[64];
02923 } u;
02924
02925 for( i=0; i < lstr.length()/64; ++i )
02926 {
02927 for( uint32 j=0; j < 64; ++j )
02928 {
02929 if( cpu.little_endian() )
02930 u.chr[j] = lstr[i*64+j];
02931 else
02932 {
02933 uint32 jr = j%4;
02934 uint32 j4 = j-jr;
02935 u.chr[j4+3-jr] = lstr[i*64+j];
02936 }
02937 }
02938 MD5_Transform( state, u.in );
02939 }
02940 for( uint32 j=0; j < 56; ++j )
02941 {
02942 if( cpu.little_endian() )
02943 u.chr[j] = lstr[i*64+j];
02944 else
02945 {
02946 uint32 jr = j%4;
02947 uint32 j4 = j-jr;
02948 u.chr[j4+3-jr] = lstr[i*64+j];
02949 }
02950 }
02951
02952 u.in[14] = (str.length()<<3)&0xffffffff;
02953 u.in[15] = (str.length()>>29)&0xffffffff;
02954
02955 MD5_Transform( state, u.in );
02956
02957 ostringstream hash;
02958 for( uint32 i=0; i < 4; ++i )
02959 hash << hex << setfill('0') << setw(8) << MD5swap(state[i]);
02960
02961 return hash.str();
02962 }
02963
02964 STATIC uint32 MD5swap( uint32 word )
02965 {
02966 DEBUG_ENTRY( "MD5swap()" );
02967
02968 union {
02969 uint32 word;
02970 unsigned char byte[4];
02971 } ui, uo;
02972
02973 ui.word = word;
02974 uo.byte[0] = ui.byte[3];
02975 uo.byte[1] = ui.byte[2];
02976 uo.byte[2] = ui.byte[1];
02977 uo.byte[3] = ui.byte[0];
02978
02979 return uo.word;
02980 }
02981
02982
02983
02984
02985
02986
02987
02988
02989
02990
02991
02992
02993
02994
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013
03014
03015
03016
03017
03018
03019
03020
03021
03022
03023
03024
03025
03026
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044
03045
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055 inline uint32 rotlFixed(uint32 x, unsigned int y)
03056 {
03057 return uint32((x<<y) | (x>>(32-y)));
03058 }
03059
03060 STATIC void MD5_Transform (uint32 *digest, const uint32 *in)
03061 {
03062 DEBUG_ENTRY( "MD5_Transform()" );
03063
03064 #define F1(x, y, z) (z ^ (x & (y ^ z)))
03065 #define F2(x, y, z) F1(z, x, y)
03066 #define F3(x, y, z) (x ^ y ^ z)
03067 #define F4(x, y, z) (y ^ (x | ~z))
03068
03069 #define MD5STEP(f, w, x, y, z, data, s) \
03070 w = rotlFixed(w + f(x, y, z) + data, s) + x
03071
03072 uint32 a, b, c, d;
03073
03074 a=digest[0];
03075 b=digest[1];
03076 c=digest[2];
03077 d=digest[3];
03078
03079 MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
03080 MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
03081 MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
03082 MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
03083 MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
03084 MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
03085 MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
03086 MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
03087 MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
03088 MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
03089 MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
03090 MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
03091 MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
03092 MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
03093 MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
03094 MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
03095
03096 MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
03097 MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
03098 MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
03099 MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
03100 MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
03101 MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
03102 MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
03103 MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
03104 MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
03105 MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
03106 MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
03107 MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
03108 MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
03109 MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
03110 MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
03111 MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
03112
03113 MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
03114 MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
03115 MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
03116 MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
03117 MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
03118 MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
03119 MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
03120 MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
03121 MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
03122 MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
03123 MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
03124 MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
03125 MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
03126 MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
03127 MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
03128 MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
03129
03130 MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
03131 MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
03132 MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
03133 MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
03134 MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
03135 MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
03136 MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
03137 MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
03138 MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
03139 MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
03140 MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
03141 MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
03142 MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
03143 MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
03144 MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
03145 MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
03146
03147 digest[0] += a;
03148 digest[1] += b;
03149 digest[2] += c;
03150 digest[3] += d;
03151 }
03152
03153
03154
03155