cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
7 inline double polevl(double x, const double coef[], int N);
8 inline double p1evl(double x, const double coef[], int N);
9 inline double chbevl(double, const double[], int);
10 
11 
12 /* the routine linfit was derived from the program slopes.f */
13 
14 /********************************************************************/
15 /* */
16 /* PROGRAM SLOPES */
17 /* */
18 /* PROGRAM TO COMPUTE THE THEORETICAL REGRESSION SLOPES */
19 /* AND UNCERTAINTIES (VIA DELTA METHOD), AND UNCERTAINTIES */
20 /* VIA BOOTSTRAP AND BIVARIATE NORMAL SIMULATION FOR A */
21 /* (X(I),Y(I)) DATA SET WITH UNKNOWN POPULATION DISTRIBUTION. */
22 /* */
23 /* WRITTEN BY ERIC FEIGELSON, PENN STATE. JUNE 1991 */
24 /* */
25 /* */
26 /* A full description of these methods can be found in: */
27 /* Isobe, T., Feigelson, E. D., Akritas, M. and Babu, G. J., */
28 /* Linear Regression in Astronomy I, Astrophys. J. 364, 104 */
29 /* (1990) */
30 /* Babu, G. J. and Feigelson, E. D., Analytical and Monte Carlo */
31 /* Comparisons of Six Different Linear Least Squares Fits, */
32 /* Communications in Statistics, Simulation & Computation, */
33 /* 21, 533 (1992) */
34 /* Feigelson, E. D. and Babu, G. J., Linear Regression in */
35 /* Astronomy II, Astrophys. J. 397, 55 (1992). */
36 /* */
37 /********************************************************************/
38 
39 /* this used to be sixlin, but only the first fit has been retained */
40 
41 /********************************************************************/
42 /************************* routine linfit ***************************/
43 /********************************************************************/
44 
45 bool linfit(
46  long n,
47  double x[], /* x[n] */
48  double y[], /* y[n] */
49  double &a,
50  double &siga,
51  double &b,
52  double &sigb
53 )
54 {
55 
56 /*
57  * linear regression
58  * written by t. isobe, g. j. babu and e. d. feigelson
59  * center for space research, m.i.t.
60  * and
61  * the pennsylvania state university
62  *
63  * rev. 1.0, september 1990
64  *
65  * this subroutine provides linear regression coefficients
66  * computed by six different methods described in isobe,
67  * feigelson, akritas, and babu 1990, astrophysical journal
68  * and babu and feigelson 1990, subm. to technometrics.
69  * the methods are ols(y/x), ols(x/y), ols bisector, orthogonal,
70  * reduced major axis, and mean-ols regressions.
71  *
72  * [Modified and translated to C/C++ by Peter van Hoof, Royal
73  * Observatory of Belgium; only the first method has been retained]
74  *
75  *
76  * input
77  * x[i] : independent variable
78  * y[i] : dependent variable
79  * n : number of data points
80  *
81  * output
82  * a : intercept coefficients
83  * b : slope coefficients
84  * siga : standard deviations of intercepts
85  * sigb : standard deviations of slopes
86  *
87  * error returns
88  * returns true when division by zero will
89  * occur (i.e. when sxy is zero)
90  */
91 
92  DEBUG_ENTRY( "linfit()" );
93 
94  /* initializations */
95  a = 0.0;
96  siga = 0.0;
97  b = 0.0;
98  sigb = 0.0;
99 
100  /* compute averages and sums */
101  double s1 = 0.0;
102  double s2 = 0.0;
103  for( long i=0; i < n; i++ )
104  {
105  s1 += x[i];
106  s2 += y[i];
107  }
108  double rn = (double)n;
109  double xavg = s1/rn;
110  double yavg = s2/rn;
111  double sxx = 0.0;
112  double sxy = 0.0;
113  for( long i=0; i < n; i++ )
114  {
115  x[i] -= xavg;
116  y[i] -= yavg;
117  sxx += pow2(x[i]);
118  sxy += x[i]*y[i];
119  }
120 
121  if( sxy == 0.0 )
122  {
123  return true;
124  }
125 
126  /* compute the slope coefficient */
127  b = sxy/sxx;
128 
129  /* compute intercept coefficient */
130  a = yavg - b*xavg;
131 
132  /* prepare for computation of variances */
133  double sum1 = 0.0;
134  for( long i=0; i < n; i++ )
135  sum1 += pow2(x[i]*(y[i] - b*x[i]));
136 
137  /* compute variance of the slope coefficient */
138  sigb = sum1/pow2(sxx);
139 
140  /* compute variance of the intercept coefficient */
141  for( long i=0; i < n; i++ )
142  siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx));
143 
144  /* convert variances to standard deviations */
145  sigb = sqrt(sigb);
146  siga = sqrt(siga)/rn;
147 
148  /* return data arrays to their original form */
149  for( long i=0; i < n; i++ )
150  {
151  x[i] += xavg;
152  y[i] += yavg;
153  }
154  return false;
155 }
156 
157 
158 /* the routines factorial and lfactorial came originally from hydro_bauman.cpp
159  * and were written by Robert Paul Bauman. lfactorial was modified by Peter van Hoof */
160 
161 /************************************************************************************************/
162 /* pre-calculated factorials */
163 /************************************************************************************************/
164 
165 static const double pre_factorial[NPRE_FACTORIAL] =
166 {
167  1.00000000000000000000e+00,
168  1.00000000000000000000e+00,
169  2.00000000000000000000e+00,
170  6.00000000000000000000e+00,
171  2.40000000000000000000e+01,
172  1.20000000000000000000e+02,
173  7.20000000000000000000e+02,
174  5.04000000000000000000e+03,
175  4.03200000000000000000e+04,
176  3.62880000000000000000e+05,
177  3.62880000000000000000e+06,
178  3.99168000000000000000e+07,
179  4.79001600000000000000e+08,
180  6.22702080000000000000e+09,
181  8.71782912000000000000e+10,
182  1.30767436800000000000e+12,
183  2.09227898880000000000e+13,
184  3.55687428096000000000e+14,
185  6.40237370572800000000e+15,
186  1.21645100408832000000e+17,
187  2.43290200817664000000e+18,
188  5.10909421717094400000e+19,
189  1.12400072777760768000e+21,
190  2.58520167388849766400e+22,
191  6.20448401733239439360e+23,
192  1.55112100433309859840e+25,
193  4.03291461126605635592e+26,
194  1.08888694504183521614e+28,
195  3.04888344611713860511e+29,
196  8.84176199373970195470e+30,
197  2.65252859812191058647e+32,
198  8.22283865417792281807e+33,
199  2.63130836933693530178e+35,
200  8.68331761881188649615e+36,
201  2.95232799039604140861e+38,
202  1.03331479663861449300e+40,
203  3.71993326789901217463e+41,
204  1.37637530912263450457e+43,
205  5.23022617466601111726e+44,
206  2.03978820811974433568e+46,
207  8.15915283247897734264e+47,
208  3.34525266131638071044e+49,
209  1.40500611775287989839e+51,
210  6.04152630633738356321e+52,
211  2.65827157478844876773e+54,
212  1.19622220865480194551e+56,
213  5.50262215981208894940e+57,
214  2.58623241511168180614e+59,
215  1.24139155925360726691e+61,
216  6.08281864034267560801e+62,
217  3.04140932017133780398e+64,
218  1.55111875328738227999e+66,
219  8.06581751709438785585e+67,
220  4.27488328406002556374e+69,
221  2.30843697339241380441e+71,
222  1.26964033536582759243e+73,
223  7.10998587804863451749e+74,
224  4.05269195048772167487e+76,
225  2.35056133128287857145e+78,
226  1.38683118545689835713e+80,
227  8.32098711274139014271e+81,
228  5.07580213877224798711e+83,
229  3.14699732603879375200e+85,
230  1.98260831540444006372e+87,
231  1.26886932185884164078e+89,
232  8.24765059208247066512e+90,
233  5.44344939077443063905e+92,
234  3.64711109181886852801e+94,
235  2.48003554243683059915e+96,
236  1.71122452428141311337e+98,
237  1.19785716699698917933e+100,
238  8.50478588567862317347e+101,
239  6.12344583768860868500e+103,
240  4.47011546151268434004e+105,
241  3.30788544151938641157e+107,
242  2.48091408113953980872e+109,
243  1.88549470166605025458e+111,
244  1.45183092028285869606e+113,
245  1.13242811782062978295e+115,
246  8.94618213078297528506e+116,
247  7.15694570462638022794e+118,
248  5.79712602074736798470e+120,
249  4.75364333701284174746e+122,
250  3.94552396972065865030e+124,
251  3.31424013456535326627e+126,
252  2.81710411438055027626e+128,
253  2.42270953836727323750e+130,
254  2.10775729837952771662e+132,
255  1.85482642257398439069e+134,
256  1.65079551609084610774e+136,
257  1.48571596448176149700e+138,
258  1.35200152767840296226e+140,
259  1.24384140546413072522e+142,
260  1.15677250708164157442e+144,
261  1.08736615665674307994e+146,
262  1.03299784882390592592e+148,
263  9.91677934870949688836e+149,
264  9.61927596824821198159e+151,
265  9.42689044888324774164e+153,
266  9.33262154439441526381e+155,
267  9.33262154439441526388e+157,
268  9.42594775983835941673e+159,
269  9.61446671503512660515e+161,
270  9.90290071648618040340e+163,
271  1.02990167451456276198e+166,
272  1.08139675824029090008e+168,
273  1.14628056373470835406e+170,
274  1.22652020319613793888e+172,
275  1.32464181945182897396e+174,
276  1.44385958320249358163e+176,
277  1.58824554152274293982e+178,
278  1.76295255109024466316e+180,
279  1.97450685722107402277e+182,
280  2.23119274865981364576e+184,
281  2.54355973347218755612e+186,
282  2.92509369349301568964e+188,
283  3.39310868445189820004e+190,
284  3.96993716080872089396e+192,
285  4.68452584975429065488e+194,
286  5.57458576120760587943e+196,
287  6.68950291344912705515e+198,
288  8.09429852527344373681e+200,
289  9.87504420083360135884e+202,
290  1.21463043670253296712e+205,
291  1.50614174151114087918e+207,
292  1.88267717688892609901e+209,
293  2.37217324288004688470e+211,
294  3.01266001845765954361e+213,
295  3.85620482362580421582e+215,
296  4.97450422247728743840e+217,
297  6.46685548922047366972e+219,
298  8.47158069087882050755e+221,
299  1.11824865119600430699e+224,
300  1.48727070609068572828e+226,
301  1.99294274616151887582e+228,
302  2.69047270731805048244e+230,
303  3.65904288195254865604e+232,
304  5.01288874827499165889e+234,
305  6.91778647261948848943e+236,
306  9.61572319694108900019e+238,
307  1.34620124757175246000e+241,
308  1.89814375907617096864e+243,
309  2.69536413788816277557e+245,
310  3.85437071718007276916e+247,
311  5.55029383273930478744e+249,
312  8.04792605747199194159e+251,
313  1.17499720439091082343e+254,
314  1.72724589045463891049e+256,
315  2.55632391787286558753e+258,
316  3.80892263763056972532e+260,
317  5.71338395644585458806e+262,
318  8.62720977423324042775e+264,
319  1.31133588568345254503e+267,
320  2.00634390509568239384e+269,
321  3.08976961384735088657e+271,
322  4.78914290146339387432e+273,
323  7.47106292628289444390e+275,
324  1.17295687942641442768e+278,
325  1.85327186949373479574e+280,
326  2.94670227249503832518e+282,
327  4.71472363599206132029e+284,
328  7.59070505394721872577e+286,
329  1.22969421873944943358e+289,
330  2.00440157654530257674e+291,
331  3.28721858553429622598e+293,
332  5.42391066613158877297e+295,
333  9.00369170577843736335e+297,
334  1.50361651486499903974e+300,
335  2.52607574497319838672e+302,
336  4.26906800900470527345e+304,
337  7.25741561530799896496e+306
338 };
339 
340 double factorial( long n )
341 {
342  DEBUG_ENTRY( "factorial()" );
343 
344  if( n < 0 || n >= NPRE_FACTORIAL )
345  {
346  fprintf( ioQQQ, "factorial: domain error\n" );
347  cdEXIT(EXIT_FAILURE);
348  }
349 
350  return pre_factorial[n];
351 }
352 
353 /* NB - this implementation is not thread-safe! */
354 
355 class t_lfact : public Singleton<t_lfact>
356 {
357  friend class Singleton<t_lfact>;
358 protected:
360  {
361  p_lf.reserve( 512 );
362  p_lf.push_back( 0. ); /* log10( 0! ) */
363  p_lf.push_back( 0. ); /* log10( 1! ) */
364  }
365 private:
366  vector<double> p_lf;
367 public:
368  double get_lfact( unsigned long n )
369  {
370  if( n < p_lf.size() )
371  {
372  return p_lf[n];
373  }
374  else
375  {
376  for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ )
377  p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) );
378  return p_lf[n];
379  }
380  }
381 };
382 
383 double lfactorial( long n )
384 {
385  /******************************/
386  /* n */
387  /* --- */
388  /* log( n! ) = > log(j) */
389  /* --- */
390  /* j=1 */
391  /******************************/
392 
393  DEBUG_ENTRY( "lfactorial()" );
394 
395  if( n < 0 )
396  {
397  fprintf( ioQQQ, "lfactorial: domain error\n" );
398  cdEXIT(EXIT_FAILURE);
399  }
400 
401  return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
402 }
403 
404 /* complex Gamma function in double precision */
405 /* this routine is a slightly modified version of the one found
406  * at http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html
407  * The following copyright applies:
408  * Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
409  * You may use, copy, modify this code for any purpose and
410  * without fee. You may distribute this ORIGINAL package. */
411 complex<double> cdgamma(complex<double> x)
412 {
413  double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
414 
415  DEBUG_ENTRY( "cdgamma()" );
416 
417  xr = x.real();
418  xi = x.imag();
419  if(xr < 0)
420  {
421  wr = 1. - xr;
422  wi = -xi;
423  }
424  else
425  {
426  wr = xr;
427  wi = xi;
428  }
429  ur = wr + 6.00009857740312429;
430  vr = ur * (wr + 4.99999857982434025) - wi * wi;
431  vi = wi * (wr + 4.99999857982434025) + ur * wi;
432  yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
433  0.293729529320536228;
434  yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
435  ur = vr * (wr + 4.00000003016801681) - vi * wi;
436  ui = vi * (wr + 4.00000003016801681) + vr * wi;
437  vr = ur * (wr + 2.99999999944915534) - ui * wi;
438  vi = ui * (wr + 2.99999999944915534) + ur * wi;
439  yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
440  yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
441  ur = vr * (wr + 2.00000000000603851) - vi * wi;
442  ui = vi * (wr + 2.00000000000603851) + vr * wi;
443  vr = ur * (wr + 0.999999999999975753) - ui * wi;
444  vi = ui * (wr + 0.999999999999975753) + ur * wi;
445  yr += ur * 10.5400280458730808 + vr;
446  yi += ui * 10.5400280458730808 + vi;
447  ur = vr * wr - vi * wi;
448  ui = vi * wr + vr * wi;
449  t = ur * ur + ui * ui;
450  vr = yr * ur + yi * ui + t * 0.0327673720261526849;
451  vi = yi * ur - yr * ui;
452  yr = wr + 7.31790632447016203;
453  ur = log(yr * yr + wi * wi) * 0.5 - 1;
454  ui = atan2(wi, yr);
455  yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
456  yi = ui * (wr - 0.5) + ur * wi;
457  ur = yr * cos(yi);
458  ui = yr * sin(yi);
459  yr = ur * vr - ui * vi;
460  yi = ui * vr + ur * vi;
461  if(xr < 0)
462  {
463  wr = xr * 3.14159265358979324;
464  wi = exp(xi * 3.14159265358979324);
465  vi = 1 / wi;
466  ur = (vi + wi) * sin(wr);
467  ui = (vi - wi) * cos(wr);
468  vr = ur * yr + ui * yi;
469  vi = ui * yr - ur * yi;
470  ur = 6.2831853071795862 / (vr * vr + vi * vi);
471  yr = ur * vr;
472  yi = ur * vi;
473  }
474  return complex<double>( yr, yi );
475 }
476 
477 /*====================================================================
478  *
479  * Below are routines from the Cephes library.
480  *
481  * This is the copyright statement included with the library:
482  *
483  * Some software in this archive may be from the book _Methods and
484  * Programs for Mathematical Functions_ (Prentice-Hall, 1989) or
485  * from the Cephes Mathematical Library, a commercial product. In
486  * either event, it is copyrighted by the author. What you see here
487  * may be used freely but it comes with no support or guarantee.
488  *
489  * The two known misprints in the book are repaired here in the
490  * source listings for the gamma function and the incomplete beta
491  * integral.
492  *
493  * Stephen L. Moshier
494  * moshier@world.std.com
495  *
496  *====================================================================*/
497 
498 /* j0.c
499  *
500  * Bessel function of order zero
501  *
502  *
503  *
504  * SYNOPSIS:
505  *
506  * double x, y, j0();
507  *
508  * y = j0( x );
509  *
510  *
511  *
512  * DESCRIPTION:
513  *
514  * Returns Bessel function of order zero of the argument.
515  *
516  * The domain is divided into the intervals [0, 5] and
517  * (5, infinity). In the first interval the following rational
518  * approximation is used:
519  *
520  *
521  * 2 2
522  * (w - r ) (w - r ) P (w) / Q (w)
523  * 1 2 3 8
524  *
525  * 2
526  * where w = x and the two r's are zeros of the function.
527  *
528  * In the second interval, the Hankel asymptotic expansion
529  * is employed with two rational functions of degree 6/6
530  * and 7/7.
531  *
532  *
533  *
534  * ACCURACY:
535  *
536  * Absolute error:
537  * arithmetic domain # trials peak rms
538  * DEC 0, 30 10000 4.4e-17 6.3e-18
539  * IEEE 0, 30 60000 4.2e-16 1.1e-16
540  *
541  */
542 /* y0.c
543  *
544  * Bessel function of the second kind, order zero
545  *
546  *
547  *
548  * SYNOPSIS:
549  *
550  * double x, y, y0();
551  *
552  * y = y0( x );
553  *
554  *
555  *
556  * DESCRIPTION:
557  *
558  * Returns Bessel function of the second kind, of order
559  * zero, of the argument.
560  *
561  * The domain is divided into the intervals [0, 5] and
562  * (5, infinity). In the first interval a rational approximation
563  * R(x) is employed to compute
564  * y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
565  * Thus a call to j0() is required.
566  *
567  * In the second interval, the Hankel asymptotic expansion
568  * is employed with two rational functions of degree 6/6
569  * and 7/7.
570  *
571  *
572  *
573  * ACCURACY:
574  *
575  * Absolute error, when y0(x) < 1; else relative error:
576  *
577  * arithmetic domain # trials peak rms
578  * DEC 0, 30 9400 7.0e-17 7.9e-18
579  * IEEE 0, 30 30000 1.3e-15 1.6e-16
580  *
581  */
582 
583 /*
584 Cephes Math Library Release 2.8: June, 2000
585 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
586 */
587 
588 /* Note: all coefficients satisfy the relative error criterion
589  * except YP, YQ which are designed for absolute error. */
590 
591 static const double b0_PP[7] = {
592  7.96936729297347051624e-4,
593  8.28352392107440799803e-2,
594  1.23953371646414299388e0,
595  5.44725003058768775090e0,
596  8.74716500199817011941e0,
597  5.30324038235394892183e0,
598  9.99999999999999997821e-1,
599 };
600 
601 static const double b0_PQ[7] = {
602  9.24408810558863637013e-4,
603  8.56288474354474431428e-2,
604  1.25352743901058953537e0,
605  5.47097740330417105182e0,
606  8.76190883237069594232e0,
607  5.30605288235394617618e0,
608  1.00000000000000000218e0,
609 };
610 
611 static const double b0_QP[8] = {
612  -1.13663838898469149931e-2,
613  -1.28252718670509318512e0,
614  -1.95539544257735972385e1,
615  -9.32060152123768231369e1,
616  -1.77681167980488050595e2,
617  -1.47077505154951170175e2,
618  -5.14105326766599330220e1,
619  -6.05014350600728481186e0,
620 };
621 
622 static const double b0_QQ[7] = {
623  /* 1.00000000000000000000e0,*/
624  6.43178256118178023184e1,
625  8.56430025976980587198e2,
626  3.88240183605401609683e3,
627  7.24046774195652478189e3,
628  5.93072701187316984827e3,
629  2.06209331660327847417e3,
630  2.42005740240291393179e2,
631 };
632 
633 static const double b0_YP[8] = {
634  1.55924367855235737965e4,
635  -1.46639295903971606143e7,
636  5.43526477051876500413e9,
637  -9.82136065717911466409e11,
638  8.75906394395366999549e13,
639  -3.46628303384729719441e15,
640  4.42733268572569800351e16,
641  -1.84950800436986690637e16,
642 };
643 
644 static const double b0_YQ[7] = {
645  /* 1.00000000000000000000e0,*/
646  1.04128353664259848412e3,
647  6.26107330137134956842e5,
648  2.68919633393814121987e8,
649  8.64002487103935000337e10,
650  2.02979612750105546709e13,
651  3.17157752842975028269e15,
652  2.50596256172653059228e17,
653 };
654 
655 /* 5.783185962946784521175995758455807035071 */
656 static const double DR1 = 5.78318596294678452118e0;
657 /* 30.47126234366208639907816317502275584842 */
658 static const double DR2 = 3.04712623436620863991e1;
659 
660 static double b0_RP[4] = {
661  -4.79443220978201773821e9,
662  1.95617491946556577543e12,
663  -2.49248344360967716204e14,
664  9.70862251047306323952e15,
665 };
666 
667 static double b0_RQ[8] = {
668  /* 1.00000000000000000000e0,*/
669  4.99563147152651017219e2,
670  1.73785401676374683123e5,
671  4.84409658339962045305e7,
672  1.11855537045356834862e10,
673  2.11277520115489217587e12,
674  3.10518229857422583814e14,
675  3.18121955943204943306e16,
676  1.71086294081043136091e18,
677 };
678 
679 static const double TWOOPI = 2./PI;
680 static const double SQ2OPI = sqrt(2./PI);
681 static const double PIO4 = PI/4.;
682 
683 double bessel_j0(double x)
684 {
685  double w, z, p, q, xn;
686 
687  DEBUG_ENTRY( "bessel_j0()" );
688 
689  if( x < 0 )
690  x = -x;
691 
692  if( x <= 5.0 )
693  {
694  z = x * x;
695  if( x < 1.0e-5 )
696  return 1.0 - z/4.0;
697 
698  p = (z - DR1) * (z - DR2);
699  p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
700  return p;
701  }
702 
703  w = 5.0/x;
704  q = 25.0/(x*x);
705  p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
706  q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
707  xn = x - PIO4;
708  p = p * cos(xn) - w * q * sin(xn);
709  return p * SQ2OPI / sqrt(x);
710 }
711 
712 /* y0() 2 */
713 /* Bessel function of second kind, order zero */
714 
715 /* Rational approximation coefficients YP[], YQ[] are used here.
716  * The function computed is y0(x) - 2 * log(x) * j0(x) / PI,
717  * whose value at x = 0 is 2 * ( log(0.5) + EULER ) / PI
718  * = 0.073804295108687225.
719  */
720 
721 double bessel_y0(double x)
722 {
723  double w, z, p, q, xn;
724 
725  DEBUG_ENTRY( "bessel_y0()" );
726 
727  if( x <= 5.0 )
728  {
729  if( x <= 0.0 )
730  {
731  fprintf( ioQQQ, "bessel_y0: domain error\n" );
732  cdEXIT(EXIT_FAILURE);
733  }
734  z = x * x;
735  w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
736  w += TWOOPI * log(x) * bessel_j0(x);
737  return w;
738  }
739 
740  w = 5.0/x;
741  z = 25.0 / (x * x);
742  p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
743  q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
744  xn = x - PIO4;
745  p = p * sin(xn) + w * q * cos(xn);
746  return p * SQ2OPI / sqrt(x);
747 }
748 
749 /* j1.c
750  *
751  * Bessel function of order one
752  *
753  *
754  *
755  * SYNOPSIS:
756  *
757  * double x, y, j1();
758  *
759  * y = j1( x );
760  *
761  *
762  *
763  * DESCRIPTION:
764  *
765  * Returns Bessel function of order one of the argument.
766  *
767  * The domain is divided into the intervals [0, 8] and
768  * (8, infinity). In the first interval a 24 term Chebyshev
769  * expansion is used. In the second, the asymptotic
770  * trigonometric representation is employed using two
771  * rational functions of degree 5/5.
772  *
773  *
774  *
775  * ACCURACY:
776  *
777  * Absolute error:
778  * arithmetic domain # trials peak rms
779  * DEC 0, 30 10000 4.0e-17 1.1e-17
780  * IEEE 0, 30 30000 2.6e-16 1.1e-16
781  *
782  *
783  */
784 /* y1.c
785  *
786  * Bessel function of second kind of order one
787  *
788  *
789  *
790  * SYNOPSIS:
791  *
792  * double x, y, y1();
793  *
794  * y = y1( x );
795  *
796  *
797  *
798  * DESCRIPTION:
799  *
800  * Returns Bessel function of the second kind of order one
801  * of the argument.
802  *
803  * The domain is divided into the intervals [0, 8] and
804  * (8, infinity). In the first interval a 25 term Chebyshev
805  * expansion is used, and a call to j1() is required.
806  * In the second, the asymptotic trigonometric representation
807  * is employed using two rational functions of degree 5/5.
808  *
809  *
810  *
811  * ACCURACY:
812  *
813  * Absolute error:
814  * arithmetic domain # trials peak rms
815  * DEC 0, 30 10000 8.6e-17 1.3e-17
816  * IEEE 0, 30 30000 1.0e-15 1.3e-16
817  *
818  * (error criterion relative when |y1| > 1).
819  *
820  */
821 
822 /*
823 Cephes Math Library Release 2.8: June, 2000
824 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
825 */
826 
827 static const double b1_RP[4] = {
828  -8.99971225705559398224e8,
829  4.52228297998194034323e11,
830  -7.27494245221818276015e13,
831  3.68295732863852883286e15,
832 };
833 
834 static const double b1_RQ[8] = {
835  /* 1.00000000000000000000E0,*/
836  6.20836478118054335476e2,
837  2.56987256757748830383e5,
838  8.35146791431949253037e7,
839  2.21511595479792499675e10,
840  4.74914122079991414898e12,
841  7.84369607876235854894e14,
842  8.95222336184627338078e16,
843  5.32278620332680085395e18,
844 };
845 
846 static const double b1_PP[7] = {
847  7.62125616208173112003e-4,
848  7.31397056940917570436e-2,
849  1.12719608129684925192e0,
850  5.11207951146807644818e0,
851  8.42404590141772420927e0,
852  5.21451598682361504063e0,
853  1.00000000000000000254e0,
854 };
855 
856 static const double b1_PQ[7] = {
857  5.71323128072548699714e-4,
858  6.88455908754495404082e-2,
859  1.10514232634061696926e0,
860  5.07386386128601488557e0,
861  8.39985554327604159757e0,
862  5.20982848682361821619e0,
863  9.99999999999999997461e-1,
864 };
865 
866 static const double b1_QP[8] = {
867  5.10862594750176621635e-2,
868  4.98213872951233449420e0,
869  7.58238284132545283818e1,
870  3.66779609360150777800e2,
871  7.10856304998926107277e2,
872  5.97489612400613639965e2,
873  2.11688757100572135698e2,
874  2.52070205858023719784e1,
875 };
876 
877 static const double b1_QQ[7] = {
878  /* 1.00000000000000000000e0,*/
879  7.42373277035675149943e1,
880  1.05644886038262816351e3,
881  4.98641058337653607651e3,
882  9.56231892404756170795e3,
883  7.99704160447350683650e3,
884  2.82619278517639096600e3,
885  3.36093607810698293419e2,
886 };
887 
888 static const double b1_YP[6] = {
889  1.26320474790178026440e9,
890  -6.47355876379160291031e11,
891  1.14509511541823727583e14,
892  -8.12770255501325109621e15,
893  2.02439475713594898196e17,
894  -7.78877196265950026825e17,
895 };
896 
897 static const double b1_YQ[8] = {
898  /* 1.00000000000000000000E0,*/
899  5.94301592346128195359E2,
900  2.35564092943068577943E5,
901  7.34811944459721705660E7,
902  1.87601316108706159478E10,
903  3.88231277496238566008E12,
904  6.20557727146953693363E14,
905  6.87141087355300489866E16,
906  3.97270608116560655612E18,
907 };
908 
909 static const double Z1 = 1.46819706421238932572E1;
910 static const double Z2 = 4.92184563216946036703E1;
911 
912 static const double THPIO4 = 3.*PI/4.;
913 
914 double bessel_j1(double x)
915 {
916  double w, z, p, q, xn;
917 
918  DEBUG_ENTRY( "bessel_j1()" );
919 
920  w = x;
921  if( x < 0 )
922  w = -x;
923 
924  if( w <= 5.0 )
925  {
926  z = x * x;
927  w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
928  w = w * x * (z - Z1) * (z - Z2);
929  return w;
930  }
931 
932  w = 5.0/x;
933  z = w * w;
934  p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
935  q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
936  xn = x - THPIO4;
937  p = p * cos(xn) - w * q * sin(xn);
938  return p * SQ2OPI / sqrt(x);
939 }
940 
941 double bessel_y1(double x)
942 {
943  double w, z, p, q, xn;
944 
945  DEBUG_ENTRY( "bessel_y1()" );
946 
947  if( x <= 5.0 )
948  {
949  if( x <= 0.0 )
950  {
951  fprintf( ioQQQ, "bessel_y1: domain error\n" );
952  cdEXIT(EXIT_FAILURE);
953  }
954  z = x * x;
955  w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
956  w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
957  return w;
958  }
959 
960  w = 5.0/x;
961  z = w * w;
962  p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
963  q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
964  xn = x - THPIO4;
965  p = p * sin(xn) + w * q * cos(xn);
966  return p * SQ2OPI / sqrt(x);
967 }
968 
969 /* jn.c
970  *
971  * Bessel function of integer order
972  *
973  *
974  *
975  * SYNOPSIS:
976  *
977  * int n;
978  * double x, y, jn();
979  *
980  * y = jn( n, x );
981  *
982  *
983  *
984  * DESCRIPTION:
985  *
986  * Returns Bessel function of order n, where n is a
987  * (possibly negative) integer.
988  *
989  * The ratio of jn(x) to j0(x) is computed by backward
990  * recurrence. First the ratio jn/jn-1 is found by a
991  * continued fraction expansion. Then the recurrence
992  * relating successive orders is applied until j0 or j1 is
993  * reached.
994  *
995  * If n = 0 or 1 the routine for j0 or j1 is called
996  * directly.
997  *
998  *
999  *
1000  * ACCURACY:
1001  *
1002  * Absolute error:
1003  * arithmetic range # trials peak rms
1004  * DEC 0, 30 5500 6.9e-17 9.3e-18
1005  * IEEE 0, 30 5000 4.4e-16 7.9e-17
1006  *
1007  *
1008  * Not suitable for large n or x. Use jv() instead.
1009  *
1010  */
1011 
1012 /* jn.c
1013 Cephes Math Library Release 2.8: June, 2000
1014 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1015 */
1016 
1017 double bessel_jn(int n, double x)
1018 {
1019  double pkm2, pkm1, pk, xk, r, ans;
1020  int k, sign;
1021 
1022  DEBUG_ENTRY( "bessel_jn()" );
1023 
1024  if( n < 0 )
1025  {
1026  n = -n;
1027  if( (n & 1) == 0 ) /* -1**n */
1028  sign = 1;
1029  else
1030  sign = -1;
1031  }
1032  else
1033  sign = 1;
1034 
1035  if( x < 0.0 )
1036  {
1037  if( n & 1 )
1038  sign = -sign;
1039  x = -x;
1040  }
1041 
1042  if( n == 0 )
1043  {
1044  return sign * bessel_j0(x);
1045  }
1046  if( n == 1 )
1047  {
1048  return sign * bessel_j1(x);
1049  }
1050  if( n == 2 )
1051  {
1052  return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
1053  }
1054 
1055  if( x < DBL_EPSILON )
1056  {
1057  return 0.0;
1058  }
1059 
1060  /* continued fraction */
1061  k = 52;
1062 
1063  pk = 2 * (n + k);
1064  ans = pk;
1065  xk = x * x;
1066 
1067  do
1068  {
1069  pk -= 2.0;
1070  ans = pk - (xk/ans);
1071  }
1072  while( --k > 0 );
1073  ans = x/ans;
1074 
1075  /* backward recurrence */
1076  pk = 1.0;
1077  pkm1 = 1.0/ans;
1078  k = n-1;
1079  r = 2 * k;
1080 
1081  do
1082  {
1083  pkm2 = (pkm1 * r - pk * x) / x;
1084  pk = pkm1;
1085  pkm1 = pkm2;
1086  r -= 2.0;
1087  }
1088  while( --k > 0 );
1089 
1090  if( fabs(pk) > fabs(pkm1) )
1091  ans = bessel_j1(x)/pk;
1092  else
1093  ans = bessel_j0(x)/pkm1;
1094  return sign * ans;
1095 }
1096 
1097 /* yn.c
1098  *
1099  * Bessel function of second kind of integer order
1100  *
1101  *
1102  *
1103  * SYNOPSIS:
1104  *
1105  * double x, y, yn();
1106  * int n;
1107  *
1108  * y = yn( n, x );
1109  *
1110  *
1111  *
1112  * DESCRIPTION:
1113  *
1114  * Returns Bessel function of order n, where n is a
1115  * (possibly negative) integer.
1116  *
1117  * The function is evaluated by forward recurrence on
1118  * n, starting with values computed by the routines
1119  * y0() and y1().
1120  *
1121  * If n = 0 or 1 the routine for y0 or y1 is called
1122  * directly.
1123  *
1124  *
1125  *
1126  * ACCURACY:
1127  *
1128  *
1129  * Absolute error, except relative
1130  * when y > 1:
1131  * arithmetic domain # trials peak rms
1132  * DEC 0, 30 2200 2.9e-16 5.3e-17
1133  * IEEE 0, 30 30000 3.4e-15 4.3e-16
1134  *
1135  *
1136  * ERROR MESSAGES:
1137  *
1138  * message condition value returned
1139  * yn singularity x = 0 MAXNUM
1140  * yn overflow MAXNUM
1141  *
1142  * Spot checked against tables for x, n between 0 and 100.
1143  *
1144  */
1145 
1146 /*
1147 Cephes Math Library Release 2.8: June, 2000
1148 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1149 */
1150 
1151 double bessel_yn(int n, double x)
1152 {
1153  double an, anm1, anm2, r;
1154  int k, sign;
1155 
1156  DEBUG_ENTRY( "bessel_yn()" );
1157 
1158  if( n < 0 )
1159  {
1160  n = -n;
1161  if( (n & 1) == 0 ) /* -1**n */
1162  sign = 1;
1163  else
1164  sign = -1;
1165  }
1166  else
1167  sign = 1;
1168 
1169  if( n == 0 )
1170  {
1171  return sign * bessel_y0(x);
1172  }
1173  if( n == 1 )
1174  {
1175  return sign * bessel_y1(x);
1176  }
1177 
1178  /* test for overflow */
1179  if( x <= 0.0 )
1180  {
1181  fprintf( ioQQQ, "bessel_yn: domain error\n" );
1182  cdEXIT(EXIT_FAILURE);
1183  }
1184 
1185  /* forward recurrence on n */
1186  anm2 = bessel_y0(x);
1187  anm1 = bessel_y1(x);
1188  k = 1;
1189  r = 2 * k;
1190  do
1191  {
1192  an = r * anm1 / x - anm2;
1193  anm2 = anm1;
1194  anm1 = an;
1195  r += 2.0;
1196  ++k;
1197  }
1198  while( k < n );
1199  return sign * an;
1200 }
1201 
1202 /* k0.c
1203  *
1204  * Modified Bessel function, third kind, order zero
1205  *
1206  *
1207  *
1208  * SYNOPSIS:
1209  *
1210  * double x, y, k0();
1211  *
1212  * y = k0( x );
1213  *
1214  *
1215  *
1216  * DESCRIPTION:
1217  *
1218  * Returns modified Bessel function of the third kind
1219  * of order zero of the argument.
1220  *
1221  * The range is partitioned into the two intervals [0,8] and
1222  * (8, infinity). Chebyshev polynomial expansions are employed
1223  * in each interval.
1224  *
1225  *
1226  *
1227  * ACCURACY:
1228  *
1229  * Tested at 2000 random points between 0 and 8. Peak absolute
1230  * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
1231  * Relative error:
1232  * arithmetic domain # trials peak rms
1233  * DEC 0, 30 3100 1.3e-16 2.1e-17
1234  * IEEE 0, 30 30000 1.2e-15 1.6e-16
1235  *
1236  * ERROR MESSAGES:
1237  *
1238  * message condition value returned
1239  * K0 domain x <= 0 MAXNUM
1240  *
1241  */
1242 /* k0e()
1243  *
1244  * Modified Bessel function, third kind, order zero,
1245  * exponentially scaled
1246  *
1247  *
1248  *
1249  * SYNOPSIS:
1250  *
1251  * double x, y, k0e();
1252  *
1253  * y = k0e( x );
1254  *
1255  *
1256  *
1257  * DESCRIPTION:
1258  *
1259  * Returns exponentially scaled modified Bessel function
1260  * of the third kind of order zero of the argument.
1261  *
1262  *
1263  *
1264  * ACCURACY:
1265  *
1266  * Relative error:
1267  * arithmetic domain # trials peak rms
1268  * IEEE 0, 30 30000 1.4e-15 1.4e-16
1269  * See k0().
1270  *
1271  */
1272 
1273 /*
1274 Cephes Math Library Release 2.8: June, 2000
1275 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1276 */
1277 
1278 /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
1279  * in the interval [0,2]. The odd order coefficients are all
1280  * zero; only the even order coefficients are listed.
1281  *
1282  * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EULER.
1283  */
1284 
1285 static const double k0_A[] =
1286 {
1287  1.37446543561352307156e-16,
1288  4.25981614279661018399e-14,
1289  1.03496952576338420167e-11,
1290  1.90451637722020886025e-9,
1291  2.53479107902614945675e-7,
1292  2.28621210311945178607e-5,
1293  1.26461541144692592338e-3,
1294  3.59799365153615016266e-2,
1295  3.44289899924628486886e-1,
1296  -5.35327393233902768720e-1
1297 };
1298 
1299 /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
1300  * in the inverted interval [2,infinity].
1301  *
1302  * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
1303  */
1304 
1305 static const double k0_B[] = {
1306  5.30043377268626276149e-18,
1307  -1.64758043015242134646e-17,
1308  5.21039150503902756861e-17,
1309  -1.67823109680541210385e-16,
1310  5.51205597852431940784e-16,
1311  -1.84859337734377901440e-15,
1312  6.34007647740507060557e-15,
1313  -2.22751332699166985548e-14,
1314  8.03289077536357521100e-14,
1315  -2.98009692317273043925e-13,
1316  1.14034058820847496303e-12,
1317  -4.51459788337394416547e-12,
1318  1.85594911495471785253e-11,
1319  -7.95748924447710747776e-11,
1320  3.57739728140030116597e-10,
1321  -1.69753450938905987466e-9,
1322  8.57403401741422608519e-9,
1323  -4.66048989768794782956e-8,
1324  2.76681363944501510342e-7,
1325  -1.83175552271911948767e-6,
1326  1.39498137188764993662e-5,
1327  -1.28495495816278026384e-4,
1328  1.56988388573005337491e-3,
1329  -3.14481013119645005427e-2,
1330  2.44030308206595545468e0
1331 };
1332 
1333 double bessel_k0(double x)
1334 {
1335  double y, z;
1336 
1337  DEBUG_ENTRY( "bessel_k0()" );
1338 
1339  if( x <= 0.0 )
1340  {
1341  fprintf( ioQQQ, "bessel_k0: domain error\n" );
1342  cdEXIT(EXIT_FAILURE);
1343  }
1344 
1345  if( x <= 2.0 )
1346  {
1347  y = x * x - 2.0;
1348  y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
1349  return y;
1350  }
1351  z = 8.0/x - 2.0;
1352  y = exp(-x) * chbevl( z, k0_B, 25 ) / sqrt(x);
1353  return y;
1354 }
1355 
1356 double bessel_k0_scaled(double x)
1357 {
1358  double y;
1359 
1360  DEBUG_ENTRY( "bessel_k0_scaled()" );
1361 
1362  if( x <= 0.0 )
1363  {
1364  fprintf( ioQQQ, "bessel_k0_scaled: domain error\n" );
1365  cdEXIT(EXIT_FAILURE);
1366  }
1367 
1368  if( x <= 2.0 )
1369  {
1370  y = x * x - 2.0;
1371  y = chbevl( y, k0_A, 10 ) - log( 0.5 * x ) * bessel_i0(x);
1372  return y * exp(x);
1373  }
1374  return chbevl( 8.0/x - 2.0, k0_B, 25 ) / sqrt(x);
1375 }
1376 
1377 /* k1.c
1378  *
1379  * Modified Bessel function, third kind, order one
1380  *
1381  *
1382  *
1383  * SYNOPSIS:
1384  *
1385  * double x, y, k1();
1386  *
1387  * y = k1( x );
1388  *
1389  *
1390  *
1391  * DESCRIPTION:
1392  *
1393  * Computes the modified Bessel function of the third kind
1394  * of order one of the argument.
1395  *
1396  * The range is partitioned into the two intervals [0,2] and
1397  * (2, infinity). Chebyshev polynomial expansions are employed
1398  * in each interval.
1399  *
1400  *
1401  *
1402  * ACCURACY:
1403  *
1404  * Relative error:
1405  * arithmetic domain # trials peak rms
1406  * DEC 0, 30 3300 8.9e-17 2.2e-17
1407  * IEEE 0, 30 30000 1.2e-15 1.6e-16
1408  *
1409  * ERROR MESSAGES:
1410  *
1411  * message condition value returned
1412  * k1 domain x <= 0 MAXNUM
1413  *
1414  */
1415 /* k1e.c
1416  *
1417  * Modified Bessel function, third kind, order one,
1418  * exponentially scaled
1419  *
1420  *
1421  *
1422  * SYNOPSIS:
1423  *
1424  * double x, y, k1e();
1425  *
1426  * y = k1e( x );
1427  *
1428  *
1429  *
1430  * DESCRIPTION:
1431  *
1432  * Returns exponentially scaled modified Bessel function
1433  * of the third kind of order one of the argument:
1434  *
1435  * k1e(x) = exp(x) * k1(x).
1436  *
1437  *
1438  *
1439  * ACCURACY:
1440  *
1441  * Relative error:
1442  * arithmetic domain # trials peak rms
1443  * IEEE 0, 30 30000 7.8e-16 1.2e-16
1444  * See k1().
1445  *
1446  */
1447 
1448 /*
1449 Cephes Math Library Release 2.8: June, 2000
1450 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1451 */
1452 
1453 /* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
1454  * in the interval [0,2].
1455  *
1456  * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
1457  */
1458 
1459 static const double k1_A[] =
1460 {
1461  -7.02386347938628759343e-18,
1462  -2.42744985051936593393e-15,
1463  -6.66690169419932900609e-13,
1464  -1.41148839263352776110e-10,
1465  -2.21338763073472585583e-8,
1466  -2.43340614156596823496e-6,
1467  -1.73028895751305206302e-4,
1468  -6.97572385963986435018e-3,
1469  -1.22611180822657148235e-1,
1470  -3.53155960776544875667e-1,
1471  1.52530022733894777053e0
1472 };
1473 
1474 /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
1475  * in the interval [2,infinity].
1476  *
1477  * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
1478  */
1479 
1480 static const double k1_B[] =
1481 {
1482  -5.75674448366501715755e-18,
1483  1.79405087314755922667e-17,
1484  -5.68946255844285935196e-17,
1485  1.83809354436663880070e-16,
1486  -6.05704724837331885336e-16,
1487  2.03870316562433424052e-15,
1488  -7.01983709041831346144e-15,
1489  2.47715442448130437068e-14,
1490  -8.97670518232499435011e-14,
1491  3.34841966607842919884e-13,
1492  -1.28917396095102890680e-12,
1493  5.13963967348173025100e-12,
1494  -2.12996783842756842877e-11,
1495  9.21831518760500529508e-11,
1496  -4.19035475934189648750e-10,
1497  2.01504975519703286596e-9,
1498  -1.03457624656780970260e-8,
1499  5.74108412545004946722e-8,
1500  -3.50196060308781257119e-7,
1501  2.40648494783721712015e-6,
1502  -1.93619797416608296024e-5,
1503  1.95215518471351631108e-4,
1504  -2.85781685962277938680e-3,
1505  1.03923736576817238437e-1,
1506  2.72062619048444266945e0
1507 };
1508 
1509 double bessel_k1(double x)
1510 {
1511  double y, z;
1512 
1513  DEBUG_ENTRY( "bessel_k1()" );
1514 
1515  z = 0.5 * x;
1516  if( z <= 0.0 )
1517  {
1518  fprintf( ioQQQ, "bessel_k1: domain error\n" );
1519  cdEXIT(EXIT_FAILURE);
1520  }
1521 
1522  if( x <= 2.0 )
1523  {
1524  y = x * x - 2.0;
1525  y = log(z) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
1526  return y;
1527  }
1528  return exp(-x) * chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
1529 }
1530 
1531 double bessel_k1_scaled(double x)
1532 {
1533  double y;
1534 
1535  DEBUG_ENTRY( "bessel_k1_scaled()" );
1536 
1537  if( x <= 0.0 )
1538  {
1539  fprintf( ioQQQ, "bessel_k1_scaled: domain error\n" );
1540  cdEXIT(EXIT_FAILURE);
1541  }
1542 
1543  if( x <= 2.0 )
1544  {
1545  y = x * x - 2.0;
1546  y = log( 0.5 * x ) * bessel_i1(x) + chbevl( y, k1_A, 11 ) / x;
1547  return y * exp(x);
1548  }
1549  return chbevl( 8.0/x - 2.0, k1_B, 25 ) / sqrt(x);
1550 }
1551 
1552 /* i0.c
1553  *
1554  * Modified Bessel function of order zero
1555  *
1556  *
1557  *
1558  * SYNOPSIS:
1559  *
1560  * double x, y, i0();
1561  *
1562  * y = i0( x );
1563  *
1564  *
1565  *
1566  * DESCRIPTION:
1567  *
1568  * Returns modified Bessel function of order zero of the
1569  * argument.
1570  *
1571  * The function is defined as i0(x) = j0( ix ).
1572  *
1573  * The range is partitioned into the two intervals [0,8] and
1574  * (8, infinity). Chebyshev polynomial expansions are employed
1575  * in each interval.
1576  *
1577  *
1578  *
1579  * ACCURACY:
1580  *
1581  * Relative error:
1582  * arithmetic domain # trials peak rms
1583  * DEC 0,30 6000 8.2e-17 1.9e-17
1584  * IEEE 0,30 30000 5.8e-16 1.4e-16
1585  *
1586  */
1587 /* i0e.c
1588  *
1589  * Modified Bessel function of order zero,
1590  * exponentially scaled
1591  *
1592  *
1593  *
1594  * SYNOPSIS:
1595  *
1596  * double x, y, i0e();
1597  *
1598  * y = i0e( x );
1599  *
1600  *
1601  *
1602  * DESCRIPTION:
1603  *
1604  * Returns exponentially scaled modified Bessel function
1605  * of order zero of the argument.
1606  *
1607  * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
1608  *
1609  *
1610  *
1611  * ACCURACY:
1612  *
1613  * Relative error:
1614  * arithmetic domain # trials peak rms
1615  * IEEE 0,30 30000 5.4e-16 1.2e-16
1616  * See i0().
1617  *
1618  */
1619 
1620 /*
1621 Cephes Math Library Release 2.8: June, 2000
1622 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1623 */
1624 
1625 /* Chebyshev coefficients for exp(-x) I0(x)
1626  * in the interval [0,8].
1627  *
1628  * lim(x->0){ exp(-x) I0(x) } = 1.
1629  */
1630 
1631 static const double i0_A[] =
1632 {
1633  -4.41534164647933937950e-18,
1634  3.33079451882223809783e-17,
1635  -2.43127984654795469359e-16,
1636  1.71539128555513303061e-15,
1637  -1.16853328779934516808e-14,
1638  7.67618549860493561688e-14,
1639  -4.85644678311192946090e-13,
1640  2.95505266312963983461e-12,
1641  -1.72682629144155570723e-11,
1642  9.67580903537323691224e-11,
1643  -5.18979560163526290666e-10,
1644  2.65982372468238665035e-9,
1645  -1.30002500998624804212e-8,
1646  6.04699502254191894932e-8,
1647  -2.67079385394061173391e-7,
1648  1.11738753912010371815e-6,
1649  -4.41673835845875056359e-6,
1650  1.64484480707288970893e-5,
1651  -5.75419501008210370398e-5,
1652  1.88502885095841655729e-4,
1653  -5.76375574538582365885e-4,
1654  1.63947561694133579842e-3,
1655  -4.32430999505057594430e-3,
1656  1.05464603945949983183e-2,
1657  -2.37374148058994688156e-2,
1658  4.93052842396707084878e-2,
1659  -9.49010970480476444210e-2,
1660  1.71620901522208775349e-1,
1661  -3.04682672343198398683e-1,
1662  6.76795274409476084995e-1
1663 };
1664 
1665 /* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
1666  * in the inverted interval [8,infinity].
1667  *
1668  * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
1669  */
1670 
1671 static const double i0_B[] =
1672 {
1673  -7.23318048787475395456e-18,
1674  -4.83050448594418207126e-18,
1675  4.46562142029675999901e-17,
1676  3.46122286769746109310e-17,
1677  -2.82762398051658348494e-16,
1678  -3.42548561967721913462e-16,
1679  1.77256013305652638360e-15,
1680  3.81168066935262242075e-15,
1681  -9.55484669882830764870e-15,
1682  -4.15056934728722208663e-14,
1683  1.54008621752140982691e-14,
1684  3.85277838274214270114e-13,
1685  7.18012445138366623367e-13,
1686  -1.79417853150680611778e-12,
1687  -1.32158118404477131188e-11,
1688  -3.14991652796324136454e-11,
1689  1.18891471078464383424e-11,
1690  4.94060238822496958910e-10,
1691  3.39623202570838634515e-9,
1692  2.26666899049817806459e-8,
1693  2.04891858946906374183e-7,
1694  2.89137052083475648297e-6,
1695  6.88975834691682398426e-5,
1696  3.36911647825569408990e-3,
1697  8.04490411014108831608e-1
1698 };
1699 
1700 double bessel_i0(double x)
1701 {
1702  double y;
1703 
1704  DEBUG_ENTRY( "bessel_i0()" );
1705 
1706  if( x < 0 )
1707  x = -x;
1708 
1709  if( x <= 8.0 )
1710  {
1711  y = (x/2.0) - 2.0;
1712  return exp(x) * chbevl( y, i0_A, 30 );
1713  }
1714  return exp(x) * chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
1715 }
1716 
1717 double bessel_i0_scaled(double x)
1718 {
1719  double y;
1720 
1721  DEBUG_ENTRY( "bessel_i0_scaled()" );
1722 
1723  if( x < 0 )
1724  x = -x;
1725 
1726  if( x <= 8.0 )
1727  {
1728  y = (x/2.0) - 2.0;
1729  return chbevl( y, i0_A, 30 );
1730  }
1731  return chbevl( 32.0/x - 2.0, i0_B, 25 ) / sqrt(x);
1732 }
1733 
1734 /* i1.c
1735  *
1736  * Modified Bessel function of order one
1737  *
1738  *
1739  *
1740  * SYNOPSIS:
1741  *
1742  * double x, y, i1();
1743  *
1744  * y = i1( x );
1745  *
1746  *
1747  *
1748  * DESCRIPTION:
1749  *
1750  * Returns modified Bessel function of order one of the
1751  * argument.
1752  *
1753  * The function is defined as i1(x) = -i j1( ix ).
1754  *
1755  * The range is partitioned into the two intervals [0,8] and
1756  * (8, infinity). Chebyshev polynomial expansions are employed
1757  * in each interval.
1758  *
1759  *
1760  *
1761  * ACCURACY:
1762  *
1763  * Relative error:
1764  * arithmetic domain # trials peak rms
1765  * DEC 0, 30 3400 1.2e-16 2.3e-17
1766  * IEEE 0, 30 30000 1.9e-15 2.1e-16
1767  *
1768  *
1769  */
1770 /* i1e.c
1771  *
1772  * Modified Bessel function of order one,
1773  * exponentially scaled
1774  *
1775  *
1776  *
1777  * SYNOPSIS:
1778  *
1779  * double x, y, i1e();
1780  *
1781  * y = i1e( x );
1782  *
1783  *
1784  *
1785  * DESCRIPTION:
1786  *
1787  * Returns exponentially scaled modified Bessel function
1788  * of order one of the argument.
1789  *
1790  * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
1791  *
1792  *
1793  *
1794  * ACCURACY:
1795  *
1796  * Relative error:
1797  * arithmetic domain # trials peak rms
1798  * IEEE 0, 30 30000 2.0e-15 2.0e-16
1799  * See i1().
1800  *
1801  */
1802 
1803 /*
1804 Cephes Math Library Release 2.8: June, 2000
1805 Copyright 1985, 1987, 2000 by Stephen L. Moshier
1806 */
1807 
1808 /* Chebyshev coefficients for exp(-x) I1(x) / x
1809  * in the interval [0,8].
1810  *
1811  * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
1812  */
1813 
1814 static double i1_A[] =
1815 {
1816  2.77791411276104639959e-18,
1817  -2.11142121435816608115e-17,
1818  1.55363195773620046921e-16,
1819  -1.10559694773538630805e-15,
1820  7.60068429473540693410e-15,
1821  -5.04218550472791168711e-14,
1822  3.22379336594557470981e-13,
1823  -1.98397439776494371520e-12,
1824  1.17361862988909016308e-11,
1825  -6.66348972350202774223e-11,
1826  3.62559028155211703701e-10,
1827  -1.88724975172282928790e-9,
1828  9.38153738649577178388e-9,
1829  -4.44505912879632808065e-8,
1830  2.00329475355213526229e-7,
1831  -8.56872026469545474066e-7,
1832  3.47025130813767847674e-6,
1833  -1.32731636560394358279e-5,
1834  4.78156510755005422638e-5,
1835  -1.61760815825896745588e-4,
1836  5.12285956168575772895e-4,
1837  -1.51357245063125314899e-3,
1838  4.15642294431288815669e-3,
1839  -1.05640848946261981558e-2,
1840  2.47264490306265168283e-2,
1841  -5.29459812080949914269e-2,
1842  1.02643658689847095384e-1,
1843  -1.76416518357834055153e-1,
1844  2.52587186443633654823e-1
1845 };
1846 
1847 /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
1848  * in the inverted interval [8,infinity].
1849  *
1850  * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
1851  */
1852 
1853 static double i1_B[] =
1854 {
1855  7.51729631084210481353e-18,
1856  4.41434832307170791151e-18,
1857  -4.65030536848935832153e-17,
1858  -3.20952592199342395980e-17,
1859  2.96262899764595013876e-16,
1860  3.30820231092092828324e-16,
1861  -1.88035477551078244854e-15,
1862  -3.81440307243700780478e-15,
1863  1.04202769841288027642e-14,
1864  4.27244001671195135429e-14,
1865  -2.10154184277266431302e-14,
1866  -4.08355111109219731823e-13,
1867  -7.19855177624590851209e-13,
1868  2.03562854414708950722e-12,
1869  1.41258074366137813316e-11,
1870  3.25260358301548823856e-11,
1871  -1.89749581235054123450e-11,
1872  -5.58974346219658380687e-10,
1873  -3.83538038596423702205e-9,
1874  -2.63146884688951950684e-8,
1875  -2.51223623787020892529e-7,
1876  -3.88256480887769039346e-6,
1877  -1.10588938762623716291e-4,
1878  -9.76109749136146840777e-3,
1879  7.78576235018280120474e-1
1880 };
1881 
1882 double bessel_i1(double x)
1883 {
1884  double y, z;
1885 
1886  DEBUG_ENTRY( "bessel_i1()" );
1887 
1888  z = fabs(x);
1889  if( z <= 8.0 )
1890  {
1891  y = (z/2.0) - 2.0;
1892  z = chbevl( y, i1_A, 29 ) * z * exp(z);
1893  }
1894  else
1895  {
1896  z = exp(z) * chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
1897  }
1898  if( x < 0.0 )
1899  z = -z;
1900  return z;
1901 }
1902 
1903 double bessel_i1_scaled(double x)
1904 {
1905  double y, z;
1906 
1907  DEBUG_ENTRY( "bessel_i1_scaled()" );
1908 
1909  z = fabs(x);
1910  if( z <= 8.0 )
1911  {
1912  y = (z/2.0) - 2.0;
1913  z = chbevl( y, i1_A, 29 ) * z;
1914  }
1915  else
1916  {
1917  z = chbevl( 32.0/z - 2.0, i1_B, 25 ) / sqrt(z);
1918  }
1919  if( x < 0.0 )
1920  z = -z;
1921  return z;
1922 }
1923 
1924 /* ellpk.c
1925  *
1926  * Complete elliptic integral of the first kind
1927  *
1928  *
1929  *
1930  * SYNOPSIS:
1931  *
1932  * double m1, y, ellpk();
1933  *
1934  * y = ellpk( m1 );
1935  *
1936  *
1937  *
1938  * DESCRIPTION:
1939  *
1940  * Approximates the integral
1941  *
1942  *
1943  *
1944  * pi/2
1945  * -
1946  * | |
1947  * | dt
1948  * K(m) = | ------------------
1949  * | 2
1950  * | | sqrt( 1 - m sin t )
1951  * -
1952  * 0
1953  *
1954  * where m = 1 - m1, using the approximation
1955  *
1956  * P(x) - log x Q(x).
1957  *
1958  * The argument m1 is used rather than m so that the logarithmic
1959  * singularity at m = 1 will be shifted to the origin; this
1960  * preserves maximum accuracy.
1961  *
1962  * K(0) = pi/2.
1963  *
1964  * ACCURACY:
1965  *
1966  * Relative error:
1967  * arithmetic domain # trials peak rms
1968  * DEC 0,1 16000 3.5e-17 1.1e-17
1969  * IEEE 0,1 30000 2.5e-16 6.8e-17
1970  *
1971  * ERROR MESSAGES:
1972  *
1973  * message condition value returned
1974  * ellpk domain x<0, x>1 0.0
1975  *
1976  */
1977 
1978 /*
1979 Cephes Math Library, Release 2.8: June, 2000
1980 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1981 */
1982 
1983 static const double elk_P[] =
1984 {
1985  1.37982864606273237150e-4,
1986  2.28025724005875567385e-3,
1987  7.97404013220415179367e-3,
1988  9.85821379021226008714e-3,
1989  6.87489687449949877925e-3,
1990  6.18901033637687613229e-3,
1991  8.79078273952743772254e-3,
1992  1.49380448916805252718e-2,
1993  3.08851465246711995998e-2,
1994  9.65735902811690126535e-2,
1995  1.38629436111989062502e0
1996 };
1997 
1998 static const double elk_Q[] =
1999 {
2000  2.94078955048598507511e-5,
2001  9.14184723865917226571e-4,
2002  5.94058303753167793257e-3,
2003  1.54850516649762399335e-2,
2004  2.39089602715924892727e-2,
2005  3.01204715227604046988e-2,
2006  3.73774314173823228969e-2,
2007  4.88280347570998239232e-2,
2008  7.03124996963957469739e-2,
2009  1.24999999999870820058e-1,
2010  4.99999999999999999821e-1
2011 };
2012 
2013 static const double C1 = 1.3862943611198906188e0; /* log(4) */
2014 
2015 double ellpk(double x)
2016 {
2017  DEBUG_ENTRY( "ellpk()" );
2018 
2019  if( x < 0.0 || x > 1.0 )
2020  {
2021  fprintf( ioQQQ, "ellpk: domain error\n" );
2022  cdEXIT(EXIT_FAILURE);
2023  }
2024 
2025  if( x > DBL_EPSILON )
2026  {
2027  return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
2028  }
2029  else
2030  {
2031  if( x == 0.0 )
2032  {
2033  fprintf( ioQQQ, "ellpk: domain error\n" );
2034  cdEXIT(EXIT_FAILURE);
2035  }
2036  else
2037  {
2038  return C1 - 0.5 * log(x);
2039  }
2040  }
2041 }
2042 
2043 /* expn.c
2044  *
2045  * Exponential integral En
2046  *
2047  *
2048  *
2049  * SYNOPSIS:
2050  *
2051  * int n;
2052  * double x, y, expn();
2053  *
2054  * y = expn( n, x );
2055  *
2056  *
2057  *
2058  * DESCRIPTION:
2059  *
2060  * Evaluates the exponential integral
2061  *
2062  * inf.
2063  * -
2064  * | | -xt
2065  * | e
2066  * E (x) = | ---- dt.
2067  * n | n
2068  * | | t
2069  * -
2070  * 1
2071  *
2072  *
2073  * Both n and x must be nonnegative.
2074  *
2075  * The routine employs either a power series, a continued
2076  * fraction, or an asymptotic formula depending on the
2077  * relative values of n and x.
2078  *
2079  * ACCURACY:
2080  *
2081  * Relative error:
2082  * arithmetic domain # trials peak rms
2083  * DEC 0, 30 5000 2.0e-16 4.6e-17
2084  * IEEE 0, 30 10000 1.7e-15 3.6e-16
2085  *
2086  */
2087 
2088 /* Cephes Math Library Release 2.8: June, 2000
2089  Copyright 1985, 2000 by Stephen L. Moshier */
2090 
2091 static const double MAXLOG = log(DBL_MAX);
2092 static const double BIG = 1.44115188075855872E+17; /* 2^57 */
2093 
2094 double expn(int n, double x)
2095 {
2096  double ans, r, t, yk, xk;
2097  double pk, pkm1, pkm2, qk, qkm1, qkm2;
2098  double psi, z;
2099  int i, k;
2100 
2101  DEBUG_ENTRY( "expn()" );
2102 
2103  if( n < 0 || x < 0. )
2104  {
2105  fprintf( ioQQQ, "expn: domain error\n" );
2106  cdEXIT(EXIT_FAILURE);
2107  }
2108 
2109  if( x > MAXLOG )
2110  {
2111  return 0.0;
2112  }
2113 
2114  if( x == 0.0 )
2115  {
2116  if( n < 2 )
2117  {
2118  fprintf( ioQQQ, "expn: domain error\n" );
2119  cdEXIT(EXIT_FAILURE);
2120  }
2121  else
2122  {
2123  return 1.0/((double)n-1.0);
2124  }
2125  }
2126 
2127  if( n == 0 )
2128  {
2129  return exp(-x)/x;
2130  }
2131 
2132  /* Expansion for large n */
2133  if( n > 5000 )
2134  {
2135  xk = x + n;
2136  yk = 1.0 / (xk * xk);
2137  t = n;
2138  ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
2139  ans = yk * (ans + t * (t - 2.0 * x));
2140  ans = yk * (ans + t);
2141  ans = (ans + 1.0) * exp( -x ) / xk;
2142  return ans;
2143  }
2144 
2145  if( x <= 1.0 )
2146  {
2147  /* Power series expansion */
2148  psi = -EULER - log(x);
2149  for( i=1; i < n; i++ )
2150  psi = psi + 1.0/i;
2151 
2152  z = -x;
2153  xk = 0.0;
2154  yk = 1.0;
2155  pk = 1.0 - n;
2156  if( n == 1 )
2157  ans = 0.0;
2158  else
2159  ans = 1.0/pk;
2160  do
2161  {
2162  xk += 1.0;
2163  yk *= z/xk;
2164  pk += 1.0;
2165  if( pk != 0.0 )
2166  {
2167  ans += yk/pk;
2168  }
2169  if( ans != 0.0 )
2170  t = fabs(yk/ans);
2171  else
2172  t = 1.0;
2173  }
2174  while( t > DBL_EPSILON );
2175  ans = powi(z,n-1)*psi/factorial(n-1) - ans;
2176  return ans;
2177  }
2178  else
2179  {
2180  /* continued fraction */
2181  k = 1;
2182  pkm2 = 1.0;
2183  qkm2 = x;
2184  pkm1 = 1.0;
2185  qkm1 = x + n;
2186  ans = pkm1/qkm1;
2187 
2188  do
2189  {
2190  k += 1;
2191  if( is_odd(k) )
2192  {
2193  yk = 1.0;
2194  xk = static_cast<double>(n + (k-1)/2);
2195  }
2196  else
2197  {
2198  yk = x;
2199  xk = static_cast<double>(k/2);
2200  }
2201  pk = pkm1 * yk + pkm2 * xk;
2202  qk = qkm1 * yk + qkm2 * xk;
2203  if( qk != 0 )
2204  {
2205  r = pk/qk;
2206  t = fabs( (ans - r)/r );
2207  ans = r;
2208  }
2209  else
2210  t = 1.0;
2211  pkm2 = pkm1;
2212  pkm1 = pk;
2213  qkm2 = qkm1;
2214  qkm1 = qk;
2215  if( fabs(pk) > BIG )
2216  {
2217  pkm2 /= BIG;
2218  pkm1 /= BIG;
2219  qkm2 /= BIG;
2220  qkm1 /= BIG;
2221  }
2222  }
2223  while( t > DBL_EPSILON );
2224 
2225  ans *= exp( -x );
2226  return ans;
2227  }
2228 }
2229 
2230 /* polevl.c
2231  * p1evl.c
2232  *
2233  * Evaluate polynomial
2234  *
2235  *
2236  *
2237  * SYNOPSIS:
2238  *
2239  * int N;
2240  * double x, y, coef[N+1], polevl[];
2241  *
2242  * y = polevl( x, coef, N );
2243  *
2244  *
2245  *
2246  * DESCRIPTION:
2247  *
2248  * Evaluates polynomial of degree N:
2249  *
2250  * 2 N
2251  * y = C + C x + C x +...+ C x
2252  * 0 1 2 N
2253  *
2254  * Coefficients are stored in reverse order:
2255  *
2256  * coef[0] = C , ..., coef[N] = C .
2257  * N 0
2258  *
2259  * The function p1evl() assumes that coef[N] = 1.0 and is
2260  * omitted from the array. Its calling arguments are
2261  * otherwise the same as polevl().
2262  *
2263  *
2264  * SPEED:
2265  *
2266  * In the interest of speed, there are no checks for out
2267  * of bounds arithmetic. This routine is used by most of
2268  * the functions in the library. Depending on available
2269  * equipment features, the user may wish to rewrite the
2270  * program in microcode or assembly language.
2271  *
2272  */
2273 
2274 /*
2275 Cephes Math Library Release 2.1: December, 1988
2276 Copyright 1984, 1987, 1988 by Stephen L. Moshier
2277 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2278 */
2279 
2280 inline double polevl(double x, const double coef[], int N)
2281 {
2282  double ans;
2283  int i;
2284  const double *p = coef;
2285 
2286  ans = *p++;
2287  i = N;
2288 
2289  do
2290  ans = ans * x + *p++;
2291  while( --i );
2292 
2293  return ans;
2294 }
2295 
2296 /* p1evl() */
2297 /* N
2298  * Evaluate polynomial when coefficient of x is 1.0.
2299  * Otherwise same as polevl.
2300  */
2301 
2302 inline double p1evl(double x, const double coef[], int N)
2303 {
2304  double ans;
2305  const double *p = coef;
2306  int i;
2307 
2308  ans = x + *p++;
2309  i = N-1;
2310 
2311  do
2312  ans = ans * x + *p++;
2313  while( --i );
2314 
2315  return ans;
2316 }
2317 
2318 /* chbevl.c
2319  *
2320  * Evaluate Chebyshev series
2321  *
2322  *
2323  *
2324  * SYNOPSIS:
2325  *
2326  * int N;
2327  * double x, y, coef[N], chebevl();
2328  *
2329  * y = chbevl( x, coef, N );
2330  *
2331  *
2332  *
2333  * DESCRIPTION:
2334  *
2335  * Evaluates the series
2336  *
2337  * N-1
2338  * - '
2339  * y = > coef[i] T (x/2)
2340  * - i
2341  * i=0
2342  *
2343  * of Chebyshev polynomials Ti at argument x/2.
2344  *
2345  * Coefficients are stored in reverse order, i.e. the zero
2346  * order term is last in the array. Note N is the number of
2347  * coefficients, not the order.
2348  *
2349  * If coefficients are for the interval a to b, x must
2350  * have been transformed to x -> 2(2x - b - a)/(b-a) before
2351  * entering the routine. This maps x from (a, b) to (-1, 1),
2352  * over which the Chebyshev polynomials are defined.
2353  *
2354  * If the coefficients are for the inverted interval, in
2355  * which (a, b) is mapped to (1/b, 1/a), the transformation
2356  * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
2357  * this becomes x -> 4a/x - 1.
2358  *
2359  *
2360  *
2361  * SPEED:
2362  *
2363  * Taking advantage of the recurrence properties of the
2364  * Chebyshev polynomials, the routine requires one more
2365  * addition per loop than evaluating a nested polynomial of
2366  * the same degree.
2367  *
2368  */
2369 
2370 /*
2371 Cephes Math Library Release 2.0: April, 1987
2372 Copyright 1985, 1987 by Stephen L. Moshier
2373 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2374 */
2375 
2376 inline double chbevl(double x, const double array[], int n)
2377 {
2378  double b0, b1, b2;
2379  const double *p = array;
2380  int i;
2381 
2382  b0 = *p++;
2383  b1 = 0.0;
2384  i = n - 1;
2385 
2386  do
2387  {
2388  b2 = b1;
2389  b1 = b0;
2390  b0 = x * b1 - b2 + *p++;
2391  }
2392  while( --i );
2393 
2394  return 0.5*(b0-b2);
2395 }
2396 
2397 
2398 /* >>refer Mersenne Twister Matsumoto, M. & Nishimura, T. 1998, ACM Transactions on Modeling
2399  * >>refercon and Computer Simulation (TOMACS), 8, 1998 */
2400 
2401 /********************************************************************
2402  * This copyright notice must accompany the following block of code *
2403  *******************************************************************/
2404 
2405 /*
2406  A C-program for MT19937, with initialization improved 2002/2/10.
2407  Coded by Takuji Nishimura and Makoto Matsumoto.
2408  This is a faster version by taking Shawn Cokus's optimization,
2409  Matthe Bellew's simplification, Isaku Wada's real version.
2410 
2411  Before using, initialize the state by using init_genrand(seed)
2412  or init_by_array(init_key, key_length).
2413 
2414  Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
2415  All rights reserved.
2416 
2417  Redistribution and use in source and binary forms, with or without
2418  modification, are permitted provided that the following conditions
2419  are met:
2420 
2421  1. Redistributions of source code must retain the above copyright
2422  notice, this list of conditions and the following disclaimer.
2423 
2424  2. Redistributions in binary form must reproduce the above copyright
2425  notice, this list of conditions and the following disclaimer in the
2426  documentation and/or other materials provided with the distribution.
2427 
2428  3. The names of its contributors may not be used to endorse or promote
2429  products derived from this software without specific prior written
2430  permission.
2431 
2432  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
2433  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
2434  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
2435  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
2436  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
2437  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
2438  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
2439  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
2440  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
2441  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
2442  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
2443 
2444 
2445  Any feedback is very welcome.
2446  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
2447  email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
2448 */
2449 
2450 /* Period parameters */
2451 #define N 624
2452 #define M 397
2453 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
2454 #define UMASK 0x80000000UL /* most significant w-r bits */
2455 #define LMASK 0x7fffffffUL /* least significant r bits */
2456 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
2457 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
2458 
2459 static unsigned long state[N]; /* the array for the state vector */
2460 static int nleft = 1;
2461 static int initf = 0;
2462 static unsigned long *next;
2463 
2464 /* initializes state[N] with a seed */
2465 void init_genrand(unsigned long s)
2466 {
2467  int j;
2468  state[0]= s & 0xffffffffUL;
2469  for (j=1; j<N; j++) {
2470  state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
2471  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
2472  /* In the previous versions, MSBs of the seed affect */
2473  /* only MSBs of the array state[]. */
2474  /* 2002/01/09 modified by Makoto Matsumoto */
2475  state[j] &= 0xffffffffUL; /* for >32 bit machines */
2476  }
2477  nleft = 1; initf = 1;
2478 }
2479 
2480 /* initialize by an array with array-length */
2481 /* init_key is the array for initializing keys */
2482 /* key_length is its length */
2483 /* slight change for C++, 2004/2/26 */
2484 void init_by_array(unsigned long init_key[], int key_length)
2485 {
2486  int i, j, k;
2487  init_genrand(19650218UL);
2488  i=1; j=0;
2489  k = (N>key_length ? N : key_length);
2490  for (; k; k--) {
2491  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
2492  + init_key[j] + j; /* non linear */
2493  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
2494  i++; j++;
2495  if (i>=N) { state[0] = state[N-1]; i=1; }
2496  if (j>=key_length) j=0;
2497  }
2498  for (k=N-1; k; k--) {
2499  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
2500  - i; /* non linear */
2501  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
2502  i++;
2503  if (i>=N) { state[0] = state[N-1]; i=1; }
2504  }
2505 
2506  state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
2507  nleft = 1; initf = 1;
2508 }
2509 
2510 static void next_state(void)
2511 {
2512  unsigned long *p=state;
2513  int j;
2514 
2515  /* if init_genrand() has not been called, */
2516  /* a default initial seed is used */
2517  if (initf==0) init_genrand(5489UL);
2518 
2519  nleft = N;
2520  next = state;
2521 
2522  for (j=N-M+1; --j; p++)
2523  *p = p[M] ^ TWIST(p[0], p[1]);
2524 
2525  for (j=M; --j; p++)
2526  *p = p[M-N] ^ TWIST(p[0], p[1]);
2527 
2528  *p = p[M-N] ^ TWIST(p[0], state[0]);
2529 }
2530 
2531 /* generates a random number on [0,0xffffffff]-interval */
2532 unsigned long genrand_int32(void)
2533 {
2534  unsigned long y;
2535 
2536  if (--nleft == 0) next_state();
2537  y = *next++;
2538 
2539  /* Tempering */
2540  y ^= (y >> 11);
2541  y ^= (y << 7) & 0x9d2c5680UL;
2542  y ^= (y << 15) & 0xefc60000UL;
2543  y ^= (y >> 18);
2544 
2545  return y;
2546 }
2547 
2548 /* generates a random number on [0,0x7fffffff]-interval */
2549 long genrand_int31(void)
2550 {
2551  unsigned long y;
2552 
2553  if (--nleft == 0) next_state();
2554  y = *next++;
2555 
2556  /* Tempering */
2557  y ^= (y >> 11);
2558  y ^= (y << 7) & 0x9d2c5680UL;
2559  y ^= (y << 15) & 0xefc60000UL;
2560  y ^= (y >> 18);
2561 
2562  return (long)(y>>1);
2563 }
2564 
2565 /* generates a random number on [0,1]-real-interval */
2566 double genrand_real1(void)
2567 {
2568  unsigned long y;
2569 
2570  if (--nleft == 0) next_state();
2571  y = *next++;
2572 
2573  /* Tempering */
2574  y ^= (y >> 11);
2575  y ^= (y << 7) & 0x9d2c5680UL;
2576  y ^= (y << 15) & 0xefc60000UL;
2577  y ^= (y >> 18);
2578 
2579  return (double)y * (1.0/4294967295.0);
2580  /* divided by 2^32-1 */
2581 }
2582 
2583 /* generates a random number on [0,1)-real-interval */
2584 double genrand_real2(void)
2585 {
2586  unsigned long y;
2587 
2588  if (--nleft == 0) next_state();
2589  y = *next++;
2590 
2591  /* Tempering */
2592  y ^= (y >> 11);
2593  y ^= (y << 7) & 0x9d2c5680UL;
2594  y ^= (y << 15) & 0xefc60000UL;
2595  y ^= (y >> 18);
2596 
2597  return (double)y * (1.0/4294967296.0);
2598  /* divided by 2^32 */
2599 }
2600 
2601 /* generates a random number on (0,1)-real-interval */
2602 double genrand_real3(void)
2603 {
2604  unsigned long y;
2605 
2606  if (--nleft == 0) next_state();
2607  y = *next++;
2608 
2609  /* Tempering */
2610  y ^= (y >> 11);
2611  y ^= (y << 7) & 0x9d2c5680UL;
2612  y ^= (y << 15) & 0xefc60000UL;
2613  y ^= (y >> 18);
2614 
2615  return ((double)y + 0.5) * (1.0/4294967296.0);
2616  /* divided by 2^32 */
2617 }
2618 
2619 /* generates a random number on [0,1) with 53-bit resolution*/
2620 double genrand_res53(void)
2621 {
2622  unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
2623  return(a*67108864.0+b)*(1.0/9007199254740992.0);
2624 }
2625 
2626 /* These real versions are due to Isaku Wada, 2002/01/09 added */
2627 
2628 /************************************************************************
2629  * This marks the end of the block of code from Matsumoto and Nishimura *
2630  ************************************************************************/

Generated for cloudy by doxygen 1.8.1.2