7 inline double polevl(
double x,
const double coef[],
int N);
8 inline double p1evl(
double x,
const double coef[],
int N);
9 inline double chbevl(
double,
const double[],
int);
103 for(
long i=0; i < n; i++ )
108 double rn = (double)n;
113 for(
long i=0; i < n; i++ )
134 for(
long i=0; i < n; i++ )
135 sum1 +=
pow2(x[i]*(y[i] - b*x[i]));
138 sigb = sum1/
pow2(sxx);
141 for(
long i=0; i < n; i++ )
142 siga +=
pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
146 siga = sqrt(siga)/rn;
149 for(
long i=0; i < n; i++ )
167 1.00000000000000000000e+00,
168 1.00000000000000000000e+00,
169 2.00000000000000000000e+00,
170 6.00000000000000000000e+00,
171 2.40000000000000000000e+01,
172 1.20000000000000000000e+02,
173 7.20000000000000000000e+02,
174 5.04000000000000000000e+03,
175 4.03200000000000000000e+04,
176 3.62880000000000000000e+05,
177 3.62880000000000000000e+06,
178 3.99168000000000000000e+07,
179 4.79001600000000000000e+08,
180 6.22702080000000000000e+09,
181 8.71782912000000000000e+10,
182 1.30767436800000000000e+12,
183 2.09227898880000000000e+13,
184 3.55687428096000000000e+14,
185 6.40237370572800000000e+15,
186 1.21645100408832000000e+17,
187 2.43290200817664000000e+18,
188 5.10909421717094400000e+19,
189 1.12400072777760768000e+21,
190 2.58520167388849766400e+22,
191 6.20448401733239439360e+23,
192 1.55112100433309859840e+25,
193 4.03291461126605635592e+26,
194 1.08888694504183521614e+28,
195 3.04888344611713860511e+29,
196 8.84176199373970195470e+30,
197 2.65252859812191058647e+32,
198 8.22283865417792281807e+33,
199 2.63130836933693530178e+35,
200 8.68331761881188649615e+36,
201 2.95232799039604140861e+38,
202 1.03331479663861449300e+40,
203 3.71993326789901217463e+41,
204 1.37637530912263450457e+43,
205 5.23022617466601111726e+44,
206 2.03978820811974433568e+46,
207 8.15915283247897734264e+47,
208 3.34525266131638071044e+49,
209 1.40500611775287989839e+51,
210 6.04152630633738356321e+52,
211 2.65827157478844876773e+54,
212 1.19622220865480194551e+56,
213 5.50262215981208894940e+57,
214 2.58623241511168180614e+59,
215 1.24139155925360726691e+61,
216 6.08281864034267560801e+62,
217 3.04140932017133780398e+64,
218 1.55111875328738227999e+66,
219 8.06581751709438785585e+67,
220 4.27488328406002556374e+69,
221 2.30843697339241380441e+71,
222 1.26964033536582759243e+73,
223 7.10998587804863451749e+74,
224 4.05269195048772167487e+76,
225 2.35056133128287857145e+78,
226 1.38683118545689835713e+80,
227 8.32098711274139014271e+81,
228 5.07580213877224798711e+83,
229 3.14699732603879375200e+85,
230 1.98260831540444006372e+87,
231 1.26886932185884164078e+89,
232 8.24765059208247066512e+90,
233 5.44344939077443063905e+92,
234 3.64711109181886852801e+94,
235 2.48003554243683059915e+96,
236 1.71122452428141311337e+98,
237 1.19785716699698917933e+100,
238 8.50478588567862317347e+101,
239 6.12344583768860868500e+103,
240 4.47011546151268434004e+105,
241 3.30788544151938641157e+107,
242 2.48091408113953980872e+109,
243 1.88549470166605025458e+111,
244 1.45183092028285869606e+113,
245 1.13242811782062978295e+115,
246 8.94618213078297528506e+116,
247 7.15694570462638022794e+118,
248 5.79712602074736798470e+120,
249 4.75364333701284174746e+122,
250 3.94552396972065865030e+124,
251 3.31424013456535326627e+126,
252 2.81710411438055027626e+128,
253 2.42270953836727323750e+130,
254 2.10775729837952771662e+132,
255 1.85482642257398439069e+134,
256 1.65079551609084610774e+136,
257 1.48571596448176149700e+138,
258 1.35200152767840296226e+140,
259 1.24384140546413072522e+142,
260 1.15677250708164157442e+144,
261 1.08736615665674307994e+146,
262 1.03299784882390592592e+148,
263 9.91677934870949688836e+149,
264 9.61927596824821198159e+151,
265 9.42689044888324774164e+153,
266 9.33262154439441526381e+155,
267 9.33262154439441526388e+157,
268 9.42594775983835941673e+159,
269 9.61446671503512660515e+161,
270 9.90290071648618040340e+163,
271 1.02990167451456276198e+166,
272 1.08139675824029090008e+168,
273 1.14628056373470835406e+170,
274 1.22652020319613793888e+172,
275 1.32464181945182897396e+174,
276 1.44385958320249358163e+176,
277 1.58824554152274293982e+178,
278 1.76295255109024466316e+180,
279 1.97450685722107402277e+182,
280 2.23119274865981364576e+184,
281 2.54355973347218755612e+186,
282 2.92509369349301568964e+188,
283 3.39310868445189820004e+190,
284 3.96993716080872089396e+192,
285 4.68452584975429065488e+194,
286 5.57458576120760587943e+196,
287 6.68950291344912705515e+198,
288 8.09429852527344373681e+200,
289 9.87504420083360135884e+202,
290 1.21463043670253296712e+205,
291 1.50614174151114087918e+207,
292 1.88267717688892609901e+209,
293 2.37217324288004688470e+211,
294 3.01266001845765954361e+213,
295 3.85620482362580421582e+215,
296 4.97450422247728743840e+217,
297 6.46685548922047366972e+219,
298 8.47158069087882050755e+221,
299 1.11824865119600430699e+224,
300 1.48727070609068572828e+226,
301 1.99294274616151887582e+228,
302 2.69047270731805048244e+230,
303 3.65904288195254865604e+232,
304 5.01288874827499165889e+234,
305 6.91778647261948848943e+236,
306 9.61572319694108900019e+238,
307 1.34620124757175246000e+241,
308 1.89814375907617096864e+243,
309 2.69536413788816277557e+245,
310 3.85437071718007276916e+247,
311 5.55029383273930478744e+249,
312 8.04792605747199194159e+251,
313 1.17499720439091082343e+254,
314 1.72724589045463891049e+256,
315 2.55632391787286558753e+258,
316 3.80892263763056972532e+260,
317 5.71338395644585458806e+262,
318 8.62720977423324042775e+264,
319 1.31133588568345254503e+267,
320 2.00634390509568239384e+269,
321 3.08976961384735088657e+271,
322 4.78914290146339387432e+273,
323 7.47106292628289444390e+275,
324 1.17295687942641442768e+278,
325 1.85327186949373479574e+280,
326 2.94670227249503832518e+282,
327 4.71472363599206132029e+284,
328 7.59070505394721872577e+286,
329 1.22969421873944943358e+289,
330 2.00440157654530257674e+291,
331 3.28721858553429622598e+293,
332 5.42391066613158877297e+295,
333 9.00369170577843736335e+297,
334 1.50361651486499903974e+300,
335 2.52607574497319838672e+302,
336 4.26906800900470527345e+304,
337 7.25741561530799896496e+306
346 fprintf(
ioQQQ,
"factorial: domain error\n" );
362 p_lf.push_back( 0. );
363 p_lf.push_back( 0. );
370 if( n <
p_lf.size() )
376 for(
unsigned long i=(
unsigned long)
p_lf.size(); i <= n; i++ )
377 p_lf.push_back(
p_lf[i-1] + log10(static_cast<double>(i)) );
397 fprintf(
ioQQQ,
"lfactorial: domain error\n" );
413 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
429 ur = wr + 6.00009857740312429;
430 vr = ur * (wr + 4.99999857982434025) - wi * wi;
431 vi = wi * (wr + 4.99999857982434025) + ur * wi;
432 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
433 0.293729529320536228;
434 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
435 ur = vr * (wr + 4.00000003016801681) - vi * wi;
436 ui = vi * (wr + 4.00000003016801681) + vr * wi;
437 vr = ur * (wr + 2.99999999944915534) - ui * wi;
438 vi = ui * (wr + 2.99999999944915534) + ur * wi;
439 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
440 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
441 ur = vr * (wr + 2.00000000000603851) - vi * wi;
442 ui = vi * (wr + 2.00000000000603851) + vr * wi;
443 vr = ur * (wr + 0.999999999999975753) - ui * wi;
444 vi = ui * (wr + 0.999999999999975753) + ur * wi;
445 yr += ur * 10.5400280458730808 + vr;
446 yi += ui * 10.5400280458730808 + vi;
447 ur = vr * wr - vi * wi;
448 ui = vi * wr + vr * wi;
449 t = ur * ur + ui * ui;
450 vr = yr * ur + yi * ui + t * 0.0327673720261526849;
451 vi = yi * ur - yr * ui;
452 yr = wr + 7.31790632447016203;
453 ur = log(yr * yr + wi * wi) * 0.5 - 1;
455 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
456 yi = ui * (wr - 0.5) + ur * wi;
459 yr = ur * vr - ui * vi;
460 yi = ui * vr + ur * vi;
463 wr = xr * 3.14159265358979324;
464 wi = exp(xi * 3.14159265358979324);
466 ur = (vi + wi) * sin(wr);
467 ui = (vi - wi) * cos(wr);
468 vr = ur * yr + ui * yi;
469 vi = ui * yr - ur * yi;
470 ur = 6.2831853071795862 / (vr * vr + vi * vi);
474 return complex<double>( yr, yi );
592 7.96936729297347051624e-4,
593 8.28352392107440799803e-2,
594 1.23953371646414299388e0,
595 5.44725003058768775090e0,
596 8.74716500199817011941e0,
597 5.30324038235394892183e0,
598 9.99999999999999997821e-1,
602 9.24408810558863637013e-4,
603 8.56288474354474431428e-2,
604 1.25352743901058953537e0,
605 5.47097740330417105182e0,
606 8.76190883237069594232e0,
607 5.30605288235394617618e0,
608 1.00000000000000000218e0,
612 -1.13663838898469149931e-2,
613 -1.28252718670509318512e0,
614 -1.95539544257735972385e1,
615 -9.32060152123768231369e1,
616 -1.77681167980488050595e2,
617 -1.47077505154951170175e2,
618 -5.14105326766599330220e1,
619 -6.05014350600728481186e0,
624 6.43178256118178023184e1,
625 8.56430025976980587198e2,
626 3.88240183605401609683e3,
627 7.24046774195652478189e3,
628 5.93072701187316984827e3,
629 2.06209331660327847417e3,
630 2.42005740240291393179e2,
634 1.55924367855235737965e4,
635 -1.46639295903971606143e7,
636 5.43526477051876500413e9,
637 -9.82136065717911466409e11,
638 8.75906394395366999549e13,
639 -3.46628303384729719441e15,
640 4.42733268572569800351e16,
641 -1.84950800436986690637e16,
646 1.04128353664259848412e3,
647 6.26107330137134956842e5,
648 2.68919633393814121987e8,
649 8.64002487103935000337e10,
650 2.02979612750105546709e13,
651 3.17157752842975028269e15,
652 2.50596256172653059228e17,
656 static const double DR1 = 5.78318596294678452118e0;
658 static const double DR2 = 3.04712623436620863991e1;
661 -4.79443220978201773821e9,
662 1.95617491946556577543e12,
663 -2.49248344360967716204e14,
664 9.70862251047306323952e15,
669 4.99563147152651017219e2,
670 1.73785401676374683123e5,
671 4.84409658339962045305e7,
672 1.11855537045356834862e10,
673 2.11277520115489217587e12,
674 3.10518229857422583814e14,
675 3.18121955943204943306e16,
676 1.71086294081043136091e18,
685 double w, z, p, q, xn;
698 p = (z -
DR1) * (z -
DR2);
708 p = p * cos(xn) - w * q * sin(xn);
709 return p *
SQ2OPI / sqrt(x);
723 double w, z, p, q, xn;
731 fprintf(
ioQQQ,
"bessel_y0: domain error\n" );
745 p = p * sin(xn) + w * q * cos(xn);
746 return p *
SQ2OPI / sqrt(x);
828 -8.99971225705559398224e8,
829 4.52228297998194034323e11,
830 -7.27494245221818276015e13,
831 3.68295732863852883286e15,
836 6.20836478118054335476e2,
837 2.56987256757748830383e5,
838 8.35146791431949253037e7,
839 2.21511595479792499675e10,
840 4.74914122079991414898e12,
841 7.84369607876235854894e14,
842 8.95222336184627338078e16,
843 5.32278620332680085395e18,
847 7.62125616208173112003e-4,
848 7.31397056940917570436e-2,
849 1.12719608129684925192e0,
850 5.11207951146807644818e0,
851 8.42404590141772420927e0,
852 5.21451598682361504063e0,
853 1.00000000000000000254e0,
857 5.71323128072548699714e-4,
858 6.88455908754495404082e-2,
859 1.10514232634061696926e0,
860 5.07386386128601488557e0,
861 8.39985554327604159757e0,
862 5.20982848682361821619e0,
863 9.99999999999999997461e-1,
867 5.10862594750176621635e-2,
868 4.98213872951233449420e0,
869 7.58238284132545283818e1,
870 3.66779609360150777800e2,
871 7.10856304998926107277e2,
872 5.97489612400613639965e2,
873 2.11688757100572135698e2,
874 2.52070205858023719784e1,
879 7.42373277035675149943e1,
880 1.05644886038262816351e3,
881 4.98641058337653607651e3,
882 9.56231892404756170795e3,
883 7.99704160447350683650e3,
884 2.82619278517639096600e3,
885 3.36093607810698293419e2,
889 1.26320474790178026440e9,
890 -6.47355876379160291031e11,
891 1.14509511541823727583e14,
892 -8.12770255501325109621e15,
893 2.02439475713594898196e17,
894 -7.78877196265950026825e17,
899 5.94301592346128195359E2,
900 2.35564092943068577943E5,
901 7.34811944459721705660E7,
902 1.87601316108706159478E10,
903 3.88231277496238566008E12,
904 6.20557727146953693363E14,
905 6.87141087355300489866E16,
906 3.97270608116560655612E18,
909 static const double Z1 = 1.46819706421238932572E1;
910 static const double Z2 = 4.92184563216946036703E1;
916 double w, z, p, q, xn;
928 w = w * x * (z -
Z1) * (z -
Z2);
937 p = p * cos(xn) - w * q * sin(xn);
938 return p *
SQ2OPI / sqrt(x);
943 double w, z, p, q, xn;
951 fprintf(
ioQQQ,
"bessel_y1: domain error\n" );
965 p = p * sin(xn) + w * q * cos(xn);
966 return p *
SQ2OPI / sqrt(x);
1019 double pkm2, pkm1, pk, xk, r, ans;
1055 if( x < DBL_EPSILON )
1070 ans = pk - (xk/ans);
1083 pkm2 = (pkm1 * r - pk * x) / x;
1090 if( fabs(pk) > fabs(pkm1) )
1153 double an, anm1, anm2, r;
1181 fprintf(
ioQQQ,
"bessel_yn: domain error\n" );
1192 an = r * anm1 / x - anm2;
1287 1.37446543561352307156e-16,
1288 4.25981614279661018399e-14,
1289 1.03496952576338420167e-11,
1290 1.90451637722020886025e-9,
1291 2.53479107902614945675e-7,
1292 2.28621210311945178607e-5,
1293 1.26461541144692592338e-3,
1294 3.59799365153615016266e-2,
1295 3.44289899924628486886e-1,
1296 -5.35327393233902768720e-1
1306 5.30043377268626276149e-18,
1307 -1.64758043015242134646e-17,
1308 5.21039150503902756861e-17,
1309 -1.67823109680541210385e-16,
1310 5.51205597852431940784e-16,
1311 -1.84859337734377901440e-15,
1312 6.34007647740507060557e-15,
1313 -2.22751332699166985548e-14,
1314 8.03289077536357521100e-14,
1315 -2.98009692317273043925e-13,
1316 1.14034058820847496303e-12,
1317 -4.51459788337394416547e-12,
1318 1.85594911495471785253e-11,
1319 -7.95748924447710747776e-11,
1320 3.57739728140030116597e-10,
1321 -1.69753450938905987466e-9,
1322 8.57403401741422608519e-9,
1323 -4.66048989768794782956e-8,
1324 2.76681363944501510342e-7,
1325 -1.83175552271911948767e-6,
1326 1.39498137188764993662e-5,
1327 -1.28495495816278026384e-4,
1328 1.56988388573005337491e-3,
1329 -3.14481013119645005427e-2,
1330 2.44030308206595545468e0
1341 fprintf(
ioQQQ,
"bessel_k0: domain error\n" );
1352 y = exp(-x) *
chbevl( z,
k0_B, 25 ) / sqrt(x);
1364 fprintf(
ioQQQ,
"bessel_k0_scaled: domain error\n" );
1374 return chbevl( 8.0/x - 2.0,
k0_B, 25 ) / sqrt(x);
1461 -7.02386347938628759343e-18,
1462 -2.42744985051936593393e-15,
1463 -6.66690169419932900609e-13,
1464 -1.41148839263352776110e-10,
1465 -2.21338763073472585583e-8,
1466 -2.43340614156596823496e-6,
1467 -1.73028895751305206302e-4,
1468 -6.97572385963986435018e-3,
1469 -1.22611180822657148235e-1,
1470 -3.53155960776544875667e-1,
1471 1.52530022733894777053e0
1482 -5.75674448366501715755e-18,
1483 1.79405087314755922667e-17,
1484 -5.68946255844285935196e-17,
1485 1.83809354436663880070e-16,
1486 -6.05704724837331885336e-16,
1487 2.03870316562433424052e-15,
1488 -7.01983709041831346144e-15,
1489 2.47715442448130437068e-14,
1490 -8.97670518232499435011e-14,
1491 3.34841966607842919884e-13,
1492 -1.28917396095102890680e-12,
1493 5.13963967348173025100e-12,
1494 -2.12996783842756842877e-11,
1495 9.21831518760500529508e-11,
1496 -4.19035475934189648750e-10,
1497 2.01504975519703286596e-9,
1498 -1.03457624656780970260e-8,
1499 5.74108412545004946722e-8,
1500 -3.50196060308781257119e-7,
1501 2.40648494783721712015e-6,
1502 -1.93619797416608296024e-5,
1503 1.95215518471351631108e-4,
1504 -2.85781685962277938680e-3,
1505 1.03923736576817238437e-1,
1506 2.72062619048444266945e0
1518 fprintf(
ioQQQ,
"bessel_k1: domain error\n" );
1528 return exp(-x) *
chbevl( 8.0/x - 2.0,
k1_B, 25 ) / sqrt(x);
1539 fprintf(
ioQQQ,
"bessel_k1_scaled: domain error\n" );
1549 return chbevl( 8.0/x - 2.0,
k1_B, 25 ) / sqrt(x);
1633 -4.41534164647933937950e-18,
1634 3.33079451882223809783e-17,
1635 -2.43127984654795469359e-16,
1636 1.71539128555513303061e-15,
1637 -1.16853328779934516808e-14,
1638 7.67618549860493561688e-14,
1639 -4.85644678311192946090e-13,
1640 2.95505266312963983461e-12,
1641 -1.72682629144155570723e-11,
1642 9.67580903537323691224e-11,
1643 -5.18979560163526290666e-10,
1644 2.65982372468238665035e-9,
1645 -1.30002500998624804212e-8,
1646 6.04699502254191894932e-8,
1647 -2.67079385394061173391e-7,
1648 1.11738753912010371815e-6,
1649 -4.41673835845875056359e-6,
1650 1.64484480707288970893e-5,
1651 -5.75419501008210370398e-5,
1652 1.88502885095841655729e-4,
1653 -5.76375574538582365885e-4,
1654 1.63947561694133579842e-3,
1655 -4.32430999505057594430e-3,
1656 1.05464603945949983183e-2,
1657 -2.37374148058994688156e-2,
1658 4.93052842396707084878e-2,
1659 -9.49010970480476444210e-2,
1660 1.71620901522208775349e-1,
1661 -3.04682672343198398683e-1,
1662 6.76795274409476084995e-1
1673 -7.23318048787475395456e-18,
1674 -4.83050448594418207126e-18,
1675 4.46562142029675999901e-17,
1676 3.46122286769746109310e-17,
1677 -2.82762398051658348494e-16,
1678 -3.42548561967721913462e-16,
1679 1.77256013305652638360e-15,
1680 3.81168066935262242075e-15,
1681 -9.55484669882830764870e-15,
1682 -4.15056934728722208663e-14,
1683 1.54008621752140982691e-14,
1684 3.85277838274214270114e-13,
1685 7.18012445138366623367e-13,
1686 -1.79417853150680611778e-12,
1687 -1.32158118404477131188e-11,
1688 -3.14991652796324136454e-11,
1689 1.18891471078464383424e-11,
1690 4.94060238822496958910e-10,
1691 3.39623202570838634515e-9,
1692 2.26666899049817806459e-8,
1693 2.04891858946906374183e-7,
1694 2.89137052083475648297e-6,
1695 6.88975834691682398426e-5,
1696 3.36911647825569408990e-3,
1697 8.04490411014108831608e-1
1714 return exp(x) *
chbevl( 32.0/x - 2.0,
i0_B, 25 ) / sqrt(x);
1731 return chbevl( 32.0/x - 2.0,
i0_B, 25 ) / sqrt(x);
1816 2.77791411276104639959e-18,
1817 -2.11142121435816608115e-17,
1818 1.55363195773620046921e-16,
1819 -1.10559694773538630805e-15,
1820 7.60068429473540693410e-15,
1821 -5.04218550472791168711e-14,
1822 3.22379336594557470981e-13,
1823 -1.98397439776494371520e-12,
1824 1.17361862988909016308e-11,
1825 -6.66348972350202774223e-11,
1826 3.62559028155211703701e-10,
1827 -1.88724975172282928790e-9,
1828 9.38153738649577178388e-9,
1829 -4.44505912879632808065e-8,
1830 2.00329475355213526229e-7,
1831 -8.56872026469545474066e-7,
1832 3.47025130813767847674e-6,
1833 -1.32731636560394358279e-5,
1834 4.78156510755005422638e-5,
1835 -1.61760815825896745588e-4,
1836 5.12285956168575772895e-4,
1837 -1.51357245063125314899e-3,
1838 4.15642294431288815669e-3,
1839 -1.05640848946261981558e-2,
1840 2.47264490306265168283e-2,
1841 -5.29459812080949914269e-2,
1842 1.02643658689847095384e-1,
1843 -1.76416518357834055153e-1,
1844 2.52587186443633654823e-1
1855 7.51729631084210481353e-18,
1856 4.41434832307170791151e-18,
1857 -4.65030536848935832153e-17,
1858 -3.20952592199342395980e-17,
1859 2.96262899764595013876e-16,
1860 3.30820231092092828324e-16,
1861 -1.88035477551078244854e-15,
1862 -3.81440307243700780478e-15,
1863 1.04202769841288027642e-14,
1864 4.27244001671195135429e-14,
1865 -2.10154184277266431302e-14,
1866 -4.08355111109219731823e-13,
1867 -7.19855177624590851209e-13,
1868 2.03562854414708950722e-12,
1869 1.41258074366137813316e-11,
1870 3.25260358301548823856e-11,
1871 -1.89749581235054123450e-11,
1872 -5.58974346219658380687e-10,
1873 -3.83538038596423702205e-9,
1874 -2.63146884688951950684e-8,
1875 -2.51223623787020892529e-7,
1876 -3.88256480887769039346e-6,
1877 -1.10588938762623716291e-4,
1878 -9.76109749136146840777e-3,
1879 7.78576235018280120474e-1
1896 z = exp(z) *
chbevl( 32.0/z - 2.0,
i1_B, 25 ) / sqrt(z);
1917 z =
chbevl( 32.0/z - 2.0,
i1_B, 25 ) / sqrt(z);
1985 1.37982864606273237150e-4,
1986 2.28025724005875567385e-3,
1987 7.97404013220415179367e-3,
1988 9.85821379021226008714e-3,
1989 6.87489687449949877925e-3,
1990 6.18901033637687613229e-3,
1991 8.79078273952743772254e-3,
1992 1.49380448916805252718e-2,
1993 3.08851465246711995998e-2,
1994 9.65735902811690126535e-2,
1995 1.38629436111989062502e0
2000 2.94078955048598507511e-5,
2001 9.14184723865917226571e-4,
2002 5.94058303753167793257e-3,
2003 1.54850516649762399335e-2,
2004 2.39089602715924892727e-2,
2005 3.01204715227604046988e-2,
2006 3.73774314173823228969e-2,
2007 4.88280347570998239232e-2,
2008 7.03124996963957469739e-2,
2009 1.24999999999870820058e-1,
2010 4.99999999999999999821e-1
2013 static const double C1 = 1.3862943611198906188e0;
2019 if( x < 0.0 || x > 1.0 )
2021 fprintf(
ioQQQ,
"ellpk: domain error\n" );
2025 if( x > DBL_EPSILON )
2033 fprintf(
ioQQQ,
"ellpk: domain error\n" );
2038 return C1 - 0.5 * log(x);
2092 static const double BIG = 1.44115188075855872E+17;
2096 double ans, r, t, yk, xk;
2097 double pk, pkm1, pkm2, qk, qkm1, qkm2;
2103 if( n < 0 || x < 0. )
2105 fprintf(
ioQQQ,
"expn: domain error\n" );
2118 fprintf(
ioQQQ,
"expn: domain error\n" );
2123 return 1.0/((double)n-1.0);
2136 yk = 1.0 / (xk * xk);
2138 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
2139 ans = yk * (ans + t * (t - 2.0 * x));
2140 ans = yk * (ans + t);
2141 ans = (ans + 1.0) * exp( -x ) / xk;
2148 psi = -
EULER - log(x);
2149 for( i=1; i < n; i++ )
2174 while( t > DBL_EPSILON );
2194 xk =
static_cast<double>(n + (k-1)/2);
2199 xk =
static_cast<double>(k/2);
2201 pk = pkm1 * yk + pkm2 * xk;
2202 qk = qkm1 * yk + qkm2 * xk;
2206 t = fabs( (ans - r)/r );
2215 if( fabs(pk) >
BIG )
2223 while( t > DBL_EPSILON );
2280 inline double polevl(
double x,
const double coef[],
int N)
2284 const double *p = coef;
2290 ans = ans * x + *p++;
2302 inline double p1evl(
double x,
const double coef[],
int N)
2305 const double *p = coef;
2312 ans = ans * x + *p++;
2376 inline double chbevl(
double x,
const double array[],
int n)
2379 const double *p = array;
2390 b0 = x * b1 - b2 + *p++;
2453 #define MATRIX_A 0x9908b0dfUL
2454 #define UMASK 0x80000000UL
2455 #define LMASK 0x7fffffffUL
2456 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
2457 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
2468 state[0]= s & 0xffffffffUL;
2469 for (j=1; j<
N; j++) {
2475 state[j] &= 0xffffffffUL;
2489 k = (
N>key_length ?
N : key_length);
2493 state[i] &= 0xffffffffUL;
2496 if (j>=key_length) j=0;
2498 for (k=
N-1; k; k--) {
2501 state[i] &= 0xffffffffUL;
2506 state[0] = 0x80000000UL;
2512 unsigned long *p=
state;
2522 for (j=
N-
M+1; --j; p++)
2523 *p = p[
M] ^
TWIST(p[0], p[1]);
2526 *p = p[
M-
N] ^
TWIST(p[0], p[1]);
2541 y ^= (y << 7) & 0x9d2c5680UL;
2542 y ^= (y << 15) & 0xefc60000UL;
2558 y ^= (y << 7) & 0x9d2c5680UL;
2559 y ^= (y << 15) & 0xefc60000UL;
2562 return (
long)(y>>1);
2575 y ^= (y << 7) & 0x9d2c5680UL;
2576 y ^= (y << 15) & 0xefc60000UL;
2579 return (
double)y * (1.0/4294967295.0);
2593 y ^= (y << 7) & 0x9d2c5680UL;
2594 y ^= (y << 15) & 0xefc60000UL;
2597 return (
double)y * (1.0/4294967296.0);
2611 y ^= (y << 7) & 0x9d2c5680UL;
2612 y ^= (y << 15) & 0xefc60000UL;
2615 return ((
double)y + 0.5) * (1.0/4294967296.0);
2623 return(a*67108864.0+b)*(1.0/9007199254740992.0);