cloudy  trunk
 All Data Structures Namespaces 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 "vectorize.h"
5 #include "thirdparty.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 inline 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 
46 bool 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 
181 static 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 
356 double 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 static const double dsf[MXDSF] =
370 {
371  1.00000000000000000e+000,
372  6.25000000000000000e-002,
373  7.81250000000000000e-003,
374  1.46484375000000000e-003,
375  3.66210937500000000e-004,
376  1.14440917968750000e-004,
377  4.29153442382812500e-005,
378  1.87754631042480469e-005,
379  9.38773155212402344e-006,
380  5.28059899806976318e-006,
381  3.30037437379360199e-006,
382  2.26900738198310137e-006,
383  1.70175553648732603e-006,
384  1.38267637339595240e-006,
385  1.20984182672145835e-006,
386  1.13422671255136720e-006,
387  1.13422671255136720e-006,
388  1.20511588208582765e-006,
389  1.35575536734655611e-006,
390  1.60995949872403538e-006,
391  2.01244937340504422e-006,
392  2.64133980259412054e-006,
393  3.63184222856691574e-006,
394  5.22077320356494170e-006,
395  7.83115980534741170e-006,
396  1.22361871958553314e-005,
397  1.98838041932649142e-005,
398  3.35539195761345408e-005,
399  5.87193592582354430e-005,
400  1.06428838655551734e-004,
401  1.99554072479159509e-004,
402  3.86636015428371569e-004,
403  7.73272030856743137e-004,
404  1.59487356364203269e-003,
405  3.38910632273931945e-003,
406  7.41367008099226132e-003,
407  1.66807576822325873e-002,
408  3.85742521401628569e-002,
409  9.16138488328867850e-002,
410  0.22330875653016155,
411  0.55827189132540389,
412  1.4305717215213474,
413  3.7552507689935370,
414  10.092236441670131,
415  27.753650214592859,
416  78.057141228542420,
417  224.41428103205945,
418  659.21695053167468,
419  1977.6508515950241,
420  6056.5557330097608,
421  18926.736665655502,
422  60328.973121776915,
423  196069.16264577498,
424  649479.10126412963,
425  2191991.9667664375,
426  7534972.3857596293,
427  26372403.350158703,
428  93951686.934940383,
429  340574865.13915890,
430  1255869815.2006485,
431  4709511807.0024319,
432  17955013764.196770,
433  69575678336.262482,
434  273954233449.03351,
435  1095816933796.1340,
436  4451756293546.7949,
437  18363494710880.527,
438  76897134101812.203,
439  326812819932701.88,
440  1409380285959776.8,
441  6166038751074023.0,
442  27361796957890976.,
443  1.23128086310509392e+017,
444  5.61771893791699072e+017,
445  2.59819500878660813e+018,
446  1.21790391036872253e+019,
447  5.78504357425143235e+019,
448  2.78405222010850181e+020,
449  1.35722545730289454e+021,
450  6.70130069543304194e+021,
451  3.35065034771652076e+022,
452  1.69626673853148848e+023,
453  8.69336703497387857e+023,
454  4.50968414939269935e+024,
455  2.36758417843116729e+025,
456  1.25777909479155770e+026,
457  6.76056263450462257e+026,
458  3.67605593251188834e+027,
459  2.02183076288153845e+028,
460  1.12464336185285569e+029,
461  6.32611891042231308e+029,
462  3.59798013030269071e+030,
463  2.06883857492404700e+031,
464  1.20251242167460229e+032,
465  7.06476047733828837e+032,
466  4.19470153341960844e+033,
467  2.51682092005176529e+034,
468  1.52582268278138270e+035,
469  9.34566393203596895e+035,
470  5.78262955794725579e+036,
471  3.61414347371703511e+037,
472  2.28142806778387824e+038,
473  1.45441039321222244e+039,
474  9.36276690630368255e+039,
475  6.08579848909739317e+040,
476  3.99380525847016458e+041,
477  2.64589598373648418e+042,
478  1.76944293912377371e+043,
479  1.19437398390854730e+044,
480  8.13667276537697798e+044,
481  5.59396252619667246e+045,
482  3.88081150254894130e+046,
483  2.71656805178425881e+047,
484  1.91857618657263275e+048,
485  1.36698553293300096e+049,
486  9.82520851795594440e+049,
487  7.12327617551805990e+050,
488  5.20889570334758147e+051,
489  3.84156058121884122e+052,
490  2.85716068228151327e+053,
491  2.14287051171113493e+054,
492  1.62054582448154580e+055,
493  1.23566619116717877e+056,
494  9.49918384459768667e+056,
495  7.36186747956320756e+057,
496  5.75145896840875570e+058,
497  4.52927393762189499e+059,
498  3.59511118798737905e+060,
499  2.87608895038990324e+061,
500  2.31884671625185959e+062,
501  1.88406295695463606e+063,
502  1.54257654600660820e+064,
503  1.27262565045545178e+065,
504  1.05787007194109428e+066,
505  8.85966185250666444e+066,
506  7.47533968805249835e+067,
507  6.35403873484462330e+068,
508  5.44064566671070888e+069,
509  4.69255688753798633e+070,
510  4.07665879604862585e+071,
511  3.56707644654254751e+072,
512  3.14348611851561972e+073,
513  2.78984393018261231e+074,
514  2.49342301260070962e+075,
515  2.24408071134063849e+076,
516  2.03369814465245358e+077,
517  1.85574955699536371e+078,
518  1.70496990548949040e+079,
519  1.57709716257777866e+080,
520  1.46867173265055644e+081,
521  1.37687974935989664e+082,
522  1.29943026345840251e+083,
523  1.23445875028548248e+084,
524  1.18045117996049257e+085,
525  1.13618426071197409e+086,
526  1.10067850256472499e+087,
527  1.07316154000060691e+088,
528  1.05303976112559556e+089,
529  1.03987676411152556e+090,
530  1.03337753433582857e+091,
531  1.03337753433582857e+092,
532  1.03983614392542748e+093,
533  1.05283409572449530e+094,
534  1.07257473501932960e+095,
535  1.09938910339481284e+096,
536  1.13374501287590075e+097,
537  1.17626045085874706e+098,
538  1.22772184558381721e+099,
539  1.28910793786300799e+100,
540  1.36162025936780216e+101,
541  1.44672152557828969e+102,
542  1.54618363046179715e+103,
543  1.66214740274643181e+104,
544  1.79719687921957945e+105,
545  1.95445160615129251e+106,
546  2.13768144422797627e+107,
547  2.35144958865077380e+108,
548  2.60129110744491837e+109,
549  2.89393635703247155e+110,
550  3.23759129943007744e+111,
551  3.64229021185883726e+112,
552  4.12034080216530982e+113,
553  4.68688766246304013e+114,
554  5.36062776394210187e+115,
555  6.16472192853341706e+116,
556  7.12795972986676387e+117,
557  8.28625318597011219e+118,
558  9.68455841110256888e+119,
559  1.13793561330455178e+121,
560  1.34418644321600176e+122,
561  1.59622140131900197e+123,
562  1.90548929782455860e+124,
563  2.28658715738947042e+125,
564  2.75819575860104856e+126,
565  3.34431235730377152e+127,
566  4.07588068546397120e+128,
567  4.99295383969336454e+129,
568  6.14757441512245555e+130,
569  7.60762333871403831e+131,
570  9.46198152752558548e+132,
571  1.18274769094069828e+134,
572  1.48582678674425232e+135,
573  1.87585631826461846e+136,
574  2.37999270379823447e+137,
575  3.03449069734274889e+138,
576  3.88794120597039710e+139,
577  5.00572430268688588e+140,
578  6.47615581660115797e+141,
579  8.41900256158150577e+142,
580  1.09973220960658419e+144,
581  1.44339852510864179e+145,
582  1.90348180498702145e+146,
583  2.52211339160780346e+147,
584  3.35756345257788818e+148,
585  4.49074111782292584e+149,
586  6.03443337707455679e+150,
587  8.14648505905065221e+151,
588  1.10486703613374480e+153,
589  1.50538133673222729e+154,
590  2.06049070465223602e+155,
591  2.83317471889682437e+156,
592  3.91332258047623873e+157,
593  5.42973508041078113e+158,
594  7.56769326832252665e+159,
595  1.05947705756515379e+161,
596  1.48988961220099760e+162,
597  2.10446907723390916e+163,
598  2.98571550332560859e+164,
599  4.25464459223899213e+165,
600  6.08946007264205694e+166,
601  8.75359885442295723e+167,
602  1.26380083460731439e+169,
603  1.83251121018060591e+170,
604  2.66859444982550725e+171,
605  3.90281938286980469e+172,
606  5.73226596859002516e+173,
607  8.45509230367028708e+174,
608  1.25241054748116123e+176,
609  1.86296068937822735e+177,
610  2.78279752975872722e+178,
611  4.17419629463809079e+179,
612  6.28738316879862387e+180,
613  9.50966704280791803e+181,
614  1.44428068212645250e+183,
615  2.20252804024284002e+184,
616  3.37262106162184907e+185,
617  5.18540488224359247e+186,
618  8.00496878696354631e+187,
619  1.24077016197934964e+189,
620  1.93094856458036277e+190,
621  3.01710713215681661e+191,
622  4.73308681357100637e+192,
623  7.45461173137433533e+193,
624  1.17876048002356679e+195,
625  1.87128226203741209e+196,
626  2.98235610512212568e+197,
627  4.77176976819540109e+198,
628  7.66465519016386263e+199,
629  1.23592564941392296e+201,
630  2.00065464498878787e+202,
631  3.25106379810678013e+203,
632  5.30329782066168502e+204,
633  8.68415018133350951e+205,
634  1.42745718605669573e+207,
635  2.35530435699354811e+208,
636  3.90097284127056435e+209,
637  6.48536734861231298e+210,
638  1.08224567629967973e+212,
639  1.81276150780196369e+213,
640  3.04770528499205135e+214,
641  5.14300266842408688e+215,
642  8.71096076964329754e+216,
643  1.48086333083936058e+218,
644  2.52672305824465885e+219,
645  4.32701323724397813e+220,
646  7.43705400151308679e+221,
647  1.28289181526100749e+223,
648  2.22100645517061936e+224,
649  3.85899871585895109e+225,
650  6.72912901077904602e+226,
651  1.17759757688633299e+228,
652  2.06815574440662224e+229,
653  3.64512449951667173e+230,
654  6.44731395852011310e+231,
655  1.14439822763732000e+233,
656  2.03845934297897628e+234,
657  3.64374607557492029e+235,
658  6.53596952306251304e+236,
659  1.17647451415125240e+238,
660  2.12500709118569980e+239,
661  3.85157535277408106e+240,
662  7.00505267285786018e+241,
663  1.27842211279655955e+243,
664  2.34111049405869977e+244,
665  4.30179053283286084e+245,
666  7.93142629491058743e+246,
667  1.46731386455845874e+248,
668  2.72370136108663892e+249,
669  5.07289378502386487e+250,
670  9.47997026076334780e+251,
671  1.77749442389312786e+253,
672  3.34391138494894698e+254,
673  6.31163273909113711e+255,
674  1.19526544996538408e+257,
675  2.27100435493422958e+258,
676  4.32910205159337539e+259,
677  8.27940767367233019e+260,
678  1.58861134738587850e+262,
679  3.05807684371781629e+263,
680  5.90591090443003272e+264,
681  1.14427023773331886e+266,
682  2.22417527459413850e+267,
683  4.33714178545856976e+268,
684  8.48453361780332716e+269,
685  1.66508972249390287e+271,
686  3.27814539115987142e+272,
687  6.47433714754074587e+273,
688  1.28272804735651034e+275,
689  2.54942199412106407e+276,
690  5.08291010077887125e+277,
691  1.01658202015577428e+279,
692  2.03951767793752200e+280,
693  4.10452932684926289e+281,
694  8.28601857857694926e+282,
695  1.67791876216183238e+284,
696  3.40827248564122204e+285,
697  6.94435518949399006e+286,
698  1.41925259185283426e+288,
699  2.90946781329831038e+289,
700  5.98259319109465094e+290,
701  1.23390984566327173e+292,
702  2.55265099321589344e+293,
703  5.29675081092297942e+294,
704  1.10238626252334501e+296,
705  2.30123132301748272e+297,
706  4.81820308256785444e+298,
707  1.01182264733924947e+300,
708  2.13115145095829414e+301,
709  4.50205744014939638e+302,
710  9.53873420131653322e+303
711 };
712 
713 // helper function for sjs calculating a scaled factorial
714 // fc2_scl(2*n+1) = (2*n)!! / 32^n = n! / 16^n
715 inline double fc2_scl( long n )
716 {
717  if( n < 1 || n >= 2*MXDSF+1 || (n%2) != 1 )
718  {
719  fprintf( ioQQQ, "fc2_scl: invalid argument\n" );
721  }
722  n >>= 1; // calculate (n-1)/2
723  return dsf[n];
724 }
725 
726 /* NB - this implementation is not thread-safe! */
727 
728 class t_lfact : public Singleton<t_lfact>
729 {
730  friend class Singleton<t_lfact>;
731 protected:
733  {
734  p_lf.reserve( 512 );
735  p_lf.push_back( 0. ); /* log10( 0! ) */
736  p_lf.push_back( 0. ); /* log10( 1! ) */
737  }
738 private:
739  vector<double> p_lf;
740 public:
741  double get_lfact( unsigned long n )
742  {
743  if( n < p_lf.size() )
744  {
745  return p_lf[n];
746  }
747  else
748  {
749  for( unsigned long i=(unsigned long)p_lf.size(); i <= n; i++ )
750  p_lf.push_back( p_lf[i-1] + log10(static_cast<double>(i)) );
751  return p_lf[n];
752  }
753  }
754 };
755 
756 double lfactorial( long n )
757 {
758  /******************************/
759  /* n */
760  /* --- */
761  /* log( n! ) = > log(j) */
762  /* --- */
763  /* j=1 */
764  /******************************/
765 
766  DEBUG_ENTRY( "lfactorial()" );
767 
768  if( n < 0 )
769  {
770  fprintf( ioQQQ, "lfactorial: domain error\n" );
772  }
773 
774  return t_lfact::Inst().get_lfact( static_cast<unsigned long>(n) );
775 }
776 
777 /*******************************************************************
778  * This marks the end of the block of code from Robert Paul Bauman *
779  *******************************************************************/
780 
781 
782 /* complex Gamma function in double precision */
783 /* this routine is a slightly modified version of the one found
784  * at http://momonga.t.u-tokyo.ac.jp/~ooura/gamerf.html
785  * The following copyright applies:
786  * Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
787  * You may use, copy, modify this code for any purpose and
788  * without fee. You may distribute this ORIGINAL package. */
789 complex<double> cdgamma(complex<double> x)
790 {
791  double xr, xi, wr, wi, ur, ui, vr, vi, yr, yi, t;
792 
793  DEBUG_ENTRY( "cdgamma()" );
794 
795  xr = x.real();
796  xi = x.imag();
797  if(xr < 0)
798  {
799  wr = 1. - xr;
800  wi = -xi;
801  }
802  else
803  {
804  wr = xr;
805  wi = xi;
806  }
807  ur = wr + 6.00009857740312429;
808  vr = ur * (wr + 4.99999857982434025) - wi * wi;
809  vi = wi * (wr + 4.99999857982434025) + ur * wi;
810  yr = ur * 13.2280130755055088 + vr * 66.2756400966213521 +
811  0.293729529320536228;
812  yi = wi * 13.2280130755055088 + vi * 66.2756400966213521;
813  ur = vr * (wr + 4.00000003016801681) - vi * wi;
814  ui = vi * (wr + 4.00000003016801681) + vr * wi;
815  vr = ur * (wr + 2.99999999944915534) - ui * wi;
816  vi = ui * (wr + 2.99999999944915534) + ur * wi;
817  yr += ur * 91.1395751189899762 + vr * 47.3821439163096063;
818  yi += ui * 91.1395751189899762 + vi * 47.3821439163096063;
819  ur = vr * (wr + 2.00000000000603851) - vi * wi;
820  ui = vi * (wr + 2.00000000000603851) + vr * wi;
821  vr = ur * (wr + 0.999999999999975753) - ui * wi;
822  vi = ui * (wr + 0.999999999999975753) + ur * wi;
823  yr += ur * 10.5400280458730808 + vr;
824  yi += ui * 10.5400280458730808 + vi;
825  ur = vr * wr - vi * wi;
826  ui = vi * wr + vr * wi;
827  t = ur * ur + ui * ui;
828  vr = yr * ur + yi * ui + t * 0.0327673720261526849;
829  vi = yi * ur - yr * ui;
830  yr = wr + 7.31790632447016203;
831  ur = log(yr * yr + wi * wi) * 0.5 - 1;
832  ui = atan2(wi, yr);
833  yr = exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t;
834  yi = ui * (wr - 0.5) + ur * wi;
835  ur = yr * cos(yi);
836  ui = yr * sin(yi);
837  yr = ur * vr - ui * vi;
838  yi = ui * vr + ur * vi;
839  if(xr < 0)
840  {
841  wr = xr * 3.14159265358979324;
842  wi = exp(xi * 3.14159265358979324);
843  vi = 1 / wi;
844  ur = (vi + wi) * sin(wr);
845  ui = (vi - wi) * cos(wr);
846  vr = ur * yr + ui * yi;
847  vi = ui * yr - ur * yi;
848  ur = 6.2831853071795862 / (vr * vr + vi * vi);
849  yr = ur * vr;
850  yi = ur * vi;
851  }
852  return complex<double>( yr, yi );
853 }
854 
855 /*************************************************************
856  * This marks the end of the block of code from Takuya OOURA *
857  *************************************************************/
858 
859 /*====================================================================
860  *
861  * Below are routines from the Cephes library.
862  *
863  * This is the copyright statement included with the library:
864  *
865  * Some software in this archive may be from the book _Methods and
866  * Programs for Mathematical Functions_ (Prentice-Hall, 1989) or
867  * from the Cephes Mathematical Library, a commercial product. In
868  * either event, it is copyrighted by the author. What you see here
869  * may be used freely but it comes with no support or guarantee.
870  *
871  * The two known misprints in the book are repaired here in the
872  * source listings for the gamma function and the incomplete beta
873  * integral.
874  *
875  * Stephen L. Moshier
876  * moshier@world.std.com
877  *
878  *====================================================================*/
879 
880 /* j0.c
881  *
882  * Bessel function of order zero
883  *
884  *
885  *
886  * SYNOPSIS:
887  *
888  * double x, y, j0();
889  *
890  * y = j0( x );
891  *
892  *
893  *
894  * DESCRIPTION:
895  *
896  * Returns Bessel function of order zero of the argument.
897  *
898  * The domain is divided into the intervals [0, 5] and
899  * (5, infinity). In the first interval the following rational
900  * approximation is used:
901  *
902  *
903  * 2 2
904  * (w - r ) (w - r ) P (w) / Q (w)
905  * 1 2 3 8
906  *
907  * 2
908  * where w = x and the two r's are zeros of the function.
909  *
910  * In the second interval, the Hankel asymptotic expansion
911  * is employed with two rational functions of degree 6/6
912  * and 7/7.
913  *
914  *
915  *
916  * ACCURACY:
917  *
918  * Absolute error:
919  * arithmetic domain # trials peak rms
920  * DEC 0, 30 10000 4.4e-17 6.3e-18
921  * IEEE 0, 30 60000 4.2e-16 1.1e-16
922  *
923  */
924 /* y0.c
925  *
926  * Bessel function of the second kind, order zero
927  *
928  *
929  *
930  * SYNOPSIS:
931  *
932  * double x, y, y0();
933  *
934  * y = y0( x );
935  *
936  *
937  *
938  * DESCRIPTION:
939  *
940  * Returns Bessel function of the second kind, of order
941  * zero, of the argument.
942  *
943  * The domain is divided into the intervals [0, 5] and
944  * (5, infinity). In the first interval a rational approximation
945  * R(x) is employed to compute
946  * y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
947  * Thus a call to j0() is required.
948  *
949  * In the second interval, the Hankel asymptotic expansion
950  * is employed with two rational functions of degree 6/6
951  * and 7/7.
952  *
953  *
954  *
955  * ACCURACY:
956  *
957  * Absolute error, when y0(x) < 1; else relative error:
958  *
959  * arithmetic domain # trials peak rms
960  * DEC 0, 30 9400 7.0e-17 7.9e-18
961  * IEEE 0, 30 30000 1.3e-15 1.6e-16
962  *
963  */
964 
965 /*
966 Cephes Math Library Release 2.8: June, 2000
967 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
968 */
969 
970 /* Note: all coefficients satisfy the relative error criterion
971  * except YP, YQ which are designed for absolute error. */
972 
973 static const double b0_PP[7] = {
974  7.96936729297347051624e-4,
975  8.28352392107440799803e-2,
976  1.23953371646414299388e0,
977  5.44725003058768775090e0,
978  8.74716500199817011941e0,
979  5.30324038235394892183e0,
980  9.99999999999999997821e-1,
981 };
982 
983 static const double b0_PQ[7] = {
984  9.24408810558863637013e-4,
985  8.56288474354474431428e-2,
986  1.25352743901058953537e0,
987  5.47097740330417105182e0,
988  8.76190883237069594232e0,
989  5.30605288235394617618e0,
990  1.00000000000000000218e0,
991 };
992 
993 static const double b0_QP[8] = {
994  -1.13663838898469149931e-2,
995  -1.28252718670509318512e0,
996  -1.95539544257735972385e1,
997  -9.32060152123768231369e1,
998  -1.77681167980488050595e2,
999  -1.47077505154951170175e2,
1000  -5.14105326766599330220e1,
1001  -6.05014350600728481186e0,
1002 };
1003 
1004 static const double b0_QQ[7] = {
1005  /* 1.00000000000000000000e0,*/
1006  6.43178256118178023184e1,
1007  8.56430025976980587198e2,
1008  3.88240183605401609683e3,
1009  7.24046774195652478189e3,
1010  5.93072701187316984827e3,
1011  2.06209331660327847417e3,
1012  2.42005740240291393179e2,
1013 };
1014 
1015 static const double b0_YP[8] = {
1016  1.55924367855235737965e4,
1017  -1.46639295903971606143e7,
1018  5.43526477051876500413e9,
1019  -9.82136065717911466409e11,
1020  8.75906394395366999549e13,
1021  -3.46628303384729719441e15,
1022  4.42733268572569800351e16,
1023  -1.84950800436986690637e16,
1024 };
1025 
1026 static const double b0_YQ[7] = {
1027  /* 1.00000000000000000000e0,*/
1028  1.04128353664259848412e3,
1029  6.26107330137134956842e5,
1030  2.68919633393814121987e8,
1031  8.64002487103935000337e10,
1032  2.02979612750105546709e13,
1033  3.17157752842975028269e15,
1034  2.50596256172653059228e17,
1035 };
1036 
1037 /* 5.783185962946784521175995758455807035071 */
1038 static const double DR1 = 5.78318596294678452118e0;
1039 /* 30.47126234366208639907816317502275584842 */
1040 static const double DR2 = 3.04712623436620863991e1;
1041 
1042 static double b0_RP[4] = {
1043  -4.79443220978201773821e9,
1044  1.95617491946556577543e12,
1045  -2.49248344360967716204e14,
1046  9.70862251047306323952e15,
1047 };
1048 
1049 static double b0_RQ[8] = {
1050  /* 1.00000000000000000000e0,*/
1051  4.99563147152651017219e2,
1052  1.73785401676374683123e5,
1053  4.84409658339962045305e7,
1054  1.11855537045356834862e10,
1055  2.11277520115489217587e12,
1056  3.10518229857422583814e14,
1057  3.18121955943204943306e16,
1058  1.71086294081043136091e18,
1059 };
1060 
1061 static const double TWOOPI = 2./PI;
1062 static const double SQ2OPI = sqrt(2./PI);
1063 static const double PIO4 = PI/4.;
1064 
1065 double bessel_j0(double x)
1066 {
1067  double w, z, p, q, xn;
1068 
1069  DEBUG_ENTRY( "bessel_j0()" );
1070 
1071  if( x < 0 )
1072  x = -x;
1073 
1074  if( x <= 5.0 )
1075  {
1076  z = x * x;
1077  if( x < 1.0e-5 )
1078  return 1.0 - z/4.0;
1079 
1080  p = (z - DR1) * (z - DR2);
1081  p = p * polevl( z, b0_RP, 3)/p1evl( z, b0_RQ, 8 );
1082  return p;
1083  }
1084 
1085  w = 5.0/x;
1086  q = 25.0/(x*x);
1087  p = polevl( q, b0_PP, 6)/polevl( q, b0_PQ, 6 );
1088  q = polevl( q, b0_QP, 7)/p1evl( q, b0_QQ, 7 );
1089  xn = x - PIO4;
1090  p = p * cos(xn) - w * q * sin(xn);
1091  return p * SQ2OPI / sqrt(x);
1092 }
1093 
1094 /* y0() 2 */
1095 /* Bessel function of second kind, order zero */
1096 
1097 /* Rational approximation coefficients YP[], YQ[] are used here.
1098  * The function computed is y0(x) - 2 * log(x) * j0(x) / PI,
1099  * whose value at x = 0 is 2 * ( log(0.5) + EULER ) / PI
1100  * = 0.073804295108687225.
1101  */
1102 
1103 double bessel_y0(double x)
1104 {
1105  double w, z, p, q, xn;
1106 
1107  DEBUG_ENTRY( "bessel_y0()" );
1108 
1109  if( x <= 5.0 )
1110  {
1111  if( x <= 0.0 )
1112  {
1113  fprintf( ioQQQ, "bessel_y0: domain error\n" );
1115  }
1116  z = x * x;
1117  w = polevl( z, b0_YP, 7 ) / p1evl( z, b0_YQ, 7 );
1118  w += TWOOPI * log(x) * bessel_j0(x);
1119  return w;
1120  }
1121 
1122  w = 5.0/x;
1123  z = 25.0 / (x * x);
1124  p = polevl( z, b0_PP, 6)/polevl( z, b0_PQ, 6 );
1125  q = polevl( z, b0_QP, 7)/p1evl( z, b0_QQ, 7 );
1126  xn = x - PIO4;
1127  p = p * sin(xn) + w * q * cos(xn);
1128  return p * SQ2OPI / sqrt(x);
1129 }
1130 
1131 /* j1.c
1132  *
1133  * Bessel function of order one
1134  *
1135  *
1136  *
1137  * SYNOPSIS:
1138  *
1139  * double x, y, j1();
1140  *
1141  * y = j1( x );
1142  *
1143  *
1144  *
1145  * DESCRIPTION:
1146  *
1147  * Returns Bessel function of order one of the argument.
1148  *
1149  * The domain is divided into the intervals [0, 8] and
1150  * (8, infinity). In the first interval a 24 term Chebyshev
1151  * expansion is used. In the second, the asymptotic
1152  * trigonometric representation is employed using two
1153  * rational functions of degree 5/5.
1154  *
1155  *
1156  *
1157  * ACCURACY:
1158  *
1159  * Absolute error:
1160  * arithmetic domain # trials peak rms
1161  * DEC 0, 30 10000 4.0e-17 1.1e-17
1162  * IEEE 0, 30 30000 2.6e-16 1.1e-16
1163  *
1164  *
1165  */
1166 /* y1.c
1167  *
1168  * Bessel function of second kind of order one
1169  *
1170  *
1171  *
1172  * SYNOPSIS:
1173  *
1174  * double x, y, y1();
1175  *
1176  * y = y1( x );
1177  *
1178  *
1179  *
1180  * DESCRIPTION:
1181  *
1182  * Returns Bessel function of the second kind of order one
1183  * of the argument.
1184  *
1185  * The domain is divided into the intervals [0, 8] and
1186  * (8, infinity). In the first interval a 25 term Chebyshev
1187  * expansion is used, and a call to j1() is required.
1188  * In the second, the asymptotic trigonometric representation
1189  * is employed using two rational functions of degree 5/5.
1190  *
1191  *
1192  *
1193  * ACCURACY:
1194  *
1195  * Absolute error:
1196  * arithmetic domain # trials peak rms
1197  * DEC 0, 30 10000 8.6e-17 1.3e-17
1198  * IEEE 0, 30 30000 1.0e-15 1.3e-16
1199  *
1200  * (error criterion relative when |y1| > 1).
1201  *
1202  */
1203 
1204 /*
1205 Cephes Math Library Release 2.8: June, 2000
1206 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
1207 */
1208 
1209 static const double b1_RP[4] = {
1210  -8.99971225705559398224e8,
1211  4.52228297998194034323e11,
1212  -7.27494245221818276015e13,
1213  3.68295732863852883286e15,
1214 };
1215 
1216 static const double b1_RQ[8] = {
1217  /* 1.00000000000000000000E0,*/
1218  6.20836478118054335476e2,
1219  2.56987256757748830383e5,
1220  8.35146791431949253037e7,
1221  2.21511595479792499675e10,
1222  4.74914122079991414898e12,
1223  7.84369607876235854894e14,
1224  8.95222336184627338078e16,
1225  5.32278620332680085395e18,
1226 };
1227 
1228 static const double b1_PP[7] = {
1229  7.62125616208173112003e-4,
1230  7.31397056940917570436e-2,
1231  1.12719608129684925192e0,
1232  5.11207951146807644818e0,
1233  8.42404590141772420927e0,
1234  5.21451598682361504063e0,
1235  1.00000000000000000254e0,
1236 };
1237 
1238 static const double b1_PQ[7] = {
1239  5.71323128072548699714e-4,
1240  6.88455908754495404082e-2,
1241  1.10514232634061696926e0,
1242  5.07386386128601488557e0,
1243  8.39985554327604159757e0,
1244  5.20982848682361821619e0,
1245  9.99999999999999997461e-1,
1246 };
1247 
1248 static const double b1_QP[8] = {
1249  5.10862594750176621635e-2,
1250  4.98213872951233449420e0,
1251  7.58238284132545283818e1,
1252  3.66779609360150777800e2,
1253  7.10856304998926107277e2,
1254  5.97489612400613639965e2,
1255  2.11688757100572135698e2,
1256  2.52070205858023719784e1,
1257 };
1258 
1259 static const double b1_QQ[7] = {
1260  /* 1.00000000000000000000e0,*/
1261  7.42373277035675149943e1,
1262  1.05644886038262816351e3,
1263  4.98641058337653607651e3,
1264  9.56231892404756170795e3,
1265  7.99704160447350683650e3,
1266  2.82619278517639096600e3,
1267  3.36093607810698293419e2,
1268 };
1269 
1270 static const double b1_YP[6] = {
1271  1.26320474790178026440e9,
1272  -6.47355876379160291031e11,
1273  1.14509511541823727583e14,
1274  -8.12770255501325109621e15,
1275  2.02439475713594898196e17,
1276  -7.78877196265950026825e17,
1277 };
1278 
1279 static const double b1_YQ[8] = {
1280  /* 1.00000000000000000000E0,*/
1281  5.94301592346128195359E2,
1282  2.35564092943068577943E5,
1283  7.34811944459721705660E7,
1284  1.87601316108706159478E10,
1285  3.88231277496238566008E12,
1286  6.20557727146953693363E14,
1287  6.87141087355300489866E16,
1288  3.97270608116560655612E18,
1289 };
1290 
1291 static const double Z1 = 1.46819706421238932572E1;
1292 static const double Z2 = 4.92184563216946036703E1;
1293 
1294 static const double THPIO4 = 3.*PI/4.;
1295 
1296 double bessel_j1(double x)
1297 {
1298  double w, z, p, q, xn;
1299 
1300  DEBUG_ENTRY( "bessel_j1()" );
1301 
1302  w = x;
1303  if( x < 0 )
1304  w = -x;
1305 
1306  if( w <= 5.0 )
1307  {
1308  z = x * x;
1309  w = polevl( z, b1_RP, 3 ) / p1evl( z, b1_RQ, 8 );
1310  w = w * x * (z - Z1) * (z - Z2);
1311  return w;
1312  }
1313 
1314  w = 5.0/x;
1315  z = w * w;
1316  p = polevl( z, b1_PP, 6)/polevl( z, b1_PQ, 6 );
1317  q = polevl( z, b1_QP, 7)/p1evl( z, b1_QQ, 7 );
1318  xn = x - THPIO4;
1319  p = p * cos(xn) - w * q * sin(xn);
1320  return p * SQ2OPI / sqrt(x);
1321 }
1322 
1323 double bessel_y1(double x)
1324 {
1325  double w, z, p, q, xn;
1326 
1327  DEBUG_ENTRY( "bessel_y1()" );
1328 
1329  if( x <= 5.0 )
1330  {
1331  if( x <= 0.0 )
1332  {
1333  fprintf( ioQQQ, "bessel_y1: domain error\n" );
1335  }
1336  z = x * x;
1337  w = x * (polevl( z, b1_YP, 5 ) / p1evl( z, b1_YQ, 8 ));
1338  w += TWOOPI * ( bessel_j1(x) * log(x) - 1.0/x );
1339  return w;
1340  }
1341 
1342  w = 5.0/x;
1343  z = w * w;
1344  p = polevl( z, b1_PP, 6 )/polevl( z, b1_PQ, 6 );
1345  q = polevl( z, b1_QP, 7 )/p1evl( z, b1_QQ, 7 );
1346  xn = x - THPIO4;
1347  p = p * sin(xn) + w * q * cos(xn);
1348  return p * SQ2OPI / sqrt(x);
1349 }
1350 
1351 /* jn.c
1352  *
1353  * Bessel function of integer order
1354  *
1355  *
1356  *
1357  * SYNOPSIS:
1358  *
1359  * int n;
1360  * double x, y, jn();
1361  *
1362  * y = jn( n, x );
1363  *
1364  *
1365  *
1366  * DESCRIPTION:
1367  *
1368  * Returns Bessel function of order n, where n is a
1369  * (possibly negative) integer.
1370  *
1371  * The ratio of jn(x) to j0(x) is computed by backward
1372  * recurrence. First the ratio jn/jn-1 is found by a
1373  * continued fraction expansion. Then the recurrence
1374  * relating successive orders is applied until j0 or j1 is
1375  * reached.
1376  *
1377  * If n = 0 or 1 the routine for j0 or j1 is called
1378  * directly.
1379  *
1380  *
1381  *
1382  * ACCURACY:
1383  *
1384  * Absolute error:
1385  * arithmetic range # trials peak rms
1386  * DEC 0, 30 5500 6.9e-17 9.3e-18
1387  * IEEE 0, 30 5000 4.4e-16 7.9e-17
1388  *
1389  *
1390  * Not suitable for large n or x. Use jv() instead.
1391  *
1392  */
1393 
1394 /* jn.c
1395 Cephes Math Library Release 2.8: June, 2000
1396 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1397 */
1398 
1399 double bessel_jn(int n, double x)
1400 {
1401  double pkm2, pkm1, pk, xk, r, ans;
1402  int k, sign;
1403 
1404  DEBUG_ENTRY( "bessel_jn()" );
1405 
1406  if( n < 0 )
1407  {
1408  n = -n;
1409  if( (n & 1) == 0 ) /* -1**n */
1410  sign = 1;
1411  else
1412  sign = -1;
1413  }
1414  else
1415  sign = 1;
1416 
1417  if( x < 0.0 )
1418  {
1419  if( n & 1 )
1420  sign = -sign;
1421  x = -x;
1422  }
1423 
1424  if( x < DBL_EPSILON )
1425  {
1426  return sign * powi(x/2.,n)/factorial(n);
1427  }
1428 
1429  if( n == 0 )
1430  {
1431  return sign * bessel_j0(x);
1432  }
1433  if( n == 1 )
1434  {
1435  return sign * bessel_j1(x);
1436  }
1437  // avoid cancellation error for very small x
1438  if( n == 2 && x > 0.1 )
1439  {
1440  return sign * (2.0 * bessel_j1(x) / x - bessel_j0(x));
1441  }
1442 
1443  /* continued fraction */
1444  k = 52;
1445 
1446  pk = 2 * (n + k);
1447  ans = pk;
1448  xk = x * x;
1449 
1450  do
1451  {
1452  pk -= 2.0;
1453  ans = pk - (xk/ans);
1454  }
1455  while( --k > 0 );
1456  ans = x/ans;
1457 
1458  /* backward recurrence */
1459  pk = 1.0;
1460  pkm1 = 1.0/ans;
1461  k = n-1;
1462  r = 2 * k;
1463 
1464  do
1465  {
1466  pkm2 = (pkm1 * r - pk * x) / x;
1467  pk = pkm1;
1468  pkm1 = pkm2;
1469  r -= 2.0;
1470  }
1471  while( --k > 0 );
1472 
1473  if( fabs(pk) > fabs(pkm1) )
1474  ans = bessel_j1(x)/pk;
1475  else
1476  ans = bessel_j0(x)/pkm1;
1477  return sign * ans;
1478 }
1479 
1480 /* yn.c
1481  *
1482  * Bessel function of second kind of integer order
1483  *
1484  *
1485  *
1486  * SYNOPSIS:
1487  *
1488  * double x, y, yn();
1489  * int n;
1490  *
1491  * y = yn( n, x );
1492  *
1493  *
1494  *
1495  * DESCRIPTION:
1496  *
1497  * Returns Bessel function of order n, where n is a
1498  * (possibly negative) integer.
1499  *
1500  * The function is evaluated by forward recurrence on
1501  * n, starting with values computed by the routines
1502  * y0() and y1().
1503  *
1504  * If n = 0 or 1 the routine for y0 or y1 is called
1505  * directly.
1506  *
1507  *
1508  *
1509  * ACCURACY:
1510  *
1511  *
1512  * Absolute error, except relative
1513  * when y > 1:
1514  * arithmetic domain # trials peak rms
1515  * DEC 0, 30 2200 2.9e-16 5.3e-17
1516  * IEEE 0, 30 30000 3.4e-15 4.3e-16
1517  *
1518  *
1519  * ERROR MESSAGES:
1520  *
1521  * message condition value returned
1522  * yn singularity x = 0 MAXNUM
1523  * yn overflow MAXNUM
1524  *
1525  * Spot checked against tables for x, n between 0 and 100.
1526  *
1527  */
1528 
1529 /*
1530 Cephes Math Library Release 2.8: June, 2000
1531 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1532 */
1533 
1534 double bessel_yn(int n, double x)
1535 {
1536  double an, anm1, anm2, r;
1537  int k, sign;
1538 
1539  DEBUG_ENTRY( "bessel_yn()" );
1540 
1541  if( n < 0 )
1542  {
1543  n = -n;
1544  if( (n & 1) == 0 ) /* -1**n */
1545  sign = 1;
1546  else
1547  sign = -1;
1548  }
1549  else
1550  sign = 1;
1551 
1552  if( n == 0 )
1553  {
1554  return sign * bessel_y0(x);
1555  }
1556  if( n == 1 )
1557  {
1558  return sign * bessel_y1(x);
1559  }
1560 
1561  /* test for overflow */
1562  if( x <= 0.0 )
1563  {
1564  fprintf( ioQQQ, "bessel_yn: domain error\n" );
1566  }
1567 
1568  /* forward recurrence on n */
1569  anm2 = bessel_y0(x);
1570  anm1 = bessel_y1(x);
1571  k = 1;
1572  r = 2.0;
1573  do
1574  {
1575  an = r * anm1 / x - anm2;
1576  anm2 = anm1;
1577  anm1 = an;
1578  r += 2.0;
1579  ++k;
1580  }
1581  while( k < n );
1582  return sign * an;
1583 }
1584 
1585 /* ellpk.c
1586  *
1587  * Complete elliptic integral of the first kind
1588  *
1589  *
1590  *
1591  * SYNOPSIS:
1592  *
1593  * double m1, y, ellpk();
1594  *
1595  * y = ellpk( m1 );
1596  *
1597  *
1598  *
1599  * DESCRIPTION:
1600  *
1601  * Approximates the integral
1602  *
1603  *
1604  *
1605  * pi/2
1606  * -
1607  * | |
1608  * | dt
1609  * K(m) = | ------------------
1610  * | 2
1611  * | | sqrt( 1 - m sin t )
1612  * -
1613  * 0
1614  *
1615  * where m = 1 - m1, using the approximation
1616  *
1617  * P(x) - log x Q(x).
1618  *
1619  * The argument m1 is used rather than m so that the logarithmic
1620  * singularity at m = 1 will be shifted to the origin; this
1621  * preserves maximum accuracy.
1622  *
1623  * K(0) = pi/2.
1624  *
1625  * ACCURACY:
1626  *
1627  * Relative error:
1628  * arithmetic domain # trials peak rms
1629  * DEC 0,1 16000 3.5e-17 1.1e-17
1630  * IEEE 0,1 30000 2.5e-16 6.8e-17
1631  *
1632  * ERROR MESSAGES:
1633  *
1634  * message condition value returned
1635  * ellpk domain x<0, x>1 0.0
1636  *
1637  */
1638 
1639 /*
1640 Cephes Math Library, Release 2.8: June, 2000
1641 Copyright 1984, 1987, 2000 by Stephen L. Moshier
1642 */
1643 
1644 static const double elk_P[] =
1645 {
1646  1.37982864606273237150e-4,
1647  2.28025724005875567385e-3,
1648  7.97404013220415179367e-3,
1649  9.85821379021226008714e-3,
1650  6.87489687449949877925e-3,
1651  6.18901033637687613229e-3,
1652  8.79078273952743772254e-3,
1653  1.49380448916805252718e-2,
1654  3.08851465246711995998e-2,
1655  9.65735902811690126535e-2,
1656  1.38629436111989062502e0
1657 };
1658 
1659 static const double elk_Q[] =
1660 {
1661  2.94078955048598507511e-5,
1662  9.14184723865917226571e-4,
1663  5.94058303753167793257e-3,
1664  1.54850516649762399335e-2,
1665  2.39089602715924892727e-2,
1666  3.01204715227604046988e-2,
1667  3.73774314173823228969e-2,
1668  4.88280347570998239232e-2,
1669  7.03124996963957469739e-2,
1670  1.24999999999870820058e-1,
1671  4.99999999999999999821e-1
1672 };
1673 
1674 static const double C1 = 1.3862943611198906188e0; /* log(4) */
1675 
1676 double ellpk(double x)
1677 {
1678  DEBUG_ENTRY( "ellpk()" );
1679 
1680  if( x < 0.0 || x > 1.0 )
1681  {
1682  fprintf( ioQQQ, "ellpk: domain error\n" );
1684  }
1685 
1686  if( x > DBL_EPSILON )
1687  {
1688  return polevl(x,elk_P,10) - log(x) * polevl(x,elk_Q,10);
1689  }
1690  else
1691  {
1692  if( x == 0.0 )
1693  {
1694  fprintf( ioQQQ, "ellpk: domain error\n" );
1696  }
1697  else
1698  {
1699  return C1 - 0.5 * log(x);
1700  }
1701  }
1702 }
1703 
1704 /* expn.c
1705  *
1706  * Exponential integral En
1707  *
1708  *
1709  *
1710  * SYNOPSIS:
1711  *
1712  * int n;
1713  * double x, y, expn();
1714  *
1715  * y = expn( n, x );
1716  *
1717  *
1718  *
1719  * DESCRIPTION:
1720  *
1721  * Evaluates the exponential integral
1722  *
1723  * inf.
1724  * -
1725  * | | -xt
1726  * | e
1727  * E (x) = | ---- dt.
1728  * n | n
1729  * | | t
1730  * -
1731  * 1
1732  *
1733  *
1734  * Both n and x must be nonnegative.
1735  *
1736  * The routine employs either a power series, a continued
1737  * fraction, or an asymptotic formula depending on the
1738  * relative values of n and x.
1739  *
1740  * ACCURACY:
1741  *
1742  * Relative error:
1743  * arithmetic domain # trials peak rms
1744  * DEC 0, 30 5000 2.0e-16 4.6e-17
1745  * IEEE 0, 30 10000 1.7e-15 3.6e-16
1746  *
1747  */
1748 
1749 /* Cephes Math Library Release 2.8: June, 2000
1750  Copyright 1985, 2000 by Stephen L. Moshier */
1751 
1752 static const double MAXLOG = log(DBL_MAX);
1753 static const double BIG = 1.44115188075855872E+17; /* 2^57 */
1754 
1755 /*expn exponential intergral for any n */
1756 double expn(int n, double x)
1757 {
1758  double ans, r, t, yk, xk;
1759  double pk, pkm1, pkm2, qk, qkm1, qkm2;
1760  double psi, z;
1761  int i, k;
1762 
1763  DEBUG_ENTRY( "expn()" );
1764 
1765  if( n < 0 || x < 0. )
1766  {
1767  fprintf( ioQQQ, "expn: domain error\n" );
1769  }
1770 
1771  if( x > MAXLOG )
1772  {
1773  return 0.0;
1774  }
1775 
1776  if( x == 0.0 )
1777  {
1778  if( n < 2 )
1779  {
1780  fprintf( ioQQQ, "expn: domain error\n" );
1782  }
1783  else
1784  {
1785  return 1.0/((double)n-1.0);
1786  }
1787  }
1788 
1789  if( n == 0 )
1790  {
1791  return exp(-x)/x;
1792  }
1793 
1794  /* Expansion for large n */
1795  if( n > 5000 )
1796  {
1797  xk = x + n;
1798  yk = 1.0 / (xk * xk);
1799  t = n;
1800  ans = yk * t * (6.0 * x * x - 8.0 * t * x + t * t);
1801  ans = yk * (ans + t * (t - 2.0 * x));
1802  ans = yk * (ans + t);
1803  ans = (ans + 1.0) * exp( -x ) / xk;
1804  return ans;
1805  }
1806 
1807  if( x <= 1.0 )
1808  {
1809  /* Power series expansion */
1810  psi = -EULER - log(x);
1811  for( i=1; i < n; i++ )
1812  psi = psi + 1.0/i;
1813 
1814  z = -x;
1815  xk = 0.0;
1816  yk = 1.0;
1817  pk = 1.0 - n;
1818  if( n == 1 )
1819  ans = 0.0;
1820  else
1821  ans = 1.0/pk;
1822  do
1823  {
1824  xk += 1.0;
1825  yk *= z/xk;
1826  pk += 1.0;
1827  if( pk != 0.0 )
1828  {
1829  ans += yk/pk;
1830  }
1831  if( ans != 0.0 )
1832  t = fabs(yk/ans);
1833  else
1834  t = 1.0;
1835  }
1836  while( t > DBL_EPSILON );
1837  ans = powi(z,n-1)*psi/factorial(n-1) - ans;
1838  return ans;
1839  }
1840  else
1841  {
1842  /* continued fraction */
1843  k = 1;
1844  pkm2 = 1.0;
1845  qkm2 = x;
1846  pkm1 = 1.0;
1847  qkm1 = x + n;
1848  ans = pkm1/qkm1;
1849 
1850  do
1851  {
1852  k += 1;
1853  if( is_odd(k) )
1854  {
1855  yk = 1.0;
1856  xk = static_cast<double>(n + (k-1)/2);
1857  }
1858  else
1859  {
1860  yk = x;
1861  xk = static_cast<double>(k/2);
1862  }
1863  pk = pkm1 * yk + pkm2 * xk;
1864  qk = qkm1 * yk + qkm2 * xk;
1865  if( qk != 0 )
1866  {
1867  r = pk/qk;
1868  t = fabs( (ans - r)/r );
1869  ans = r;
1870  }
1871  else
1872  t = 1.0;
1873  pkm2 = pkm1;
1874  pkm1 = pk;
1875  qkm2 = qkm1;
1876  qkm1 = qk;
1877  if( fabs(pk) > BIG )
1878  {
1879  pkm2 /= BIG;
1880  pkm1 /= BIG;
1881  qkm2 /= BIG;
1882  qkm1 /= BIG;
1883  }
1884  }
1885  while( t > DBL_EPSILON );
1886 
1887  ans *= exp( -x );
1888  return ans;
1889  }
1890 }
1891 
1892 /* erf.c
1893  *
1894  * Error function
1895  *
1896  *
1897  *
1898  * SYNOPSIS:
1899  *
1900  * double x, y, erf();
1901  *
1902  * y = erf( x );
1903  *
1904  *
1905  *
1906  * DESCRIPTION:
1907  *
1908  * The integral is
1909  *
1910  * x
1911  * -
1912  * 2 | | 2
1913  * erf(x) = -------- | exp( - t ) dt.
1914  * sqrt(pi) | |
1915  * -
1916  * 0
1917  *
1918  * The magnitude of x is limited to 9.231948545 for DEC
1919  * arithmetic; 1 or -1 is returned outside this range.
1920  *
1921  * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
1922  * erf(x) = 1 - erfc(x).
1923  *
1924  *
1925  *
1926  * ACCURACY:
1927  *
1928  * Relative error:
1929  * arithmetic domain # trials peak rms
1930  * DEC 0,1 14000 4.7e-17 1.5e-17
1931  * IEEE 0,1 30000 3.7e-16 1.0e-16
1932  *
1933  */
1934 /* erfc.c
1935  *
1936  * Complementary error function
1937  *
1938  *
1939  *
1940  * SYNOPSIS:
1941  *
1942  * double x, y, erfc();
1943  *
1944  * y = erfc( x );
1945  *
1946  *
1947  *
1948  * DESCRIPTION:
1949  *
1950  *
1951  * 1 - erf(x) =
1952  *
1953  * inf.
1954  * -
1955  * 2 | | 2
1956  * erfc(x) = -------- | exp( - t ) dt
1957  * sqrt(pi) | |
1958  * -
1959  * x
1960  *
1961  *
1962  * For small x, erfc(x) = 1 - erf(x); otherwise rational
1963  * approximations are computed.
1964  *
1965  * A special function expx2.c is used to suppress error amplification
1966  * in computing exp(-x^2).
1967  *
1968  *
1969  * ACCURACY:
1970  *
1971  * Relative error:
1972  * arithmetic domain # trials peak rms
1973  * IEEE 0,26.6417 30000 1.3e-15 2.2e-16
1974  *
1975  *
1976  * ERROR MESSAGES:
1977  *
1978  * message condition value returned
1979  * erfc underflow x > 9.231948545 (DEC) 0.0
1980  *
1981  *
1982  */
1983 
1984 /*
1985 Cephes Math Library Release 2.9: November, 2000
1986 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
1987 */
1988 
1989 static double erf_P[] = {
1990  2.46196981473530512524e-10,
1991  5.64189564831068821977e-1,
1992  7.46321056442269912687e0,
1993  4.86371970985681366614e1,
1994  1.96520832956077098242e2,
1995  5.26445194995477358631e2,
1996  9.34528527171957607540e2,
1997  1.02755188689515710272e3,
1998  5.57535335369399327526e2
1999 };
2000 static double erf_Q[] = {
2001 /* 1.00000000000000000000e0,*/
2002  1.32281951154744992508e1,
2003  8.67072140885989742329e1,
2004  3.54937778887819891062e2,
2005  9.75708501743205489753e2,
2006  1.82390916687909736289e3,
2007  2.24633760818710981792e3,
2008  1.65666309194161350182e3,
2009  5.57535340817727675546e2
2010 };
2011 static double erf_R[] = {
2012  5.64189583547755073984e-1,
2013  1.27536670759978104416e0,
2014  5.01905042251180477414e0,
2015  6.16021097993053585195e0,
2016  7.40974269950448939160e0,
2017  2.97886665372100240670e0
2018 };
2019 static double erf_S[] = {
2020 /* 1.00000000000000000000e0,*/
2021  2.26052863220117276590e0,
2022  9.39603524938001434673e0,
2023  1.20489539808096656605e1,
2024  1.70814450747565897222e1,
2025  9.60896809063285878198e0,
2026  3.36907645100081516050e0
2027 };
2028 
2029 
2030 #ifndef HAVE_ERFC
2031 
2032 STATIC double expx2(double x, int sign);
2033 
2034 /* Define this macro to suppress error propagation in exp(x^2)
2035  by using the expx2 function. The tradeoff is that doing so
2036  generates two calls to the exponential function instead of one. */
2037 const bool lgUSE_EXPXSQ = true;
2038 
2039 double erfc(double a)
2040 {
2041  double p,q,x,y,z;
2042 
2043  DEBUG_ENTRY( "erfc()" );
2044 
2045  x = abs(a);
2046 
2047  if( x < 1.0 )
2048  return 1.0 - erf(a);
2049 
2050  z = -a * a;
2051 
2052  if( z < -MAXLOG )
2053  return ( a < 0.0 ) ? 2.0 : 0.0;
2054 
2055  if( lgUSE_EXPXSQ )
2056  /* Compute z = exp(z). */
2057  z = expx2(a, -1);
2058  else
2059  z = exp(z);
2060 
2061  if( x < 8.0 )
2062  {
2063  p = polevl( x, erf_P, 8 );
2064  q = p1evl( x, erf_Q, 8 );
2065  }
2066  else
2067  {
2068  p = polevl( x, erf_R, 5 );
2069  q = p1evl( x, erf_S, 6 );
2070  }
2071  y = (z * p)/q;
2072 
2073  if( a < 0 )
2074  y = 2.0 - y;
2075 
2076  if( y == 0.0 )
2077  return ( a < 0. ) ? 2.0 : 0.0;
2078 
2079  return y;
2080 }
2081 
2082 #endif
2083 
2084 
2085 /* Exponentially scaled erfc function
2086  exp(x^2) erfc(x)
2087  valid for x > 1.
2088  Use with ndtr and expx2. */
2089 double erfce(double x)
2090 {
2091  double p,q;
2092 
2093  DEBUG_ENTRY( "erfce()" );
2094 
2095  if( x < 8.0 )
2096  {
2097  p = polevl( x, erf_P, 8 );
2098  q = p1evl( x, erf_Q, 8 );
2099  }
2100  else
2101  {
2102  p = polevl( x, erf_R, 5 );
2103  q = p1evl( x, erf_S, 6 );
2104  }
2105  return p/q;
2106 }
2107 
2108 
2109 #ifndef HAVE_ERF
2110 
2111 static double erf_T[] = {
2112  9.60497373987051638749e0,
2113  9.00260197203842689217e1,
2114  2.23200534594684319226e3,
2115  7.00332514112805075473e3,
2116  5.55923013010394962768e4
2117 };
2118 static double erf_U[] = {
2119 /* 1.00000000000000000000e0,*/
2120  3.35617141647503099647e1,
2121  5.21357949780152679795e2,
2122  4.59432382970980127987e3,
2123  2.26290000613890934246e4,
2124  4.92673942608635921086e4
2125 };
2126 
2127 double erf(double x)
2128 {
2129  double y, z;
2130 
2131  DEBUG_ENTRY( "erf()" );
2132 
2133  if( abs(x) > 1.0 )
2134  return 1.0 - erfc(x);
2135  z = x * x;
2136  y = x * polevl( z, erf_T, 4 ) / p1evl( z, erf_U, 5 );
2137  return y;
2138 }
2139 
2140 #endif
2141 
2142 
2143 /* expx2.c
2144  *
2145  * Exponential of squared argument
2146  *
2147  *
2148  *
2149  * SYNOPSIS:
2150  *
2151  * double x, y, expx2();
2152  * int sign;
2153  *
2154  * y = expx2( x, sign );
2155  *
2156  *
2157  *
2158  * DESCRIPTION:
2159  *
2160  * Computes y = exp(x*x) while suppressing error amplification
2161  * that would ordinarily arise from the inexactness of the
2162  * exponential argument x*x.
2163  *
2164  * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
2165  *
2166  *
2167  * ACCURACY:
2168  *
2169  * Relative error:
2170  * arithmetic domain # trials peak rms
2171  * IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17
2172  *
2173  */
2174 
2175 /*
2176 Cephes Math Library Release 2.9: June, 2000
2177 Copyright 2000 by Stephen L. Moshier
2178 */
2179 
2180 
2181 #ifndef HAVE_ERFC
2182 
2183 const double exp_M = 128.0;
2184 const double exp_MINV = .0078125;
2185 
2186 STATIC double expx2(double x, int sign)
2187 {
2188  double u, u1, m, f;
2189 
2190  DEBUG_ENTRY( "expx2()" );
2191 
2192  x = abs(x);
2193  if( sign < 0 )
2194  x = -x;
2195 
2196  /* Represent x as an exact multiple of exp_M plus a residual.
2197  exp_M is a power of 2 chosen so that exp(m * m) does not overflow
2198  or underflow and so that |x - m| is small. */
2199  m = exp_MINV * floor(exp_M * x + 0.5);
2200  f = x - m;
2201 
2202  /* x^2 = m^2 + 2mf + f^2 */
2203  u = m * m;
2204  u1 = 2 * m * f + f * f;
2205 
2206  if( sign < 0 )
2207  {
2208  u = -u;
2209  u1 = -u1;
2210  }
2211 
2212  if( (u+u1) > MAXLOG )
2213  return DBL_MAX;
2214 
2215  /* u is exact, u1 is small. */
2216  u = exp(u) * exp(u1);
2217  return u;
2218 }
2219 
2220 #endif
2221 
2222 
2223 /* igam.c
2224  *
2225  * Incomplete gamma integral
2226  *
2227  *
2228  *
2229  * SYNOPSIS:
2230  *
2231  * double a, x, y, igam();
2232  *
2233  * y = igam( a, x );
2234  *
2235  * DESCRIPTION:
2236  *
2237  * The function is defined by
2238  *
2239  * x
2240  * -
2241  * 1 | | -t a-1
2242  * igam(a,x) = ----- | e t dt.
2243  * - | |
2244  * | (a) -
2245  * 0
2246  *
2247  *
2248  * In this implementation both arguments must be positive.
2249  * The integral is evaluated by either a power series or
2250  * continued fraction expansion, depending on the relative
2251  * values of a and x.
2252  *
2253  * ACCURACY:
2254  *
2255  * Relative error:
2256  * arithmetic domain # trials peak rms
2257  * IEEE 0,30 200000 3.6e-14 2.9e-15
2258  * IEEE 0,100 300000 9.9e-14 1.5e-14
2259  */
2260 /* igamc()
2261  *
2262  * Complemented incomplete gamma integral
2263  *
2264  *
2265  *
2266  * SYNOPSIS:
2267  *
2268  * double a, x, y, igamc();
2269  *
2270  * y = igamc( a, x );
2271  *
2272  * DESCRIPTION:
2273  *
2274  * The function is defined by
2275  *
2276  *
2277  * igamc(a,x) = 1 - igam(a,x)
2278  *
2279  * inf.
2280  * -
2281  * 1 | | -t a-1
2282  * = ----- | e t dt.
2283  * - | |
2284  * | (a) -
2285  * x
2286  *
2287  *
2288  * In this implementation both arguments must be positive.
2289  * The integral is evaluated by either a power series or
2290  * continued fraction expansion, depending on the relative
2291  * values of a and x.
2292  *
2293  * ACCURACY:
2294  *
2295  * Tested at random a, x.
2296  * a x Relative error:
2297  * arithmetic domain domain # trials peak rms
2298  * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
2299  * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
2300  */
2301 /*
2302 Cephes Math Library Release 2.8: June, 2000
2303 Copyright 1985, 1987, 2000 by Stephen L. Moshier
2304 */
2305 
2306 STATIC double igamc_fraction( double a, double x );
2307 
2308 static const double igam_big = 4.503599627370496e15;
2309 static const double igam_biginv = 2.22044604925031308085e-16;
2310 
2311 double igamc( double a, double x )
2312 {
2313  DEBUG_ENTRY( "igamc()" );
2314 
2315  ASSERT( a > 0. && x > 0. );
2316 
2317  if( x < 1.0 || x < a )
2318  return 1.0 - igam(a,x);
2319 
2320  double ax = a * log(x) - x - lgamma(a);
2321  if( ax < -MAXLOG )
2322  return 0.0;
2323  ax = exp(ax);
2324 
2325  return ax*igamc_fraction(a,x);
2326 }
2327 
2328 /* calculate igamc(a,x)*exp(x) which prevents underflow for very large x */
2329 double igamc_scaled( double a, double x )
2330 {
2331  DEBUG_ENTRY( "igamc()" );
2332 
2333  ASSERT( a > 0. && x > 0. );
2334 
2335  if( x < 1.0 || x < a )
2336  return (1.0 - igam(a,x))*exp(x);
2337 
2338  double ax = a * log(x) - lgamma(a);
2339  if( ax < -MAXLOG )
2340  return 0.0;
2341  ax = exp(ax);
2342 
2343  return ax*igamc_fraction(a,x);
2344 }
2345 
2346 /* left tail of incomplete gamma function:
2347  *
2348  * inf. k
2349  * a -x - x
2350  * x e > ----------
2351  * - -
2352  * k=0 | (a+k+1)
2353  *
2354  */
2355 
2356 double igam( double a, double x )
2357 {
2358  DEBUG_ENTRY( "igam()" );
2359 
2360  ASSERT( a > 0. && x > 0. );
2361 
2362  double ans, ax, c, r;
2363 
2364  if( x > 1.0 && x > a )
2365  return 1.0 - igamc(a,x);
2366 
2367  /* Compute x**a * exp(-x) / gamma(a) */
2368  ax = a * log(x) - x - lgamma(a);
2369  if( ax < -MAXLOG )
2370  return 0.0;
2371  ax = exp(ax);
2372 
2373  /* power series */
2374  r = a;
2375  c = 1.0;
2376  ans = 1.0;
2377 
2378  do
2379  {
2380  r += 1.0;
2381  c *= x/r;
2382  ans += c;
2383  }
2384  while( c/ans > DBL_EPSILON );
2385 
2386  return ans * ax/a;
2387 }
2388 
2389 STATIC double igamc_fraction( double a, double x )
2390 {
2391  double ans, c, yc, r, t, y, z;
2392  double pk, pkm1, pkm2, qk, qkm1, qkm2;
2393 
2394  /* continued fraction */
2395  y = 1.0 - a;
2396  z = x + y + 1.0;
2397  c = 0.0;
2398  pkm2 = 1.0;
2399  qkm2 = x;
2400  pkm1 = x + 1.0;
2401  qkm1 = z * x;
2402  ans = pkm1/qkm1;
2403 
2404  do
2405  {
2406  c += 1.0;
2407  y += 1.0;
2408  z += 2.0;
2409  yc = y * c;
2410  pk = pkm1 * z - pkm2 * yc;
2411  qk = qkm1 * z - qkm2 * yc;
2412  if( qk != 0 )
2413  {
2414  r = pk/qk;
2415  t = fabs( (ans - r)/r );
2416  ans = r;
2417  }
2418  else
2419  t = 1.0;
2420  pkm2 = pkm1;
2421  pkm1 = pk;
2422  qkm2 = qkm1;
2423  qkm1 = qk;
2424  if( fabs(pk) > igam_big )
2425  {
2426  pkm2 *= igam_biginv;
2427  pkm1 *= igam_biginv;
2428  qkm2 *= igam_biginv;
2429  qkm1 *= igam_biginv;
2430  }
2431  }
2432  while( t > DBL_EPSILON );
2433 
2434  return ans;
2435 }
2436 
2437 
2438 /* polevl.c
2439  * p1evl.c
2440  *
2441  * Evaluate polynomial
2442  *
2443  *
2444  *
2445  * SYNOPSIS:
2446  *
2447  * int N;
2448  * double x, y, coef[N+1], polevl[];
2449  *
2450  * y = polevl( x, coef, N );
2451  *
2452  *
2453  *
2454  * DESCRIPTION:
2455  *
2456  * Evaluates polynomial of degree N:
2457  *
2458  * 2 N
2459  * y = C + C x + C x +...+ C x
2460  * 0 1 2 N
2461  *
2462  * Coefficients are stored in reverse order:
2463  *
2464  * coef[0] = C , ..., coef[N] = C .
2465  * N 0
2466  *
2467  * The function p1evl() assumes that coef[N] = 1.0 and is
2468  * omitted from the array. Its calling arguments are
2469  * otherwise the same as polevl().
2470  *
2471  *
2472  * SPEED:
2473  *
2474  * In the interest of speed, there are no checks for out
2475  * of bounds arithmetic. This routine is used by most of
2476  * the functions in the library. Depending on available
2477  * equipment features, the user may wish to rewrite the
2478  * program in microcode or assembly language.
2479  *
2480  */
2481 
2482 /*
2483 Cephes Math Library Release 2.1: December, 1988
2484 Copyright 1984, 1987, 1988 by Stephen L. Moshier
2485 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2486 */
2487 
2488 inline double polevl(double x, const double coef[], int N)
2489 {
2490  double ans;
2491  int i;
2492  const double *p = coef;
2493 
2494  ans = *p++;
2495  i = N;
2496 
2497  do
2498  ans = ans * x + *p++;
2499  while( --i );
2500 
2501  return ans;
2502 }
2503 
2504 /* p1evl() */
2505 /* N
2506  * Evaluate polynomial when coefficient of x is 1.0.
2507  * Otherwise same as polevl.
2508  */
2509 
2510 inline double p1evl(double x, const double coef[], int N)
2511 {
2512  double ans;
2513  const double *p = coef;
2514  int i;
2515 
2516  ans = x + *p++;
2517  i = N-1;
2518 
2519  do
2520  ans = ans * x + *p++;
2521  while( --i );
2522 
2523  return ans;
2524 }
2525 
2526 /* chbevl.c
2527  *
2528  * Evaluate Chebyshev series
2529  *
2530  *
2531  *
2532  * SYNOPSIS:
2533  *
2534  * int N;
2535  * double x, y, coef[N], chebevl();
2536  *
2537  * y = chbevl( x, coef, N );
2538  *
2539  *
2540  *
2541  * DESCRIPTION:
2542  *
2543  * Evaluates the series
2544  *
2545  * N-1
2546  * - '
2547  * y = > coef[i] T (x/2)
2548  * - i
2549  * i=0
2550  *
2551  * of Chebyshev polynomials Ti at argument x/2.
2552  *
2553  * Coefficients are stored in reverse order, i.e. the zero
2554  * order term is last in the array. Note N is the number of
2555  * coefficients, not the order.
2556  *
2557  * If coefficients are for the interval a to b, x must
2558  * have been transformed to x -> 2(2x - b - a)/(b-a) before
2559  * entering the routine. This maps x from (a, b) to (-1, 1),
2560  * over which the Chebyshev polynomials are defined.
2561  *
2562  * If the coefficients are for the inverted interval, in
2563  * which (a, b) is mapped to (1/b, 1/a), the transformation
2564  * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
2565  * this becomes x -> 4a/x - 1.
2566  *
2567  *
2568  *
2569  * SPEED:
2570  *
2571  * Taking advantage of the recurrence properties of the
2572  * Chebyshev polynomials, the routine requires one more
2573  * addition per loop than evaluating a nested polynomial of
2574  * the same degree.
2575  *
2576  */
2577 
2578 /*
2579 Cephes Math Library Release 2.0: April, 1987
2580 Copyright 1985, 1987 by Stephen L. Moshier
2581 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
2582 */
2583 
2584 inline double chbevl(double x, const double array[], int n)
2585 {
2586  double b0, b1, b2;
2587  const double *p = array;
2588  int i;
2589 
2590  b0 = *p++;
2591  b1 = 0.0;
2592  i = n - 1;
2593 
2594  do
2595  {
2596  b2 = b1;
2597  b1 = b0;
2598  b0 = x * b1 - b2 + *p++;
2599  }
2600  while( --i );
2601 
2602  return 0.5*(b0-b2);
2603 }
2604 
2605 /*******************************************************************
2606  * This marks the end of the block of code from the Cephes library *
2607  *******************************************************************/
2608 
2609 /****************************************************************************
2610  The code below is adapted from the Boost library version 1.62.0
2611  obtained from http://www.boost.org and subject to the following license:
2612 
2613 Boost Software License - Version 1.0 - August 17th, 2003
2614 
2615 Permission is hereby granted, free of charge, to any person or organization
2616 obtaining a copy of the software and accompanying documentation covered by
2617 this license (the "Software") to use, reproduce, display, distribute,
2618 execute, and transmit the Software, and to prepare derivative works of the
2619 Software, and to permit third-parties to whom the Software is furnished to
2620 do so, all subject to the following:
2621 
2622 The copyright notices in the Software and this entire statement, including
2623 the above license grant, this restriction and the following disclaimer,
2624 must be included in all copies of the Software, in whole or in part, and
2625 all derivative works of the Software, unless such copies or derivative
2626 works are solely in the form of machine-executable object code generated by
2627 a source language processor.
2628 
2629 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
2630 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
2631 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
2632 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
2633 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
2634 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
2635 DEALINGS IN THE SOFTWARE.
2636 *****************************************************************************/
2637 
2638 inline double ratevl_5_4(double x, const double a[5], const double b[4]);
2639 inline double ratevl_6_3(double x, const double a[6], const double b[3]);
2640 inline double ratevl_6_4(double x, const double a[6], const double b[4]);
2641 inline double ratevl_7_8(double x, const double a[7], const double b[8]);
2642 inline double ratevl_8_7(double x, const double a[8], const double b[7]);
2643 inline double ratevl_10_11(double x, const double a[10], const double b[11]);
2644 inline double ratevl_11_10(double x, const double a[11], const double b[10]);
2645 inline double ratevl_15_6(double x, const double a[15], const double b[6]);
2646 
2647 static const double BESSEL_K0_P1[] = {
2648  2.4708152720399552679e+03, 5.9169059852270512312e+03, 4.6850901201934832188e+02,
2649  1.1999463724910714109e+01, 1.3166052564989571850e-01, 5.8599221412826100000e-04
2650 };
2651 static const double BESSEL_K0_Q1[] = {
2652  2.1312714303849120380e+04, -2.4994418972832303646e+02, 1.0
2653 };
2654 static const double BESSEL_K0_P2[] = {
2655  -1.6128136304458193998e+06, -3.7333769444840079748e+05, -1.7984434409411765813e+04,
2656  -2.9501657892958843865e+02, -1.6414452837299064100e+00
2657 };
2658 static const double BESSEL_K0_Q2[] = {
2659  -1.6128136304458193998e+06, 2.9865713163054025489e+04, -2.5064972445877992730e+02, 1.0
2660 };
2661 static const double BESSEL_K0_P3[] = {
2662  1.1600249425076035558e+02, 2.3444738764199315021e+03, 1.8321525870183537725e+04,
2663  7.1557062783764037541e+04, 1.5097646353289914539e+05, 1.7398867902565686251e+05,
2664  1.0577068948034021957e+05, 3.1075408980684392399e+04, 3.6832589957340267940e+03,
2665  1.1394980557384778174e+02
2666 };
2667 static const double BESSEL_K0_Q3[] = {
2668  9.2556599177304839811e+01, 1.8821890840982713696e+03, 1.4847228371802360957e+04,
2669  5.8824616785857027752e+04, 1.2689839587977598727e+05, 1.5144644673520157801e+05,
2670  9.7418829762268075784e+04, 3.1474655750295278825e+04, 4.4329628889746408858e+03,
2671  2.0013443064949242491e+02, 1.0
2672 };
2673 
2674 static const double BESSEL_K1_P1[] = {
2675  -2.2149374878243304548e+06, 7.1938920065420586101e+05, 1.7733324035147015630e+05,
2676  7.1885382604084798576e+03, 9.9991373567429309922e+01, 4.8127070456878442310e-01
2677 };
2678 static const double BESSEL_K1_Q1[] = {
2679  -2.2149374878243304548e+06, 3.7264298672067697862e+04, -2.8143915754538725829e+02, 1.0
2680 };
2681 static const double BESSEL_K1_P2[] = {
2682  0.0, -1.3531161492785421328e+06, -1.4758069205414222471e+05,
2683  -4.5051623763436087023e+03, -5.3103913335180275253e+01, -2.2795590826955002390e-01
2684 };
2685 static const double BESSEL_K1_Q2[] = {
2686  -2.7062322985570842656e+06, 4.3117653211351080007e+04, -3.0507151578787595807e+02, 1.0
2687 };
2688 static const double BESSEL_K1_P3[] = {
2689  2.2196792496874548962e+00, 4.4137176114230414036e+01, 3.4122953486801312910e+02,
2690  1.3319486433183221990e+03, 2.8590657697910288226e+03, 3.4540675585544584407e+03,
2691  2.3123742209168871550e+03, 8.1094256146537402173e+02, 1.3182609918569941308e+02,
2692  7.5584584631176030810e+00, 6.4257745859173138767e-02
2693 };
2694 static const double BESSEL_K1_Q3[] = {
2695  1.7710478032601086579e+00, 3.4552228452758912848e+01, 2.5951223655579051357e+02,
2696  9.6929165726802648634e+02, 1.9448440788918006154e+03, 2.1181000487171943810e+03,
2697  1.2082692316002348638e+03, 3.3031020088765390854e+02, 3.6001069306861518855e+01, 1.0
2698 };
2699 
2700 double bessel_k0(double x)
2701 {
2702  DEBUG_ENTRY( "bessel_k0()" );
2703 
2704  if( x <= 0. )
2705  throw domain_error( "bessel_k0" );
2706  if( x <= 1. )
2707  {
2708  double y = x*x;
2709  double r1 = ratevl_6_3(y, BESSEL_K0_P1, BESSEL_K0_Q1);
2710  double r2 = ratevl_5_4(y, BESSEL_K0_P2, BESSEL_K0_Q2);
2711  double factor = log(x);
2712  return r1 - factor * r2;
2713  }
2714  else
2715  {
2716  double y = 1./x;
2717  double r = ratevl_10_11(y, BESSEL_K0_P3, BESSEL_K0_Q3);
2718  double factor = exp(-x) / sqrt(x);
2719  return factor * r;
2720  }
2721 }
2722 
2723 double bessel_k0_scaled(double x)
2724 {
2725  DEBUG_ENTRY( "bessel_k0_scaled()" );
2726 
2727  if( x <= 0. )
2728  throw domain_error( "bessel_k0_scaled" );
2729  if( x <= 1. )
2730  {
2731  double y = x*x;
2732  double r1 = ratevl_6_3(y, BESSEL_K0_P1, BESSEL_K0_Q1);
2733  double r2 = ratevl_5_4(y, BESSEL_K0_P2, BESSEL_K0_Q2);
2734  double factor = log(x);
2735  return exp(x) * (r1 - factor * r2);
2736  }
2737  else
2738  {
2739  double y = 1./x;
2740  double r = ratevl_10_11(y, BESSEL_K0_P3, BESSEL_K0_Q3);
2741  double factor = 1. / sqrt(x);
2742  return factor * r;
2743  }
2744 }
2745 
2746 double bessel_k1(double x)
2747 {
2748  DEBUG_ENTRY( "bessel_k1()" );
2749 
2750  if( x <= 0. )
2751  throw domain_error( "bessel_k1" );
2752  if( x <= 1. )
2753  {
2754  double y = x*x;
2755  double r1 = ratevl_6_4(y, BESSEL_K1_P1, BESSEL_K1_Q1);
2756  double r2 = ratevl_6_4(y, BESSEL_K1_P2, BESSEL_K1_Q2);
2757  double factor = log(x);
2758  return (r1 + factor * r2) / x;
2759  }
2760  else
2761  {
2762  double y = 1./x;
2763  double r1 = ratevl_11_10(y, BESSEL_K1_P3, BESSEL_K1_Q3);
2764  double factor = exp(-x) / sqrt(x);
2765  return factor * r1;
2766  }
2767 }
2768 
2769 double bessel_k1_scaled(double x)
2770 {
2771  DEBUG_ENTRY( "bessel_k1_scaled()" );
2772 
2773  if( x <= 0. )
2774  throw domain_error( "bessel_k1_scaled" );
2775  if( x <= 1. )
2776  {
2777  double y = x*x;
2778  double r1 = ratevl_6_4(y, BESSEL_K1_P1, BESSEL_K1_Q1);
2779  double r2 = ratevl_6_4(y, BESSEL_K1_P2, BESSEL_K1_Q2);
2780  double factor = log(x);
2781  return exp(x) * (r1 + factor * r2) / x;
2782  }
2783  else
2784  {
2785  double y = 1./x;
2786  double r1 = ratevl_11_10(y, BESSEL_K1_P3, BESSEL_K1_Q3);
2787  double factor = 1. / sqrt(x);
2788  return factor * r1;
2789  }
2790 }
2791 
2792 void bessel_k0_k1(double x, double* k0val, double* k1val)
2793 {
2794  DEBUG_ENTRY( "bessel_k0_k1()" );
2795 
2796  if( x <= 0. )
2797  throw domain_error( "bessel_k0_k1" );
2798  if( x <= 1. )
2799  {
2800  double y = x*x;
2801  double r1_0 = ratevl_6_3(y, BESSEL_K0_P1, BESSEL_K0_Q1);
2802  double r2_0 = ratevl_5_4(y, BESSEL_K0_P2, BESSEL_K0_Q2);
2803  double r1_1 = ratevl_6_4(y, BESSEL_K1_P1, BESSEL_K1_Q1);
2804  double r2_1 = ratevl_6_4(y, BESSEL_K1_P2, BESSEL_K1_Q2);
2805  double factor = log(x);
2806  *k0val = r1_0 - factor * r2_0;
2807  *k1val = (r1_1 + factor * r2_1) / x;
2808  }
2809  else
2810  {
2811  double y = 1./x;
2812  double r0 = ratevl_10_11(y, BESSEL_K0_P3, BESSEL_K0_Q3);
2813  double r1 = ratevl_11_10(y, BESSEL_K1_P3, BESSEL_K1_Q3);
2814  double factor = exp(-x) / sqrt(x);
2815  *k0val = factor * r0;
2816  *k1val = factor * r1;
2817  }
2818 }
2819 
2820 void bessel_k0_k1_scaled(double x, double* k0val, double* k1val)
2821 {
2822  DEBUG_ENTRY( "bessel_k0_k1_scaled()" );
2823 
2824  if( x <= 0. )
2825  throw domain_error( "bessel_k0_k1_scaled" );
2826  if( x <= 1. )
2827  {
2828  double y = x*x;
2829  double r1_0 = ratevl_6_3(y, BESSEL_K0_P1, BESSEL_K0_Q1);
2830  double r2_0 = ratevl_5_4(y, BESSEL_K0_P2, BESSEL_K0_Q2);
2831  double r1_1 = ratevl_6_4(y, BESSEL_K1_P1, BESSEL_K1_Q1);
2832  double r2_1 = ratevl_6_4(y, BESSEL_K1_P2, BESSEL_K1_Q2);
2833  double f1 = log(x);
2834  double f2 = exp(x);
2835  *k0val = f2 * (r1_0 - f1 * r2_0);
2836  *k1val = f2 * (r1_1 + f1 * r2_1) / x;
2837  }
2838  else
2839  {
2840  double y = 1./x;
2841  double r0 = ratevl_10_11(y, BESSEL_K0_P3, BESSEL_K0_Q3);
2842  double r1 = ratevl_11_10(y, BESSEL_K1_P3, BESSEL_K1_Q3);
2843  double factor = 1. / sqrt(x);
2844  *k0val = factor * r0;
2845  *k1val = factor * r1;
2846  }
2847 }
2848 
2849 static const double BESSEL_I0_P1[] = {
2850  -2.2335582639474375249e+15, -5.5050369673018427753e+14, -3.2940087627407749166e+13,
2851  -8.4925101247114157499e+11, -1.1912746104985237192e+10, -1.0313066708737980747e+08,
2852  -5.9545626019847898221e+05, -2.4125195876041896775e+03, -7.0935347449210549190e+00,
2853  -1.5453977791786851041e-02, -2.5172644670688975051e-05, -3.0517226450451067446e-08,
2854  -2.6843448573468483278e-11, -1.5982226675653184646e-14, -5.2487866627945699800e-18
2855 };
2856 static const double BESSEL_I0_Q1[] = {
2857  -2.2335582639474375245e+15, 7.8858692566751002988e+12, -1.2207067397808979846e+10,
2858  1.0377081058062166144e+07, -4.8527560179962773045e+03, 1.0
2859 };
2860 static const double BESSEL_I0_P2[] = {
2861  -2.2210262233306573296e-04, 1.3067392038106924055e-02, -4.4700805721174453923e-01,
2862  5.5674518371240761397e+00, -2.3517945679239481621e+01, 3.1611322818701131207e+01,
2863  -9.6090021968656180000e+00
2864 };
2865 static const double BESSEL_I0_Q2[] = {
2866  -5.5194330231005480228e-04, 3.2547697594819615062e-02, -1.1151759188741312645e+00,
2867  1.3982595353892851542e+01, -6.0228002066743340583e+01, 8.5539563258012929600e+01,
2868  -3.1446690275135491500e+01, 1.0
2869 };
2870 
2871 static const double BESSEL_I1_P1[] = {
2872  -1.4577180278143463643e+15, -1.7732037840791591320e+14, -6.9876779648010090070e+12,
2873  -1.3357437682275493024e+11, -1.4828267606612366099e+09, -1.0588550724769347106e+07,
2874  -5.1894091982308017540e+04, -1.8225946631657315931e+02, -4.7207090827310162436e-01,
2875  -9.1746443287817501309e-04, -1.3466829827635152875e-06, -1.4831904935994647675e-09,
2876  -1.1928788903603238754e-12, -6.5245515583151902910e-16, -1.9705291802535139930e-19
2877 };
2878 static const double BESSEL_I1_Q1[] = {
2879  -2.9154360556286927285e+15, 9.7887501377547640438e+12, -1.4386907088588283434e+10,
2880  1.1594225856856884006e+07, -5.1326864679904189920e+03, 1.0
2881 };
2882 static const double BESSEL_I1_P2[] = {
2883  1.4582087408985668208e-05, -8.9359825138577646443e-04, 2.9204895411257790122e-02,
2884  -3.4198728018058047439e-01, 1.3960118277609544334e+00, -1.9746376087200685843e+00,
2885  8.5591872901933459000e-01, -6.0437159056137599999e-02
2886 };
2887 static const double BESSEL_I1_Q2[] = {
2888  3.7510433111922824643e-05, -2.2835624489492512649e-03, 7.4212010813186530069e-02,
2889  -8.5017476463217924408e-01, 3.2593714889036996297e+00, -3.8806586721556593450e+00, 1.0
2890 };
2891 
2892 double bessel_i0(double x)
2893 {
2894  DEBUG_ENTRY( "bessel_i0()" );
2895 
2896  x = fabs(x);
2897  if( x == 0. )
2898  return 1.;
2899  if( x <= 15. )
2900  {
2901  double y = x*x;
2902  return ratevl_15_6(y, BESSEL_I0_P1, BESSEL_I0_Q1);
2903  }
2904  else
2905  {
2906  double y = 1./x - 1./15.;
2907  double r = ratevl_7_8(y, BESSEL_I0_P2, BESSEL_I0_Q2);
2908  double factor = exp(x) / sqrt(x);
2909  return factor * r;
2910  }
2911 }
2912 
2913 double bessel_i0_scaled(double x)
2914 {
2915  DEBUG_ENTRY( "bessel_i0_scaled()" );
2916 
2917  x = fabs(x);
2918  if( x == 0. )
2919  return 1.;
2920  if( x <= 15. )
2921  {
2922  double y = x*x;
2923  return exp(-x) * ratevl_15_6(y, BESSEL_I0_P1, BESSEL_I0_Q1);
2924  }
2925  else
2926  {
2927  double y = 1./x - 1./15.;
2928  double r = ratevl_7_8(y, BESSEL_I0_P2, BESSEL_I0_Q2);
2929  double factor = 1. / sqrt(x);
2930  return factor * r;
2931  }
2932 }
2933 
2934 double bessel_i1(double x)
2935 {
2936  DEBUG_ENTRY( "bessel_i1()" );
2937 
2938  double w = fabs(x);
2939  if( w == 0. )
2940  return 0.;
2941  if( w <= 15. )
2942  {
2943  double y = x*x;
2944  return x * ratevl_15_6(y, BESSEL_I1_P1, BESSEL_I1_Q1);
2945  }
2946  else
2947  {
2948  double y = 1./w - 1./15.;
2949  double r = ratevl_8_7(y, BESSEL_I1_P2, BESSEL_I1_Q2);
2950  double factor = exp(w) / sqrt(w);
2951  return ( x < 0. ) ? -factor*r : factor*r;
2952  }
2953 }
2954 
2955 double bessel_i1_scaled(double x)
2956 {
2957  DEBUG_ENTRY( "bessel_i1_scaled()" );
2958 
2959  double w = fabs(x);
2960  if( w == 0. )
2961  return 0.;
2962  if( w <= 15. )
2963  {
2964  double y = x*x;
2965  return x * exp(-w) * ratevl_15_6(y, BESSEL_I1_P1, BESSEL_I1_Q1);
2966  }
2967  else
2968  {
2969  double y = 1./w - 1./15.;
2970  double r = ratevl_8_7(y, BESSEL_I1_P2, BESSEL_I1_Q2);
2971  double factor = 1. / sqrt(w);
2972  return ( x < 0. ) ? -factor*r : factor*r;
2973  }
2974 }
2975 
2976 void bessel_i0_i1(double x, double* i0val, double* i1val)
2977 {
2978  DEBUG_ENTRY( "bessel_i0_i1()" );
2979 
2980  double w = fabs(x);
2981  if( w == 0. )
2982  {
2983  *i0val = 1.;
2984  *i1val = 0.;
2985  }
2986  if( w <= 15. )
2987  {
2988  double y = x*x;
2989  *i0val = ratevl_15_6(y, BESSEL_I0_P1, BESSEL_I0_Q1);
2990  *i1val = x * ratevl_15_6(y, BESSEL_I1_P1, BESSEL_I1_Q1);
2991  }
2992  else
2993  {
2994  double y = 1./w - 1./15.;
2995  double r0 = ratevl_7_8(y, BESSEL_I0_P2, BESSEL_I0_Q2);
2996  double r1 = ratevl_8_7(y, BESSEL_I1_P2, BESSEL_I1_Q2);
2997  double factor = exp(w) / sqrt(w);
2998  *i0val = factor * r0;
2999  *i1val = ( x < 0. ) ? -factor*r1 : factor*r1;
3000  }
3001 }
3002 
3003 void bessel_i0_i1_scaled(double x, double* i0val, double* i1val)
3004 {
3005  DEBUG_ENTRY( "bessel_i0_i1_scaled()" );
3006 
3007  double w = fabs(x);
3008  if( w == 0. )
3009  {
3010  *i0val = 1.;
3011  *i1val = 0.;
3012  }
3013  if( w <= 15. )
3014  {
3015  double y = x*x;
3016  double r0 = ratevl_15_6(y, BESSEL_I0_P1, BESSEL_I0_Q1);
3017  double r1 = x * ratevl_15_6(y, BESSEL_I1_P1, BESSEL_I1_Q1);
3018  double factor = exp(-w);
3019  *i0val = factor * r0;
3020  *i1val = factor * r1;
3021  }
3022  else
3023  {
3024  double y = 1./w - 1./15.;
3025  double r0 = ratevl_7_8(y, BESSEL_I0_P2, BESSEL_I0_Q2);
3026  double r1 = ratevl_8_7(y, BESSEL_I1_P2, BESSEL_I1_Q2);
3027  double factor = 1. / sqrt(w);
3028  *i0val = factor * r0;
3029  *i1val = ( x < 0. ) ? -factor*r1 : factor*r1;
3030  }
3031 }
3032 
3033 inline double ratevl_5_4(double x, const double a[5], const double b[4])
3034 {
3035  double x2 = x * x;
3036  double t[4];
3037  t[0] = a[4] * x2 + a[2];
3038  t[1] = a[3] * x2 + a[1];
3039  t[2] = b[3] * x2 + b[1];
3040  t[3] = b[2] * x2 + b[0];
3041  t[0] *= x2;
3042  t[1] *= x;
3043  t[2] *= x;
3044  t[0] += a[0];
3045  return (t[0] + t[1]) / (t[2] + t[3]);
3046 }
3047 
3048 inline double ratevl_6_3(double x, const double a[6], const double b[3])
3049 {
3050  double x2 = x * x;
3051  double t[4];
3052  t[0] = a[5] * x2 + a[3];
3053  t[1] = a[4] * x2 + a[2];
3054  t[2] = b[2] * x2 + b[0];
3055  t[0] *= x2;
3056  t[1] *= x2;
3057  t[3] = b[1];
3058  t[0] += a[1];
3059  t[1] += a[0];
3060  t[3] *= x;
3061  t[0] *= x;
3062  return (t[0] + t[1]) / (t[2] + t[3]);
3063 }
3064 
3065 inline double ratevl_6_4(double x, const double a[6], const double b[4])
3066 {
3067  double x2 = x * x;
3068  double t[4];
3069  t[0] = a[5] * x2 + a[3];
3070  t[1] = a[4] * x2 + a[2];
3071  t[2] = b[3] * x2 + b[1];
3072  t[3] = b[2] * x2 + b[0];
3073  t[0] *= x2;
3074  t[1] *= x2;
3075  t[2] *= x;
3076  t[0] += a[1];
3077  t[1] += a[0];
3078  t[0] *= x;
3079  return (t[0] + t[1]) / (t[2] + t[3]);
3080 }
3081 
3082 inline double ratevl_7_8(double x, const double a[7], const double b[8])
3083 {
3084  double x2 = x * x;
3085  double t[4];
3086  t[0] = a[6] * x2 + a[4];
3087  t[1] = a[5] * x2 + a[3];
3088  t[2] = b[7] * x2 + b[5];
3089  t[3] = b[6] * x2 + b[4];
3090  t[0] *= x2;
3091  t[1] *= x2;
3092  t[2] *= x2;
3093  t[3] *= x2;
3094  t[0] += a[2];
3095  t[1] += a[1];
3096  t[2] += b[3];
3097  t[3] += b[2];
3098  t[0] *= x2;
3099  t[1] *= x;
3100  t[2] *= x2;
3101  t[3] *= x2;
3102  t[0] += a[0];
3103  t[2] += b[1];
3104  t[3] += b[0];
3105  t[2] *= x;
3106  return (t[0] + t[1]) / (t[2] + t[3]);
3107 }
3108 
3109 inline double ratevl_8_7(double x, const double a[8], const double b[7])
3110 {
3111  double x2 = x * x;
3112  double t[4];
3113  t[0] = a[7] * x2 + a[5];
3114  t[1] = a[6] * x2 + a[4];
3115  t[2] = b[6] * x2 + b[4];
3116  t[3] = b[5] * x2 + b[3];
3117  t[0] *= x2;
3118  t[1] *= x2;
3119  t[2] *= x2;
3120  t[3] *= x2;
3121  t[0] += a[3];
3122  t[1] += a[2];
3123  t[2] += b[2];
3124  t[3] += b[1];
3125  t[0] *= x2;
3126  t[1] *= x2;
3127  t[2] *= x2;
3128  t[3] *= x;
3129  t[0] += a[1];
3130  t[1] += a[0];
3131  t[2] += b[0];
3132  t[0] *= x;
3133  return (t[0] + t[1]) / (t[2] + t[3]);
3134 }
3135 
3136 inline double ratevl_10_11(double x, const double a[10], const double b[11])
3137 {
3138  double x2 = x * x;
3139  double t[4];
3140  t[0] = a[9] * x2 + a[7];
3141  t[1] = a[8] * x2 + a[6];
3142  t[2] = b[10] * x2 + b[8];
3143  t[3] = b[9] * x2 + b[7];
3144  t[0] *= x2;
3145  t[1] *= x2;
3146  t[2] *= x2;
3147  t[3] *= x2;
3148  t[0] += a[5];
3149  t[1] += a[4];
3150  t[2] += b[6];
3151  t[3] += b[5];
3152  t[0] *= x2;
3153  t[1] *= x2;
3154  t[2] *= x2;
3155  t[3] *= x2;
3156  t[0] += a[3];
3157  t[1] += a[2];
3158  t[2] += b[4];
3159  t[3] += b[3];
3160  t[0] *= x2;
3161  t[1] *= x2;
3162  t[2] *= x2;
3163  t[3] *= x2;
3164  t[0] += a[1];
3165  t[1] += a[0];
3166  t[2] += b[2];
3167  t[3] += b[1];
3168  t[2] *= x2;
3169  t[0] *= x;
3170  t[3] *= x;
3171  t[2] += b[0];
3172  return (t[0] + t[1]) / (t[2] + t[3]);
3173 }
3174 
3175 inline double ratevl_11_10(double x, const double a[11], const double b[10])
3176 {
3177  double x2 = x * x;
3178  double t[4];
3179  t[0] = a[10] * x2 + a[8];
3180  t[1] = a[9] * x2 + a[7];
3181  t[2] = b[9] * x2 + b[7];
3182  t[3] = b[8] * x2 + b[6];
3183  t[0] *= x2;
3184  t[1] *= x2;
3185  t[2] *= x2;
3186  t[3] *= x2;
3187  t[0] += a[6];
3188  t[1] += a[5];
3189  t[2] += b[5];
3190  t[3] += b[4];
3191  t[0] *= x2;
3192  t[1] *= x2;
3193  t[2] *= x2;
3194  t[3] *= x2;
3195  t[0] += a[4];
3196  t[1] += a[3];
3197  t[2] += b[3];
3198  t[3] += b[2];
3199  t[0] *= x2;
3200  t[1] *= x2;
3201  t[2] *= x2;
3202  t[3] *= x2;
3203  t[0] += a[2];
3204  t[1] += a[1];
3205  t[2] += b[1];
3206  t[3] += b[0];
3207  t[0] *= x2;
3208  t[1] *= x;
3209  t[2] *= x;
3210  t[0] += a[0];
3211  return (t[0] + t[1]) / (t[2] + t[3]);
3212 }
3213 
3214 inline double ratevl_15_6(double x, const double a[15], const double b[6])
3215 {
3216  double x2 = x * x;
3217  double t[4];
3218  t[0] = a[14] * x2 + a[12];
3219  t[1] = a[13] * x2 + a[11];
3220  t[2] = b[5] * x2 + b[3];
3221  t[3] = b[4] * x2 + b[2];
3222  t[0] *= x2;
3223  t[1] *= x2;
3224  t[2] *= x2;
3225  t[3] *= x2;
3226  t[0] += a[10];
3227  t[1] += a[9];
3228  t[2] += b[1];
3229  t[3] += b[0];
3230  t[0] *= x2;
3231  t[1] *= x2;
3232  t[2] *= x;
3233  t[0] += a[8];
3234  t[1] += a[7];
3235  t[0] *= x2;
3236  t[1] *= x2;
3237  t[0] += a[6];
3238  t[1] += a[5];
3239  t[0] *= x2;
3240  t[1] *= x2;
3241  t[0] += a[4];
3242  t[1] += a[3];
3243  t[0] *= x2;
3244  t[1] *= x2;
3245  t[0] += a[2];
3246  t[1] += a[1];
3247  t[0] *= x2;
3248  t[1] *= x;
3249  t[0] += a[0];
3250  return (t[0] + t[1]) / (t[2] + t[3]);
3251 }
3252 
3253 /*******************************************************************
3254  * This marks the end of the block of code from the Boost library *
3255  *******************************************************************/
3256 
3257 /*e1 first exponential integral */
3258 double e1(double x)
3259 {
3260  DEBUG_ENTRY( "e1()" );
3261 
3262  /* computes the exponential integral E1(x),
3263  * from Abramowitz and Stegun
3264  * stops with error condition for non-positive argument,
3265  * returns zero in large x limit
3266  * */
3267 
3268  /* error - does not accept non-positive arguments */
3269  if( x <= 0. )
3270  {
3271  fprintf( ioQQQ, " DISASTER invalid argument x=%g for function e1, requires x > 0\n", x );
3273  }
3274  else if( x < 1. )
3275  {
3276  // Abramowitz and Stegun Sect. 5.1.53
3277  // abs. accuracy better than 2e-7 for the polynomial excluding the log() term
3278  const double a[6]= {-.57721566,.99999193,-.24991055,
3279  .05519968,-.00976004,.00107857};
3280  return ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
3281  }
3282  else
3283  {
3284  // Abramowitz and Stegun Sect. 5.1.56
3285  // for the rational function fit to x*exp(x)*E1(x) this holds:
3286  // abs. accuracy better than 2e-8
3287  const double a[4]= {8.5733287401,18.0590169730,8.6347608925,.2677737343};
3288  const double b[4]= {9.5733223454,25.6329561486,21.0996530827,3.9584969228};
3289  double top = (((x + a[0])*x + a[1])*x + a[2])*x + a[3];
3290  double bot = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
3291  return top/bot/x*exp(-x);
3292  }
3293 }
3294 
3295 /*e1_scaled is exp(x)*e1(x) */
3296 double e1_scaled(double x)
3297 {
3298  DEBUG_ENTRY( "e1_scaled()" );
3299 
3300  /* computes the function exp(x)*E1(x),
3301  * stops with error condition for non-positive argument,
3302  * */
3303 
3304  /* error - does not accept non-positive arguments */
3305  if( x <= 0. )
3306  {
3307  fprintf( ioQQQ, " DISASTER invalid argument x=%g for function e1_scaled, requires x > 0\n", x );
3309  }
3310  else if( x < 1. )
3311  {
3312  return exp(x)*e1(x);
3313  }
3314  else
3315  {
3316  // Abramowitz and Stegun Sect. 5.1.56
3317  // for the rational function fit to x*exp(x)*E1(x) this holds:
3318  // abs. accuracy better than 2e-8
3319  const double a[4]= {8.5733287401,18.0590169730,8.6347608925,.2677737343};
3320  const double b[4]= {9.5733223454,25.6329561486,21.0996530827,3.9584969228};
3321  double top = (((x + a[0])*x + a[1])*x + a[2])*x + a[3];
3322  double bot = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
3323  return top/bot/x;
3324  }
3325 }
3326 
3327 /*e2 second exponential integral */
3328 double e2(double x)
3329 {
3330  DEBUG_ENTRY( "e2()" );
3331 
3332  // the rel. accuracy of e2() is better than 4e-9 everywhere it doesn't underflow to zero
3333  if( x < 0. )
3334  {
3335  fprintf( ioQQQ, " DISASTER invalid argument x=%g for function e2, requires x >= 0\n", x );
3337  }
3338  else if( x == 0. )
3339  {
3340  return 1.;
3341  }
3342  else if( x < 0.45 )
3343  {
3344  // for the polynomial fit with the x*log(x) term excluded this holds:
3345  // rel. accuracy better than 3.7e-9, abs. accuracy better than 2.6e-9
3346  const double a[6] = { 1.0000000000, -0.4227844628, -0.4999962685,
3347  0.0832964152, -0.0137254444, 0.0017460335 };
3348  return ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] + x*log(x);
3349  }
3350  else if( x < 1. )
3351  {
3352  // for the polynomial fit with the x*log(x) term excluded this holds:
3353  // rel. accuracy better than 1.4e-9, abs. accuracy better than 4.3e-10
3354  const double a[7] = { 1.0000006037, -0.4227920037, -0.4999593292, 0.0832162574,
3355  -0.0136902608, 0.0018824593, -0.0001622201 };
3356  return (((((a[6]*x + a[5])*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] + x*log(x);
3357  }
3358  else if( x < 2.6666667 )
3359  {
3360  // for the rational function fit to x*exp(x)*E2(x) this holds:
3361  // rel. accuracy better than 4.0e-9, abs. accuracy better than 2.5e-9
3362  const double a[4] = { 8.5977992972, 17.4635101763, 7.5697246936, 0.0512652659 };
3363  const double b[4] = { 10.5974966324, 32.6715286491, 33.0501272716, 8.6019987326 };
3364  double y = 1./x;
3365  double top = (((a[3]*y + a[2])*y + a[1])*y + a[0])*y + 1.;
3366  double bot = (((b[3]*y + b[2])*y + b[1])*y + b[0])*y + 1.;
3367  return top/bot*y*exp(-x);
3368  }
3369  else if( x < 21. )
3370  {
3371  // for the rational function fit to x*exp(x)*E2(x) this holds:
3372  // rel. accuracy better than 3.4e-9, abs. accuracy better than 3.2e-9
3373  const double a[4] = { 13.3703104742, 46.8560764900, 39.5905932936, 0.9608837426 };
3374  const double b[4] = { 15.3703105373, 71.5967206495, 114.5581563870, 49.5883983926 };
3375  double y = 1./x;
3376  double top = (((a[3]*y + a[2])*y + a[1])*y + a[0])*y + 1.;
3377  double bot = (((b[3]*y + b[2])*y + b[1])*y + b[0])*y + 1.;
3378  return top/bot*y*exp(-x);
3379  }
3380  else if( x < 752. )
3381  {
3382  // for the rational function fit to x*exp(x)*E2(x) this holds:
3383  // rel. accuracy better than 2.3e-9, abs. accuracy better than 2.1e-9
3384  const double a[3] = { -1.7217469460, -40.8941038520, -13.3258180489 };
3385  const double b[3] = { 0.2782514490, -46.3371070179, -83.7227541235 };
3386  double y = 1./x;
3387  double top = ((a[2]*y + a[1])*y + a[0])*y + 1.;
3388  double bot = ((b[2]*y + b[1])*y + b[0])*y + 1.;
3389  return top/bot*y*exp(-x);
3390  }
3391  else
3392  {
3393  // result would underflow to zero anyway...
3394  return 0.;
3395  }
3396 }
3397 
3398 inline double expn2_scaled(double x)
3399 {
3400  ASSERT(x >= 0.0);
3401  if (x == 0.)
3402  return 1.;
3403 
3404  double z = 1.0/x;
3405  double val = expn(2,z);
3406  if (val > 0.0)
3407  {
3408  val *= exp(z)*z;
3409  }
3410  else
3411  {
3412  val = 1.-2.0*x*(1-3.0*x*(1-4.0*x));
3413  }
3414  return val;
3415 }
3416 
3417 void chbfit(double a, double b, vector<double>& c, double (*func)(double))
3418 {
3419  const size_t n = c.size();
3420  vector<double> fv(n);
3421  const double fa = PI/n;
3422  for (size_t k=0; k<n; ++k)
3423  {
3424  double frac = 0.5*(1.+cos(fa*(k+0.5)));
3425  fv[k] = func(a+(b-a)*frac);
3426  }
3427  for (size_t j=0; j<n; ++j)
3428  {
3429  double s = 0.0;
3430  for (size_t k=0; k<n; ++k)
3431  {
3432  s += fv[k]*cos(fa*(k+0.5)*j);
3433  }
3434  c[n-1-j] = s*2./n;
3435  }
3436 }
3437 
3439 {
3440  vector<double> c(50);
3441  double zmax=200.0;
3442  chbfit(0.,1.,c,expn2_scaled);
3443  fprintf(ioQQQ,"Chebyshev coefficients:");
3444  for (size_t i=0; i<c.size(); ++i)
3445  {
3446  fprintf(ioQQQ,"%.12e,",c[c.size()-1-i]);
3447  }
3448  fprintf(ioQQQ,"\n");
3449  if (0)
3450  {
3451  int nterm = 10; // Number of coefficients to use in Chebyshev fit
3452  for (double z=1.0; z<zmax; z *= 1.1)
3453  {
3454  double chfrac = 4/z-2.;
3455  double expnval = expn2_scaled(1./z);
3456  double chbval = chbevl(chfrac,&c[0]+(c.size()-nterm),nterm);
3457  fprintf(ioQQQ,"%.12e e2 %.12e cheb %.12e error %.4e\n",z,expnval,
3458  chbval,2*(chbval-expnval)/(chbval+expnval));
3459  }
3460  }
3461 }
3462 
3463 // Evaluate the Gegenbauer (aka ultraspherical) polynomial C_n^(alpha)(x)
3464 double gegenbauer(long n, double al, double x)
3465 {
3466  DEBUG_ENTRY( "gegenbauer()" );
3467 
3468  // using recurrence relation
3469  double cpp = 1.0;
3470  if (n == 0)
3471  return cpp;
3472  double c = 2.*al*x;
3473  double afac = 2*(al-1.);
3474  for (long nn = 2; nn<=n; ++nn)
3475  {
3476  double cp = c;
3477  c = (x*(2*nn+afac)*cp-(nn+afac)*cpp)/double(nn);
3478  cpp = cp;
3479  }
3480  return c;
3481 }
3482 
3483 inline double sg(long S)
3484 {
3485  // S is in 2*J notation, this routine returns pow(-1.0,S/2)
3486  if( abs(S%2) == 1 )
3487  TotalInsanity();
3488  return ( S%4 == 0 ) ? 1.0 : -1.0;
3489 }
3490 
3491 // Wigner 6J symbols evaluation, original routine Fortran 6j for autostructure
3492 /* the six quantum number arguments have twice their physical value; */
3493 /* scaled factorials must be supplied by fc2_scl(2*n+1) = n! / 16**n */
3494 double sjs( long j1, long j2, long j3, long l1, long l2, long l3 )
3495 {
3496  DEBUG_ENTRY( "sjs()" );
3497 
3498  if( !Triangle2( j1, j2, j3 ) || !Triangle2( l1, l2, j3 ) ||
3499  !Triangle2( j1, l2, l3 ) || !Triangle2( j2, l1, l3 ) )
3500  return 0.;
3501 
3502  long ij0 = j1+j2+j3+2;
3503  long ij1 = j1+l2+l3+2;
3504  long ij2 = l1+j2+l3+2;
3505  long ij3 = l1+l2+j3+2;
3506 
3507  /*some corrections have been done here to translate from original fortran routine
3508  * as fc2_scl ranks from 0 to num-1 */
3509  long iwmin=max(max(max(ij0,ij1),ij2),ij3)+1;
3510 
3511  long id1 = ij0+ij1-j1-j1;
3512  long id2 = ij0+ij2-j2-j2;
3513  long id3 = ij0+ij3-j3-j3;
3514 
3515  long iwmax = min(min(id1,id2),id3)-1;
3516 
3517  double omega=0.;
3518  if( iwmax >= iwmin )
3519  {
3520  for ( long iw=iwmin; iw <= iwmax; iw += 2 )
3521  {
3522  omega += sg(iw+1)*fc2_scl(iw)/
3523  (fc2_scl(id1-iw)*fc2_scl(id2-iw)*fc2_scl(id3-iw)*
3524  fc2_scl(iw-ij0)*fc2_scl(iw-ij1)*fc2_scl(iw-ij2)*fc2_scl(iw-ij3));
3525  }
3526 
3527  ij0++;
3528  ij1++;
3529  ij2++;
3530  ij3++;
3531 
3532  omega *= sqrt((fc2_scl(id1-ij0)*fc2_scl(id2-ij0)*fc2_scl(id3-ij0)/fc2_scl(ij0))*
3533  (fc2_scl(id1-ij1)*fc2_scl(id2-ij1)*fc2_scl(id3-ij1)/fc2_scl(ij1))*
3534  (fc2_scl(id1-ij2)*fc2_scl(id2-ij2)*fc2_scl(id3-ij2)/fc2_scl(ij2))*
3535  (fc2_scl(id1-ij3)*fc2_scl(id2-ij3)*fc2_scl(id3-ij3)/fc2_scl(ij3)))/16;
3536  }
3537  return omega;
3538 }
3539 
3540 inline double fc2(long n2)
3541 {
3542  // return n!, where n2 = 2*n
3543 
3544  if( abs(n2%2) == 1 )
3545  TotalInsanity();
3546  return factorial(n2/2);
3547 }
3548 
3549 inline double Delta(long j1, long j2, long j3)
3550 {
3551  // this is the triangle coefficient taken from Eq. 5 of
3552  // Latha K.V.P., Angom D., Das B.P., 2008, arXiv:0805.2723v1
3553  // (note that one factorial sign is missing in the paper)
3554 
3555  return fc2(j1+j2-j3)*fc2(j1-j2+j3)*fc2(-j1+j2+j3)/fc2(j1+j2+j3+2);
3556 }
3557 
3558 double SixJFull(long j1, long j2, long j3, long j4, long j5, long j6)
3559 {
3560  DEBUG_ENTRY( "SixJFull()" );
3561 
3562  // evaluate a 6j symbol using the Racah formula, this version is valid
3563  // for arbitrary large arguments, but is slower... it may also suffer
3564  // from cancellation errors for very large arguments.
3565 
3566  // Written by Peter van Hoof.
3567 
3568  // The expression below is taken from Eq. 4 of
3569  // Latha K.V.P., Angom D., Das B.P., 2008, arXiv:0805.2723v1
3570 
3571  if( !Triangle2( j1, j2, j3 ) || !Triangle2( j4, j5, j3 ) ||
3572  !Triangle2( j1, j5, j6 ) || !Triangle2( j2, j4, j6 ) )
3573  return 0.;
3574 
3575  long tlo = max(max(max(j1+j2+j3,j1+j5+j6),j4+j2+j6),j4+j5+j3);
3576  long thi = min(min(j1+j2+j4+j5,j2+j3+j5+j6),j3+j1+j6+j4);
3577 
3578  double sum = 0.;
3579  double term = sg(tlo)*fc2(tlo+2)/
3580  (fc2(tlo-j1-j2-j3)*fc2(tlo-j1-j5-j6)*fc2(tlo-j4-j2-j6)*fc2(tlo-j4-j5-j3)*
3581  fc2(j1+j2+j4+j5-tlo)*fc2(j2+j3+j5+j6-tlo)*fc2(j3+j1+j6+j4-tlo));
3582  for( long t=tlo; t <= thi; t += 2 )
3583  {
3584  if (t != tlo)
3585  term *= -(t+2)
3586  /double((t-j1-j2-j3)*(t-j1-j5-j6)*(t-j4-j2-j6)*(t-j4-j5-j3))
3587  *double((j1+j2+j4+j5+2-t)*(j2+j3+j5+j6+2-t)*(j1+j3+j4+j6+2-t));
3588  sum += term;
3589  }
3590 
3591  if( sum != 0. )
3592  sum *= sqrt( Delta(j1,j2,j3)*Delta(j1,j5,j6)*Delta(j4,j2,j6)*Delta(j4,j5,j3) );
3593 
3594  return sum;
3595 }
3596 
3597 inline double frac(double d)
3598 {
3599  return d-floor(d);
3600 }
3601 
3602 /* Recursive routine for SixJ coefficients, based on the routine of
3603  * >>refer K. Schulten & R.G. Gordon, Comput. Phys. Commun. 11(1976)269
3604  *
3605  * Included with permission of the author, please cite the original
3606  * publications if used independently of Cloudy
3607  */
3608 void rec6j(double *sixcof, double l2, double l3,
3609  double l4, double l5, double l6, double *l1min,
3610  double *l1max, double *lmatch, long ndim, long *ier)
3611 {
3612  DEBUG_ENTRY( "rec6j()" );
3613 
3614  double eps = .01;
3615  double tiny = 1e-300;
3616  double srtiny = sqrt(tiny);
3617  double huge = 1e300;
3618  double srhuge = sqrt(huge);
3619 
3620 /* J1-RECURSION OF 6J-COEFFICIENTS */
3621 
3622 /* */
3623 /* ROUTINE TO GENERATE THE SET OF 6J-COEFFICIENTS L1 L2 L3 */
3624 /* L4 L5 L6 */
3625 /* BY RECURSION FROM L1MIN = MAX0(/L2-L3/,/L5-L6/) */
3626 /* TO L1MAX = MIN0( L2+L3 , L5+L6 ) */
3627 /* */
3628 /* THE RESULTING 6J-COEFFICIENTS ARE STORED AS SIXCOF(L1-L1MIN). */
3629 /* */
3630 /* FOR A DISCUSSION OF THE RECURSION EQUATION USED SEE K. SCHULTEN */
3631 /* AND R.G. GORDON, J. MATH. PHYS. 16, 1961-1970 (1975), IBID. 16, */
3632 /* 1971-1988 (1975) */
3633 /* FOR THE SAKE OF NUMERICAL STABILITY THE RECURSION WILL PROCEED */
3634 /* SIMULTANEOUSLY FORWARD AND BACKWARDS, STARTING AT L1MIN AND */
3635 /* L1MAX, RESPECTIVELY. */
3636 /* LMATCH IS THE L1-VALUE AT WHICH FORWARD AND BACKWARD RECURSION */
3637 /* ARE MATCHED. ITS VALUE WILL BE RETURNED, THOUGH IT IS NOT OF */
3638 /* ACTUAL USE. */
3639 /* */
3640 /* NDIM IS THE LENGTH OF THE ARRAY SIXCOF. */
3641 /* */
3642 /* IER IS SET TO -1 IF EITHER NO 6J-COEFFICIENT SATISFIES TRIANGULAR */
3643 /* CONDITION OR IF L2+L3+L5+L6 OR L2+L4+L6 NOT */
3644 /* INTEGER, IN WHICH CASE ALL 6J-COEFFICIENTS */
3645 /* VANISH. */
3646 /* IER IS SET TO -2 IF NUMBER OF POSSIBLE 6J-COEFFICIENTS EXCEEDS NDIM */
3647 /* */
3648 /* TINY SHOULD BE SET CLOSE TO THE SMALLEST POSITIVE FLOATING POINT */
3649 /* NUMBER WHICH IS REPRESENTABLE ON THE COMPUTER. SRTINY IS SQUARE */
3650 /* ROOT OF TINY . */
3651 /* */
3652 /* */
3653 /* HUGE SHOULD BE SET CLOSE TO LARGEST POSITIVE FLOATING POINT */
3654 /* NUMBER WHICH IS REPRESENTABLE ON THE COMPUTER. SRHUGE IS */
3655 /* SQUARE ROOT OF HUGE . */
3656 /* */
3657  *lmatch = 0.;
3658 
3659  /* CHECK IF 6J-COEFFICIENTS OBEY SELECTION RULES */
3660 
3661  /* LIMITS FOR L1 */
3662  *l1min = max(fabs(l2 - l3),fabs(l5 - l6));
3663  *l1max = min(l2 + l3,l5 + l6);
3664 
3665  if ( (frac(l2 + l3 + l5 + l6 + eps) >= 2*eps) ||
3666  (frac(l4 + l2 + l6 + eps) >= 2*eps) ||
3667  (l4 + l2 - l6 < 0.) ||
3668  (l4 - l2 + l6 < 0.) ||
3669  (-(l4) + l2 + l6 < 0.) ||
3670  (l4 + l3 - l5 < 0.) ||
3671  (l4 - l3 + l5 < 0.) ||
3672  (-(l4) + l3 + l5 < 0.) ||
3673  (*l1min >= *l1max + eps)
3674  ) {
3675  /* THIS IS REACHED IF TRIANGULAR CONDITION NOT SATISFIED OR IF
3676  * L2+L3+L5+L6 OR L2+L4+L6 NOT INTEGER */
3677  (*ier) = -1;
3678  if (0)
3679  fprintf(ioQQQ," 6J-COEFFICIENTS L1 %7.1f%7.1f "
3680  "DO NOT SATISFY TRIANGULAR CONDITIONS OR\n"
3681  " %7.1f%7.1f%7.1f "
3682  "L2+L3+L5+L6 OR L2+L4+L6 NOT INTEGER\n",
3683  l2,l3,l4,l5,l6);
3684  return;
3685  }
3686 
3687  int nfin = int (*l1max - *l1min + 1. + eps);
3688  if (ndim < nfin) {
3689  /* THIS IS REACHED IF ARRAY SIXCOF NOT LARGE ENOUGH TO HOLD ALL */
3690  /* 6J - COEFFICIENTS REQUIRED */
3691  (*ier) = -2;
3692  if (0)
3693  fprintf(ioQQQ," 6J-COEFFICIENTS L1%7.1f%7.1f\n"
3694  " %7.1f%7.1f%7.1f"
3695  " EXCEED STORAGE PROVIDED (%4d,%4ld)\n",
3696  l2,l3,l4,l5,l6,nfin,ndim);
3697  return;
3698  }
3699 
3700  int sign2 = 0x1 & int (l2 + l3 + l5 + l6 + eps) ? -1 : 1;
3701  if (nfin == 1) {
3702  /* THIS IS REACHED IN CASE THAT L1 CAN TAKE ONLY ONE VALUE */
3703  (*ier) = 0;
3704  sixcof[0] = sign2 / sqrt((2*(*l1min) + 1.) * ( 2*l4 + 1.));
3705  return;
3706  }
3707 
3708  /* START OF FORWARD RECURSION */
3709  double l1 = *l1min;
3710  sixcof[0] = srtiny;
3711  double sum1 = (2*l1 + 1.) * tiny;
3712  (*ier) = 0;
3713 
3714  int i, n;
3715  double c1old=0., oldfac=0., sumfor, x;
3716 
3717  int lstep;
3718  for (lstep=1; ; ++lstep) {
3719  l1 += 1.;
3720  double a1 = (l1 + l2 + l3 + 1.) * (l1 - l2 + l3) *
3721  (l1 + l2 - l3) * (-l1 + l2 + l3 + 1.);
3722  double a2 = (l1 + l5 + l6 + 1.) * (l1 - l5 + l6) *
3723  (l1 + l5 - l6) * (-l1 + l5 + l6 + 1.);
3724  double newfac = sqrt(a1 * a2);
3725  double c1, denom=0.;
3726  if (l1 >= 1. + eps) {
3727  double dv = 2. * (l2 * (l2 + 1.) * l5 * (l5 + 1.) +
3728  l3 * (l3 + 1.) * l6 * (l6 + 1.) -
3729  l1 * (l1 - 1.) * l4 * (l4 + 1.)) -
3730  (l2 * (l2 + 1.) + l3 * (l3 + 1.) - l1 * (l1 - 1.)) *
3731  (l5 * (l5 + 1.) + l6 * (l6 + 1.) - l1 * (l1 - 1.));
3732  denom = (l1 - 1.) * newfac;
3733  c1 = -(2*l1 - 1.) * dv / denom;
3734  }
3735  else
3736  {
3737  /* IF L1 = 1 (L1 - 1) HAS TO BE FACTORED OUT OF DV, HENCE */
3738  c1 = -2. * (l2 * (l2 + 1.) + l5 * (l5 + 1.) - l4 * (l4 + 1.)) /
3739  newfac;
3740  }
3741 
3742  if (lstep == 1) {
3743  /* IF L1 = L1MIN + 1 THE THIRD TERM IN RECURSION EQUATION
3744  * VANISHES */
3745  x = srtiny * c1;
3746  sixcof[1] = x;
3747  sum1 += tiny * (2*l1 + 1.) * c1 * c1;
3748  if (lstep+1 == nfin) {
3749  break;
3750  }
3751  }
3752  else
3753  {
3754  /* RECURSION TO THE NEXT 6J - COEFFICIENT X */
3755  double c2 = -l1 * oldfac / denom;
3756  x = c1 * sixcof[lstep-1] + c2 * sixcof[lstep - 2];
3757  sixcof[lstep] = x;
3758  sumfor = sum1;
3759  sum1 += (2*l1 + 1.) * x * x;
3760 
3761  /* AS LONG AS THE COEFFICIENT /C1/ IS DECREASING THE
3762  * RECURSION PROCEEDS TOWARDS INCREASING 6J-VALUES AND,
3763  * HENCE, IS NUMERICALLY STABLE. ONCE AN INCREASE OF /C1/
3764  * IS DETECTED, THE RECURSION DIRECTION IS REVERSED. */
3765  if (lstep+1 == nfin || c1old <= fabs(c1)) {
3766  break;
3767  }
3768 
3769  /* SEE IF LAST UNNORMALIZED 6J-COEFFICIENT EXCEEDS SRHUGE */
3770  if (fabs(x) >= srhuge) {
3771  /* THIS IS REACHED IF LAST 6J-COEFFICIENT LARGER THAN
3772  * SRHUGE SO THAT THE RECURSION SERIES SIXCOF(0),
3773  * ... ,SIXCOF(LSTEP) HAS TO BE RESCALED TO PREVENT
3774  * OVERFLOW */
3775  ++(*ier);
3776  for (i = 0; i < lstep+1; ++i) {
3777  sixcof[i] /= srhuge;
3778  }
3779  sum1 /= huge;
3780  sumfor /= huge;
3781  x /= srhuge;
3782  }
3783  }
3784  c1old = fabs(c1);
3785  oldfac = newfac;
3786  }
3787 
3788  double sumuni;
3789  if (lstep+1 == nfin) {
3790  sumuni = sum1;
3791  }
3792  else
3793  {
3794  /* KEEP THREE 6J-COEFFICIENTS AROUND LMATCH FOR COMPARISION LATER
3795  * WITH BACKWARD RECURSION. */
3796  *lmatch = l1 - 1;
3797  double x1 = x;
3798  double x2 = sixcof[lstep-1];
3799  double x3 = sixcof[lstep-2];
3800  /* STARTING BACKWARD RECURSION FROM L1MAX TAKING NSTEP1+1 STEPS, SO
3801  * THAT FORWARD AND BACKWARD RECURSION OVERLAP AT THE THREE POINTS
3802  * L1 = LMATCH+1, LMATCH, LMATCH-1. */
3803  int nlim = lstep-1;
3804  l1 = *l1max;
3805  sixcof[nfin-1] = srtiny;
3806  double sum2 = (2*l1 + 1.) * tiny;
3807  l1 += 2.;
3808  double sumbac=0.;
3809  double y;
3810  for(lstep=1;;++lstep) {
3811  l1 -= 1.;
3812  double a1s = (l1 + l2 + l3) * (l1 - l2 + l3 - 1.) *
3813  (l1 + l2 - l3 - 1.) * (-l1 + l2 + l3 + 2.);
3814  double a2s = (l1 + l5 + l6) * (l1 - l5 + l6 - 1.) *
3815  (l1 + l5 - l6 - 1.) * (-l1 + l5 + l6 + 2.);
3816  double newfac = sqrt(a1s * a2s);
3817  double dv = 2. * (l2 * (l2 + 1.) * l5 * (l5 + 1.) + l3 * (l3 + 1.) *
3818  l6 * (l6 + 1.) - l1 * (l1 - 1.) * l4 * (l4 + 1.)) -
3819  (l2 * (l2 + 1.) + l3 * (l3 + 1.) - l1 * (l1 - 1.)) *
3820  (l5 * (l5 + 1.) + l6 * (l6 + 1.) - l1 * (l1 - 1.));
3821  double denom = l1 * newfac;
3822  double c1 = -(2*l1 - 1.) * dv / denom;
3823  if (lstep == 1) {
3824  /* IF L1 = L1MAX + 1 THE THIRD TERM IN THE RECURSION
3825  * EQUATION VANISHES */
3826  y = srtiny * c1;
3827  sixcof[nfin-2] = y;
3828  if (lstep == nfin-nlim) {
3829  break;
3830  }
3831  sumbac = sum2;
3832  sum2 += (2*l1 - 3.) * c1 * c1 * tiny;
3833  }
3834  else
3835  {
3836  double c2 = -(l1 - 1.) * oldfac / denom;
3837  /* RECURSION TO THE NEXT 6J - COEFFICIENT Y */
3838  y = c1 * sixcof[nfin-lstep] + c2 * sixcof[nfin-lstep+1];
3839  if (lstep == nfin-nlim) {
3840  break;
3841  }
3842  sixcof[nfin-lstep-1] = y;
3843  sumbac = sum2;
3844  sum2 += (2*l1 - 3.) * y * y;
3845  /* SEE IF LAST UNNORMALIZED 6J-COEFFICIENT EXCEEDS
3846  * SRHUGE */
3847  if (fabs(y) >= srhuge) {
3848  /* THIS IS REACHED IF LAST 6J-COEFFICIENT LARGER THAN
3849  * SRHUGE SO THAT THE RECURSION SERIES SIXCOF(NFIN-1),
3850  * ... ,SIXCOF(NFIN-LSTEP-1) HAS TO BE RESCALED TO
3851  * PREVENT OVERFLOW */
3852  ++(*ier);
3853  for (i = 0; i < lstep+1; ++i) {
3854  int index = nfin-i-1;
3855  sixcof[index] /= srhuge;
3856  }
3857  sumbac /= huge;
3858  sum2 /= huge;
3859  }
3860  }
3861  oldfac = newfac;
3862  }
3863 
3864  /* THE FORWARD RECURSION 6J-COEFFICIENTS X1, X2, X3 ARE TO BE
3865  * MATCHED WITH THE CORRESPONDING BACKWARD RECURSION VALUES Y1, Y2,
3866  * Y3. */
3867  double y3 = y;
3868  double y2 = sixcof[nfin-lstep];
3869  double y1 = sixcof[nfin-lstep+1];
3870  /* DETERMINE NOW RATIO SUCH THAT YI = RATIO * XI (I=1,2,3) HOLDS
3871  * WITH MINIMAL ERROR. */
3872  double ratio = (x1 * y1 + x2 * y2 + x3 * y3) /
3873  (x1 * x1 + x2 * x2 + x3 * x3);
3874  if (fabs(ratio) >= 1.) {
3875  for (n = 0; n < nlim; ++n) {
3876  sixcof[n] = ratio * sixcof[n];
3877  }
3878  sumuni = ratio * ratio * sumfor + sumbac;
3879  }
3880  else
3881  {
3882  ratio = 1. / ratio;
3883  for (n = nlim; n < nfin; ++n) {
3884  sixcof[n] = ratio * sixcof[n];
3885  }
3886  sumuni = sumfor + ratio * ratio * sumbac;
3887  }
3888  }
3889 
3890  /* NORMALIZE 6J-COEFFICIENTS */
3891  /* SIGN CONVENTION FOR LAST 6J-COEFF. DETERMINES OVERALL PHASE */
3892  int sign1 = sixcof[nfin-1] >= 0. ? 1 : -1;
3893  double cnorm = sign1*sign2 / sqrt((2*l4 + 1.) * sumuni);
3894 
3895  for (n = 0; n < nfin; ++n) {
3896  sixcof[n] = cnorm * sixcof[n];
3897  }
3898 }
3899 /* End of rec6j */
3900 
3901 /* >>refer Mersenne Twister Matsumoto, M. & Nishimura, T. 1998, ACM Transactions on Modeling
3902  * >>refercon and Computer Simulation (TOMACS), 8, 1998 */
3903 
3904 /********************************************************************
3905  * This copyright notice must accompany the following block of code *
3906  *******************************************************************/
3907 
3908 /*
3909  A C-program for MT19937, with initialization improved 2002/2/10.
3910  Coded by Takuji Nishimura and Makoto Matsumoto.
3911  This is a faster version by taking Shawn Cokus's optimization,
3912  Matthe Bellew's simplification, Isaku Wada's real version.
3913 
3914  Before using, initialize the state by using init_genrand(seed)
3915  or init_by_array(init_key, key_length).
3916 
3917  Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
3918  All rights reserved.
3919 
3920  Redistribution and use in source and binary forms, with or without
3921  modification, are permitted provided that the following conditions
3922  are met:
3923 
3924  1. Redistributions of source code must retain the above copyright
3925  notice, this list of conditions and the following disclaimer.
3926 
3927  2. Redistributions in binary form must reproduce the above copyright
3928  notice, this list of conditions and the following disclaimer in the
3929  documentation and/or other materials provided with the distribution.
3930 
3931  3. The names of its contributors may not be used to endorse or promote
3932  products derived from this software without specific prior written
3933  permission.
3934 
3935  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
3936  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
3937  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
3938  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
3939  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
3940  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
3941  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
3942  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
3943  LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
3944  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
3945  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
3946 
3947 
3948  Any feedback is very welcome.
3949  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
3950  email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
3951 */
3952 
3953 /* Period parameters */
3954 static const int N = 624;
3955 static const int M = 397;
3956 static const unsigned long MATRIX_A = 0x9908b0dfUL; /* constant vector a */
3957 static const unsigned long UMASK = 0x80000000UL; /* most significant w-r bits */
3958 static const unsigned long LMASK = 0x7fffffffUL; /* least significant r bits */
3959 inline unsigned long MIXBITS(unsigned long u, unsigned long v)
3960 {
3961  return (u & UMASK) | (v & LMASK);
3962 }
3963 inline unsigned long TWIST(unsigned long u, unsigned long v)
3964 {
3965  return (MIXBITS(u,v) >> 1) ^ (v&1UL ? MATRIX_A : 0UL);
3966 }
3967 
3968 static unsigned long state[N]; /* the array for the state vector */
3969 static int nleft = 1;
3970 static int initf = 0;
3971 static unsigned long *nexxt;
3972 
3973 /* initializes state[N] with a seed */
3974 void init_genrand(unsigned long s)
3975 {
3976  int j;
3977  state[0]= s & 0xffffffffUL;
3978  for (j=1; j<N; j++) {
3979  state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
3980  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
3981  /* In the previous versions, MSBs of the seed affect */
3982  /* only MSBs of the array state[]. */
3983  /* 2002/01/09 modified by Makoto Matsumoto */
3984  state[j] &= 0xffffffffUL; /* for >32 bit machines */
3985  }
3986  nleft = 1; initf = 1;
3987 }
3988 
3989 /* initialize by an array with array-length */
3990 /* init_key is the array for initializing keys */
3991 /* key_length is its length */
3992 /* slight change for C++, 2004/2/26 */
3993 void init_by_array(unsigned long init_key[], int key_length)
3994 {
3995  int i, j, k;
3996  init_genrand(19650218UL);
3997  i=1; j=0;
3998  k = (N>key_length ? N : key_length);
3999  for (; k; k--) {
4000  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
4001  + init_key[j] + j; /* non linear */
4002  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
4003  i++; j++;
4004  if (i>=N) { state[0] = state[N-1]; i=1; }
4005  if (j>=key_length) j=0;
4006  }
4007  for (k=N-1; k; k--) {
4008  state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
4009  - i; /* non linear */
4010  state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
4011  i++;
4012  if (i>=N) { state[0] = state[N-1]; i=1; }
4013  }
4014 
4015  state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
4016  nleft = 1; initf = 1;
4017 }
4018 
4019 static void next_state()
4020 {
4021  unsigned long *p=state;
4022  int j;
4023 
4024  /* if init_genrand() has not been called, */
4025  /* a default initial seed is used */
4026  if (initf==0) init_genrand(5489UL);
4027 
4028  nleft = N;
4029  nexxt = state;
4030 
4031  for (j=N-M+1; --j; p++)
4032  *p = p[M] ^ TWIST(p[0], p[1]);
4033 
4034  for (j=M; --j; p++)
4035  *p = p[M-N] ^ TWIST(p[0], p[1]);
4036 
4037  *p = p[M-N] ^ TWIST(p[0], state[0]);
4038 }
4039 
4040 /* generates a random number on [0,0xffffffff]-interval */
4041 unsigned long genrand_int32()
4042 {
4043  unsigned long y;
4044 
4045  if (--nleft == 0) next_state();
4046  y = *nexxt++;
4047 
4048  /* Tempering */
4049  y ^= (y >> 11);
4050  y ^= (y << 7) & 0x9d2c5680UL;
4051  y ^= (y << 15) & 0xefc60000UL;
4052  y ^= (y >> 18);
4053 
4054  return y;
4055 }
4056 
4057 /* generates a random number on [0,0x7fffffff]-interval */
4059 {
4060  unsigned long y;
4061 
4062  if (--nleft == 0) next_state();
4063  y = *nexxt++;
4064 
4065  /* Tempering */
4066  y ^= (y >> 11);
4067  y ^= (y << 7) & 0x9d2c5680UL;
4068  y ^= (y << 15) & 0xefc60000UL;
4069  y ^= (y >> 18);
4070 
4071  return (long)(y>>1);
4072 }
4073 
4074 /* generates a random number on [0,1]-real-interval */
4076 {
4077  unsigned long y;
4078 
4079  if (--nleft == 0) next_state();
4080  y = *nexxt++;
4081 
4082  /* Tempering */
4083  y ^= (y >> 11);
4084  y ^= (y << 7) & 0x9d2c5680UL;
4085  y ^= (y << 15) & 0xefc60000UL;
4086  y ^= (y >> 18);
4087 
4088  return (double)y * (1.0/4294967295.0);
4089  /* divided by 2^32-1 */
4090 }
4091 
4092 /* generates a random number on [0,1)-real-interval */
4094 {
4095  unsigned long y;
4096 
4097  if (--nleft == 0) next_state();
4098  y = *nexxt++;
4099 
4100  /* Tempering */
4101  y ^= (y >> 11);
4102  y ^= (y << 7) & 0x9d2c5680UL;
4103  y ^= (y << 15) & 0xefc60000UL;
4104  y ^= (y >> 18);
4105 
4106  return (double)y * (1.0/4294967296.0);
4107  /* divided by 2^32 */
4108 }
4109 
4110 /* generates a random number on (0,1)-real-interval */
4112 {
4113  unsigned long y;
4114 
4115  if (--nleft == 0) next_state();
4116  y = *nexxt++;
4117 
4118  /* Tempering */
4119  y ^= (y >> 11);
4120  y ^= (y << 7) & 0x9d2c5680UL;
4121  y ^= (y << 15) & 0xefc60000UL;
4122  y ^= (y >> 18);
4123 
4124  return ((double)y + 0.5) * (1.0/4294967296.0);
4125  /* divided by 2^32 */
4126 }
4127 
4128 /* generates a random number on [0,1) with 53-bit resolution*/
4129 double genrand_res53()
4130 {
4131  unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
4132  return (a*67108864.0+b)*(1.0/9007199254740992.0);
4133 }
4134 
4135 /* These real versions are due to Isaku Wada, 2002/01/09 added */
4136 
4137 /************************************************************************
4138  * This marks the end of the block of code from Matsumoto and Nishimura *
4139  ************************************************************************/
4140 
4141 
4142 const int N_DAWSON = 100;
4143 
4144 // tabulated function values of Dawson's function
4145 // calculated with xmaxima 5.22.1
4146 static const double tbl_dawson[N_DAWSON+1] = {
4147  0.000000000000000000e+0,
4148  9.933599239785286114e-2,
4149  1.947510333680280496e-1,
4150  2.826316650213119286e-1,
4151  3.599434819348881042e-1,
4152  4.244363835020222959e-1,
4153  4.747632036629779306e-1,
4154  5.105040575592317787e-1,
4155  5.321017070563654290e-1,
4156  5.407243187262986750e-1,
4157  5.380794985696822201e-1,
4158  5.262066799705525356e-1,
4159  5.072734964077396141e-1,
4160  4.833975173848241360e-1,
4161  4.565072375268972572e-1,
4162  4.282490710853986254e-1,
4163  3.999398943230814126e-1,
4164  3.725593489740788250e-1,
4165  3.467727691148722451e-1,
4166  3.229743193228178382e-1,
4167  3.013403889237919660e-1,
4168  2.818849389255277938e-1,
4169  2.645107599508319576e-1,
4170  2.490529568377666955e-1,
4171  2.353130556638425762e-1,
4172  2.230837221674354811e-1,
4173  2.121651242424990041e-1,
4174  2.023745109105139857e-1,
4175  1.935507238593667923e-1,
4176  1.855552345354997718e-1,
4177  1.782710306105582873e-1,
4178  1.716003557161446928e-1,
4179  1.654619998786752031e-1,
4180  1.597885804744950500e-1,
4181  1.545240577369634525e-1,
4182  1.496215930807564847e-1,
4183  1.450417730540888593e-1,
4184  1.407511741154154125e-1,
4185  1.367212216746364963e-1,
4186  1.329272910810892667e-1,
4187  1.293480012360051155e-1,
4188  1.259646584343461329e-1,
4189  1.227608160065229225e-1,
4190  1.197219228088390280e-1,
4191  1.168350399532972540e-1,
4192  1.140886102268249801e-1,
4193  1.114722685321253072e-1,
4194  1.089766845891902214e-1,
4195  1.065934312832810744e-1,
4196  1.043148736220704706e-1,
4197  1.021340744242768354e-1,
4198  1.000447137201476355e-1,
4199  9.804101948507806734e-2,
4200  9.611770781195023073e-2,
4201  9.426993099823683440e-2,
4202  9.249323231075475996e-2,
4203  9.078350641561113352e-2,
4204  8.913696463869524546e-2,
4205  8.755010436413927499e-2,
4206  8.601968199264808016e-2,
4207  8.454268897454385223e-2,
4208  8.311633050835148859e-2,
4209  8.173800655824702378e-2,
4210  8.040529489538834118e-2,
4211  7.911593591113373223e-2,
4212  7.786781898606987138e-2,
4213  7.665897022891428948e-2,
4214  7.548754142476270211e-2,
4215  7.435180005364808165e-2,
4216  7.325012025863538754e-2,
4217  7.218097465823629202e-2,
4218  7.114292691123472774e-2,
4219  7.013462495342931107e-2,
4220  6.915479483562112754e-2,
4221  6.820223510065167384e-2,
4222  6.727581164463061598e-2,
4223  6.637445301385693290e-2,
4224  6.549714609447248675e-2,
4225  6.464293215671364666e-2,
4226  6.381090321984490301e-2,
4227  6.300019870755338791e-2,
4228  6.221000236682679049e-2,
4229  6.143953942619040819e-2,
4230  6.068807397169402555e-2,
4231  5.995490652126037542e-2,
4232  5.923937177997213955e-2,
4233  5.854083656061641403e-2,
4234  5.785869785535237086e-2,
4235  5.719238104574369009e-2,
4236  5.654133823962313062e-2,
4237  5.590504672435046070e-2,
4238  5.528300752700259690e-2,
4239  5.467474407290986634e-2,
4240  5.407980093473671650e-2,
4241  5.349774266500934015e-2,
4242  5.292815270562564644e-2,
4243  5.237063236845275277e-2,
4244  5.182479988163068060e-2,
4245  5.129028949666435102e-2,
4246  5.076675065180469937e-2,
4247  5.025384718759852810e-2
4248 };
4249 
4250 //
4251 // this routine calculates Dawson's function:
4252 //
4253 // x
4254 // -
4255 // | | 2 2
4256 // | (y - x )
4257 // | e dy
4258 // | |
4259 // -
4260 // 0
4261 //
4262 // the precomputed values are stored in the array tbl_dawson above.
4263 // tbl_dawson[i] is the value of the integral for x = double(i)/10.
4264 //
4265 // here we do 1st or 3rd order interpolation on this array using
4266 // the fact that the grid has equidistant spacing of 0.1.
4267 //
4268 // the accuracy of 3rd order interpolation is much better, but to
4269 // keep the routine fast we sometimes revert to 1st order when
4270 // the higher accuracy of 3rd order is not required.
4271 //
4272 // The Laurent series for this function is given in Mihalas Eq. 9-59.
4273 // It is not implemented here.
4274 //
4275 // dawson has been written by Peter van Hoof (ROB).
4276 //
4277 inline double dawson(double x, int order)
4278 {
4279  double x10 = x*10.;
4280  if( order == 1 )
4281  {
4282  int ind = min(max(int(x10),0),N_DAWSON-1);
4283  double p = x10 - double(ind);
4284  return tbl_dawson[ind] + p*(tbl_dawson[ind+1]-tbl_dawson[ind]);
4285  }
4286  else if( order == 3 )
4287  {
4288  int ind = min(max(int(x10-1.),0),N_DAWSON-3);
4289  double p = x10 - double(ind+1);
4290  double pm1 = p-1.;
4291  double pm2 = p-2.;
4292  double pp1 = p+1.;
4293  // Lagrange 4-point interpolation
4294  return p*pm1*(pp1*tbl_dawson[ind+3]-pm2*tbl_dawson[ind])/6. +
4295  pm2*pp1*(pm1*tbl_dawson[ind+1]-p*tbl_dawson[ind+2])/2.;
4296  }
4297  else
4298  {
4299  TotalInsanity();
4300  }
4301 }
4302 
4303 
4304 //
4305 // this is a fast version of the Voigt function that is only valid for small a.
4306 // The theory is described in:
4307 // >>refer rt Voigt Hjerting F., 1938, ApJ 88, 508
4308 //
4309 // FastVoigtH has been written by Peter van Hoof (ROB).
4310 //
4312 {
4313  DEBUG_ENTRY( "FastVoigtH()" );
4314 
4315  ASSERT( a <= 0.101f );
4316 
4317  //
4318  // This routine is guaranteed to give results to better than 0.25% relative accuracy
4319  // over its entire range of validity. The largest discrepancies occur between 1 < v < 10,
4320  // peaking around v = 2 - 7, as shown in the table below:
4321  //
4322  // a = 1e-10 : v = 5.6234e+00 dH/H = 7.6e-06
4323  // a = 1e-5 : v = 7.0795e+00 dH/H = 7.5e-06
4324  // a = 1e-4 : v = 3.7584e+00 dH/H = 1.3e-05
4325  // a = 3e-4 : v = 3.5481e+00 dH/H = 1.6e-05
4326  // a = 1e-3 : v = 3.3497e+00 dH/H = 1.9e-05
4327  // a = 3e-3 : v = 3.1623e+00 dH/H = 2.2e-05
4328  // a = 0.01 : v = 3.1623e+00 dH/H = 1.0e-05
4329  // a = 0.03 : v = 2.8184e+00 dH/H = 1.8e-04
4330  // a = 0.1 : v = 2.6607e+00 dH/H = 2.4e-03
4331  //
4332  // to get a guaranteed relative accuracy <= 1.e-4, use a < 0.0235
4333  // to get a guaranteed relative accuracy <= 1.e-3, use a < 0.066
4334  //
4335  // for a > 0.1 the series expansions used in this routine quickly start breaking down
4336  // and the algorithm becomes useless, so never use this routine for a > 0.1 !
4337  //
4338 
4339  v = abs(v); // the profile is symmetrical in v
4340 
4341  if( v > 9.f )
4342  {
4343  // use Laurent series; this is Eq. 7 of Hjerting with one higher order term added
4344  //
4345  // The complete series is:
4346  //
4347  // inf
4348  // ----
4349  // a \ (2 n + 1)!!
4350  // H(a,v) = ----------- | -----------
4351  // 2 / n 2n
4352  // sqrt(pi) v ---- 2 v
4353  // n=0
4354  //
4355  realnum vm2 = 1.f/pow2(v);
4356  return a*vm2/realnum(SQRTPI)*(((105.f/8.f*vm2 + 15.f/4.f)*vm2 + 3.f/2.f)*vm2 + 1.f);
4357  }
4358  else
4359  {
4360  realnum v2 = pow2(v);
4361  // NOTE: it is important to use dsexp here so that we avoid expf(). The
4362  // latter can be significantly slower, at least on non-Suse Linux platforms!
4363  // see: https://bugzilla.redhat.com/show_bug.cgi?id=521190
4364  // revert this statement to: emv2 = exp(-v2) once this is solved.
4365  realnum emv2 = realnum(dsexp(v2));
4366  int order = ( a > 0.003f || v > 1.5f ) ? 3 : 1;
4367  // this is Eq. 3 of Hjerting with an additional term:
4368  //
4369  // the last term in Eq. 4 of Hjerting can be expanded to lowest order in a as:
4370  //
4371  // inf
4372  // -
4373  // | | 2 2 2
4374  // 1 | (a x) - x 2 2 - v
4375  // -------- | ------ exp(----) cos(v x) dx = - a (2 v - 1) e
4376  // sqrt(pi) | | 2 4
4377  // -
4378  // 0
4379  //
4380  return emv2*(1.f-pow2(a)*(2.f*v2-1.f)) +
4381  2.f*a/realnum(SQRTPI)*(-1.f + 2.f*v*realnum(dawson(v,order)));
4382  }
4383 }
4384 
4385 // vector version of the above
4386 void FastVoigtH(realnum a, const realnum v[], realnum y[], size_t n)
4387 {
4388  DEBUG_ENTRY( "FastVoigtH()" );
4389 
4390  ASSERT( a <= 0.101f );
4391 
4392  // check that v[i] is strictly monotonically increasing!
4393  // we don't want to do this every time to save CPU cycles
4394  //for( size_t i=1; i < n; ++i )
4395  // ASSERT( v[i] > v[i-1] );
4396 
4397  // find the core of the profile
4398  size_t ilo = 0, ihi = n;
4399  for( size_t i=0; i < n; ++i )
4400  {
4401  if( v[i] < -9.f )
4402  ++ilo;
4403  else if( v[i] > 9.f )
4404  {
4405  ihi = i;
4406  break;
4407  }
4408  }
4409 
4410  for( size_t i=0; i < ilo; ++i )
4411  {
4412  realnum vm2 = 1.f/pow2(v[i]);
4413  y[i] = a*vm2/realnum(SQRTPI)*(((105.f/8.f*vm2 + 15.f/4.f)*vm2 + 3.f/2.f)*vm2 + 1.f);
4414  }
4415 
4416  avx_ptr<realnum> arg(ilo,ihi), expval(ilo,ihi);
4417  for( size_t i=ilo; i < ihi; ++i )
4418  arg[i] = -pow2(v[i]);
4419  vexp( arg.ptr0(), expval.ptr0(), ilo, ihi );
4420  for( size_t i=ilo; i < ihi; ++i )
4421  {
4422  realnum vv = abs(v[i]);
4423  int order = ( a > 0.003f || vv > 1.5f ) ? 3 : 1;
4424  y[i] = expval[i]*(1.f-pow2(a)*(2.f*pow2(vv)-1.f)) +
4425  2.f*a/realnum(SQRTPI)*(-1.f + 2.f*vv*realnum(dawson(vv,order)));
4426  }
4427 
4428  for( size_t i=ihi; i < n; ++i )
4429  {
4430  realnum vm2 = 1.f/pow2(v[i]);
4431  y[i] = a*vm2/realnum(SQRTPI)*(((105.f/8.f*vm2 + 15.f/4.f)*vm2 + 3.f/2.f)*vm2 + 1.f);
4432  }
4433 }
4434 
4435 /*
4436  Calculate the Voigt profile aka Faddeeva function with relative
4437  error less than 10^(-R).
4438 
4439  R0=1.51*EXP(1.144*R) and R1=1.60*EXP(0.554*R) can be set by the the
4440  user subject to the constraints 14.88<R0<460.4 and 4.85<R1<25.5
4441 
4442  Translated from a Fortran version by R.J. Wells,
4443  see http://www.atm.ox.ac.uk/user/wells/voigt.html
4444 
4445  >>refer rt Voigt Wells, R.J. 1999, JQSRT, 62, 29
4446 */
4447 void humlik(int n, const realnum x[], realnum y, realnum k[])
4448 {
4449  DEBUG_ENTRY( "humlik()" );
4450 
4451  /* n IN Number of points
4452  x IN Input x array
4453  y IN Input y value >=0.0
4454  k OUT (Voigt) array */
4455 
4456  // use doubles internally to avoid overflow for very large y values (above roughly 5000)
4457  // these very high values can occur in X-ray lines; the end result is converted back to realnum
4458 
4459  /* Initialized data */
4460  static const double c[6] = { 1.0117281, -.75197147, .012557727, .010022008, -2.4206814e-4, 5.0084806e-7 };
4461  static const double s[6] = { 1.393237, .23115241, -.15535147, .0062183662, 9.1908299e-5, -6.2752596e-7 };
4462  static const double t[6] = { .31424038, .94778839, 1.5976826, 2.2795071, 3.020637, 3.8897249 };
4463 
4464  static const double R0 = 146.7, R1 = 14.67;
4465 
4466  /* Local variables */
4467  double d, ki;
4468  int i, j;
4469  double a0, d0, e0, d2, e2, h0, e4, h2, h4, h6, p0,
4470  p2, p4, p6, p8, z0, z2, z4, z6, z8, mf[6], pf[6],
4471  mq[6], yf, pq[6], xm[6], ym[6], xp[6], xq, yq, yp[6];
4472  bool rg1, rg2, rg3;
4473  double abx, ypy0, xlim0, xlim1, xlim2, xlim3, xlim4, ypy0q, yrrtpi;
4474 
4475  rg1 = rg2 = rg3 = true;
4476  // Initialization to quiet warnings
4477  z0 = z2 = z4 = z6 = z8 = 0.;
4478  p0 = p2 = p4 = p6 = p8 = 0.;
4479  h0 = h2 = h4 = h6 = 0.;
4480  a0 = d0 = d2 = e0 = e2 = e4 = 0.;
4481 
4482  yq = y * y;
4483  yrrtpi = y * .56418958; // 1/SQRT(pi)
4484  /* Region boundaries */
4485  xlim0 = R0 - y;
4486  xlim1 = R1 - y;
4487  xlim3 = y * 3.097 - .45;
4488  xlim2 = 6.8 - y;
4489  xlim4 = y * 18.1 + 1.65;
4490  if (y <= 1e-6)
4491  { /* avoid W4 algorithm */
4492  xlim1 = xlim0;
4493  xlim2 = xlim0;
4494  }
4495 
4496  for (i = 0; i < n; ++i)
4497  {
4498  abx = fabs(x[i]);
4499  xq = abx * abx;
4500  if (abx > xlim0)
4501  { /* Region 0 algorithm */
4502  k[i] = realnum(yrrtpi / (xq + yq));
4503  }
4504  else if (abx > xlim1)
4505  { /* Humlicek W4 Region 1 */
4506  if (rg1)
4507  {
4508  /* First point in Region 1 */
4509  rg1 = false;
4510  a0 = yq + .5;
4511  d0 = a0 * a0;
4512  d2 = yq + yq - 1.;
4513  }
4514  d = .56418958 / (d0 + xq * (d2 + xq));
4515  k[i] = realnum(d * y * (a0 + xq));
4516  }
4517  else if (abx > xlim2)
4518  { /* Humlicek W4 Region 2 */
4519  if (rg2)
4520  { /* First point in Region 2 */
4521  rg2 = false;
4522  h0 = yq * (yq * (yq * (yq + 6.) + 10.5) + 4.5) + .5625;
4523  h2 = yq * (yq * (yq * 4. + 6.) + 9.) - 4.5;
4524  h4 = 10.5 - yq * (6. - yq * 6.);
4525  h6 = yq * 4. - 6.;
4526  e0 = yq * (yq * (yq + 5.5) + 8.25) + 1.875;
4527  e2 = yq * (yq * 3. + 1.) + 5.25;
4528  e4 = h6 * .75;
4529  }
4530  d = .56418958 / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
4531  k[i] = realnum(d * y * (e0 + xq * (e2 + xq * (e4 + xq))));
4532  }
4533  else if (abx < xlim3)
4534  { /* Humlicek W4 Region 3 */
4535  if (rg3)
4536  {
4537  /* First point in Region 3 */
4538  rg3 = false;
4539  z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 13.3988) +
4540  88.26741) + 369.1989) + 1074.409) + 2256.981) + 3447.629) +
4541  3764.966) + 2802.87) + 1280.829) + 272.1014;
4542  z2 = y * (y * (y * (y * (y * (y * (y * (y * 5. + 53.59518) +
4543  266.2987) + 793.4273) + 1549.675) + 2037.31) + 1758.336) +
4544  902.3066) + 211.678;
4545  z4 = y * (y * (y * (y * (y * (y * 10. + 80.39278) + 269.2916) +
4546  479.2576) + 497.3014) + 308.1852) + 78.86585;
4547  z6 = y * (y * (y * (y * 10. + 53.59518) + 92.75679) + 55.02933) +
4548  22.03523;
4549  z8 = y * (y * 5. + 13.3988) + 1.49646;
4550  p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * .3183291 +
4551  4.264678) + 27.93941) + 115.3772) + 328.2151) + 662.8097) +
4552  946.897) + 919.4955) + 549.3954) + 153.5168;
4553  p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163 + 12.79458) +
4554  56.81652) + 139.4665) + 189.773) + 124.5975) - 1.322256) -
4555  34.16955;
4556  p4 = y * (y * (y * (y * (y * 1.9099744 + 12.79568) + 29.81482) +
4557  24.01655) + 10.46332) + 2.584042;
4558  p6 = y * (y * (y * 1.273316 + 4.266322) + .9377051) - .07272979;
4559  p8 = y * .3183291 + 5.480304e-4;
4560  }
4561  d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
4562  k[i] = realnum(d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8)))));
4563  }
4564  else
4565  { /* Humlicek CPF12 algorithm */
4566  ypy0 = y + 1.5;
4567  ypy0q = ypy0 * ypy0;
4568  ki = 0.;
4569  for (j = 0; j <= 5; ++j)
4570  {
4571  d = x[i] - t[j];
4572  mq[j] = d * d;
4573  mf[j] = 1. / (mq[j] + ypy0q);
4574  xm[j] = mf[j] * d;
4575  ym[j] = mf[j] * ypy0;
4576  d = x[i] + t[j];
4577  pq[j] = d * d;
4578  pf[j] = 1. / (pq[j] + ypy0q);
4579  xp[j] = pf[j] * d;
4580  yp[j] = pf[j] * ypy0;
4581  }
4582  if (abx <= xlim4)
4583  { /* Humlicek CPF12 Region I */
4584  for (j = 0; j <= 5; ++j)
4585  {
4586  ki = ki + c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
4587  }
4588  }
4589  else
4590  { /* Humlicek CPF12 Region II */
4591  yf = y + 3.;
4592  for (j = 0; j <= 5; ++j)
4593  {
4594  ki = ki + (c[j] * (mq[j] * mf[j] - ym[j] * 1.5) + s[j] * yf * xm[j]) / (mq[j] + 2.25) +
4595  (c[j] * (pq[j] * pf[j] - yp[j] * 1.5) - s[j] * yf * xp[j]) / (pq[j] + 2.25);
4596  }
4597  ki = y * ki + exp(-xq);
4598  }
4599  k[i] = realnum(ki);
4600  }
4601  }
4602 }
4603 
4604 /******************************************************************
4605  * This marks the end of the block of code for the Voigt function *
4606  ******************************************************************/
4607 
4608 STATIC uint32 MD5swap( uint32 word );
4609 STATIC void MD5_Transform (uint32 *digest, const uint32 *in);
4610 
4611 //
4612 // The routines MD5file(), MD5datafile(), MD5datastream(), MD5string() and MD5swap() were written by Peter van Hoof
4613 //
4614 // this version returns the md5sum of a file and is identical to the well known md5sum algorithm
4615 string MD5file(const char* fnam, access_scheme scheme)
4616 {
4617  DEBUG_ENTRY( "MD5file()" );
4618 
4619  fstream ioFile;
4620  open_data( ioFile, fnam, mode_r, scheme );
4621 
4622  char c;
4623  string content;
4624  while( ioFile.get(c) )
4625  content += c;
4626 
4627  return MD5string( content );
4628 }
4629 
4630 
4631 // this version returns the md5sum of a text file. It filters out the eol characters,
4632 // which makes it incompatible with the md5sum algorithm, but it is OS agnostic...
4633 // comment lines that start with the hash symbol are also skipped
4634 string MD5datafile(const char* fnam, access_scheme scheme)
4635 {
4636  DEBUG_ENTRY( "MD5datafile()" );
4637 
4638  fstream ioFile;
4639  open_data( ioFile, fnam, mode_r, scheme );
4640  return MD5datastream(ioFile);
4641 }
4642 
4643 
4644 string MD5datastream(fstream& ioFile)
4645 {
4646  DEBUG_ENTRY( "MD5datastream()" );
4647 
4648  // get file size
4649  ioFile.seekg( 0, ios::end );
4650  streampos fsize = ioFile.tellg();
4651  ioFile.seekg( 0, ios::beg );
4652 
4653  // get strings with suitable sizes
4654  string line, content;
4655  line.reserve(INPUT_LINE_LENGTH);
4656  content.reserve(fsize);
4657 
4658  while( getline( ioFile, line ) )
4659  if( line[0] != '#' )
4660  content += line;
4661 
4662  return MD5string( content );
4663 }
4664 
4665 
4666 // calculate the md5sum of an arbitrary string
4667 string MD5string(const string& str)
4668 {
4669  DEBUG_ENTRY( "MD5string()" );
4670 
4671  uint32 state[4];
4672 
4673  state[0] = 0x67452301L;
4674  state[1] = 0xefcdab89L;
4675  state[2] = 0x98badcfeL;
4676  state[3] = 0x10325476L;
4677 
4678  string lstr = str;
4679 
4680  // pad the string following RFC 1321 3.1 Step 1.
4681  uint32 bytes = str.length()%64;
4682  uint32 padlen = ( bytes >= 56 ) ? 64 + 56 - bytes : 56 - bytes;
4683  lstr += '\x80';
4684  for( uint32 i=1; i < padlen; ++i )
4685  lstr += '\0';
4686 
4687  ASSERT( lstr.length()%64 == 56 );
4688 
4689  uint32 i;
4690  union {
4691  uint32 in[16];
4692  unsigned char chr[64];
4693  } u;
4694 
4695  for( i=0; i < lstr.length()/64; ++i )
4696  {
4697  for( uint32 j=0; j < 64; ++j )
4698  {
4699  if( cpu.i().little_endian() )
4700  u.chr[j] = lstr[i*64+j];
4701  else
4702  {
4703  uint32 jr = j%4;
4704  uint32 j4 = j-jr;
4705  u.chr[j4+3-jr] = lstr[i*64+j];
4706  }
4707  }
4708  MD5_Transform( state, u.in );
4709  }
4710  for( uint32 j=0; j < 56; ++j )
4711  {
4712  if( cpu.i().little_endian() )
4713  u.chr[j] = lstr[i*64+j];
4714  else
4715  {
4716  uint32 jr = j%4;
4717  uint32 j4 = j-jr;
4718  u.chr[j4+3-jr] = lstr[i*64+j];
4719  }
4720  }
4721  // append the length of the string in _bits_ following RFC 1321 3.1 Step 2.
4722  u.in[14] = (str.length()<<3)&0xffffffff;
4723  u.in[15] = (str.length()>>29)&0xffffffff;
4724 
4725  MD5_Transform( state, u.in );
4726 
4727  ostringstream hash;
4728  for( uint32 i=0; i < 4; ++i )
4729  hash << hex << setfill('0') << setw(8) << MD5swap(state[i]);
4730 
4731  return hash.str();
4732 }
4733 
4734 STATIC uint32 MD5swap( uint32 word )
4735 {
4736  DEBUG_ENTRY( "MD5swap()" );
4737 
4738  union {
4739  uint32 word;
4740  unsigned char byte[4];
4741  } ui, uo;
4742 
4743  ui.word = word;
4744  uo.byte[0] = ui.byte[3];
4745  uo.byte[1] = ui.byte[2];
4746  uo.byte[2] = ui.byte[1];
4747  uo.byte[3] = ui.byte[0];
4748 
4749  return uo.word;
4750 }
4751 
4752 //
4753 // The following implementation of the MD5 algorithm was taken from the
4754 // Crypto++ library (http://www.cryptopp.com/) and is subject to the
4755 // following license:
4756 //
4757 //
4758 // Compilation Copyright (c) 1995-2010 by Wei Dai. All rights reserved.
4759 // This copyright applies only to this software distribution package
4760 // as a compilation, and does not imply a copyright on any particular
4761 // file in the package.
4762 //
4763 // All individual files in this compilation are placed in the public domain by
4764 // Wei Dai and other contributors.
4765 //
4766 // I would like to thank the following authors for placing their works into
4767 // the public domain:
4768 //
4769 // Joan Daemen - 3way.cpp
4770 // Leonard Janke - cast.cpp, seal.cpp
4771 // Steve Reid - cast.cpp
4772 // Phil Karn - des.cpp
4773 // Andrew M. Kuchling - md2.cpp, md4.cpp
4774 // Colin Plumb - md5.cpp
4775 // Seal Woods - rc6.cpp
4776 // Chris Morgan - rijndael.cpp
4777 // Paulo Baretto - rijndael.cpp, skipjack.cpp, square.cpp
4778 // Richard De Moliner - safer.cpp
4779 // Matthew Skala - twofish.cpp
4780 // Kevin Springle - camellia.cpp, shacal2.cpp, ttmac.cpp, whrlpool.cpp, ripemd.cpp
4781 //
4782 // Permission to use, copy, modify, and distribute this compilation for
4783 // any purpose, including commercial applications, is hereby granted
4784 // without fee, subject to the following restrictions:
4785 //
4786 // 1. Any copy or modification of this compilation in any form, except
4787 // in object code form as part of an application software, must include
4788 // the above copyright notice and this license.
4789 //
4790 // 2. Users of this software agree that any modification or extension
4791 // they provide to Wei Dai will be considered public domain and not
4792 // copyrighted unless it includes an explicit copyright notice.
4793 //
4794 // 3. Wei Dai makes no warranty or representation that the operation of the
4795 // software in this compilation will be error-free, and Wei Dai is under no
4796 // obligation to provide any services, by way of maintenance, update, or
4797 // otherwise. THE SOFTWARE AND ANY DOCUMENTATION ARE PROVIDED "AS IS"
4798 // WITHOUT EXPRESS OR IMPLIED WARRANTY INCLUDING, BUT NOT LIMITED TO,
4799 // THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
4800 // PURPOSE. IN NO EVENT WILL WEI DAI OR ANY OTHER CONTRIBUTOR BE LIABLE FOR
4801 // DIRECT, INCIDENTAL OR CONSEQUENTIAL DAMAGES, EVEN IF
4802 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
4803 //
4804 // 4. Users will not use Wei Dai or any other contributor's name in any
4805 // publicity or advertising, without prior written consent in each case.
4806 //
4807 // 5. Export of this software from the United States may require a
4808 // specific license from the United States Government. It is the
4809 // responsibility of any person or organization contemplating export
4810 // to obtain such a license before exporting.
4811 //
4812 // 6. Certain parts of this software may be protected by patents. It
4813 // is the users' responsibility to obtain the appropriate
4814 // licenses before using those parts.
4815 //
4816 // If this compilation is used in object code form in an application
4817 // software, acknowledgement of the author is not required but would be
4818 // appreciated. The contribution of any useful modifications or extensions
4819 // to Wei Dai is not required but would also be appreciated.
4820 //
4821 
4822 // md5.cpp - modified by Wei Dai from Colin Plumb's public domain md5.c
4823 // any modifications are placed in the public domain
4824 
4825 inline uint32 rotlFixed(uint32 x, unsigned int y)
4826 {
4827  return uint32((x<<y) | (x>>(32-y)));
4828 }
4829 
4830 STATIC void MD5_Transform (uint32 *digest, const uint32 *in)
4831 {
4832  DEBUG_ENTRY( "MD5_Transform()" );
4833 
4834 #define F1(x, y, z) (z ^ (x & (y ^ z)))
4835 #define F2(x, y, z) F1(z, x, y)
4836 #define F3(x, y, z) (x ^ y ^ z)
4837 #define F4(x, y, z) (y ^ (x | ~z))
4838 
4839 #define MD5STEP(f, w, x, y, z, data, s) \
4840  w = rotlFixed(w + f(x, y, z) + data, s) + x
4841 
4842  uint32 a, b, c, d;
4843 
4844  a=digest[0];
4845  b=digest[1];
4846  c=digest[2];
4847  d=digest[3];
4848 
4849  MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
4850  MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
4851  MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
4852  MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
4853  MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
4854  MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
4855  MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
4856  MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
4857  MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
4858  MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
4859  MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
4860  MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
4861  MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
4862  MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
4863  MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
4864  MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
4865 
4866  MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
4867  MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
4868  MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
4869  MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
4870  MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
4871  MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
4872  MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
4873  MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
4874  MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
4875  MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
4876  MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
4877  MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
4878  MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
4879  MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
4880  MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
4881  MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
4882 
4883  MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
4884  MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
4885  MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
4886  MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
4887  MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
4888  MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
4889  MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
4890  MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
4891  MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
4892  MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
4893  MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
4894  MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
4895  MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
4896  MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
4897  MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
4898  MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
4899 
4900  MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
4901  MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
4902  MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
4903  MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
4904  MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
4905  MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
4906  MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
4907  MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
4908  MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
4909  MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
4910  MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
4911  MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
4912  MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
4913  MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
4914  MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
4915  MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
4916 
4917  digest[0] += a;
4918  digest[1] += b;
4919  digest[2] += c;
4920  digest[3] += d;
4921 }
4922 
4923 /***************************************************************
4924  * This marks the end of the block of code from Wei Dai et al. *
4925  ***************************************************************/
4926 
4927 /* From http://www.netlib.org/fmm/svd.f */
4928 /* Downloaded 11/Feb/2014 */
4929 void svd(const int nm, const int m, const int n, double *a,
4930  double *w, bool matu, double *u, bool matv,
4931  double *v, int *ierr, double *rv1)
4932 {
4933  double c, f, g, h;
4934  int i, j, k, l=0;
4935  double s, x, y, z;
4936  int i1, k1, l1;
4937  double scale, anorm;
4938 
4939 
4940 
4941 /* this subroutine is a translation of the algol procedure svd, */
4942 /* num. math. 14, 403-420(1970) by golub and reinsch. */
4943 /* handbook for auto. comp., vol ii-linear algebra, 134-151(1971). */
4944 
4945 /* this subroutine determines the singular value decomposition */
4946 /* t */
4947 /* a=usv of a real m by n rectangular matrix. householder */
4948 /* bidiagonalization and a variant of the qr algorithm are used. */
4949 
4950 /* on input. */
4951 
4952 /* nm must be set to the row dimension of two-dimensional */
4953 /* array parameters as declared in the calling program */
4954 /* dimension statement. note that nm must be at least */
4955 /* as large as the maximum of m and n. */
4956 
4957 /* m is the number of rows of a (and u). */
4958 
4959 /* n is the number of columns of a (and u) and the order of v. */
4960 
4961 /* a contains the rectangular input matrix to be decomposed. */
4962 
4963 /* matu should be set to .true. if the u matrix in the */
4964 /* decomposition is desired, and to .false. otherwise. */
4965 
4966 /* matv should be set to .true. if the v matrix in the */
4967 /* decomposition is desired, and to .false. otherwise. */
4968 
4969 /* on output. */
4970 
4971 /* a is unaltered (unless overwritten by u or v). */
4972 
4973 /* w contains the n (non-negative) singular values of a (the */
4974 /* diagonal elements of s). they are unordered. if an */
4975 /* error exit is made, the singular values should be correct */
4976 /* for indices ierr+1,ierr+2,...,n. */
4977 
4978 /* u contains the matrix u (orthogonal column vectors) of the */
4979 /* decomposition if matu has been set to .true. otherwise */
4980 /* u is used as a temporary array. u may coincide with a. */
4981 /* if an error exit is made, the columns of u corresponding */
4982 /* to indices of correct singular values should be correct. */
4983 
4984 /* v contains the matrix v (orthogonal) of the decomposition if */
4985 /* matv has been set to .true. otherwise v is not referenced. */
4986 /* v may also coincide with a if u is not needed. if an error */
4987 /* exit is made, the columns of v corresponding to indices of */
4988 /* correct singular values should be correct. */
4989 
4990 /* ierr is set to */
4991 /* zero for normal return, */
4992 /* k if the k-th singular value has not been */
4993 /* determined after 30 iterations. */
4994 
4995 /* rv1 is a temporary storage array. */
4996 
4997 /* this is a modified version of a routine from the eispack */
4998 /* collection by the nats project */
4999 
5000 /* modified to eliminate machep */
5001 
5002  /* Parameter adjustments */
5003  --rv1;
5004  const int v_dim1 = nm;
5005  const int v_offset = 1 + v_dim1;
5006  v -= v_offset;
5007  const int u_dim1 = nm;
5008  const int u_offset = 1 + u_dim1;
5009  u -= u_offset;
5010  --w;
5011  const int a_dim1 = nm;
5012  const int a_offset = 1 + a_dim1;
5013  a -= a_offset;
5014 
5015  *ierr = 0;
5016 
5017  for (i = 1; i <= m; ++i)
5018  {
5019  for (j = 1; j <= n; ++j)
5020  {
5021  u[i + j * u_dim1] = a[i + j * a_dim1];
5022  }
5023  }
5024  /* householder reduction to bidiagonal form */
5025  g = 0.;
5026  scale = 0.;
5027  anorm = 0.;
5028 
5029  for (i = 1; i <= n; ++i)
5030  {
5031  l = i + 1;
5032  rv1[i] = scale * g;
5033  g = 0.;
5034  s = 0.;
5035  scale = 0.;
5036  if (i <= m) {
5037  for (k = i; k <= m; ++k)
5038  {
5039  scale += fabs(u[k + i * u_dim1]);
5040  }
5041 
5042  if (scale != 0.)
5043  {
5044  for (k = i; k <= m; ++k)
5045  {
5046  u[k + i * u_dim1] /= scale;
5047  double d__1 = u[k + i * u_dim1];
5048  s += d__1 * d__1;
5049  }
5050 
5051  f = u[i + i * u_dim1];
5052  g = -copysign(sqrt(s), f);
5053  h = f * g - s;
5054  u[i + i * u_dim1] = f - g;
5055  if (i != n)
5056  {
5057  for (j = l; j <= n; ++j)
5058  {
5059  s = 0.;
5060  for (k = i; k <= m; ++k)
5061  {
5062  s += u[k + i * u_dim1] * u[k + j * u_dim1];
5063  }
5064  f = s / h;
5065 
5066  for (k = i; k <= m; ++k)
5067  {
5068  u[k + j * u_dim1] += f * u[k + i * u_dim1];
5069  }
5070  }
5071  }
5072  for (k = i; k <= m; ++k)
5073  {
5074  u[k + i * u_dim1] = scale * u[k + i * u_dim1];
5075  }
5076  }
5077  }
5078  w[i] = scale * g;
5079  g = 0.;
5080  s = 0.;
5081  scale = 0.;
5082  if (i <= m && i != n)
5083  {
5084  for (k = l; k <= n; ++k)
5085  {
5086  scale += fabs(u[i + k * u_dim1]);
5087  }
5088 
5089  if (scale != 0.)
5090  {
5091  for (k = l; k <= n; ++k)
5092  {
5093  u[i + k * u_dim1] /= scale;
5094  double d__1 = u[i + k * u_dim1];
5095  s += d__1 * d__1;
5096  }
5097 
5098  f = u[i + l * u_dim1];
5099  g = -copysign(sqrt(s), f);
5100  h = f * g - s;
5101  u[i + l * u_dim1] = f - g;
5102 
5103  for (k = l; k <= n; ++k)
5104  {
5105  rv1[k] = u[i + k * u_dim1] / h;
5106  }
5107 
5108  if (i != m)
5109  {
5110  for (j = l; j <= m; ++j)
5111  {
5112  s = 0.;
5113  for (k = l; k <= n; ++k)
5114  {
5115  s += u[j + k * u_dim1] * u[i + k * u_dim1];
5116  }
5117  for (k = l; k <= n; ++k)
5118  {
5119  u[j + k * u_dim1] += s * rv1[k];
5120  }
5121  }
5122  }
5123  for (k = l; k <= n; ++k)
5124  {
5125  u[i + k * u_dim1] = scale * u[i + k * u_dim1];
5126  }
5127  }
5128  }
5129  anorm = fmax(anorm,fabs(w[i]) + fabs(rv1[i]));
5130  }
5131 
5132  /* accumulation of right-hand transformations */
5133  if (matv)
5134  {
5135  /* for i=n step -1 until 1 do -- */
5136  for (int ii = 1; ii <= n; ++ii)
5137  {
5138  i = n + 1 - ii;
5139  if (i != n)
5140  {
5141  if (g != 0.)
5142  {
5143 
5144  for (j = l; j <= n; ++j)
5145  {
5146  /* double division avoids possible underflow */
5147  v[j + i * v_dim1] = u[i + j * u_dim1] / u[i + l * u_dim1] / g;
5148  }
5149 
5150  for (j = l; j <= n; ++j)
5151  {
5152  s = 0.;
5153  for (k = l; k <= n; ++k)
5154  {
5155  s += u[i + k * u_dim1] * v[k + j * v_dim1];
5156  }
5157 
5158  for (k = l; k <= n; ++k)
5159  {
5160  v[k + j * v_dim1] += s * v[k + i * v_dim1];
5161  }
5162  }
5163  }
5164  for (j = l; j <= n; ++j) {
5165  v[i + j * v_dim1] = 0.;
5166  v[j + i * v_dim1] = 0.;
5167  }
5168  }
5169  v[i + i * v_dim1] = 1.;
5170  g = rv1[i];
5171  l = i;
5172  }
5173  }
5174 
5175  /* accumulation of left-hand transformations */
5176  if (matu)
5177  {
5178  /* for i=min(m,n) step -1 until 1 do -- */
5179  int mn = n;
5180  if (m < n)
5181  {
5182  mn = m;
5183  }
5184 
5185  for (int ii = 1; ii <= mn; ++ii)
5186  {
5187  i = mn + 1 - ii;
5188  l = i + 1;
5189  g = w[i];
5190  if (i != n)
5191  {
5192  for (j = l; j <= n; ++j)
5193  {
5194  u[i + j * u_dim1] = 0.;
5195  }
5196  }
5197 
5198  if (g != 0.)
5199  {
5200  if (i != mn)
5201  {
5202  for (j = l; j <= n; ++j)
5203  {
5204  s = 0.;
5205  for (k = l; k <= m; ++k)
5206  {
5207  s += u[k + i * u_dim1] * u[k + j * u_dim1];
5208  }
5209  /* double division avoids possible underflow */
5210  f = s / u[i + i * u_dim1] / g;
5211 
5212  for (k = i; k <= m; ++k)
5213  {
5214  u[k + j * u_dim1] += f * u[k + i * u_dim1];
5215  }
5216  }
5217  }
5218 
5219  for (j = i; j <= m; ++j)
5220  {
5221  u[j + i * u_dim1] /= g;
5222  }
5223  }
5224  else
5225  {
5226  for (j = i; j <= m; ++j)
5227  {
5228  u[j + i * u_dim1] = 0.;
5229  }
5230  }
5231  u[i + i * u_dim1] += 1.;
5232  }
5233  }
5234  /* diagonalization of the bidiagonal form */
5235  /* for k=n step -1 until 1 do -- */
5236  for (int kk = 1; kk <= n; ++kk)
5237  {
5238  k1 = n - kk;
5239  k = k1 + 1;
5240  int its = 0;
5241  while (1)
5242  {
5243  /* test for splitting. */
5244  /* for l=k step -1 until 1 do -- */
5245  for (int ll = 1; ll <= k; ++ll)
5246  {
5247  l1 = k - ll;
5248  l = l1 + 1;
5249  /* rv1(1) is always zero, so there is no exit */
5250  /* through the bottom of the loop */
5251  if (fabs(rv1[l]) + anorm == anorm)
5252  {
5253  break;
5254  }
5255  if (fabs(w[l1]) + anorm == anorm)
5256  {
5257  /* cancellation of rv1(l) if l greater than 1 */
5258  c = 0.;
5259  s = 1.;
5260 
5261  for (i = l; i <= k; ++i)
5262  {
5263  f = s * rv1[i];
5264  rv1[i] = c * rv1[i];
5265  if (fabs(f) + anorm == anorm)
5266  {
5267  break;
5268  }
5269  g = w[i];
5270  h = sqrt(f * f + g * g);
5271  w[i] = h;
5272  c = g / h;
5273  s = -f / h;
5274  if (matu)
5275  {
5276  for (j = 1; j <= m; ++j)
5277  {
5278  y = u[j + l1 * u_dim1];
5279  z = u[j + i * u_dim1];
5280  u[j + l1 * u_dim1] = y * c + z * s;
5281  u[j + i * u_dim1] = -y * s + z * c;
5282  }
5283  }
5284  }
5285  break;
5286  }
5287  }
5288 
5289  /* test for convergence */
5290  z = w[k];
5291  if (l == k)
5292  {
5293  break;
5294  }
5295  if (its == 30)
5296  {
5297  /* set error -- no convergence to a singular value after
5298  30 iterations
5299  */
5300  *ierr = k;
5301  return;
5302  }
5303  ++its;
5304  /* shift from bottom 2 by 2 minor */
5305  x = w[l];
5306  y = w[k1];
5307  g = rv1[k1];
5308  h = rv1[k];
5309  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (h * 2. * y);
5310  g = sqrt(f * f + 1.);
5311  f = ((x - z) * (x + z) + h * (y / (f + copysign(g, f)) - h)) / x;
5312  /* next qr transformation */
5313  c = 1.;
5314  s = 1.;
5315 
5316  for (i1 = l; i1 <= k1; ++i1)
5317  {
5318  i = i1 + 1;
5319  g = rv1[i];
5320  y = w[i];
5321  h = s * g;
5322  g = c * g;
5323  z = sqrt(f * f + h * h);
5324  rv1[i1] = z;
5325  c = f / z;
5326  s = h / z;
5327  f = x * c + g * s;
5328  g = -x * s + g * c;
5329  h = y * s;
5330  y *= c;
5331  if (matv)
5332  {
5333  for (j = 1; j <= n; ++j)
5334  {
5335  x = v[j + i1 * v_dim1];
5336  z = v[j + i * v_dim1];
5337  v[j + i1 * v_dim1] = x * c + z * s;
5338  v[j + i * v_dim1] = -x * s + z * c;
5339  }
5340  }
5341 
5342  z = sqrt(f * f + h * h);
5343  w[i1] = z;
5344  /* rotation can be arbitrary if z is zero */
5345  if (z != 0.)
5346  {
5347  c = f / z;
5348  s = h / z;
5349  }
5350  f = c * g + s * y;
5351  x = -s * g + c * y;
5352  if (matu)
5353  {
5354  for (j = 1; j <= m; ++j)
5355  {
5356  y = u[j + i1 * u_dim1];
5357  z = u[j + i * u_dim1];
5358  u[j + i1 * u_dim1] = y * c + z * s;
5359  u[j + i * u_dim1] = -y * s + z * c;
5360  }
5361  }
5362  }
5363 
5364  rv1[l] = 0.;
5365  rv1[k] = f;
5366  w[k] = x;
5367  }
5368  /* convergence */
5369  if (z < 0.)
5370  {
5371  /* w(k) is made non-negative */
5372  w[k] = -z;
5373  if (matv)
5374  {
5375  for (j = 1; j <= n; ++j)
5376  {
5377  v[j + k * v_dim1] = -v[j + k * v_dim1];
5378  }
5379  }
5380  }
5381  }
5382  return;
5383 }
double sg(long S)
void bessel_i0_i1(double x, double *i0val, double *i1val)
double bessel_k0_scaled(double x)
static const double b1_QQ[7]
static const double BESSEL_K1_Q1[]
string MD5datastream(fstream &ioFile)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
long genrand_int31()
static int initf
static double x2[63]
static double erf_P[]
double bessel_i1_scaled(double x)
bool is_odd(int j)
Definition: cddefines.h:753
double ratevl_7_8(double x, const double a[7], const double b[8])
T * ptr0()
Definition: vectorize.h:240
static unsigned long * nexxt
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
void chbfit(double a, double b, vector< double > &c, double(*func)(double))
double genrand_real1()
static double x1[83]
static const double b0_YP[8]
t_cpu_i & i()
Definition: cpu.h:419
double ratevl_15_6(double x, const double a[15], const double b[6])
double e1(double x)
static const double SQ2OPI
void svd(const int nm, const int m, const int n, double *a, double *w, bool matu, double *u, bool matv, double *v, int *ierr, double *rv1)
static const int N
double lfactorial(long n)
Definition: thirdparty.cpp:756
double bessel_i1(double x)
double SixJFull(long j1, long j2, long j3, long j4, long j5, long j6)
realnum FastVoigtH(realnum a, realnum v)
static const double BESSEL_K1_P3[]
double p1evl(double x, const double coef[], int N)
T sign(T x, T y)
Definition: cddefines.h:842
bool Triangle2(long J1, long J2, long J3)
Definition: thirdparty.h:497
static const double BESSEL_K1_Q3[]
access_scheme
Definition: cpu.h:262
static const double b0_PP[7]
Definition: thirdparty.cpp:973
#define F4(x, y, z)
static const double BESSEL_I1_Q2[]
static const double igam_big
static const double BESSEL_I0_Q1[]
static const double DR1
uint32 rotlFixed(uint32 x, unsigned int y)
static const double TWOOPI
static const double b0_QP[8]
Definition: thirdparty.cpp:993
static const double elk_Q[]
double Delta(long j1, long j2, long j3)
unsigned long genrand_int32()
static const double b1_PQ[7]
STATIC void order(long, realnum[], long *, long *, long *)
static const double C1
FILE * ioQQQ
Definition: cddefines.cpp:7
string MD5datafile(const char *fnam, access_scheme scheme)
double sjs(long j1, long j2, long j3, long l1, long l2, long l3)
static const double BESSEL_I0_Q2[]
void bessel_i0_i1_scaled(double x, double *i0val, double *i1val)
double expn(int n, double x)
static const double BESSEL_K1_P1[]
complex< double > cdgamma(complex< double > x)
Definition: thirdparty.cpp:789
static const int MXDSF
Definition: thirdparty.h:28
double dsexp(double x)
Definition: service.cpp:1038
static double erf_Q[]
static double erf_S[]
static const double BESSEL_K0_P3[]
bool little_endian() const
Definition: cpu.h:355
static const double PIO4
double ratevl_6_3(double x, const double a[6], const double b[3])
static t_lfact & Inst()
Definition: cddefines.h:209
void bessel_k0_k1(double x, double *k0val, double *k1val)
void vexp(const double x[], double y[], long nlo, long nhi)
static const int NPRE_FACTORIAL
Definition: thirdparty.h:25
double dawson(double x, int order)
#define F1(x, y, z)
double ratevl_11_10(double x, const double a[11], const double b[10])
void init_by_array(unsigned long init_key[], int key_length)
double bessel_k1(double x)
double erfc(double)
static const unsigned long MATRIX_A
static realnum b0[83]
static const double tbl_dawson[N_DAWSON+1]
static double b2[63]
const ios_base::openmode mode_r
Definition: cpu.h:267
static double erf_R[]
double bessel_i0_scaled(double x)
#define STATIC
Definition: cddefines.h:118
static const double b0_QQ[7]
static const double DR2
double bessel_i0(double x)
double genrand_real2()
double bessel_jn(int n, double x)
static const int M
void bessel_k0_k1_scaled(double x, double *k0val, double *k1val)
static const double b0_PQ[7]
Definition: thirdparty.cpp:983
static const double BESSEL_K0_P1[]
float realnum
Definition: cddefines.h:124
double fc2_scl(long n)
Definition: thirdparty.cpp:715
#define EXIT_FAILURE
Definition: cddefines.h:168
static const double THPIO4
double bessel_y0(double x)
STATIC void MD5_Transform(uint32 *digest, const uint32 *in)
double ratevl_5_4(double x, const double a[5], const double b[4])
string MD5file(const char *fnam, access_scheme scheme)
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
static const double BESSEL_K0_Q1[]
double expn2_scaled(double x)
long max(int a, long b)
Definition: cddefines.h:817
#define cdEXIT(FAIL)
Definition: cddefines.h:482
static double a1[83]
double bessel_j0(double x)
double powi(double, long int)
Definition: service.cpp:594
static const double BIG
static realnum b1[83]
long min(int a, long b)
Definition: cddefines.h:762
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
static const double b1_PP[7]
double bessel_j1(double x)
static const double pre_factorial[NPRE_FACTORIAL]
Definition: thirdparty.cpp:181
static const double BESSEL_I1_P1[]
t_state state
Definition: state.cpp:18
double igamc(double a, double x)
double ratevl_8_7(double x, const double a[8], const double b[7])
static const double BESSEL_K0_P2[]
static const double b1_RP[4]
bool linfit(long n, const double xorg[], const double yorg[], double &a, double &siga, double &b, double &sigb)
Definition: thirdparty.cpp:46
void test_expn()
static const unsigned long UMASK
static const double BESSEL_K1_P2[]
double igam(double a, double x)
static const double b1_QP[8]
double ratevl_10_11(double x, const double a[10], const double b[11])
double factorial(long n)
Definition: thirdparty.cpp:356
STATIC double igamc_fraction(double a, double x)
#define ASSERT(exp)
Definition: cddefines.h:613
double fc2(long n2)
double bessel_yn(int n, double x)
static const double BESSEL_I1_P2[]
static void next_state()
double bessel_k1_scaled(double x)
double erfce(double x)
double chbevl(double, const double[], int)
T pow2(T a)
Definition: cddefines.h:981
unsigned long TWIST(unsigned long u, unsigned long v)
static int nleft
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define F3(x, y, z)
static const double b1_YQ[8]
unsigned long MIXBITS(unsigned long u, unsigned long v)
void humlik(int n, const realnum x[], realnum y, realnum k[])
double bessel_y1(double x)
double genrand_res53()
double erf(double)
#define MD5STEP(f, w, x, y, z, data, s)
double e1_scaled(double x)
static const unsigned long LMASK
void rec6j(double *sixcof, double l2, double l3, double l4, double l5, double l6, double *l1min, double *l1max, double *lmatch, long ndim, long *ier)
vector< double > p_lf
Definition: thirdparty.cpp:739
const int N_DAWSON
double ratevl_6_4(double x, const double a[6], const double b[4])
double ellpk(double x)
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
static const double BESSEL_K1_Q2[]
#define F2(x, y, z)
static const double BESSEL_K0_Q3[]
double igamc_scaled(double a, double x)
double polevl(double x, const double coef[], int N)
double e2(double x)
static const double Z2
#define S(I_, J_)
static const double Z1
static const double b0_YQ[7]
double get_lfact(unsigned long n)
Definition: thirdparty.cpp:741
void init_genrand(unsigned long s)
static const double BESSEL_I0_P2[]
static const double BESSEL_K0_Q2[]
double gegenbauer(long n, double al, double x)
static const double elk_P[]
static const double BESSEL_I0_P1[]
STATIC uint32 MD5swap(uint32 word)
static t_cpu cpu
Definition: cpu.h:427
static double b0_RP[4]
static const double b1_YP[6]
double bessel_k0(double x)
double genrand_real3()
double frac(double d)
string MD5string(const string &str)
static double a0[83]
static const double MAXLOG
static double b0_RQ[8]
static double a2[63]
static const double b1_RQ[8]
static const double igam_biginv
static const double dsf[MXDSF]
Definition: thirdparty.cpp:369
static const double BESSEL_I1_Q1[]