cloudy trunk
Loading...
Searching...
No Matches
thirdparty.cpp
Go to the documentation of this file.
1/* This file contains routines (perhaps in modified form) written by third parties.
2 * Use and distribution of these works are determined by their respective copyrights. */
3#include "cddefines.h"
4#include "thirdparty.h"
5#include "physconst.h"
6
7inline double polevl(double x, const double coef[], int N);
8inline double p1evl(double x, const double coef[], int N);
9inline double chbevl(double, const double[], int);
10inline double dawson(double x, int order);
11
12
13/* the routine linfit was derived from the program slopes.f */
14
15/********************************************************************/
16/* */
17/* PROGRAM SLOPES */
18/* */
19/* PROGRAM TO COMPUTE THE THEORETICAL REGRESSION SLOPES */
20/* AND UNCERTAINTIES (VIA DELTA METHOD), AND UNCERTAINTIES */
21/* VIA BOOTSTRAP AND BIVARIATE NORMAL SIMULATION FOR A */
22/* (X(I),Y(I)) DATA SET WITH UNKNOWN POPULATION DISTRIBUTION. */
23/* */
24/* WRITTEN BY ERIC FEIGELSON, PENN STATE. JUNE 1991 */
25/* */
26/* */
27/* A full description of these methods can be found in: */
28/* Isobe, T., Feigelson, E. D., Akritas, M. and Babu, G. J., */
29/* Linear Regression in Astronomy I, Astrophys. J. 364, 104 */
30/* (1990) */
31/* Babu, G. J. and Feigelson, E. D., Analytical and Monte Carlo */
32/* Comparisons of Six Different Linear Least Squares Fits, */
33/* Communications in Statistics, Simulation & Computation, */
34/* 21, 533 (1992) */
35/* Feigelson, E. D. and Babu, G. J., Linear Regression in */
36/* Astronomy II, Astrophys. J. 397, 55 (1992). */
37/* */
38/********************************************************************/
39
40/* this used to be sixlin, but only the first fit has been retained */
41
42/********************************************************************/
43/************************* routine linfit ***************************/
44/********************************************************************/
45
46bool linfit(
47 long n,
48 const double xorg[], /* x[n] */
49 const double yorg[], /* y[n] */
50 double &a,
51 double &siga,
52 double &b,
53 double &sigb
54)
55{
56
57/*
58 * linear regression
59 * written by t. isobe, g. j. babu and e. d. feigelson
60 * center for space research, m.i.t.
61 * and
62 * the pennsylvania state university
63 *
64 * rev. 1.0, september 1990
65 *
66 * this subroutine provides linear regression coefficients
67 * computed by six different methods described in isobe,
68 * feigelson, akritas, and babu 1990, astrophysical journal
69 * and babu and feigelson 1990, subm. to technometrics.
70 * the methods are ols(y/x), ols(x/y), ols bisector, orthogonal,
71 * reduced major axis, and mean-ols regressions.
72 *
73 * [Modified and translated to C/C++ by Peter van Hoof, Royal
74 * Observatory of Belgium; only the first method has been retained]
75 *
76 *
77 * input
78 * x[i] : independent variable
79 * y[i] : dependent variable
80 * n : number of data points
81 *
82 * output
83 * a : intercept coefficients
84 * b : slope coefficients
85 * siga : standard deviations of intercepts
86 * sigb : standard deviations of slopes
87 *
88 * error returns
89 * returns true when division by zero will
90 * occur (i.e. when sxy is zero)
91 */
92
93 DEBUG_ENTRY( "linfit()" );
94
95 ASSERT( n >= 2 );
96
97 valarray<double> x(n);
98 valarray<double> y(n);
99
100 for( long i=0; i < n; i++ )
101 {
102 x[i] = xorg[i];
103 y[i] = yorg[i];
104 }
105
106 /* initializations */
107 a = 0.0;
108 siga = 0.0;
109 b = 0.0;
110 sigb = 0.0;
111
112 /* compute averages and sums */
113 double s1 = 0.0;
114 double s2 = 0.0;
115 for( long i=0; i < n; i++ )
116 {
117 s1 += x[i];
118 s2 += y[i];
119 }
120 double rn = (double)n;
121 double xavg = s1/rn;
122 double yavg = s2/rn;
123 double sxx = 0.0;
124 double sxy = 0.0;
125 for( long i=0; i < n; i++ )
126 {
127 x[i] -= xavg;
128 y[i] -= yavg;
129 sxx += pow2(x[i]);
130 sxy += x[i]*y[i];
131 }
132
133 if( pow2(sxx) == 0.0 )
134 {
135 return true;
136 }
137
138 /* compute the slope coefficient */
139 b = sxy/sxx;
140
141 /* compute intercept coefficient */
142 a = yavg - b*xavg;
143
144 /* prepare for computation of variances */
145 double sum1 = 0.0;
146 for( long i=0; i < n; i++ )
147 sum1 += pow2(x[i]*(y[i] - b*x[i]));
148
149 /* compute variance of the slope coefficient */
150 sigb = sum1/pow2(sxx);
151
152 /* compute variance of the intercept coefficient */
153 for( long i=0; i < n; i++ )
154 siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
155
156 /* convert variances to standard deviations */
157 sigb = sqrt(sigb);
158 siga = sqrt(siga)/rn;
159
160 /* return data arrays to their original form */
161 for( long i=0; i < n; i++ )
162 {
163 x[i] += xavg;
164 y[i] += yavg;
165 }
166 return false;
167}
168
169/************************************************************************
170 * This marks the end of the block of code from Isobe, Babu & Feigelson *
171 ************************************************************************/
172
173
174/* the routines factorial and lfactorial came originally from hydro_bauman.cpp
175 * and were written by Robert Paul Bauman. lfactorial was modified by Peter van Hoof */
176
177/************************************************************************************************/
178/* pre-calculated factorials */
179/************************************************************************************************/
180
181static const double pre_factorial[NPRE_FACTORIAL] =
182{
183 1.00000000000000000000e+00,
184 1.00000000000000000000e+00,
185 2.00000000000000000000e+00,
186 6.00000000000000000000e+00,
187 2.40000000000000000000e+01,
188 1.20000000000000000000e+02,
189 7.20000000000000000000e+02,
190 5.04000000000000000000e+03,
191 4.03200000000000000000e+04,
192 3.62880000000000000000e+05,
193 3.62880000000000000000e+06,
194 3.99168000000000000000e+07,
195 4.79001600000000000000e+08,
196 6.22702080000000000000e+09,
197 8.71782912000000000000e+10,
198 1.30767436800000000000e+12,
199 2.09227898880000000000e+13,
200 3.55687428096000000000e+14,
201 6.40237370572800000000e+15,
202 1.21645100408832000000e+17,
203 2.43290200817664000000e+18,
204 5.10909421717094400000e+19,
205 1.12400072777760768000e+21,
206 2.58520167388849766400e+22,
207 6.20448401733239439360e+23,
208 1.55112100433309859840e+25,
209 4.03291461126605635592e+26,
210 1.08888694504183521614e+28,
211 3.04888344611713860511e+29,
212 8.84176199373970195470e+30,
213 2.65252859812191058647e+32,
214 8.22283865417792281807e+33,
215 2.63130836933693530178e+35,
216 8.68331761881188649615e+36,
217 2.95232799039604140861e+38,
218 1.03331479663861449300e+40,
219 3.71993326789901217463e+41,
220 1.37637530912263450457e+43,
221 5.23022617466601111726e+44,
222 2.03978820811974433568e+46,
223 8.15915283247897734264e+47,
224 3.34525266131638071044e+49,
225 1.40500611775287989839e+51,
226 6.04152630633738356321e+52,
227 2.65827157478844876773e+54,
228 1.19622220865480194551e+56,
229 5.50262215981208894940e+57,
230 2.58623241511168180614e+59,
231 1.24139155925360726691e+61,
232 6.08281864034267560801e+62,
233 3.04140932017133780398e+64,
234 1.55111875328738227999e+66,
235 8.06581751709438785585e+67,
236 4.27488328406002556374e+69,
237 2.30843697339241380441e+71,
238 1.26964033536582759243e+73,
239 7.10998587804863451749e+74,
240 4.05269195048772167487e+76,
241 2.35056133128287857145e+78,
242 1.38683118545689835713e+80,
243 8.32098711274139014271e+81,
244 5.07580213877224798711e+83,
245 3.14699732603879375200e+85,
246 1.98260831540444006372e+87,
247 1.26886932185884164078e+89,
248 8.24765059208247066512e+90,
249 5.44344939077443063905e+92,
250 3.64711109181886852801e+94,
251 2.48003554243683059915e+96,
252 1.71122452428141311337e+98,
253 1.19785716699698917933e+100,
254 8.50478588567862317347e+101,
255 6.12344583768860868500e+103,
256 4.47011546151268434004e+105,
257 3.30788544151938641157e+107,
258 2.48091408113953980872e+109,
259 1.88549470166605025458e+111,
260 1.45183092028285869606e+113,
261 1.13242811782062978295e+115,
262 8.94618213078297528506e+116,
263 7.15694570462638022794e+118,
264 5.79712602074736798470e+120,
265 4.75364333701284174746e+122,
266 3.94552396972065865030e+124,
267 3.31424013456535326627e+126,
268 2.81710411438055027626e+128,
269 2.42270953836727323750e+130,
270 2.10775729837952771662e+132,
271 1.85482642257398439069e+134,
272 1.65079551609084610774e+136,
273 1.48571596448176149700e+138,
274 1.35200152767840296226e+140,
275 1.24384140546413072522e+142,
276 1.15677250708164157442e+144,
277 1.08736615665674307994e+146,
278 1.03299784882390592592e+148,
279 9.91677934870949688836e+149,
280 9.61927596824821198159e+151,
281 9.42689044888324774164e+153,
282 9.33262154439441526381e+155,
283 9.33262154439441526388e+157,
284 9.42594775983835941673e+159,
285 9.61446671503512660515e+161,
286 9.90290071648618040340e+163,
287 1.02990167451456276198e+166,
288 1.08139675824029090008e+168,
289 1.14628056373470835406e+170,
290 1.22652020319613793888e+172,
291 1.32464181945182897396e+174,
292 1.44385958320249358163e+176,
293 1.58824554152274293982e+178,
294 1.76295255109024466316e+180,
295 1.97450685722107402277e+182,
296 2.23119274865981364576e+184,
297 2.54355973347218755612e+186,
298 2.92509369349301568964e+188,
299 3.39310868445189820004e+190,
300 3.96993716080872089396e+192,
301 4.68452584975429065488e+194,
302 5.57458576120760587943e+196,
303 6.68950291344912705515e+198,
304 8.09429852527344373681e+200,
305 9.87504420083360135884e+202,
306 1.21463043670253296712e+205,
307 1.50614174151114087918e+207,
308 1.88267717688892609901e+209,
309 2.37217324288004688470e+211,
310 3.01266001845765954361e+213,
311 3.85620482362580421582e+215,
312 4.97450422247728743840e+217,
313 6.46685548922047366972e+219,
314 8.47158069087882050755e+221,
315 1.11824865119600430699e+224,
316 1.48727070609068572828e+226,
317 1.99294274616151887582e+228,
318 2.69047270731805048244e+230,
319 3.65904288195254865604e+232,
320 5.01288874827499165889e+234,
321 6.91778647261948848943e+236,
322 9.61572319694108900019e+238,
323 1.34620124757175246000e+241,
324 1.89814375907617096864e+243,
325 2.69536413788816277557e+245,
326 3.85437071718007276916e+247,
327 5.55029383273930478744e+249,
328 8.04792605747199194159e+251,
329 1.17499720439091082343e+254,
330 1.72724589045463891049e+256,
331 2.55632391787286558753e+258,
332 3.80892263763056972532e+260,
333 5.71338395644585458806e+262,
334 8.62720977423324042775e+264,
335 1.31133588568345254503e+267,
336 2.00634390509568239384e+269,
337 3.08976961384735088657e+271,
338 4.78914290146339387432e+273,
339 7.47106292628289444390e+275,
340 1.17295687942641442768e+278,
341 1.85327186949373479574e+280,
342 2.94670227249503832518e+282,
343 4.71472363599206132029e+284,
344 7.59070505394721872577e+286,
345 1.22969421873944943358e+289,
346 2.00440157654530257674e+291,
347 3.28721858553429622598e+293,
348 5.42391066613158877297e+295,
349 9.00369170577843736335e+297,
350 1.50361651486499903974e+300,
351 2.52607574497319838672e+302,
352 4.26906800900470527345e+304,
353 7.25741561530799896496e+306
354};
355
356double factorial( long n )
357{
358 DEBUG_ENTRY( "factorial()" );
359
360 if( n < 0 || n >= NPRE_FACTORIAL )
361 {
362 fprintf( ioQQQ, "factorial: domain error\n" );
364 }
365
366 return pre_factorial[n];
367}
368
369/* NB - this implementation is not thread-safe! */
370
371class t_lfact : public Singleton<t_lfact>
372{
373 friend class Singleton<t_lfact>;
374protected:
376 {
377 p_lf.reserve( 512 );
378 p_lf.push_back( 0. ); /* log10( 0! ) */
379 p_lf.push_back( 0. ); /* log10( 1! ) */
380 }
381private:
382 vector<double> p_lf;
383public:
384 double get_lfact( unsigned long n )
385 {
386 if( n < p_lf.size() )
387 {
388 return p_lf[n];
389 }
390 else
391 {
392 for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ )
393 p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) );
394 return p_lf[n];
395 }
396 }
397};
398
399double lfactorial( long n )
400{
401 /******************************/
402 /* n */
403 /* --- */
404 /* log( n! ) = > log(j) */
405 /* --- */
406 /* j=1 */
407 /******************************/
408
409 DEBUG_ENTRY( "lfactorial()" );
410
411 if( n < 0 )
412 {
413 fprintf( ioQQQ, "lfactorial: domain error\n" );
415 }
416
417 return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
418}
419
420/*******************************************************************
421 * This marks the end of the block of code from Robert Paul Bauman *
422 *******************************************************************/
423
424
425/* complex Gamma function in double precision */
426/* this routine is a slightly modified version of the one found
427 * at http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html
428 * The following copyright applies:
429 * Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
430 * You may use, copy, modify this code for any purpose and
431 * without fee. You may distribute this ORIGINAL package. */
432complex<double> cdgamma(complex<double> x)
433{
434 double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
435
436 DEBUG_ENTRY( "cdgamma()" );
437
438 xr = x.real();
439 xi = x.imag();
440 if(xr < 0)
441 {
442 wr = 1. - xr;
443 wi = -xi;
444 }
445 else
446 {
447 wr = xr;
448 wi = xi;
449 }
450 ur = wr + 6.00009857740312429;
451 vr = ur * (wr + 4.99999857982434025) - wi * wi;
452 vi = wi * (wr + 4.99999857982434025) + ur * wi;
453 yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
454 0.293729529320536228;
455 yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
456 ur = vr * (wr + 4.00000003016801681) - vi * wi;
457 ui = vi * (wr + 4.00000003016801681) + vr * wi;
458 vr = ur * (wr + 2.99999999944915534) - ui * wi;
459 vi = ui * (wr + 2.99999999944915534) + ur * wi;
460 yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
461 yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
462 ur = vr * (wr + 2.00000000000603851) - vi * wi;
463 ui = vi * (wr + 2.00000000000603851) + vr * wi;
464 vr = ur * (wr + 0.999999999999975753) - ui * wi;
465 vi = ui * (wr + 0.999999999999975753) + ur * wi;
466 yr += ur * 10.5400280458730808 + vr;
467 yi += ui * 10.5400280458730808 + vi;
468 ur = vr * wr - vi * wi;
469 ui = vi * wr + vr * wi;
470 t = ur * ur + ui * ui;
471 vr = yr * ur + yi * ui + t * 0.0327673720261526849;
472 vi = yi * ur - yr * ui;
473 yr = wr + 7.31790632447016203;
474 ur = log(yr * yr + wi * wi) * 0.5 - 1;
475 ui = atan2(wi, yr);
476 yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
477 yi = ui * (wr - 0.5) + ur * wi;
478 ur = yr * cos(yi);
479 ui = yr * sin(yi);
480 yr = ur * vr - ui * vi;
481 yi = ui * vr + ur * vi;
482 if(xr < 0)
483 {
484 wr = xr * 3.14159265358979324;
485 wi = exp(xi * 3.14159265358979324);
486 vi = 1 / wi;
487 ur = (vi + wi) * sin(wr);
488 ui = (vi - wi) * cos(wr);
489 vr = ur * yr + ui * yi;
490 vi = ui * yr - ur * yi;
491 ur = 6.2831853071795862 / (vr * vr + vi * vi);
492 yr = ur * vr;
493 yi = ur * vi;
494 }
495 return complex<double>( yr, yi );
496}
497
498/*************************************************************
499 * This marks the end of the block of code from Takuya OOURA *
500 *************************************************************/
501
502/*====================================================================
503 *
504 * Below are routines from the Cephes library.
505 *
506 * This is the copyright statement included with the library:
507 *
508 * Some software in this archive may be from the book _Methods and
509 * Programs for Mathematical Functions_ (Prentice-Hall, 1989) or
510 * from the Cephes Mathematical Library, a commercial product. In
511 * either event, it is copyrighted by the author. What you see here
512 * may be used freely but it comes with no support or guarantee.
513 *
514 * The two known misprints in the book are repaired here in the
515 * source listings for the gamma function and the incomplete beta
516 * integral.
517 *
518 * Stephen L. Moshier
519 * moshier@world.std.com
520 *
521 *====================================================================*/
522
523/* j0.c
524 *
525 * Bessel function of order zero
526 *
527 *
528 *
529 * SYNOPSIS:
530 *
531 * double x, y, j0();
532 *
533 * y = j0( x );
534 *
535 *
536 *
537 * DESCRIPTION:
538 *
539 * Returns Bessel function of order zero of the argument.
540 *
541 * The domain is divided into the intervals [0, 5] and
542 * (5, infinity). In the first interval the following rational
543 * approximation is used:
544 *
545 *
546 * 2 2
547 * (w - r ) (w - r ) P (w) / Q (w)
548 * 1 2 3 8
549 *
550 * 2
551 * where w = x and the two r's are zeros of the function.
552 *
553 * In the second interval, the Hankel asymptotic expansion
554 * is employed with two rational functions of degree 6/6
555 * and 7/7.
556 *
557 *
558 *
559 * ACCURACY:
560 *
561 * Absolute error:
562 * arithmetic domain # trials peak rms
563 * DEC 0, 30 10000 4.4e-17 6.3e-18
564 * IEEE 0, 30 60000 4.2e-16 1.1e-16
565 *
566 */
567/* y0.c
568 *
569 * Bessel function of the second kind, order zero
570 *
571 *
572 *
573 * SYNOPSIS:
574 *
575 * double x, y, y0();
576 *
577 * y = y0( x );
578 *
579 *
580 *
581 * DESCRIPTION:
582 *
583 * Returns Bessel function of the second kind, of order
584 * zero, of the argument.
585 *
586 * The domain is divided into the intervals [0, 5] and
587 * (5, infinity). In the first interval a rational approximation
588 * R(x) is employed to compute
589 * y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
590 * Thus a call to j0() is required.
591 *
592 * In the second interval, the Hankel asymptotic expansion
593 * is employed with two rational functions of degree 6/6
594 * and 7/7.
595 *
596 *
597 *
598 * ACCURACY:
599 *
600 * Absolute error, when y0(x) < 1; else relative error:
601 *
602 * arithmetic domain # trials peak rms
603 * DEC 0, 30 9400 7.0e-17 7.9e-18
604 * IEEE 0, 30 30000 1.3e-15 1.6e-16
605 *
606 */
607
608/*
609Cephes Math Library Release 2.8: June, 2000
610Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
611*/
612
613/* Note: all coefficients satisfy the relative error criterion
614 * except YP, YQ which are designed for absolute error. */
615
616static const double b0_PP[7] = {
617 7.96936729297347051624e-4,
618 8.28352392107440799803e-2,
619 1.23953371646414299388e0,
620 5.44725003058768775090e0,
621 8.74716500199817011941e0,
622 5.30324038235394892183e0,
623 9.99999999999999997821e-1,
624};
625
626static const double b0_PQ[7] = {
627 9.24408810558863637013e-4,
628 8.56288474354474431428e-2,
629 1.25352743901058953537e0,
630 5.47097740330417105182e0,
631 8.76190883237069594232e0,
632 5.30605288235394617618e0,
633 1.00000000000000000218e0,
634};
635
636static const double b0_QP[8] = {
637 -1.13663838898469149931e-2,
638 -1.28252718670509318512e0,
639 -1.95539544257735972385e1,
640 -9.32060152123768231369e1,
641 -1.77681167980488050595e2,
642 -1.47077505154951170175e2,
643 -5.14105326766599330220e1,
644 -6.05014350600728481186e0,
645};
646
647static const double b0_QQ[7] = {
648 /* 1.00000000000000000000e0,*/
649 6.43178256118178023184e1,
650 8.56430025976980587198e2,
651 3.88240183605401609683e3,
652 7.24046774195652478189e3,
653 5.93072701187316984827e3,
654 2.06209331660327847417e3,
655 2.42005740240291393179e2,
656};
657
658static const double b0_YP[8] = {
659 1.55924367855235737965e4,
660 -1.46639295903971606143e7,
661 5.43526477051876500413e9,
662 -9.82136065717911466409e11,
663 8.75906394395366999549e13,
664 -3.46628303384729719441e15,
665 4.42733268572569800351e16,
666 -1.84950800436986690637e16,
667};
668
669static const double b0_YQ[7] = {
670 /* 1.00000000000000000000e0,*/
671 1.04128353664259848412e3,
672 6.26107330137134956842e5,
673 2.68919633393814121987e8,
674 8.64002487103935000337e10,
675 2.02979612750105546709e13,
676 3.17157752842975028269e15,
677 2.50596256172653059228e17,
678};
679
680/* 5.783185962946784521175995758455807035071 */
681static const double DR1 = 5.78318596294678452118e0;
682/* 30.47126234366208639907816317502275584842 */
683static const double DR2 = 3.04712623436620863991e1;
684
685static double b0_RP[4] = {
686 -4.79443220978201773821e9,
687 1.95617491946556577543e12,
688 -2.49248344360967716204e14,
689 9.70862251047306323952e15,
690};
691
692static double b0_RQ[8] = {
693 /* 1.00000000000000000000e0,*/
694 4.99563147152651017219e2,
695 1.73785401676374683123e5,
696 4.84409658339962045305e7,
697 1.11855537045356834862e10,
698 2.11277520115489217587e12,
699 3.10518229857422583814e14,
700 3.18121955943204943306e16,
701 1.71086294081043136091e18,
702};
703
704static const double TWOOPI = 2./PI;
705static const double SQ2OPI = sqrt(2./PI);
706static const double PIO4 = PI/4.;
707
708double bessel_j0(double x)
709{
710 double w, z, p, q, xn;
711
712 DEBUG_ENTRY( "bessel_j0()" );
713
714 if( x < 0 )
715 x = -x;
716
717 if( x <= 5.0 )
718 {
719 z = x * x;
720 if( x < 1.0e-5 )
721 return 1.0 - z/4.0;
722
723 p = (z - DR1) * (z - DR2);
724 p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
725 return p;
726 }
727
728 w = 5.0/x;
729 q = 25.0/(x*x);
730 p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
731 q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
732 xn = x - PIO4;
733 p = p * cos(xn) - w * q * sin(xn);
734 return p * SQ2OPI / sqrt(x);
735}
736
737/* y0() 2 */
738/* Bessel function of second kind, order zero */
739
740/* Rational approximation coefficients YP[], YQ[] are used here.
741 * The function computed is y0(x) - 2 * log(x) * j0(x) / PI,
742 * whose value at x = 0 is 2 * ( log(0.5) + EULER ) / PI
743 * = 0.073804295108687225.
744 */
745
746double bessel_y0(double x)
747{
748 double w, z, p, q, xn;
749
750 DEBUG_ENTRY( "bessel_y0()" );
751
752 if( x <= 5.0 )
753 {
754 if( x <= 0.0 )
755 {
756 fprintf( ioQQQ, "bessel_y0: domain error\n" );
758 }
759 z = x * x;
760 w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
761 w += TWOOPI * log(x) * bessel_j0(x);
762 return w;
763 }
764
765 w = 5.0/x;
766 z = 25.0 / (x * x);
767 p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
768 q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
769 xn = x - PIO4;
770 p = p * sin(xn) + w * q * cos(xn);
771 return p * SQ2OPI / sqrt(x);
772}
773
774/* j1.c
775 *
776 * Bessel function of order one
777 *
778 *
779 *
780 * SYNOPSIS:
781 *
782 * double x, y, j1();
783 *
784 * y = j1( x );
785 *
786 *
787 *
788 * DESCRIPTION:
789 *
790 * Returns Bessel function of order one of the argument.
791 *
792 * The domain is divided into the intervals [0, 8] and
793 * (8, infinity). In the first interval a 24 term Chebyshev
794 * expansion is used. In the second, the asymptotic
795 * trigonometric representation is employed using two
796 * rational functions of degree 5/5.
797 *
798 *
799 *
800 * ACCURACY:
801 *
802 * Absolute error:
803 * arithmetic domain # trials peak rms
804 * DEC 0, 30 10000 4.0e-17 1.1e-17
805 * IEEE 0, 30 30000 2.6e-16 1.1e-16
806 *
807 *
808 */
809/* y1.c
810 *
811 * Bessel function of second kind of order one
812 *
813 *
814 *
815 * SYNOPSIS:
816 *
817 * double x, y, y1();
818 *
819 * y = y1( x );
820 *
821 *
822 *
823 * DESCRIPTION:
824 *
825 * Returns Bessel function of the second kind of order one
826 * of the argument.
827 *
828 * The domain is divided into the intervals [0, 8] and
829 * (8, infinity). In the first interval a 25 term Chebyshev
830 * expansion is used, and a call to j1() is required.
831 * In the second, the asymptotic trigonometric representation
832 * is employed using two rational functions of degree 5/5.
833 *
834 *
835 *
836 * ACCURACY:
837 *
838 * Absolute error:
839 * arithmetic domain # trials peak rms
840 * DEC 0, 30 10000 8.6e-17 1.3e-17
841 * IEEE 0, 30 30000 1.0e-15 1.3e-16
842 *
843 * (error criterion relative when |y1| > 1).
844 *
845 */
846
847/*
848Cephes Math Library Release 2.8: June, 2000
849Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
850*/
851
852static const double b1_RP[4] = {
853 -8.99971225705559398224e8,
854 4.52228297998194034323e11,
855 -7.27494245221818276015e13,
856 3.68295732863852883286e15,
857};
858
859static const double b1_RQ[8] = {
860 /* 1.00000000000000000000E0,*/
861 6.20836478118054335476e2,
862 2.56987256757748830383e5,
863 8.35146791431949253037e7,
864 2.21511595479792499675e10,
865 4.74914122079991414898e12,
866 7.84369607876235854894e14,
867 8.95222336184627338078e16,
868 5.32278620332680085395e18,
869};
870
871static const double b1_PP[7] = {
872 7.62125616208173112003e-4,
873 7.31397056940917570436e-2,
874 1.12719608129684925192e0,
875 5.11207951146807644818e0,
876 8.42404590141772420927e0,
877 5.21451598682361504063e0,
878 1.00000000000000000254e0,
879};
880
881static const double b1_PQ[7] = {
882 5.71323128072548699714e-4,
883 6.88455908754495404082e-2,
884 1.10514232634061696926e0,
885 5.07386386128601488557e0,
886 8.39985554327604159757e0,
887 5.20982848682361821619e0,
888 9.99999999999999997461e-1,
889};
890
891static const double b1_QP[8] = {
892 5.10862594750176621635e-2,
893 4.98213872951233449420e0,
894 7.58238284132545283818e1,
895 3.66779609360150777800e2,
896 7.10856304998926107277e2,
897 5.97489612400613639965e2,
898 2.11688757100572135698e2,
899 2.52070205858023719784e1,
900};
901
902static const double b1_QQ[7] = {
903 /* 1.00000000000000000000e0,*/
904 7.42373277035675149943e1,
905 1.05644886038262816351e3,
906 4.98641058337653607651e3,
907 9.56231892404756170795e3,
908 7.99704160447350683650e3,
909 2.82619278517639096600e3,
910 3.36093607810698293419e2,
911};
912
913static const double b1_YP[6] = {
914 1.26320474790178026440e9,
915 -6.47355876379160291031e11,
916 1.14509511541823727583e14,
917 -8.12770255501325109621e15,
918 2.02439475713594898196e17,
919 -7.78877196265950026825e17,
920};
921
922static const double b1_YQ[8] = {
923 /* 1.00000000000000000000E0,*/
924 5.94301592346128195359E2,
925 2.35564092943068577943E5,
926 7.34811944459721705660E7,
927 1.87601316108706159478E10,
928 3.88231277496238566008E12,
929 6.20557727146953693363E14,
930 6.87141087355300489866E16,
931 3.97270608116560655612E18,
932};
933
934static const double Z1 = 1.46819706421238932572E1;
935static const double Z2 = 4.92184563216946036703E1;
936
937static const double THPIO4 = 3.*PI/4.;
938
939double bessel_j1(double x)
940{
941 double w, z, p, q, xn;
942
943 DEBUG_ENTRY( "bessel_j1()" );
944
945 w = x;
946 if( x < 0 )
947 w = -x;
948
949 if( w <= 5.0 )
950 {
951 z = x * x;
952 w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
953 w = w * x * (z - Z1) * (z - Z2);
954 return w;
955 }
956
957 w = 5.0/x;
958 z = w * w;
959 p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
960 q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
961 xn = x - THPIO4;
962 p = p * cos(xn) - w * q * sin(xn);
963 return p * SQ2OPI / sqrt(x);
964}
965
966double bessel_y1(double x)
967{
968 double w, z, p, q, xn;
969
970 DEBUG_ENTRY( "bessel_y1()" );
971
972 if( x <= 5.0 )
973 {
974 if( x <= 0.0 )
975 {
976 fprintf( ioQQQ, "bessel_y1: domain error\n" );
978 }
979 z = x * x;
980 w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
981 w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
982 return w;
983 }
984
985 w = 5.0/x;
986 z = w * w;
987 p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
988 q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
989 xn = x - THPIO4;
990 p = p * sin(xn) + w * q * cos(xn);
991 return p * SQ2OPI / sqrt(x);
992}
993
994/* jn.c
995 *
996 * Bessel function of integer order
997 *
998 *
999 *
1000 * SYNOPSIS:
1001 *
1002 * int n;
1003 * double x, y, jn();
1004 *
1005 * y = jn( n, x );
1006 *
1007 *
1008 *
1009 * DESCRIPTION:
1010 *
1011 * Returns Bessel function of order n, where n is a
1012 * (possibly negative) integer.
1013 *
1014 * The ratio of jn(x) to j0(x) is computed by backward
1015 * recurrence. First the ratio jn/jn-1 is found by a
1016 * continued fraction expansion. Then the recurrence
1017 * relating successive orders is applied until j0 or j1 is
1018 * reached.
1019 *
1020 * If n = 0 or 1 the routine for j0 or j1 is called
1021 * directly.
1022 *
1023 *
1024 *
1025 * ACCURACY:
1026 *
1027 * Absolute error:
1028 * arithmetic range # trials peak rms
1029 * DEC 0, 30 5500 6.9e-17 9.3e-18
1030 * IEEE 0, 30 5000 4.4e-16 7.9e-17
1031 *
1032 *
1033 * Not suitable for large n or x. Use jv() instead.
1034 *
1035 */
1036
1037/* jn.c
1038Cephes Math Library Release 2.8: June, 2000
1039Copyright 1984, 1987, 2000 by Stephen L. Moshier
1040*/
1041
1042double bessel_jn(int n, double x)
1043{
1044 double pkm2, pkm1, pk, xk, r, ans;
1045 int k, sign;
1046
1047 DEBUG_ENTRY( "bessel_jn()" );
1048
1049 if( n < 0 )
1050 {
1051 n = -n;
1052 if( (n & 1) == 0 ) /* -1**n */
1053 sign = 1;
1054 else
1055 sign = -1;
1056 }
1057 else
1058 sign = 1;
1059
1060 if( x < 0.0 )
1061 {
1062 if( n & 1 )
1063 sign = -sign;
1064 x = -x;
1065 }
1066
1067 if( x < DBL_EPSILON )
1068 {
1069 return sign * powi(x/2.,n)/factorial(n);
1070 }
1071
1072 if( n == 0 )
1073 {
1074 return sign * bessel_j0(x);
1075 }
1076 if( n == 1 )
1077 {
1078 return sign * bessel_j1(x);
1079 }
1080 // avoid cancellation error for very small x
1081 if( n == 2 && x > 0.1 )
1082 {
1083 return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
1084 }
1085
1086 /* continued fraction */
1087 k = 52;
1088
1089 pk = 2 * (n + k);
1090 ans = pk;
1091 xk = x * x;
1092
1093 do
1094 {
1095 pk -= 2.0;
1096 ans = pk - (xk/ans);
1097 }
1098 while( --k > 0 );
1099 ans = x/ans;
1100
1101 /* backward recurrence */
1102 pk = 1.0;
1103 pkm1 = 1.0/ans;
1104 k = n-1;
1105 r = 2 * k;
1106
1107 do
1108 {
1109 pkm2 = (pkm1 * r - pk * x) / x;
1110 pk = pkm1;
1111 pkm1 = pkm2;
1112 r -= 2.0;
1113 }
1114 while( --k > 0 );
1115
1116 if( fabs(pk) > fabs(pkm1) )
1117 ans = bessel_j1(x)/pk;
1118 else
1119 ans = bessel_j0(x)/pkm1;
1120 return sign * ans;
1121}
1122
1123/* yn.c
1124 *
1125 * Bessel function of second kind of integer order
1126 *
1127 *
1128 *
1129 * SYNOPSIS:
1130 *
1131 * double x, y, yn();
1132 * int n;
1133 *
1134 * y = yn( n, x );
1135 *
1136 *
1137 *
1138 * DESCRIPTION:
1139 *
1140 * Returns Bessel function of order n, where n is a
1141 * (possibly negative) integer.
1142 *
1143 * The function is evaluated by forward recurrence on
1144 * n, starting with values computed by the routines
1145 * y0() and y1().
1146 *
1147 * If n = 0 or 1 the routine for y0 or y1 is called
1148 * directly.
1149 *
1150 *
1151 *
1152 * ACCURACY:
1153 *
1154 *
1155 * Absolute error, except relative
1156 * when y > 1:
1157 * arithmetic domain # trials peak rms
1158 * DEC 0, 30 2200 2.9e-16 5.3e-17
1159 * IEEE 0, 30 30000 3.4e-15 4.3e-16
1160 *
1161 *
1162 * ERROR MESSAGES:
1163 *
1164 * message condition value returned
1165 * yn singularity x = 0 MAXNUM
1166 * yn overflow MAXNUM
1167 *
1168 * Spot checked against tables for x, n between 0 and 100.
1169 *
1170 */
1171
1172/*
1173Cephes Math Library Release 2.8: June, 2000
1174Copyright 1984, 1987, 2000 by Stephen L. Moshier
1175*/
1176
1177double bessel_yn(int n, double x)
1178{
1179 double an, anm1, anm2, r;
1180 int k, sign;
1181
1182 DEBUG_ENTRY( "bessel_yn()" );
1183
1184 if( n < 0 )
1185 {
1186 n = -n;
1187 if( (n & 1) == 0 ) /* -1**n */
1188 sign = 1;
1189 else
1190 sign = -1;
1191 }
1192 else
1193 sign = 1;
1194
1195 if( n == 0 )
1196 {
1197 return sign * bessel_y0(x);
1198 }
1199 if( n == 1 )
1200 {
1201 return sign * bessel_y1(x);
1202 }
1203
1204 /* test for overflow */
1205 if( x <= 0.0 )
1206 {
1207 fprintf( ioQQQ, "bessel_yn: domain error\n" );
1209 }
1210
1211 /* forward recurrence on n */
1212 anm2 = bessel_y0(x);
1213 anm1 = bessel_y1(x);
1214 k = 1;
1215 r = 2 * k;
1216 do
1217 {
1218 an = r * anm1 / x - anm2;
1219 anm2 = anm1;
1220 anm1 = an;
1221 r += 2.0;
1222 ++k;
1223 }
1224 while( k < n );
1225 return sign * an;
1226}
1227
1228/* k0.c
1229 *
1230 * Modified Bessel function, third kind, order zero
1231 *
1232 *
1233 *
1234 * SYNOPSIS:
1235 *
1236 * double x, y, k0();
1237 *
1238 * y = k0( x );
1239 *
1240 *
1241 *
1242 * DESCRIPTION:
1243 *
1244 * Returns modified Bessel function of the third kind
1245 * of order zero of the argument.
1246 *
1247 * The range is partitioned into the two intervals [0,8] and
1248 * (8, infinity). Chebyshev polynomial expansions are employed
1249 * in each interval.
1250 *
1251 *
1252 *
1253 * ACCURACY:
1254 *
1255 * Tested at 2000 random points between 0 and 8. Peak absolute
1256 * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
1257 * Relative error:
1258 * arithmetic domain # trials peak rms
1259 * DEC 0, 30 3100 1.3e-16 2.1e-17
1260 * IEEE 0, 30 30000 1.2e-15 1.6e-16
1261 *
1262 * ERROR MESSAGES:
1263 *
1264 * message condition value returned
1265 * K0 domain x <= 0 MAXNUM
1266 *
1267 */
1268/* k0e()
1269 *
1270 * Modified Bessel function, third kind, order zero,
1271 * exponentially scaled
1272 *
1273 *
1274 *
1275 * SYNOPSIS:
1276 *
1277 * double x, y, k0e();
1278 *
1279 * y = k0e( x );
1280 *
1281 *
1282 *
1283 * DESCRIPTION:
1284 *
1285 * Returns exponentially scaled modified Bessel function
1286 * of the third kind of order zero of the argument.
1287 *
1288 *
1289 *
1290 * ACCURACY:
1291 *
1292 * Relative error:
1293 * arithmetic domain # trials peak rms
1294 * IEEE 0, 30 30000 1.4e-15 1.4e-16
1295 * See k0().
1296 *
1297 */
1298
1299/*
1300Cephes Math Library Release 2.8: June, 2000
1301Copyright 1984, 1987, 2000 by Stephen L. Moshier
1302*/
1303
1304/* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
1305 * in the interval [0,2]. The odd order coefficients are all
1306 * zero; only the even order coefficients are listed.
1307 *
1308 * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EULER.
1309 */
1310
1311static const double k0_A[] =
1312{
1313 1.37446543561352307156e-16,
1314 4.25981614279661018399e-14,
1315 1.03496952576338420167e-11,
1316 1.90451637722020886025e-9,
1317 2.53479107902614945675e-7,
1318 2.28621210311945178607e-5,
1319 1.26461541144692592338e-3,
1320 3.59799365153615016266e-2,
1321 3.44289899924628486886e-1,
1322 -5.35327393233902768720e-1
1323};
1324
1325/* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
1326 * in the inverted interval [2,infinity].
1327 *
1328 * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
1329 */
1330
1331static const double k0_B[] = {
1332 5.30043377268626276149e-18,
1333 -1.64758043015242134646e-17,
1334 5.21039150503902756861e-17,
1335 -1.67823109680541210385e-16,
1336 5.51205597852431940784e-16,
1337 -1.84859337734377901440e-15,
1338 6.34007647740507060557e-15,
1339 -2.22751332699166985548e-14,
1340 8.03289077536357521100e-14,
1341 -2.98009692317273043925e-13,
1342 1.14034058820847496303e-12,
1343 -4.51459788337394416547e-12,
1344 1.85594911495471785253e-11,
1345 -7.95748924447710747776e-11,
1346 3.57739728140030116597e-10,
1347 -1.69753450938905987466e-9,
1348 8.57403401741422608519e-9,
1349 -4.66048989768794782956e-8,
1350 2.76681363944501510342e-7,
1351 -1.83175552271911948767e-6,
1352 1.39498137188764993662e-5,
1353 -1.28495495816278026384e-4,
1354 1.56988388573005337491e-3,
1355 -3.14481013119645005427e-2,
1356 2.44030308206595545468e0
1357};
1358
1359double bessel_k0(double x)
1360{
1361 double y, z;
1362
1363 DEBUG_ENTRY( "bessel_k0()" );
1364
1365 if( x <= 0.0 )
1366 {
1367 fprintf( ioQQQ, "bessel_k0: domain error\n" );
1369 }
1370
1371 if( x <= 2.0 )
1372 {
1373 y = x * x - 2.0;
1374 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
1375 return y;
1376 }
1377 z = 8.0/x - 2.0;
1378 y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x);
1379 return y;
1380}
1381
1382double bessel_k0_scaled(double x)
1383{
1384 double y;
1385
1386 DEBUG_ENTRY( "bessel_k0_scaled()" );
1387
1388 if( x <= 0.0 )
1389 {
1390 fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" );
1392 }
1393
1394 if( x <= 2.0 )
1395 {
1396 y = x * x - 2.0;
1397 y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
1398 return y * exp(x);
1399 }
1400 return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x);
1401}
1402
1403/* k1.c
1404 *
1405 * Modified Bessel function, third kind, order one
1406 *
1407 *
1408 *
1409 * SYNOPSIS:
1410 *
1411 * double x, y, k1();
1412 *
1413 * y = k1( x );
1414 *
1415 *
1416 *
1417 * DESCRIPTION:
1418 *
1419 * Computes the modified Bessel function of the third kind
1420 * of order one of the argument.
1421 *
1422 * The range is partitioned into the two intervals [0,2] and
1423 * (2, infinity). Chebyshev polynomial expansions are employed
1424 * in each interval.
1425 *
1426 *
1427 *
1428 * ACCURACY:
1429 *
1430 * Relative error:
1431 * arithmetic domain # trials peak rms
1432 * DEC 0, 30 3300 8.9e-17 2.2e-17
1433 * IEEE 0, 30 30000 1.2e-15 1.6e-16
1434 *
1435 * ERROR MESSAGES:
1436 *
1437 * message condition value returned
1438 * k1 domain x <= 0 MAXNUM
1439 *
1440 */
1441/* k1e.c
1442 *
1443 * Modified Bessel function, third kind, order one,
1444 * exponentially scaled
1445 *
1446 *
1447 *
1448 * SYNOPSIS:
1449 *
1450 * double x, y, k1e();
1451 *
1452 * y = k1e( x );
1453 *
1454 *
1455 *
1456 * DESCRIPTION:
1457 *
1458 * Returns exponentially scaled modified Bessel function
1459 * of the third kind of order one of the argument:
1460 *
1461 * k1e(x) = exp(x) * k1(x).
1462 *
1463 *
1464 *
1465 * ACCURACY:
1466 *
1467 * Relative error:
1468 * arithmetic domain # trials peak rms
1469 * IEEE 0, 30 30000 7.8e-16 1.2e-16
1470 * See k1().
1471 *
1472 */
1473
1474/*
1475Cephes Math Library Release 2.8: June, 2000
1476Copyright 1984, 1987, 2000 by Stephen L. Moshier
1477*/
1478
1479/* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
1480 * in the interval [0,2].
1481 *
1482 * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
1483 */
1484
1485static const double k1_A[] =
1486{
1487 -7.02386347938628759343e-18,
1488 -2.42744985051936593393e-15,
1489 -6.66690169419932900609e-13,
1490 -1.41148839263352776110e-10,
1491 -2.21338763073472585583e-8,
1492 -2.43340614156596823496e-6,
1493 -1.73028895751305206302e-4,
1494 -6.97572385963986435018e-3,
1495 -1.22611180822657148235e-1,
1496 -3.53155960776544875667e-1,
1497 1.52530022733894777053e0
1498};
1499
1500/* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
1501 * in the interval [2,infinity].
1502 *
1503 * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
1504 */
1505
1506static const double k1_B[] =
1507{
1508 -5.75674448366501715755e-18,
1509 1.79405087314755922667e-17,
1510 -5.68946255844285935196e-17,
1511 1.83809354436663880070e-16,
1512 -6.05704724837331885336e-16,
1513 2.03870316562433424052e-15,
1514 -7.01983709041831346144e-15,
1515 2.47715442448130437068e-14,
1516 -8.97670518232499435011e-14,
1517 3.34841966607842919884e-13,
1518 -1.28917396095102890680e-12,
1519 5.13963967348173025100e-12,
1520 -2.12996783842756842877e-11,
1521 9.21831518760500529508e-11,
1522 -4.19035475934189648750e-10,
1523 2.01504975519703286596e-9,
1524 -1.03457624656780970260e-8,
1525 5.74108412545004946722e-8,
1526 -3.50196060308781257119e-7,
1527 2.40648494783721712015e-6,
1528 -1.93619797416608296024e-5,
1529 1.95215518471351631108e-4,
1530 -2.85781685962277938680e-3,
1531 1.03923736576817238437e-1,
1532 2.72062619048444266945e0
1533};
1534
1535double bessel_k1(double x)
1536{
1537 double y, z;
1538
1539 DEBUG_ENTRY( "bessel_k1()" );
1540
1541 z = 0.5 * x;
1542 if( z <= 0.0 )
1543 {
1544 fprintf( ioQQQ, "bessel_k1: domain error\n" );
1546 }
1547
1548 if( x <= 2.0 )
1549 {
1550 y = x * x - 2.0;
1551 y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
1552 return y;
1553 }
1554 return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
1555}
1556
1557double bessel_k1_scaled(double x)
1558{
1559 double y;
1560
1561 DEBUG_ENTRY( "bessel_k1_scaled()" );
1562
1563 if( x <= 0.0 )
1564 {
1565 fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" );
1567 }
1568
1569 if( x <= 2.0 )
1570 {
1571 y = x * x - 2.0;
1572 y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
1573 return y * exp(x);
1574 }
1575 return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
1576}
1577
1578/* i0.c
1579 *
1580 * Modified Bessel function of order zero
1581 *
1582 *
1583 *
1584 * SYNOPSIS:
1585 *
1586 * double x, y, i0();
1587 *
1588 * y = i0( x );
1589 *
1590 *
1591 *
1592 * DESCRIPTION:
1593 *
1594 * Returns modified Bessel function of order zero of the
1595 * argument.
1596 *
1597 * The function is defined as i0(x) = j0( ix ).
1598 *
1599 * The range is partitioned into the two intervals [0,8] and
1600 * (8, infinity). Chebyshev polynomial expansions are employed
1601 * in each interval.
1602 *
1603 *
1604 *
1605 * ACCURACY:
1606 *
1607 * Relative error:
1608 * arithmetic domain # trials peak rms
1609 * DEC 0,30 6000 8.2e-17 1.9e-17
1610 * IEEE 0,30 30000 5.8e-16 1.4e-16
1611 *
1612 */
1613/* i0e.c
1614 *
1615 * Modified Bessel function of order zero,
1616 * exponentially scaled
1617 *
1618 *
1619 *
1620 * SYNOPSIS:
1621 *
1622 * double x, y, i0e();
1623 *
1624 * y = i0e( x );
1625 *
1626 *
1627 *
1628 * DESCRIPTION:
1629 *
1630 * Returns exponentially scaled modified Bessel function
1631 * of order zero of the argument.
1632 *
1633 * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
1634 *
1635 *
1636 *
1637 * ACCURACY:
1638 *
1639 * Relative error:
1640 * arithmetic domain # trials peak rms
1641 * IEEE 0,30 30000 5.4e-16 1.2e-16
1642 * See i0().
1643 *
1644 */
1645
1646/*
1647Cephes Math Library Release 2.8: June, 2000
1648Copyright 1984, 1987, 2000 by Stephen L. Moshier
1649*/
1650
1651/* Chebyshev coefficients for exp(-x) I0(x)
1652 * in the interval [0,8].
1653 *
1654 * lim(x->0){ exp(-x) I0(x) } = 1.
1655 */
1656
1657static const double i0_A[] =
1658{
1659 -4.41534164647933937950e-18,
1660 3.33079451882223809783e-17,
1661 -2.43127984654795469359e-16,
1662 1.71539128555513303061e-15,
1663 -1.16853328779934516808e-14,
1664 7.67618549860493561688e-14,
1665 -4.85644678311192946090e-13,
1666 2.95505266312963983461e-12,
1667 -1.72682629144155570723e-11,
1668 9.67580903537323691224e-11,
1669 -5.18979560163526290666e-10,
1670 2.65982372468238665035e-9,
1671 -1.30002500998624804212e-8,
1672 6.04699502254191894932e-8,
1673 -2.67079385394061173391e-7,
1674 1.11738753912010371815e-6,
1675 -4.41673835845875056359e-6,
1676 1.64484480707288970893e-5,
1677 -5.75419501008210370398e-5,
1678 1.88502885095841655729e-4,
1679 -5.76375574538582365885e-4,
1680 1.63947561694133579842e-3,
1681 -4.32430999505057594430e-3,
1682 1.05464603945949983183e-2,
1683 -2.37374148058994688156e-2,
1684 4.93052842396707084878e-2,
1685 -9.49010970480476444210e-2,
1686 1.71620901522208775349e-1,
1687 -3.04682672343198398683e-1,
1688 6.76795274409476084995e-1
1689};
1690
1691/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
1692 * in the inverted interval [8,infinity].
1693 *
1694 * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
1695 */
1696
1697static const double i0_B[] =
1698{
1699 -7.23318048787475395456e-18,
1700 -4.83050448594418207126e-18,
1701 4.46562142029675999901e-17,
1702 3.46122286769746109310e-17,
1703 -2.82762398051658348494e-16,
1704 -3.42548561967721913462e-16,
1705 1.77256013305652638360e-15,
1706 3.81168066935262242075e-15,
1707 -9.55484669882830764870e-15,
1708 -4.15056934728722208663e-14,
1709 1.54008621752140982691e-14,
1710 3.85277838274214270114e-13,
1711 7.18012445138366623367e-13,
1712 -1.79417853150680611778e-12,
1713 -1.32158118404477131188e-11,
1714 -3.14991652796324136454e-11,
1715 1.18891471078464383424e-11,
1716 4.94060238822496958910e-10,
1717 3.39623202570838634515e-9,
1718 2.26666899049817806459e-8,
1719 2.04891858946906374183e-7,
1720 2.89137052083475648297e-6,
1721 6.88975834691682398426e-5,
1722 3.36911647825569408990e-3,
1723 8.04490411014108831608e-1
1724};
1725
1726double bessel_i0(double x)
1727{
1728 double y;
1729
1730 DEBUG_ENTRY( "bessel_i0()" );
1731
1732 if( x < 0 )
1733 x = -x;
1734
1735 if( x <= 8.0 )
1736 {
1737 y = (x/2.0) - 2.0;
1738 return exp(x) * chbevl( y, i0_A, 30 );
1739 }
1740 return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
1741}
1742
1743double bessel_i0_scaled(double x)
1744{
1745 double y;
1746
1747 DEBUG_ENTRY( "bessel_i0_scaled()" );
1748
1749 if( x < 0 )
1750 x = -x;
1751
1752 if( x <= 8.0 )
1753 {
1754 y = (x/2.0) - 2.0;
1755 return chbevl( y, i0_A, 30 );
1756 }
1757 return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
1758}
1759
1760/* i1.c
1761 *
1762 * Modified Bessel function of order one
1763 *
1764 *
1765 *
1766 * SYNOPSIS:
1767 *
1768 * double x, y, i1();
1769 *
1770 * y = i1( x );
1771 *
1772 *
1773 *
1774 * DESCRIPTION:
1775 *
1776 * Returns modified Bessel function of order one of the
1777 * argument.
1778 *
1779 * The function is defined as i1(x) = -i j1( ix ).
1780 *
1781 * The range is partitioned into the two intervals [0,8] and
1782 * (8, infinity). Chebyshev polynomial expansions are employed
1783 * in each interval.
1784 *
1785 *
1786 *
1787 * ACCURACY:
1788 *
1789 * Relative error:
1790 * arithmetic domain # trials peak rms
1791 * DEC 0, 30 3400 1.2e-16 2.3e-17
1792 * IEEE 0, 30 30000 1.9e-15 2.1e-16
1793 *
1794 *
1795 */
1796/* i1e.c
1797 *
1798 * Modified Bessel function of order one,
1799 * exponentially scaled
1800 *
1801 *
1802 *
1803 * SYNOPSIS:
1804 *
1805 * double x, y, i1e();
1806 *
1807 * y = i1e( x );
1808 *
1809 *
1810 *
1811 * DESCRIPTION:
1812 *
1813 * Returns exponentially scaled modified Bessel function
1814 * of order one of the argument.
1815 *
1816 * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
1817 *
1818 *
1819 *
1820 * ACCURACY:
1821 *
1822 * Relative error:
1823 * arithmetic domain # trials peak rms
1824 * IEEE 0, 30 30000 2.0e-15 2.0e-16
1825 * See i1().
1826 *
1827 */
1828
1829/*
1830Cephes Math Library Release 2.8: June, 2000
1831Copyright 1985, 1987, 2000 by Stephen L. Moshier
1832*/
1833
1834/* Chebyshev coefficients for exp(-x) I1(x) / x
1835 * in the interval [0,8].
1836 *
1837 * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
1838 */
1839
1840static double i1_A[] =
1841{
1842 2.77791411276104639959e-18,
1843 -2.11142121435816608115e-17,
1844 1.55363195773620046921e-16,
1845 -1.10559694773538630805e-15,
1846 7.60068429473540693410e-15,
1847 -5.04218550472791168711e-14,
1848 3.22379336594557470981e-13,
1849 -1.98397439776494371520e-12,
1850 1.17361862988909016308e-11,
1851 -6.66348972350202774223e-11,
1852 3.62559028155211703701e-10,
1853 -1.88724975172282928790e-9,
1854 9.38153738649577178388e-9,
1855 -4.44505912879632808065e-8,
1856 2.00329475355213526229e-7,
1857 -8.56872026469545474066e-7,
1858 3.47025130813767847674e-6,
1859 -1.32731636560394358279e-5,
1860 4.78156510755005422638e-5,
1861 -1.61760815825896745588e-4,
1862 5.12285956168575772895e-4,
1863 -1.51357245063125314899e-3,
1864 4.15642294431288815669e-3,
1865 -1.05640848946261981558e-2,
1866 2.47264490306265168283e-2,
1867 -5.29459812080949914269e-2,
1868 1.02643658689847095384e-1,
1869 -1.76416518357834055153e-1,
1870 2.52587186443633654823e-1
1871};
1872
1873/* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
1874 * in the inverted interval [8,infinity].
1875 *
1876 * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
1877 */
1878
1879static double i1_B[] =
1880{
1881 7.51729631084210481353e-18,
1882 4.41434832307170791151e-18,
1883 -4.65030536848935832153e-17,
1884 -3.20952592199342395980e-17,
1885 2.96262899764595013876e-16,
1886 3.30820231092092828324e-16,
1887 -1.88035477551078244854e-15,
1888 -3.81440307243700780478e-15,
1889 1.04202769841288027642e-14,
1890 4.27244001671195135429e-14,
1891 -2.10154184277266431302e-14,
1892 -4.08355111109219731823e-13,
1893 -7.19855177624590851209e-13,
1894 2.03562854414708950722e-12,
1895 1.41258074366137813316e-11,
1896 3.25260358301548823856e-11,
1897 -1.89749581235054123450e-11,
1898 -5.58974346219658380687e-10,
1899 -3.83538038596423702205e-9,
1900 -2.63146884688951950684e-8,
1901 -2.51223623787020892529e-7,
1902 -3.88256480887769039346e-6,
1903 -1.10588938762623716291e-4,
1904 -9.76109749136146840777e-3,
1905 7.78576235018280120474e-1
1906};
1907
1908double bessel_i1(double x)
1909{
1910 double y, z;
1911
1912 DEBUG_ENTRY( "bessel_i1()" );
1913
1914 z = fabs(x);
1915 if( z <= 8.0 )
1916 {
1917 y = (z/2.0) - 2.0;
1918 z = chbevl( y, i1_A, 29 ) * z * exp(z);
1919 }
1920 else
1921 {
1922 z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
1923 }
1924 if( x < 0.0 )
1925 z = -z;
1926 return z;
1927}
1928
1929double bessel_i1_scaled(double x)
1930{
1931 double y, z;
1932
1933 DEBUG_ENTRY( "bessel_i1_scaled()" );
1934
1935 z = fabs(x);
1936 if( z <= 8.0 )
1937 {
1938 y = (z/2.0) - 2.0;
1939 z = chbevl( y, i1_A, 29 ) * z;
1940 }
1941 else
1942 {
1943 z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
1944 }
1945 if( x < 0.0 )
1946 z = -z;
1947 return z;
1948}
1949
1950/* ellpk.c
1951 *
1952 * Complete elliptic integral of the first kind
1953 *
1954 *
1955 *
1956 * SYNOPSIS:
1957 *
1958 * double m1, y, ellpk();
1959 *
1960 * y = ellpk( m1 );
1961 *
1962 *
1963 *
1964 * DESCRIPTION:
1965 *
1966 * Approximates the integral
1967 *
1968 *
1969 *
1970 * pi/2
1971 * -
1972 * | |
1973 * | dt
1974 * K(m) = | ------------------
1975 * | 2
1976 * | | sqrt( 1 - m sin t )
1977 * -
1978 * 0
1979 *
1980 * where m = 1 - m1, using the approximation
1981 *
1982 * P(x) - log x Q(x).
1983 *
1984 * The argument m1 is used rather than m so that the logarithmic
1985 * singularity at m = 1 will be shifted to the origin; this
1986 * preserves maximum accuracy.
1987 *
1988 * K(0) = pi/2.
1989 *
1990 * ACCURACY:
1991 *
1992 * Relative error:
1993 * arithmetic domain # trials peak rms
1994 * DEC 0,1 16000 3.5e-17 1.1e-17
1995 * IEEE 0,1 30000 2.5e-16 6.8e-17
1996 *
1997 * ERROR MESSAGES:
1998 *
1999 * message condition value returned
2000 * ellpk domain x<0, x>1 0.0
2001 *
2002 */
2003
2004/*
2005Cephes Math Library, Release 2.8: June, 2000
2006Copyright 1984, 1987, 2000 by Stephen L. Moshier
2007*/
2008
2009static const double elk_P[] =
2010{
2011 1.37982864606273237150e-4,
2012 2.28025724005875567385e-3,
2013 7.97404013220415179367e-3,
2014 9.85821379021226008714e-3,
2015 6.87489687449949877925e-3,
2016 6.18901033637687613229e-3,
2017 8.79078273952743772254e-3,
2018 1.49380448916805252718e-2,
2019 3.08851465246711995998e-2,
2020 9.65735902811690126535e-2,
2021 1.38629436111989062502e0
2022};
2023
2024static const double elk_Q[] =
2025{
2026 2.94078955048598507511e-5,
2027 9.14184723865917226571e-4,
2028 5.94058303753167793257e-3,
2029 1.54850516649762399335e-2,
2030 2.39089602715924892727e-2,
2031 3.01204715227604046988e-2,
2032 3.73774314173823228969e-2,
2033 4.88280347570998239232e-2,
2034 7.03124996963957469739e-2,
2035 1.24999999999870820058e-1,
2036 4.99999999999999999821e-1
2037};
2038
2039static const double C1 = 1.3862943611198906188e0; /* log(4) */
2040
2041double ellpk(double x)
2042{
2043 DEBUG_ENTRY( "ellpk()" );
2044
2045 if( x < 0.0 || x > 1.0 )
2046 {
2047 fprintf( ioQQQ, "ellpk: domain error\n" );
2049 }
2050
2051 if( x > DBL_EPSILON )
2052 {
2053 return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
2054 }
2055 else
2056 {
2057 if( x == 0.0 )
2058 {
2059 fprintf( ioQQQ, "ellpk: domain error\n" );
2061 }
2062 else
2063 {
2064 return C1 - 0.5 * log(x);
2065 }
2066 }
2067}
2068
2069/* expn.c
2070 *
2071 * Exponential integral En
2072 *
2073 *
2074 *
2075 * SYNOPSIS:
2076 *
2077 * int n;
2078 * double x, y, expn();
2079 *
2080 * y = expn( n, x );
2081 *
2082 *
2083 *
2084 * DESCRIPTION:
2085 *
2086 * Evaluates the exponential integral
2087 *
2088 * inf.
2089 * -
2090 * | | -xt
2091 * | e
2092 * E (x) = | ---- dt.
2093 * n | n
2094 * | | t
2095 * -
2096 * 1
2097 *
2098 *
2099 * Both n and x must be nonnegative.
2100 *
2101 * The routine employs either a power series, a continued
2102 * fraction, or an asymptotic formula depending on the
2103 * relative values of n and x.
2104 *
2105 * ACCURACY:
2106 *
2107 * Relative error:
2108 * arithmetic domain # trials peak rms
2109 * DEC 0, 30 5000 2.0e-16 4.6e-17
2110 * IEEE 0, 30 10000 1.7e-15 3.6e-16
2111 *
2112 */
2113
2114/* Cephes Math Library Release 2.8: June, 2000
2115 Copyright 1985, 2000 by Stephen L. Moshier */
2116
2117static const double MAXLOG = log(DBL_MAX);
2118static const double BIG = 1.44115188075855872E+17; /* 2^57 */
2119
2120/*expn exponential intergral for any n */
2121double expn(int n, double x)
2122{
2123 double ans, r, t, yk, xk;
2124 double pk, pkm1, pkm2, qk, qkm1, qkm2;
2125 double psi, z;
2126 int i, k;
2127
2128 DEBUG_ENTRY( "expn()" );
2129
2130 if( n < 0 || x < 0. )
2131 {
2132 fprintf( ioQQQ, "expn: domain error\n" );
2134 }
2135
2136 if( x > MAXLOG )
2137 {
2138 return 0.0;
2139 }
2140
2141 if( x == 0.0 )
2142 {
2143 if( n < 2 )
2144 {
2145 fprintf( ioQQQ, "expn: domain error\n" );
2147 }
2148 else
2149 {
2150 return 1.0/((double)n-1.0);
2151 }
2152 }
2153
2154 if( n == 0 )
2155 {
2156 return exp(-x)/x;
2157 }
2158
2159 /* Expansion for large n */
2160 if( n > 5000 )
2161 {
2162 xk = x + n;
2163 yk = 1.0 / (xk * xk);
2164 t = n;
2165 ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
2166 ans = yk * (ans + t * (t - 2.0 * x));
2167 ans = yk * (ans + t);
2168 ans = (ans + 1.0) * exp( -x ) / xk;
2169 return ans;
2170 }
2171
2172 if( x <= 1.0 )
2173 {
2174 /* Power series expansion */
2175 psi = -EULER - log(x);
2176 for( i=1; i < n; i++ )
2177 psi = psi + 1.0/i;
2178
2179 z = -x;
2180 xk = 0.0;
2181 yk = 1.0;
2182 pk = 1.0 - n;
2183 if( n == 1 )
2184 ans = 0.0;
2185 else
2186 ans = 1.0/pk;
2187 do
2188 {
2189 xk += 1.0;
2190 yk *= z/xk;
2191 pk += 1.0;
2192 if( pk != 0.0 )
2193 {
2194 ans += yk/pk;
2195 }
2196 if( ans != 0.0 )
2197 t = fabs(yk/ans);
2198 else
2199 t = 1.0;
2200 }
2201 while( t > DBL_EPSILON );
2202 ans = powi(z,n-1)*psi/factorial(n-1) - ans;
2203 return ans;
2204 }
2205 else
2206 {
2207 /* continued fraction */
2208 k = 1;
2209 pkm2 = 1.0;
2210 qkm2 = x;
2211 pkm1 = 1.0;
2212 qkm1 = x + n;
2213 ans = pkm1/qkm1;
2214
2215 do
2216 {
2217 k += 1;
2218 if( is_odd(k) )
2219 {
2220 yk = 1.0;
2221 xk = static_cast<double>(n + (k-1)/2);
2222 }
2223 else
2224 {
2225 yk = x;
2226 xk = static_cast<double>(k/2);
2227 }
2228 pk = pkm1 * yk + pkm2 * xk;
2229 qk = qkm1 * yk + qkm2 * xk;
2230 if( qk != 0 )
2231 {
2232 r = pk/qk;
2233 t = fabs( (ans - r)/r );
2234 ans = r;
2235 }
2236 else
2237 t = 1.0;
2238 pkm2 = pkm1;
2239 pkm1 = pk;
2240 qkm2 = qkm1;
2241 qkm1 = qk;
2242 if( fabs(pk) > BIG )
2243 {
2244 pkm2 /= BIG;
2245 pkm1 /= BIG;
2246 qkm2 /= BIG;
2247 qkm1 /= BIG;
2248 }
2249 }
2250 while( t > DBL_EPSILON );
2251
2252 ans *= exp( -x );
2253 return ans;
2254 }
2255}
2256
2257/* erf.c
2258 *
2259 * Error function
2260 *
2261 *
2262 *
2263 * SYNOPSIS:
2264 *
2265 * double x, y, erf();
2266 *
2267 * y = erf( x );
2268 *
2269 *
2270 *
2271 * DESCRIPTION:
2272 *
2273 * The integral is
2274 *
2275 * x
2276 * -
2277 * 2 | | 2
2278 * erf(x) = -------- | exp( - t ) dt.
2279 * sqrt(pi) | |
2280 * -
2281 * 0
2282 *
2283 * The magnitude of x is limited to 9.231948545 for DEC
2284 * arithmetic; 1 or -1 is returned outside this range.
2285 *
2286 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
2287 * erf(x) = 1 - erfc(x).
2288 *
2289 *
2290 *
2291 * ACCURACY:
2292 *
2293 * Relative error:
2294 * arithmetic domain # trials peak rms
2295 * DEC 0,1 14000 4.7e-17 1.5e-17
2296 * IEEE 0,1 30000 3.7e-16 1.0e-16
2297 *
2298 */
2299/* erfc.c
2300 *
2301 * Complementary error function
2302 *
2303 *
2304 *
2305 * SYNOPSIS:
2306 *
2307 * double x, y, erfc();
2308 *
2309 * y = erfc( x );
2310 *
2311 *
2312 *
2313 * DESCRIPTION:
2314 *
2315 *
2316 * 1 - erf(x) =
2317 *
2318 * inf.
2319 * -
2320 * 2 | | 2
2321 * erfc(x) = -------- | exp( - t ) dt
2322 * sqrt(pi) | |
2323 * -
2324 * x
2325 *
2326 *
2327 * For small x, erfc(x) = 1 - erf(x); otherwise rational
2328 * approximations are computed.
2329 *
2330 * A special function expx2.c is used to suppress error amplification
2331 * in computing exp(-x^2).
2332 *
2333 *
2334 * ACCURACY:
2335 *
2336 * Relative error:
2337 * arithmetic domain # trials peak rms
2338 * IEEE 0,26.6417 30000 1.3e-15 2.2e-16
2339 *
2340 *
2341 * ERROR MESSAGES:
2342 *
2343 * message condition value returned
2344 * erfc underflow x > 9.231948545 (DEC) 0.0
2345 *
2346 *
2347 */
2348
2349/*
2350Cephes Math Library Release 2.9: November, 2000
2351Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
2352*/
2353
2354static double erf_P[] = {
2355 2.46196981473530512524e-10,
2356 5.64189564831068821977e-1,
2357 7.46321056442269912687e0,
2358 4.86371970985681366614e1,
2359 1.96520832956077098242e2,
2360 5.26445194995477358631e2,
2361 9.34528527171957607540e2,
2362 1.02755188689515710272e3,
2363 5.57535335369399327526e2
2364};
2365static double erf_Q[] = {
2366/* 1.00000000000000000000e0,*/
2367 1.32281951154744992508e1,
2368 8.67072140885989742329e1,
2369 3.54937778887819891062e2,
2370 9.75708501743205489753e2,
2371 1.82390916687909736289e3,
2372 2.24633760818710981792e3,
2373 1.65666309194161350182e3,
2374 5.57535340817727675546e2
2375};
2376static double erf_R[] = {
2377 5.64189583547755073984e-1,
2378 1.27536670759978104416e0,
2379 5.01905042251180477414e0,
2380 6.16021097993053585195e0,
2381 7.40974269950448939160e0,
2382 2.97886665372100240670e0
2383};
2384static double erf_S[] = {
2385/* 1.00000000000000000000e0,*/
2386 2.26052863220117276590e0,
2387 9.39603524938001434673e0,
2388 1.20489539808096656605e1,
2389 1.70814450747565897222e1,
2390 9.60896809063285878198e0,
2391 3.36907645100081516050e0
2392};
2393
2394
2395#ifndef HAVE_ERFC
2396
2397STATIC double expx2(double x, int sign);
2398
2399/* Define this macro to suppress error propagation in exp(x^2)
2400 by using the expx2 function. The tradeoff is that doing so
2401 generates two calls to the exponential function instead of one. */
2402const bool lgUSE_EXPXSQ = true;
2403
2404double erfc(double a)
2405{
2406 double p,q,x,y,z;
2407
2408 DEBUG_ENTRY( "erfc()" );
2409
2410 x = abs(a);
2411
2412 if( x < 1.0 )
2413 return 1.0 - erf(a);
2414
2415 z = -a * a;
2416
2417 if( z < -MAXLOG )
2418 return ( a < 0.0 ) ? 2.0 : 0.0;
2419
2420 if( lgUSE_EXPXSQ )
2421 /* Compute z = exp(z). */
2422 z = expx2(a, -1);
2423 else
2424 z = exp(z);
2425
2426 if( x < 8.0 )
2427 {
2428 p = polevl( x, erf_P, 8 );
2429 q = p1evl( x, erf_Q, 8 );
2430 }
2431 else
2432 {
2433 p = polevl( x, erf_R, 5 );
2434 q = p1evl( x, erf_S, 6 );
2435 }
2436 y = (z * p)/q;
2437
2438 if( a < 0 )
2439 y = 2.0 - y;
2440
2441 if( y == 0.0 )
2442 return ( a < 0. ) ? 2.0 : 0.0;
2443
2444 return y;
2445}
2446
2447#endif
2448
2449
2450/* Exponentially scaled erfc function
2451 exp(x^2) erfc(x)
2452 valid for x > 1.
2453 Use with ndtr and expx2. */
2454double erfce(double x)
2455{
2456 double p,q;
2457
2458 DEBUG_ENTRY( "erfce()" );
2459
2460 if( x < 8.0 )
2461 {
2462 p = polevl( x, erf_P, 8 );
2463 q = p1evl( x, erf_Q, 8 );
2464 }
2465 else
2466 {
2467 p = polevl( x, erf_R, 5 );
2468 q = p1evl( x, erf_S, 6 );
2469 }
2470 return p/q;
2471}
2472
2473
2474#ifndef HAVE_ERF
2475
2476static double erf_T[] = {
2477 9.60497373987051638749e0,
2478 9.00260197203842689217e1,
2479 2.23200534594684319226e3,
2480 7.00332514112805075473e3,
2481 5.55923013010394962768e4
2482};
2483static double erf_U[] = {
2484/* 1.00000000000000000000e0,*/
2485 3.35617141647503099647e1,
2486 5.21357949780152679795e2,
2487 4.59432382970980127987e3,
2488 2.26290000613890934246e4,
2489 4.92673942608635921086e4
2490};
2491
2492double erf(double x)
2493{
2494 double y, z;
2495
2496 DEBUG_ENTRY( "erf()" );
2497
2498 if( abs(x) > 1.0 )
2499 return 1.0 - erfc(x);
2500 z = x * x;
2501 y = x * polevl( z, erf_T, 4 ) / p1evl( z, erf_U, 5 );
2502 return y;
2503}
2504
2505#endif
2506
2507
2508/* expx2.c
2509 *
2510 * Exponential of squared argument
2511 *
2512 *
2513 *
2514 * SYNOPSIS:
2515 *
2516 * double x, y, expx2();
2517 * int sign;
2518 *
2519 * y = expx2( x, sign );
2520 *
2521 *
2522 *
2523 * DESCRIPTION:
2524 *
2525 * Computes y = exp(x*x) while suppressing error amplification
2526 * that would ordinarily arise from the inexactness of the
2527 * exponential argument x*x.
2528 *
2529 * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
2530 *
2531 *
2532 * ACCURACY:
2533 *
2534 * Relative error:
2535 * arithmetic domain # trials peak rms
2536 * IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17
2537 *
2538 */
2539
2540/*
2541Cephes Math Library Release 2.9: June, 2000
2542Copyright 2000 by Stephen L. Moshier
2543*/
2544
2545
2546#ifndef HAVE_ERFC
2547
2548const double exp_M = 128.0;
2549const double exp_MINV = .0078125;
2550
2551STATIC double expx2(double x, int sign)
2552{
2553 double u, u1, m, f;
2554
2555 DEBUG_ENTRY( "expx2()" );
2556
2557 x = abs(x);
2558 if( sign < 0 )
2559 x = -x;
2560
2561 /* Represent x as an exact multiple of exp_M plus a residual.
2562 exp_M is a power of 2 chosen so that exp(m * m) does not overflow
2563 or underflow and so that |x - m| is small. */
2564 m = exp_MINV * floor(exp_M * x + 0.5);
2565 f = x - m;
2566
2567 /* x^2 = m^2 + 2mf + f^2 */
2568 u = m * m;
2569 u1 = 2 * m * f + f * f;
2570
2571 if( sign < 0 )
2572 {
2573 u = -u;
2574 u1 = -u1;
2575 }
2576
2577 if( (u+u1) > MAXLOG )
2578 return DBL_MAX;
2579
2580 /* u is exact, u1 is small. */
2581 u = exp(u) * exp(u1);
2582 return u;
2583}
2584
2585#endif
2586
2587
2588/* polevl.c
2589 * p1evl.c
2590 *
2591 * Evaluate polynomial
2592 *
2593 *
2594 *
2595 * SYNOPSIS:
2596 *
2597 * int N;
2598 * double x, y, coef[N+1], polevl[];
2599 *
2600 * y = polevl( x, coef, N );
2601 *
2602 *
2603 *
2604 * DESCRIPTION:
2605 *
2606 * Evaluates polynomial of degree N:
2607 *
2608 * 2 N
2609 * y = C + C x + C x +...+ C x
2610 * 0 1 2 N
2611 *
2612 * Coefficients are stored in reverse order:
2613 *
2614 * coef[0] = C , ..., coef[N] = C .
2615 * N 0
2616 *
2617 * The function p1evl() assumes that coef[N] = 1.0 and is
2618 * omitted from the array. Its calling arguments are
2619 * otherwise the same as polevl().
2620 *
2621 *
2622 * SPEED:
2623 *
2624 * In the interest of speed, there are no checks for out
2625 * of bounds arithmetic. This routine is used by most of
2626 * the functions in the library. Depending on available
2627 * equipment features, the user may wish to rewrite the
2628 * program in microcode or assembly language.
2629 *
2630 */
2631
2632/*
2633Cephes Math Library Release 2.1: December, 1988
2634Copyright 1984, 1987, 1988 by Stephen L. Moshier
2635Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2636*/
2637
2638inline double polevl(double x, const double coef[], int N)
2639{
2640 double ans;
2641 int i;
2642 const double *p = coef;
2643
2644 ans = *p++;
2645 i = N;
2646
2647 do
2648 ans = ans * x + *p++;
2649 while( --i );
2650
2651 return ans;
2652}
2653
2654/* p1evl() */
2655/* N
2656 * Evaluate polynomial when coefficient of x is 1.0.
2657 * Otherwise same as polevl.
2658 */
2659
2660inline double p1evl(double x, const double coef[], int N)
2661{
2662 double ans;
2663 const double *p = coef;
2664 int i;
2665
2666 ans = x + *p++;
2667 i = N-1;
2668
2669 do
2670 ans = ans * x + *p++;
2671 while( --i );
2672
2673 return ans;
2674}
2675
2676/* chbevl.c
2677 *
2678 * Evaluate Chebyshev series
2679 *
2680 *
2681 *
2682 * SYNOPSIS:
2683 *
2684 * int N;
2685 * double x, y, coef[N], chebevl();
2686 *
2687 * y = chbevl( x, coef, N );
2688 *
2689 *
2690 *
2691 * DESCRIPTION:
2692 *
2693 * Evaluates the series
2694 *
2695 * N-1
2696 * - '
2697 * y = > coef[i] T (x/2)
2698 * - i
2699 * i=0
2700 *
2701 * of Chebyshev polynomials Ti at argument x/2.
2702 *
2703 * Coefficients are stored in reverse order, i.e. the zero
2704 * order term is last in the array. Note N is the number of
2705 * coefficients, not the order.
2706 *
2707 * If coefficients are for the interval a to b, x must
2708 * have been transformed to x -> 2(2x - b - a)/(b-a) before
2709 * entering the routine. This maps x from (a, b) to (-1, 1),
2710 * over which the Chebyshev polynomials are defined.
2711 *
2712 * If the coefficients are for the inverted interval, in
2713 * which (a, b) is mapped to (1/b, 1/a), the transformation
2714 * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
2715 * this becomes x -> 4a/x - 1.
2716 *
2717 *
2718 *
2719 * SPEED:
2720 *
2721 * Taking advantage of the recurrence properties of the
2722 * Chebyshev polynomials, the routine requires one more
2723 * addition per loop than evaluating a nested polynomial of
2724 * the same degree.
2725 *
2726 */
2727
2728/*
2729Cephes Math Library Release 2.0: April, 1987
2730Copyright 1985, 1987 by Stephen L. Moshier
2731Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2732*/
2733
2734inline double chbevl(double x, const double array[], int n)
2735{
2736 double b0, b1, b2;
2737 const double *p = array;
2738 int i;
2739
2740 b0 = *p++;
2741 b1 = 0.0;
2742 i = n - 1;
2743
2744 do
2745 {
2746 b2 = b1;
2747 b1 = b0;
2748 b0 = x * b1 - b2 + *p++;
2749 }
2750 while( --i );
2751
2752 return 0.5*(b0-b2);
2753}
2754
2755/*******************************************************************
2756 * This marks the end of the block of code from the Cephes library *
2757 *******************************************************************/
2758
2759
2760
2761/* >>refer Mersenne Twister Matsumoto, M. & Nishimura, T. 1998, ACM Transactions on Modeling
2762 * >>refercon and Computer Simulation (TOMACS), 8, 1998 */
2763
2764/********************************************************************
2765 * This copyright notice must accompany the following block of code *
2766 *******************************************************************/
2767
2768/*
2769 A C-program for MT19937, with initialization improved 2002/2/10.
2770 Coded by Takuji Nishimura and Makoto Matsumoto.
2771 This is a faster version by taking Shawn Cokus's optimization,
2772 Matthe Bellew's simplification, Isaku Wada's real version.
2773
2774 Before using, initialize the state by using init_genrand(seed)
2775 or init_by_array(init_key, key_length).
2776
2777 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
2778 All rights reserved.
2779
2780 Redistribution and use in source and binary forms, with or without
2781 modification, are permitted provided that the following conditions
2782 are met:
2783
2784 1. Redistributions of source code must retain the above copyright
2785 notice, this list of conditions and the following disclaimer.
2786
2787 2. Redistributions in binary form must reproduce the above copyright
2788 notice, this list of conditions and the following disclaimer in the
2789 documentation and/or other materials provided with the distribution.
2790
2791 3. The names of its contributors may not be used to endorse or promote
2792 products derived from this software without specific prior written
2793 permission.
2794
2795 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
2796 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
2797 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
2798 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
2799 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
2800 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
2801 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
2802 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
2803 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
2804 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
2805 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
2806
2807
2808 Any feedback is very welcome.
2809 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
2810 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
2811*/
2812
2813/* Period parameters */
2814static const int N = 624;
2815static const int M = 397;
2816static const unsigned long MATRIX_A = 0x9908b0dfUL; /* constant vector a */
2817static const unsigned long UMASK = 0x80000000UL; /* most significant w-r bits */
2818static const unsigned long LMASK = 0x7fffffffUL; /* least significant r bits */
2819inline unsigned long MIXBITS(unsigned long u, unsigned long v)
2820{
2821 return (u & UMASK) | (v & LMASK);
2822}
2823inline unsigned long TWIST(unsigned long u, unsigned long v)
2824{
2825 return (MIXBITS(u,v) >> 1) ^ (v&1UL ? MATRIX_A : 0UL);
2826}
2827
2828static unsigned long state[N]; /* the array for the state vector */
2829static int nleft = 1;
2830static int initf = 0;
2831static unsigned long *nexxt;
2832
2833/* initializes state[N] with a seed */
2834void init_genrand(unsigned long s)
2835{
2836 int j;
2837 state[0]= s & 0xffffffffUL;
2838 for (j=1; j<N; j++) {
2839 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
2840 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
2841 /* In the previous versions, MSBs of the seed affect */
2842 /* only MSBs of the array state[]. */
2843 /* 2002/01/09 modified by Makoto Matsumoto */
2844 state[j] &= 0xffffffffUL; /* for >32 bit machines */
2845 }
2846 nleft = 1; initf = 1;
2847}
2848
2849/* initialize by an array with array-length */
2850/* init_key is the array for initializing keys */
2851/* key_length is its length */
2852/* slight change for C++, 2004/2/26 */
2853void init_by_array(unsigned long init_key[], int key_length)
2854{
2855 int i, j, k;
2856 init_genrand(19650218UL);
2857 i=1; j=0;
2858 k = (N>key_length ? N : key_length);
2859 for (; k; k--) {
2860 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
2861 + init_key[j] + j; /* non linear */
2862 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
2863 i++; j++;
2864 if (i>=N) { state[0] = state[N-1]; i=1; }
2865 if (j>=key_length) j=0;
2866 }
2867 for (k=N-1; k; k--) {
2868 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
2869 - i; /* non linear */
2870 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
2871 i++;
2872 if (i>=N) { state[0] = state[N-1]; i=1; }
2873 }
2874
2875 state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
2876 nleft = 1; initf = 1;
2877}
2878
2879static void next_state()
2880{
2881 unsigned long *p=state;
2882 int j;
2883
2884 /* if init_genrand() has not been called, */
2885 /* a default initial seed is used */
2886 if (initf==0) init_genrand(5489UL);
2887
2888 nleft = N;
2889 nexxt = state;
2890
2891 for (j=N-M+1; --j; p++)
2892 *p = p[M] ^ TWIST(p[0], p[1]);
2893
2894 for (j=M; --j; p++)
2895 *p = p[M-N] ^ TWIST(p[0], p[1]);
2896
2897 *p = p[M-N] ^ TWIST(p[0], state[0]);
2898}
2899
2900/* generates a random number on [0,0xffffffff]-interval */
2901unsigned long genrand_int32()
2902{
2903 unsigned long y;
2904
2905 if (--nleft == 0) next_state();
2906 y = *nexxt++;
2907
2908 /* Tempering */
2909 y ^= (y >> 11);
2910 y ^= (y << 7) & 0x9d2c5680UL;
2911 y ^= (y << 15) & 0xefc60000UL;
2912 y ^= (y >> 18);
2913
2914 return y;
2915}
2916
2917/* generates a random number on [0,0x7fffffff]-interval */
2919{
2920 unsigned long y;
2921
2922 if (--nleft == 0) next_state();
2923 y = *nexxt++;
2924
2925 /* Tempering */
2926 y ^= (y >> 11);
2927 y ^= (y << 7) & 0x9d2c5680UL;
2928 y ^= (y << 15) & 0xefc60000UL;
2929 y ^= (y >> 18);
2930
2931 return (long)(y>>1);
2932}
2933
2934/* generates a random number on [0,1]-real-interval */
2936{
2937 unsigned long y;
2938
2939 if (--nleft == 0) next_state();
2940 y = *nexxt++;
2941
2942 /* Tempering */
2943 y ^= (y >> 11);
2944 y ^= (y << 7) & 0x9d2c5680UL;
2945 y ^= (y << 15) & 0xefc60000UL;
2946 y ^= (y >> 18);
2947
2948 return (double)y * (1.0/4294967295.0);
2949 /* divided by 2^32-1 */
2950}
2951
2952/* generates a random number on [0,1)-real-interval */
2954{
2955 unsigned long y;
2956
2957 if (--nleft == 0) next_state();
2958 y = *nexxt++;
2959
2960 /* Tempering */
2961 y ^= (y >> 11);
2962 y ^= (y << 7) & 0x9d2c5680UL;
2963 y ^= (y << 15) & 0xefc60000UL;
2964 y ^= (y >> 18);
2965
2966 return (double)y * (1.0/4294967296.0);
2967 /* divided by 2^32 */
2968}
2969
2970/* generates a random number on (0,1)-real-interval */
2972{
2973 unsigned long y;
2974
2975 if (--nleft == 0) next_state();
2976 y = *nexxt++;
2977
2978 /* Tempering */
2979 y ^= (y >> 11);
2980 y ^= (y << 7) & 0x9d2c5680UL;
2981 y ^= (y << 15) & 0xefc60000UL;
2982 y ^= (y >> 18);
2983
2984 return ((double)y + 0.5) * (1.0/4294967296.0);
2985 /* divided by 2^32 */
2986}
2987
2988/* generates a random number on [0,1) with 53-bit resolution*/
2990{
2991 unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
2992 return (a*67108864.0+b)*(1.0/9007199254740992.0);
2993}
2994
2995/* These real versions are due to Isaku Wada, 2002/01/09 added */
2996
2997/************************************************************************
2998 * This marks the end of the block of code from Matsumoto and Nishimura *
2999 ************************************************************************/
3000
3001
3002const int N_DAWSON = 100;
3003
3004// tabulated function values of Dawson's function
3005// calculated with xmaxima 5.22.1
3006static const double tbl_dawson[N_DAWSON+1] = {
3007 0.000000000000000000e+0,
3008 9.933599239785286114e-2,
3009 1.947510333680280496e-1,
3010 2.826316650213119286e-1,
3011 3.599434819348881042e-1,
3012 4.244363835020222959e-1,
3013 4.747632036629779306e-1,
3014 5.105040575592317787e-1,
3015 5.321017070563654290e-1,
3016 5.407243187262986750e-1,
3017 5.380794985696822201e-1,
3018 5.262066799705525356e-1,
3019 5.072734964077396141e-1,
3020 4.833975173848241360e-1,
3021 4.565072375268972572e-1,
3022 4.282490710853986254e-1,
3023 3.999398943230814126e-1,
3024 3.725593489740788250e-1,
3025 3.467727691148722451e-1,
3026 3.229743193228178382e-1,
3027 3.013403889237919660e-1,
3028 2.818849389255277938e-1,
3029 2.645107599508319576e-1,
3030 2.490529568377666955e-1,
3031 2.353130556638425762e-1,
3032 2.230837221674354811e-1,
3033 2.121651242424990041e-1,
3034 2.023745109105139857e-1,
3035 1.935507238593667923e-1,
3036 1.855552345354997718e-1,
3037 1.782710306105582873e-1,
3038 1.716003557161446928e-1,
3039 1.654619998786752031e-1,
3040 1.597885804744950500e-1,
3041 1.545240577369634525e-1,
3042 1.496215930807564847e-1,
3043 1.450417730540888593e-1,
3044 1.407511741154154125e-1,
3045 1.367212216746364963e-1,
3046 1.329272910810892667e-1,
3047 1.293480012360051155e-1,
3048 1.259646584343461329e-1,
3049 1.227608160065229225e-1,
3050 1.197219228088390280e-1,
3051 1.168350399532972540e-1,
3052 1.140886102268249801e-1,
3053 1.114722685321253072e-1,
3054 1.089766845891902214e-1,
3055 1.065934312832810744e-1,
3056 1.043148736220704706e-1,
3057 1.021340744242768354e-1,
3058 1.000447137201476355e-1,
3059 9.804101948507806734e-2,
3060 9.611770781195023073e-2,
3061 9.426993099823683440e-2,
3062 9.249323231075475996e-2,
3063 9.078350641561113352e-2,
3064 8.913696463869524546e-2,
3065 8.755010436413927499e-2,
3066 8.601968199264808016e-2,
3067 8.454268897454385223e-2,
3068 8.311633050835148859e-2,
3069 8.173800655824702378e-2,
3070 8.040529489538834118e-2,
3071 7.911593591113373223e-2,
3072 7.786781898606987138e-2,
3073 7.665897022891428948e-2,
3074 7.548754142476270211e-2,
3075 7.435180005364808165e-2,
3076 7.325012025863538754e-2,
3077 7.218097465823629202e-2,
3078 7.114292691123472774e-2,
3079 7.013462495342931107e-2,
3080 6.915479483562112754e-2,
3081 6.820223510065167384e-2,
3082 6.727581164463061598e-2,
3083 6.637445301385693290e-2,
3084 6.549714609447248675e-2,
3085 6.464293215671364666e-2,
3086 6.381090321984490301e-2,
3087 6.300019870755338791e-2,
3088 6.221000236682679049e-2,
3089 6.143953942619040819e-2,
3090 6.068807397169402555e-2,
3091 5.995490652126037542e-2,
3092 5.923937177997213955e-2,
3093 5.854083656061641403e-2,
3094 5.785869785535237086e-2,
3095 5.719238104574369009e-2,
3096 5.654133823962313062e-2,
3097 5.590504672435046070e-2,
3098 5.528300752700259690e-2,
3099 5.467474407290986634e-2,
3100 5.407980093473671650e-2,
3101 5.349774266500934015e-2,
3102 5.292815270562564644e-2,
3103 5.237063236845275277e-2,
3104 5.182479988163068060e-2,
3105 5.129028949666435102e-2,
3106 5.076675065180469937e-2,
3107 5.025384718759852810e-2
3108};
3109
3110//
3111// this routine calculates Dawson's function:
3112//
3113// x
3114// -
3115// | | 2 2
3116// | (y - x )
3117// | e dy
3118// | |
3119// -
3120// 0
3121//
3122// the precomputed values are stored in the array tbl_dawson above.
3123// tbl_dawson[i] is the value of the integral for x = double(i)/10.
3124//
3125// here we do 1st or 3rd order interpolation on this array using
3126// the fact that the grid has equidistant spacing of 0.1.
3127//
3128// the accuracy of 3rd order interpolation is much better, but to
3129// keep the routine fast we sometimes revert to 1st order when
3130// the higher accuracy of 3rd order is not required.
3131//
3132// The Laurent series for this function is given in Mihalas Eq. 9-59.
3133// It is not implemented here.
3134//
3135// dawson has been written by Peter van Hoof (ROB).
3136//
3137inline double dawson(double x, int order)
3138{
3139 double x10 = x*10.;
3140 if( order == 1 )
3141 {
3142 int ind = min(max(int(x10),0),N_DAWSON-1);
3143 double p = x10 - double(ind);
3144 return tbl_dawson[ind] + p*(tbl_dawson[ind+1]-tbl_dawson[ind]);
3145 }
3146 else if( order == 3 )
3147 {
3148 int ind = min(max(int(x10-1.),0),N_DAWSON-3);
3149 double p = x10 - double(ind+1);
3150 double pm1 = p-1.;
3151 double pm2 = p-2.;
3152 double pp1 = p+1.;
3153 // Lagrange 4-point interpolation
3154 return p*pm1*(pp1*tbl_dawson[ind+3]-pm2*tbl_dawson[ind])/6. +
3155 pm2*pp1*(pm1*tbl_dawson[ind+1]-p*tbl_dawson[ind+2])/2.;
3156 }
3157 else
3158 {
3159 TotalInsanity();
3160 }
3161}
3162
3163
3164//
3165// this is a fast version of the Voigt function that is only valid for small a.
3166// The theory is described in:
3167// >>refer rt Voigt Hjerting F., 1938, ApJ 88, 508
3168//
3169// FastVoigtH has been written by Peter van Hoof (ROB).
3170//
3172{
3173 DEBUG_ENTRY( "FastVoigtH()" );
3174
3175 ASSERT( a <= 0.101f );
3176
3177 //
3178 // This routine is guaranteed to give results to better than 0.25% relative accuracy
3179 // over its entire range of validity. The largest discrepancies occur between 1 < v < 10,
3180 // peaking around v = 2 - 7, as shown in the table below:
3181 //
3182 // a = 1e-10 : v = 5.6234e+00 dH/H = 7.6e-06
3183 // a = 1e-5 : v = 7.0795e+00 dH/H = 7.5e-06
3184 // a = 1e-4 : v = 3.7584e+00 dH/H = 1.3e-05
3185 // a = 3e-4 : v = 3.5481e+00 dH/H = 1.6e-05
3186 // a = 1e-3 : v = 3.3497e+00 dH/H = 1.9e-05
3187 // a = 3e-3 : v = 3.1623e+00 dH/H = 2.2e-05
3188 // a = 0.01 : v = 3.1623e+00 dH/H = 1.0e-05
3189 // a = 0.03 : v = 2.8184e+00 dH/H = 1.8e-04
3190 // a = 0.1 : v = 2.6607e+00 dH/H = 2.4e-03
3191 //
3192 // to get a guaranteed relative accuracy <= 1.e-4, use a < 0.0235
3193 // to get a guaranteed relative accuracy <= 1.e-3, use a < 0.066
3194 //
3195 // for a > 0.1 the series expansions used in this routine quickly start breaking down
3196 // and the algorithm becomes useless, so never use this routine for a > 0.1 !
3197 //
3198
3199 v = abs(v); // the profile is symmetrical in v
3200
3201 if( v > 9.f )
3202 {
3203 // use Laurent series; this is Eq. 7 of Hjerting with one higher order term added
3204 //
3205 // The complete series is:
3206 //
3207 // inf
3208 // ----
3209 // a \ (2 n + 1)!!
3210 // H(a,v) = ----------- | -----------
3211 // 2 / n 2n
3212 // sqrt(pi) v ---- 2 v
3213 // n=0
3214 //
3215 realnum vm2 = 1.f/pow2(v);
3216 return a*vm2/realnum(SQRTPI)*(((105.f/8.f*vm2 + 15.f/4.f)*vm2 + 3.f/2.f)*vm2 + 1.f);
3217 }
3218 else
3219 {
3220 realnum v2 = pow2(v);
3221 // NOTE: it is important to use dsexp here so that we avoid expf(). The
3222 // latter can be significantly slower, at least on non-Suse Linux platforms!
3223 // see: https://bugzilla.redhat.com/show_bug.cgi?id=521190
3224 // revert this statement to: emv2 = exp(-v2) once this is solved.
3225 realnum emv2 = realnum(dsexp(v2));
3226 int order = ( a > 0.003f || v > 1.5f ) ? 3 : 1;
3227 // this is Eq. 3 of Hjerting with an additional term:
3228 //
3229 // the last term in Eq. 4 of Hjerting can be expanded to lowest order in a as:
3230 //
3231 // inf
3232 // -
3233 // | | 2 2 2
3234 // 1 | (a x) - x 2 2 - v
3235 // -------- | ------ exp(----) cos(v x) dx = - a (2 v - 1) e
3236 // sqrt(pi) | | 2 4
3237 // -
3238 // 0
3239 //
3240 return emv2*(1.f-pow2(a)*(2.f*v2-1.f)) +
3241 2.f*a/realnum(SQRTPI)*(-1.f + 2.f*v*realnum(dawson(v,order)));
3242 }
3243}
3244
3245/*
3246 Calculate the Voigt profile aka Faddeeva function with relative
3247 error less than 10^(-R).
3248
3249 R0=1.51*EXP(1.144*R) and R1=1.60*EXP(0.554*R) can be set by the the
3250 user subject to the constraints 14.88<R0<460.4 and 4.85<R1<25.5
3251
3252 Translated from a Fortran version by R.J. Wells,
3253 see http://www.atm.ox.ac.uk/user/wells/voigt.html
3254
3255 >>refer rt Voigt Wells, R.J. 1999, JQSRT, 62, 29
3256*/
3257void humlik(int n, const realnum x[], realnum y, realnum k[])
3258{
3259 DEBUG_ENTRY( "humlik()" );
3260
3261 /* n IN Number of points
3262 x IN Input x array
3263 y IN Input y value >=0.0
3264 k OUT (Voigt) array */
3265
3266 // use doubles internally to avoid overflow for very large y values (above roughly 5000)
3267 // these very high values can occur in X-ray lines; the end result is converted back to realnum
3268
3269 /* Initialized data */
3270 static const double c[6] = { 1.0117281, -.75197147, .012557727, .010022008, -2.4206814e-4, 5.0084806e-7 };
3271 static const double s[6] = { 1.393237, .23115241, -.15535147, .0062183662, 9.1908299e-5, -6.2752596e-7 };
3272 static const double t[6] = { .31424038, .94778839, 1.5976826, 2.2795071, 3.020637, 3.8897249 };
3273
3274 static const double R0 = 146.7, R1 = 14.67;
3275
3276 /* Local variables */
3277 double d, ki;
3278 int i, j;
3279 double a0, d0, e0, d2, e2, h0, e4, h2, h4, h6, p0,
3280 p2, p4, p6, p8, z0, z2, z4, z6, z8, mf[6], pf[6],
3281 mq[6], yf, pq[6], xm[6], ym[6], xp[6], xq, yq, yp[6];
3282 bool rg1, rg2, rg3;
3283 double abx, ypy0, xlim0, xlim1, xlim2, xlim3, xlim4, ypy0q, yrrtpi;
3284
3285 rg1 = rg2 = rg3 = true;
3286 // Initialization to quiet warnings
3287 z0 = z2 = z4 = z6 = z8 = 0.;
3288 p0 = p2 = p4 = p6 = p8 = 0.;
3289 h0 = h2 = h4 = h6 = 0.;
3290 a0 = d0 = d2 = e0 = e2 = e4 = 0.;
3291
3292 yq = y * y;
3293 yrrtpi = y * .56418958; // 1/SQRT(pi)
3294 /* Region boundaries */
3295 xlim0 = R0 - y;
3296 xlim1 = R1 - y;
3297 xlim3 = y * 3.097 - .45;
3298 xlim2 = 6.8 - y;
3299 xlim4 = y * 18.1 + 1.65;
3300 if (y <= 1e-6)
3301 { /* avoid W4 algorithm */
3302 xlim1 = xlim0;
3303 xlim2 = xlim0;
3304 }
3305
3306 for (i = 0; i < n; ++i)
3307 {
3308 abx = fabs(x[i]);
3309 xq = abx * abx;
3310 if (abx > xlim0)
3311 { /* Region 0 algorithm */
3312 k[i] = realnum(yrrtpi / (xq + yq));
3313 }
3314 else if (abx > xlim1)
3315 { /* Humlicek W4 Region 1 */
3316 if (rg1)
3317 {
3318 /* First point in Region 1 */
3319 rg1 = false;
3320 a0 = yq + .5;
3321 d0 = a0 * a0;
3322 d2 = yq + yq - 1.;
3323 }
3324 d = .56418958 / (d0 + xq * (d2 + xq));
3325 k[i] = realnum(d * y * (a0 + xq));
3326 }
3327 else if (abx > xlim2)
3328 { /* Humlicek W4 Region 2 */
3329 if (rg2)
3330 { /* First point in Region 2 */
3331 rg2 = false;
3332 h0 = yq * (yq * (yq * (yq + 6.) + 10.5) + 4.5) + .5625;
3333 h2 = yq * (yq * (yq * 4. + 6.) + 9.) - 4.5;
3334 h4 = 10.5 - yq * (6. - yq * 6.);
3335 h6 = yq * 4. - 6.;
3336 e0 = yq * (yq * (yq + 5.5) + 8.25) + 1.875;
3337 e2 = yq * (yq * 3. + 1.) + 5.25;
3338 e4 = h6 * .75;
3339 }
3340 d = .56418958 / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
3341 k[i] = realnum(d * y * (e0 + xq * (e2 + xq * (e4 + xq))));
3342 }
3343 else if (abx < xlim3)
3344 { /* Humlicek W4 Region 3 */
3345 if (rg3)
3346 {
3347 /* First point in Region 3 */
3348 rg3 = false;
3349 z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 13.3988) +
3350 88.26741) + 369.1989) + 1074.409) + 2256.981) + 3447.629) +
3351 3764.966) + 2802.87) + 1280.829) + 272.1014;
3352 z2 = y * (y * (y * (y * (y * (y * (y * (y * 5. + 53.59518) +
3353 266.2987) + 793.4273) + 1549.675) + 2037.31) + 1758.336) +
3354 902.3066) + 211.678;
3355 z4 = y * (y * (y * (y * (y * (y * 10. + 80.39278) + 269.2916) +
3356 479.2576) + 497.3014) + 308.1852) + 78.86585;
3357 z6 = y * (y * (y * (y * 10. + 53.59518) + 92.75679) + 55.02933) +
3358 22.03523;
3359 z8 = y * (y * 5. + 13.3988) + 1.49646;
3360 p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * .3183291 +
3361 4.264678) + 27.93941) + 115.3772) + 328.2151) + 662.8097) +
3362 946.897) + 919.4955) + 549.3954) + 153.5168;
3363 p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163 + 12.79458) +
3364 56.81652) + 139.4665) + 189.773) + 124.5975) - 1.322256) -
3365 34.16955;
3366 p4 = y * (y * (y * (y * (y * 1.9099744 + 12.79568) + 29.81482) +
3367 24.01655) + 10.46332) + 2.584042;
3368 p6 = y * (y * (y * 1.273316 + 4.266322) + .9377051) - .07272979;
3369 p8 = y * .3183291 + 5.480304e-4;
3370 }
3371 d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
3372 k[i] = realnum(d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8)))));
3373 }
3374 else
3375 { /* Humlicek CPF12 algorithm */
3376 ypy0 = y + 1.5;
3377 ypy0q = ypy0 * ypy0;
3378 ki = 0.;
3379 for (j = 0; j <= 5; ++j)
3380 {
3381 d = x[i] - t[j];
3382 mq[j] = d * d;
3383 mf[j] = 1. / (mq[j] + ypy0q);
3384 xm[j] = mf[j] * d;
3385 ym[j] = mf[j] * ypy0;
3386 d = x[i] + t[j];
3387 pq[j] = d * d;
3388 pf[j] = 1. / (pq[j] + ypy0q);
3389 xp[j] = pf[j] * d;
3390 yp[j] = pf[j] * ypy0;
3391 }
3392 if (abx <= xlim4)
3393 { /* Humlicek CPF12 Region I */
3394 for (j = 0; j <= 5; ++j)
3395 {
3396 ki = ki + c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
3397 }
3398 }
3399 else
3400 { /* Humlicek CPF12 Region II */
3401 yf = y + 3.;
3402 for (j = 0; j <= 5; ++j)
3403 {
3404 ki = ki + (c[j] * (mq[j] * mf[j] - ym[j] * 1.5) + s[j] * yf * xm[j]) / (mq[j] + 2.25) +
3405 (c[j] * (pq[j] * pf[j] - yp[j] * 1.5) - s[j] * yf * xp[j]) / (pq[j] + 2.25);
3406 }
3407 ki = y * ki + exp(-xq);
3408 }
3409 k[i] = realnum(ki);
3410 }
3411 }
3412}
3413
3414/******************************************************************
3415 * This marks the end of the block of code for the Voigt function *
3416 ******************************************************************/
3417
3418STATIC uint32 MD5swap( uint32 word );
3419STATIC void MD5_Transform (uint32 *digest, const uint32 *in);
3420
3421//
3422// The routines MD5file(), MD5textfile(), MD5string() and MD5swap() were written by Peter van Hoof
3423//
3424// this version returns the md5sum of a file and is identical to the well known md5sum algorithm
3425string MD5file(const char* fnam, access_scheme scheme)
3426{
3427 DEBUG_ENTRY( "MD5file()" );
3428
3429 fstream ioFile;
3430 open_data( ioFile, fnam, mode_r, scheme );
3431
3432 char c;
3433 string content;
3434 while( ioFile.get(c) )
3435 content += c;
3436
3437 return MD5string( content );
3438}
3439
3440
3441// this version returns the md5sum of a text file. It filters out the eol characters,
3442// which makes it incompatible with the md5sum algorithm, but it is OS agnostic...
3443// comment lines that start with the hash symbol are also skipped
3444string MD5datafile(const char* fnam, access_scheme scheme)
3445{
3446 DEBUG_ENTRY( "MD5datafile()" );
3447
3448 fstream ioFile;
3449 open_data( ioFile, fnam, mode_r, scheme );
3450
3451 string line, content;
3452 while( getline( ioFile, line ) )
3453 if( line[0] != '#' )
3454 content += line;
3455
3456 return MD5string( content );
3457}
3458
3459
3460// calculate the md5sum of an arbitrary string
3461string MD5string(const string& str)
3462{
3463 DEBUG_ENTRY( "MD5string()" );
3464
3465 uint32 state[4];
3466
3467 state[0] = 0x67452301L;
3468 state[1] = 0xefcdab89L;
3469 state[2] = 0x98badcfeL;
3470 state[3] = 0x10325476L;
3471
3472 string lstr = str;
3473
3474 // pad the string following RFC 1321 3.1 Step 1.
3475 uint32 bytes = str.length()%64;
3476 uint32 padlen = ( bytes >= 56 ) ? 64 + 56 - bytes : 56 - bytes;
3477 lstr += '\x80';
3478 for( uint32 i=1; i < padlen; ++i )
3479 lstr += '\0';
3480
3481 ASSERT( lstr.length()%64 == 56 );
3482
3483 uint32 i;
3484 union {
3485 uint32 in[16];
3486 unsigned char chr[64];
3487 } u;
3488
3489 for( i=0; i < lstr.length()/64; ++i )
3490 {
3491 for( uint32 j=0; j < 64; ++j )
3492 {
3493 if( cpu.i().little_endian() )
3494 u.chr[j] = lstr[i*64+j];
3495 else
3496 {
3497 uint32 jr = j%4;
3498 uint32 j4 = j-jr;
3499 u.chr[j4+3-jr] = lstr[i*64+j];
3500 }
3501 }
3502 MD5_Transform( state, u.in );
3503 }
3504 for( uint32 j=0; j < 56; ++j )
3505 {
3506 if( cpu.i().little_endian() )
3507 u.chr[j] = lstr[i*64+j];
3508 else
3509 {
3510 uint32 jr = j%4;
3511 uint32 j4 = j-jr;
3512 u.chr[j4+3-jr] = lstr[i*64+j];
3513 }
3514 }
3515 // append the length of the string in _bits_ following RFC 1321 3.1 Step 2.
3516 u.in[14] = (str.length()<<3)&0xffffffff;
3517 u.in[15] = (str.length()>>29)&0xffffffff;
3518
3519 MD5_Transform( state, u.in );
3520
3521 ostringstream hash;
3522 for( uint32 i=0; i < 4; ++i )
3523 hash << hex << setfill('0') << setw(8) << MD5swap(state[i]);
3524
3525 return hash.str();
3526}
3527
3528STATIC uint32 MD5swap( uint32 word )
3529{
3530 DEBUG_ENTRY( "MD5swap()" );
3531
3532 union {
3533 uint32 word;
3534 unsigned char byte[4];
3535 } ui, uo;
3536
3537 ui.word = word;
3538 uo.byte[0] = ui.byte[3];
3539 uo.byte[1] = ui.byte[2];
3540 uo.byte[2] = ui.byte[1];
3541 uo.byte[3] = ui.byte[0];
3542
3543 return uo.word;
3544}
3545
3546//
3547// The following implementation of the MD5 algorithm was taken from the
3548// Crypto++ library (http://www.cryptopp.com/) and is subject to the
3549// following license:
3550//
3551//
3552// Compilation Copyright (c) 1995-2010 by Wei Dai. All rights reserved.
3553// This copyright applies only to this software distribution package
3554// as a compilation, and does not imply a copyright on any particular
3555// file in the package.
3556//
3557// All individual files in this compilation are placed in the public domain by
3558// Wei Dai and other contributors.
3559//
3560// I would like to thank the following authors for placing their works into
3561// the public domain:
3562//
3563// Joan Daemen - 3way.cpp
3564// Leonard Janke - cast.cpp, seal.cpp
3565// Steve Reid - cast.cpp
3566// Phil Karn - des.cpp
3567// Andrew M. Kuchling - md2.cpp, md4.cpp
3568// Colin Plumb - md5.cpp
3569// Seal Woods - rc6.cpp
3570// Chris Morgan - rijndael.cpp
3571// Paulo Baretto - rijndael.cpp, skipjack.cpp, square.cpp
3572// Richard De Moliner - safer.cpp
3573// Matthew Skala - twofish.cpp
3574// Kevin Springle - camellia.cpp, shacal2.cpp, ttmac.cpp, whrlpool.cpp, ripemd.cpp
3575//
3576// Permission to use, copy, modify, and distribute this compilation for
3577// any purpose, including commercial applications, is hereby granted
3578// without fee, subject to the following restrictions:
3579//
3580// 1. Any copy or modification of this compilation in any form, except
3581// in object code form as part of an application software, must include
3582// the above copyright notice and this license.
3583//
3584// 2. Users of this software agree that any modification or extension
3585// they provide to Wei Dai will be considered public domain and not
3586// copyrighted unless it includes an explicit copyright notice.
3587//
3588// 3. Wei Dai makes no warranty or representation that the operation of the
3589// software in this compilation will be error-free, and Wei Dai is under no
3590// obligation to provide any services, by way of maintenance, update, or
3591// otherwise. THE SOFTWARE AND ANY DOCUMENTATION ARE PROVIDED "AS IS"
3592// WITHOUT EXPRESS OR IMPLIED WARRANTY INCLUDING, BUT NOT LIMITED TO,
3593// THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
3594// PURPOSE. IN NO EVENT WILL WEI DAI OR ANY OTHER CONTRIBUTOR BE LIABLE FOR
3595// DIRECT, INCIDENTAL OR CONSEQUENTIAL DAMAGES, EVEN IF
3596// ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
3597//
3598// 4. Users will not use Wei Dai or any other contributor's name in any
3599// publicity or advertising, without prior written consent in each case.
3600//
3601// 5. Export of this software from the United States may require a
3602// specific license from the United States Government. It is the
3603// responsibility of any person or organization contemplating export
3604// to obtain such a license before exporting.
3605//
3606// 6. Certain parts of this software may be protected by patents. It
3607// is the users' responsibility to obtain the appropriate
3608// licenses before using those parts.
3609//
3610// If this compilation is used in object code form in an application
3611// software, acknowledgement of the author is not required but would be
3612// appreciated. The contribution of any useful modifications or extensions
3613// to Wei Dai is not required but would also be appreciated.
3614//
3615
3616// md5.cpp - modified by Wei Dai from Colin Plumb's public domain md5.c
3617// any modifications are placed in the public domain
3618
3619inline uint32 rotlFixed(uint32 x, unsigned int y)
3620{
3621 return uint32((x<<y) | (x>>(32-y)));
3622}
3623
3624STATIC void MD5_Transform (uint32 *digest, const uint32 *in)
3625{
3626 DEBUG_ENTRY( "MD5_Transform()" );
3627
3628#define F1(x, y, z) (z ^ (x & (y ^ z)))
3629#define F2(x, y, z) F1(z, x, y)
3630#define F3(x, y, z) (x ^ y ^ z)
3631#define F4(x, y, z) (y ^ (x | ~z))
3632
3633#define MD5STEP(f, w, x, y, z, data, s) \
3634 w = rotlFixed(w + f(x, y, z) + data, s) + x
3635
3636 uint32 a, b, c, d;
3637
3638 a=digest[0];
3639 b=digest[1];
3640 c=digest[2];
3641 d=digest[3];
3642
3643 MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
3644 MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
3645 MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
3646 MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
3647 MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
3648 MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
3649 MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
3650 MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
3651 MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
3652 MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
3653 MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
3654 MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
3655 MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
3656 MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
3657 MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
3658 MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
3659
3660 MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
3661 MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
3662 MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
3663 MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
3664 MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
3665 MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
3666 MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
3667 MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
3668 MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
3669 MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
3670 MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
3671 MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
3672 MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
3673 MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
3674 MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
3675 MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
3676
3677 MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
3678 MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
3679 MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
3680 MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
3681 MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
3682 MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
3683 MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
3684 MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
3685 MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
3686 MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
3687 MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
3688 MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
3689 MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
3690 MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
3691 MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
3692 MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
3693
3694 MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
3695 MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
3696 MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
3697 MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
3698 MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
3699 MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
3700 MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
3701 MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
3702 MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
3703 MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
3704 MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
3705 MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
3706 MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
3707 MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
3708 MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
3709 MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
3710
3711 digest[0] += a;
3712 digest[1] += b;
3713 digest[2] += c;
3714 digest[3] += d;
3715}
3716
3717/***************************************************************
3718 * This marks the end of the block of code from Wei Dai et al. *
3719 ***************************************************************/
static realnum b1[83]
static double b2[63]
static realnum b0[83]
static double a0[83]
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
bool is_odd(int j)
Definition cddefines.h:714
#define STATIC
Definition cddefines.h:97
T pow2(T a)
Definition cddefines.h:931
T sign(T x, T y)
Definition cddefines.h:800
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
long max(int a, long b)
Definition cddefines.h:775
double e2(double t)
Definition service.cpp:299
NORETURN void TotalInsanity(void)
Definition service.cpp:886
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
double powi(double, long int)
Definition service.cpp:604
static t_lfact & Inst()
Definition cddefines.h:175
vector< double > p_lf
double get_lfact(unsigned long n)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
const ios_base::openmode mode_r
Definition cpu.h:212
access_scheme
Definition cpu.h:207
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
STATIC void order(long, realnum[], long *, long *, long *)
UNUSED const double PI
Definition physconst.h:29
UNUSED const double EULER
Definition physconst.h:26
UNUSED const double SQRTPI
Definition physconst.h:44
t_state state
Definition state.cpp:19
double genrand_real3()
STATIC uint32 MD5swap(uint32 word)
#define MD5STEP(f, w, x, y, z, data, s)
static const double tbl_dawson[N_DAWSON+1]
STATIC void MD5_Transform(uint32 *digest, const uint32 *in)
static double i1_B[]
static const double k1_A[]
double bessel_j1(double x)
const int N_DAWSON
double bessel_k1_scaled(double x)
static const double BIG
string MD5file(const char *fnam, access_scheme scheme)
static const double PIO4
double factorial(long n)
static const double MAXLOG
static double erf_Q[]
static const double b0_PQ[7]
static const double b1_YQ[8]
double bessel_k0(double x)
static const double b1_YP[6]
static double b0_RQ[8]
double ellpk(double x)
static const double b0_QP[8]
static const double k0_A[]
double bessel_i1_scaled(double x)
double bessel_yn(int n, double x)
#define F1(x, y, z)
static const double DR2
double genrand_real1()
static const double b0_QQ[7]
static const double k0_B[]
static const double b0_YQ[7]
static const unsigned long UMASK
double bessel_y0(double x)
double bessel_k0_scaled(double x)
static const unsigned long LMASK
void humlik(int n, const realnum x[], realnum y, realnum k[])
double chbevl(double, const double[], int)
void init_genrand(unsigned long s)
static int nleft
double bessel_jn(int n, double x)
#define F4(x, y, z)
string MD5string(const string &str)
static const double Z2
static const double k1_B[]
unsigned long TWIST(unsigned long u, unsigned long v)
double lfactorial(long n)
static const double i0_A[]
string MD5datafile(const char *fnam, access_scheme scheme)
double genrand_res53()
static const double elk_Q[]
complex< double > cdgamma(complex< double > x)
static double erf_P[]
#define F3(x, y, z)
static const double C1
long genrand_int31()
double bessel_i1(double x)
unsigned long genrand_int32()
static const int M
static const double b1_QP[8]
static const double DR1
double bessel_i0_scaled(double x)
realnum FastVoigtH(realnum a, realnum v)
static double i1_A[]
static const double b0_PP[7]
double erfce(double x)
double bessel_j0(double x)
static double erf_R[]
double dawson(double x, int order)
static const unsigned long MATRIX_A
static const int N
double bessel_i0(double x)
static const double b1_RQ[8]
double genrand_real2()
static const double b1_QQ[7]
bool linfit(long n, const double xorg[], const double yorg[], double &a, double &siga, double &b, double &sigb)
void init_by_array(unsigned long init_key[], int key_length)
static const double b1_RP[4]
double polevl(double x, const double coef[], int N)
double p1evl(double x, const double coef[], int N)
static const double pre_factorial[NPRE_FACTORIAL]
double bessel_k1(double x)
uint32 rotlFixed(uint32 x, unsigned int y)
static const double elk_P[]
static double erf_S[]
static double b0_RP[4]
static const double b1_PQ[7]
static const double Z1
static int initf
#define F2(x, y, z)
double expn(int n, double x)
static const double b1_PP[7]
static const double THPIO4
static const double TWOOPI
static const double i0_B[]
unsigned long MIXBITS(unsigned long u, unsigned long v)
static void next_state()
double bessel_y1(double x)
static unsigned long * nexxt
static const double b0_YP[8]
static const double SQ2OPI
double erfc(double)
static const int NPRE_FACTORIAL
Definition thirdparty.h:24
double erf(double)