Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
SpecFuncCephes.cxx
Go to the documentation of this file.
1 //
2 //
3 // gamma and related functions from Cephes library
4 // see: http://www.netlib.org/cephes
5 //
6 // Copyright 1985, 1987, 2000 by Stephen L. Moshier
7 //
8 //
9 
10 #include "SpecFuncCephes.h"
11 #include "Math/Math.h"
12 
13 
14 #include <cmath>
15 
16 #include <limits>
17 
18 
19 
20 namespace ROOT {
21 namespace Math {
22 
23 namespace Cephes {
24 
25 
26 static double kBig = 4.503599627370496e15;
27 static double kBiginv = 2.22044604925031308085e-16;
28 
29 /* log( sqrt( 2*pi ) ) */
30 static double LS2PI = 0.91893853320467274178;
31 
32 
33 // incomplete gamma function (complement integral)
34 // igamc(a,x) = 1 - igam(a,x)
35 //
36 // inf.
37 // -
38 // 1 | | -t a-1
39 // = ----- | e t dt.
40 // - | |
41 // | (a) -
42 // x
43 //
44 //
45 
46 // In this implementation both arguments must be positive.
47 // The integral is evaluated by either a power series or
48 // continued fraction expansion, depending on the relative
49 // values of a and x.
50 
51 double igamc( double a, double x )
52 {
53 
54  double ans, ax, c, yc, r, t, y, z;
55  double pk, pkm1, pkm2, qk, qkm1, qkm2;
56 
57  // LM: for negative values returns 0.0
58  // This is correct if a is a negative integer since Gamma(-n) = +/- inf
59  if (a <= 0) return 0.0;
60 
61  if (x <= 0) return 1.0;
62 
63  if( (x < 1.0) || (x < a) )
64  return( 1.0 - igam(a,x) );
65 
66  ax = a * std::log(x) - x - lgam(a);
67  if( ax < -kMAXLOG )
68  return( 0.0 );
69 
70  ax = std::exp(ax);
71 
72 /* continued fraction */
73  y = 1.0 - a;
74  z = x + y + 1.0;
75  c = 0.0;
76  pkm2 = 1.0;
77  qkm2 = x;
78  pkm1 = x + 1.0;
79  qkm1 = z * x;
80  ans = pkm1/qkm1;
81 
82  do
83  {
84  c += 1.0;
85  y += 1.0;
86  z += 2.0;
87  yc = y * c;
88  pk = pkm1 * z - pkm2 * yc;
89  qk = qkm1 * z - qkm2 * yc;
90  if(qk)
91  {
92  r = pk/qk;
93  t = std::abs( (ans - r)/r );
94  ans = r;
95  }
96  else
97  t = 1.0;
98  pkm2 = pkm1;
99  pkm1 = pk;
100  qkm2 = qkm1;
101  qkm1 = qk;
102  if( std::abs(pk) > kBig )
103  {
104  pkm2 *= kBiginv;
105  pkm1 *= kBiginv;
106  qkm2 *= kBiginv;
107  qkm1 *= kBiginv;
108  }
109  }
110  while( t > kMACHEP );
111 
112  return( ans * ax );
113 }
114 
115 
116 
117 /* left tail of incomplete gamma function:
118  *
119  * inf. k
120  * a -x - x
121  * x e > ----------
122  * - -
123  * k=0 | (a+k+1)
124  *
125  */
126 
127 double igam( double a, double x )
128 {
129  double ans, ax, c, r;
130 
131  // LM: for negative values returns 1.0 instead of zero
132  // This is correct if a is a negative integer since Gamma(-n) = +/- inf
133  if (a <= 0) return 1.0;
134 
135  if (x <= 0) return 0.0;
136 
137  if( (x > 1.0) && (x > a ) )
138  return( 1.0 - igamc(a,x) );
139 
140 /* Compute x**a * exp(-x) / gamma(a) */
141  ax = a * std::log(x) - x - lgam(a);
142  if( ax < -kMAXLOG )
143  return( 0.0 );
144 
145  ax = std::exp(ax);
146 
147 /* power series */
148  r = a;
149  c = 1.0;
150  ans = 1.0;
151 
152  do
153  {
154  r += 1.0;
155  c *= x/r;
156  ans += c;
157  }
158  while( c/ans > kMACHEP );
159 
160  return( ans * ax/a );
161 }
162 
163 /*---------------------------------------------------------------------------*/
164 
165 /* Logarithm of gamma function */
166 /* A[]: Stirling's formula expansion of log gamma
167  * B[], C[]: log gamma function between 2 and 3
168  */
169 
170 static double A[] = {
171  8.11614167470508450300E-4,
172  -5.95061904284301438324E-4,
173  7.93650340457716943945E-4,
174  -2.77777777730099687205E-3,
175  8.33333333333331927722E-2
176 };
177 
178 static double B[] = {
179  -1.37825152569120859100E3,
180  -3.88016315134637840924E4,
181  -3.31612992738871184744E5,
182  -1.16237097492762307383E6,
183  -1.72173700820839662146E6,
184  -8.53555664245765465627E5
185 };
186 
187 static double C[] = {
188 /* 1.00000000000000000000E0, */
189  -3.51815701436523470549E2,
190  -1.70642106651881159223E4,
191  -2.20528590553854454839E5,
192  -1.13933444367982507207E6,
193  -2.53252307177582951285E6,
194  -2.01889141433532773231E6
195 };
196 
197 double lgam( double x )
198 {
199  double p, q, u, w, z;
200  int i;
201 
202  int sgngam = 1;
203 
204  if (x >= std::numeric_limits<double>::infinity())
205  return(std::numeric_limits<double>::infinity());
206 
207  if( x < -34.0 )
208  {
209  q = -x;
210  w = lgam(q);
211  p = std::floor(q);
212  if( p==q )//_unur_FP_same(p,q)
213  return (std::numeric_limits<double>::infinity());
214  i = (int) p;
215  if( (i & 1) == 0 )
216  sgngam = -1;
217  else
218  sgngam = 1;
219  z = q - p;
220  if( z > 0.5 )
221  {
222  p += 1.0;
223  z = p - q;
224  }
225  z = q * std::sin( ROOT::Math::Pi() * z );
226  if( z == 0 )
227  return (std::numeric_limits<double>::infinity());
228 /* z = log(ROOT::Math::Pi()) - log( z ) - w;*/
229  z = std::log(ROOT::Math::Pi()) - std::log( z ) - w;
230  return( z );
231  }
232 
233  if( x < 13.0 )
234  {
235  z = 1.0;
236  p = 0.0;
237  u = x;
238  while( u >= 3.0 )
239  {
240  p -= 1.0;
241  u = x + p;
242  z *= u;
243  }
244  while( u < 2.0 )
245  {
246  if( u == 0 )
247  return (std::numeric_limits<double>::infinity());
248  z /= u;
249  p += 1.0;
250  u = x + p;
251  }
252  if( z < 0.0 )
253  {
254  sgngam = -1;
255  z = -z;
256  }
257  else
258  sgngam = 1;
259  if( u == 2.0 )
260  return( std::log(z) );
261  p -= 2.0;
262  x = x + p;
263  p = x * Polynomialeval(x, B, 5 ) / Polynomial1eval( x, C, 6);
264  return( std::log(z) + p );
265  }
266 
267  if( x > kMAXLGM )
268  return( sgngam * std::numeric_limits<double>::infinity() );
269 
270  q = ( x - 0.5 ) * std::log(x) - x + LS2PI;
271  if( x > 1.0e8 )
272  return( q );
273 
274  p = 1.0/(x*x);
275  if( x >= 1000.0 )
276  q += (( 7.9365079365079365079365e-4 * p
277  - 2.7777777777777777777778e-3) *p
278  + 0.0833333333333333333333) / x;
279  else
280  q += Polynomialeval( p, A, 4 ) / x;
281  return( q );
282 }
283 
284 /*---------------------------------------------------------------------------*/
285 static double P[] = {
286  1.60119522476751861407E-4,
287  1.19135147006586384913E-3,
288  1.04213797561761569935E-2,
289  4.76367800457137231464E-2,
290  2.07448227648435975150E-1,
291  4.94214826801497100753E-1,
292  9.99999999999999996796E-1
293 };
294 static double Q[] = {
295  -2.31581873324120129819E-5,
296  5.39605580493303397842E-4,
297  -4.45641913851797240494E-3,
298  1.18139785222060435552E-2,
299  3.58236398605498653373E-2,
300  -2.34591795718243348568E-1,
301  7.14304917030273074085E-2,
302  1.00000000000000000320E0
303 };
304 
305 /* Stirling's formula for the gamma function */
306 static double STIR[5] = {
307  7.87311395793093628397E-4,
308  -2.29549961613378126380E-4,
309  -2.68132617805781232825E-3,
310  3.47222221605458667310E-3,
311  8.33333333333482257126E-2,
312 };
313 
314 #define SQTPI std::sqrt(2*ROOT::Math::Pi()) /* sqrt(2*pi) */
315 /* Stirling formula for the gamma function */
316 static double stirf( double x)
317 {
318  double y, w, v;
319 
320  w = 1.0/x;
321  w = 1.0 + w * Polynomialeval( w, STIR, 4 );
322  y = exp(x);
323 
324 /* #define kMAXSTIR kMAXLOG/log(kMAXLOG) */
325 
326  if( x > kMAXSTIR )
327  { /* Avoid overflow in pow() */
328  v = pow( x, 0.5 * x - 0.25 );
329  y = v * (v / y);
330  }
331  else
332  {
333  y = pow( x, x - 0.5 ) / y;
334  }
335  y = SQTPI * y * w;
336  return( y );
337 }
338 
339 double gamma( double x )
340 {
341  double p, q, z;
342  int i;
343 
344  int sgngam = 1;
345 
346  if (x >=std::numeric_limits<double>::infinity())
347  return(x);
348 
349  q = std::abs(x);
350 
351  if( q > 33.0 )
352  {
353  if( x < 0.0 )
354  {
355  p = std::floor(q);
356  if( p == q )
357  {
358  return( sgngam * std::numeric_limits<double>::infinity());
359  }
360  i = (int) p;
361  if( (i & 1) == 0 )
362  sgngam = -1;
363  z = q - p;
364  if( z > 0.5 )
365  {
366  p += 1.0;
367  z = q - p;
368  }
369  z = q * std::sin( ROOT::Math::Pi() * z );
370  if( z == 0 )
371  {
372  return( sgngam * std::numeric_limits<double>::infinity());
373  }
374  z = std::abs(z);
375  z = ROOT::Math::Pi()/(z * stirf(q) );
376  }
377  else
378  {
379  z = stirf(x);
380  }
381  return( sgngam * z );
382  }
383 
384  z = 1.0;
385  while( x >= 3.0 )
386  {
387  x -= 1.0;
388  z *= x;
389  }
390 
391  while( x < 0.0 )
392  {
393  if( x > -1.E-9 )
394  goto small;
395  z /= x;
396  x += 1.0;
397  }
398 
399  while( x < 2.0 )
400  {
401  if( x < 1.e-9 )
402  goto small;
403  z /= x;
404  x += 1.0;
405  }
406 
407  if( x == 2.0 )
408  return(z);
409 
410  x -= 2.0;
411  p = Polynomialeval( x, P, 6 );
412  q = Polynomialeval( x, Q, 7 );
413  return( z * p / q );
414 
415 small:
416  if( x == 0 )
417  return( std::numeric_limits<double>::infinity() );
418  else
419  return( z/((1.0 + 0.5772156649015329 * x) * x) );
420 }
421 
422 /*---------------------------------------------------------------------------*/
423 
424 //#define kMAXLGM 2.556348e305
425 
426 /*---------------------------------------------------------------------------*/
427 /* implementation based on cephes log gamma */
428 double beta(double z, double w)
429 {
430  return std::exp(ROOT::Math::Cephes::lgam(z) + ROOT::Math::Cephes::lgam(w)- ROOT::Math::Cephes::lgam(z+w) );
431 }
432 
433 
434 /*---------------------------------------------------------------------------*/
435 /* inplementation of the incomplete beta function */
436 /**
437  * DESCRIPTION:
438  *
439  * Returns incomplete beta integral of the arguments, evaluated
440  * from zero to x. The function is defined as
441  *
442  * x
443  * - -
444  * | (a+b) | | a-1 b-1
445  * ----------- | t (1-t) dt.
446  * - - | |
447  * | (a) | (b) -
448  * 0
449  *
450  * The domain of definition is 0 <= x <= 1. In this
451  * implementation a and b are restricted to positive values.
452  * The integral from x to 1 may be obtained by the symmetry
453  * relation
454  *
455  * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
456  *
457  * The integral is evaluated by a continued fraction expansion
458  * or, when b*x is small, by a power series.
459  *
460  * ACCURACY:
461  *
462  * Tested at uniformly distributed random points (a,b,x) with a and b
463  * in "domain" and x between 0 and 1.
464  * Relative error
465  * arithmetic domain # trials peak rms
466  * IEEE 0,5 10000 6.9e-15 4.5e-16
467  * IEEE 0,85 250000 2.2e-13 1.7e-14
468  * IEEE 0,1000 30000 5.3e-12 6.3e-13
469  * IEEE 0,10000 250000 9.3e-11 7.1e-12
470  * IEEE 0,100000 10000 8.7e-10 4.8e-11
471  * Outputs smaller than the IEEE gradual underflow threshold
472  * were excluded from these statistics.
473  *
474  * ERROR MESSAGES:
475  * message condition value returned
476  * incbet domain x<0, x>1 0.0
477  * incbet underflow 0.0
478  *
479  * Cephes Math Library, Release 2.8: June, 2000
480  * Copyright 1984, 1995, 2000 by Stephen L. Moshier
481  */
482 
483 
484 double incbet( double aa, double bb, double xx )
485 {
486  double a, b, t, x, xc, w, y;
487  int flag;
488 
489  if( aa <= 0.0 || bb <= 0.0 )
490  return( 0.0 );
491 
492  // LM: changed: for X > 1 return 1.
493  if (xx <= 0.0) return( 0.0 );
494  if ( xx >= 1.0) return( 1.0 );
495 
496  flag = 0;
497 
498 /* - to test if that way is better for large b/ (comment out from Cephes version)
499  if( (bb * xx) <= 1.0 && xx <= 0.95)
500  {
501  t = pseries(aa, bb, xx);
502  goto done;
503  }
504 
505 **/
506  w = 1.0 - xx;
507 
508 /* Reverse a and b if x is greater than the mean. */
509 /* aa,bb > 1 -> sharp rise at x=aa/(aa+bb) */
510  if( xx > (aa/(aa+bb)) )
511  {
512  flag = 1;
513  a = bb;
514  b = aa;
515  xc = xx;
516  x = w;
517  }
518  else
519  {
520  a = aa;
521  b = bb;
522  xc = w;
523  x = xx;
524  }
525 
526  if( flag == 1 && (b * x) <= 1.0 && x <= 0.95)
527  {
528  t = pseries(a, b, x);
529  goto done;
530  }
531 
532 /* Choose expansion for better convergence. */
533  y = x * (a+b-2.0) - (a-1.0);
534  if( y < 0.0 )
535  w = incbcf( a, b, x );
536  else
537  w = incbd( a, b, x ) / xc;
538 
539 /* Multiply w by the factor
540  a b _ _ _
541  x (1-x) | (a+b) / ( a | (a) | (b) ) . */
542 
543  y = a * std::log(x);
544  t = b * std::log(xc);
545  if( (a+b) < kMAXSTIR && std::abs(y) < kMAXLOG && std::abs(t) < kMAXLOG )
546  {
547  t = pow(xc,b);
548  t *= pow(x,a);
549  t /= a;
550  t *= w;
551  t *= ROOT::Math::Cephes::gamma(a+b) / (ROOT::Math::Cephes::gamma(a) * ROOT::Math::Cephes::gamma(b));
552  goto done;
553  }
554 /* Resort to logarithms. */
555  y += t + lgam(a+b) - lgam(a) - lgam(b);
556  y += std::log(w/a);
557  if( y < kMINLOG )
558  t = 0.0;
559  else
560  t = std::exp(y);
561 
562 done:
563 
564  if( flag == 1 )
565  {
566  if( t <= kMACHEP )
567  t = 1.0 - kMACHEP;
568  else
569  t = 1.0 - t;
570  }
571  return( t );
572 }
573 /*---------------------------------------------------------------------------*/
574 
575 /*---------------------------------------------------------------------------*/
576 
577 /* Continued fraction expansion #1
578  * for incomplete beta integral
579  */
580 
581 double incbcf( double a, double b, double x )
582 {
583  double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
584  double k1, k2, k3, k4, k5, k6, k7, k8;
585  double r, t, ans, thresh;
586  int n;
587 
588  k1 = a;
589  k2 = a + b;
590  k3 = a;
591  k4 = a + 1.0;
592  k5 = 1.0;
593  k6 = b - 1.0;
594  k7 = k4;
595  k8 = a + 2.0;
596 
597  pkm2 = 0.0;
598  qkm2 = 1.0;
599  pkm1 = 1.0;
600  qkm1 = 1.0;
601  ans = 1.0;
602  r = 1.0;
603  n = 0;
604  thresh = 3.0 * kMACHEP;
605  do
606  {
607 
608  xk = -( x * k1 * k2 )/( k3 * k4 );
609  pk = pkm1 + pkm2 * xk;
610  qk = qkm1 + qkm2 * xk;
611  pkm2 = pkm1;
612  pkm1 = pk;
613  qkm2 = qkm1;
614  qkm1 = qk;
615 
616  xk = ( x * k5 * k6 )/( k7 * k8 );
617  pk = pkm1 + pkm2 * xk;
618  qk = qkm1 + qkm2 * xk;
619  pkm2 = pkm1;
620  pkm1 = pk;
621  qkm2 = qkm1;
622  qkm1 = qk;
623 
624  if( qk !=0 )
625  r = pk/qk;
626  if( r != 0 )
627  {
628  t = std::abs( (ans - r)/r );
629  ans = r;
630  }
631  else
632  t = 1.0;
633 
634  if( t < thresh )
635  goto cdone;
636 
637  k1 += 1.0;
638  k2 += 1.0;
639  k3 += 2.0;
640  k4 += 2.0;
641  k5 += 1.0;
642  k6 -= 1.0;
643  k7 += 2.0;
644  k8 += 2.0;
645 
646  if( (std::abs(qk) + std::abs(pk)) > kBig )
647  {
648  pkm2 *= kBiginv;
649  pkm1 *= kBiginv;
650  qkm2 *= kBiginv;
651  qkm1 *= kBiginv;
652  }
653  if( (std::abs(qk) < kBiginv) || (std::abs(pk) < kBiginv) )
654  {
655  pkm2 *= kBig;
656  pkm1 *= kBig;
657  qkm2 *= kBig;
658  qkm1 *= kBig;
659  }
660  }
661  while( ++n < 300 );
662 
663 cdone:
664  return(ans);
665 }
666 
667 
668 /*---------------------------------------------------------------------------*/
669 
670 /* Continued fraction expansion #2
671  * for incomplete beta integral
672  */
673 
674 double incbd( double a, double b, double x )
675 {
676  double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
677  double k1, k2, k3, k4, k5, k6, k7, k8;
678  double r, t, ans, z, thresh;
679  int n;
680 
681  k1 = a;
682  k2 = b - 1.0;
683  k3 = a;
684  k4 = a + 1.0;
685  k5 = 1.0;
686  k6 = a + b;
687  k7 = a + 1.0;;
688  k8 = a + 2.0;
689 
690  pkm2 = 0.0;
691  qkm2 = 1.0;
692  pkm1 = 1.0;
693  qkm1 = 1.0;
694  z = x / (1.0-x);
695  ans = 1.0;
696  r = 1.0;
697  n = 0;
698  thresh = 3.0 * kMACHEP;
699  do
700  {
701 
702  xk = -( z * k1 * k2 )/( k3 * k4 );
703  pk = pkm1 + pkm2 * xk;
704  qk = qkm1 + qkm2 * xk;
705  pkm2 = pkm1;
706  pkm1 = pk;
707  qkm2 = qkm1;
708  qkm1 = qk;
709 
710  xk = ( z * k5 * k6 )/( k7 * k8 );
711  pk = pkm1 + pkm2 * xk;
712  qk = qkm1 + qkm2 * xk;
713  pkm2 = pkm1;
714  pkm1 = pk;
715  qkm2 = qkm1;
716  qkm1 = qk;
717 
718  if( qk != 0 )
719  r = pk/qk;
720  if( r != 0 )
721  {
722  t = std::abs( (ans - r)/r );
723  ans = r;
724  }
725  else
726  t = 1.0;
727 
728  if( t < thresh )
729  goto cdone;
730 
731  k1 += 1.0;
732  k2 -= 1.0;
733  k3 += 2.0;
734  k4 += 2.0;
735  k5 += 1.0;
736  k6 += 1.0;
737  k7 += 2.0;
738  k8 += 2.0;
739 
740  if( (std::abs(qk) + std::abs(pk)) > kBig )
741  {
742  pkm2 *= kBiginv;
743  pkm1 *= kBiginv;
744  qkm2 *= kBiginv;
745  qkm1 *= kBiginv;
746  }
747  if( (std::abs(qk) < kBiginv) || (std::abs(pk) < kBiginv) )
748  {
749  pkm2 *= kBig;
750  pkm1 *= kBig;
751  qkm2 *= kBig;
752  qkm1 *= kBig;
753  }
754  }
755  while( ++n < 300 );
756 cdone:
757  return(ans);
758 }
759 
760 
761 /*---------------------------------------------------------------------------*/
762 
763 /* Power series for incomplete beta integral.
764  Use when b*x is small and x not too close to 1. */
765 
766 double pseries( double a, double b, double x )
767 {
768  double s, t, u, v, n, t1, z, ai;
769 
770  ai = 1.0 / a;
771  u = (1.0 - b) * x;
772  v = u / (a + 1.0);
773  t1 = v;
774  t = u;
775  n = 2.0;
776  s = 0.0;
777  z = kMACHEP * ai;
778  while( std::abs(v) > z )
779  {
780  u = (n - b) * x / n;
781  t *= u;
782  v = t / (a + n);
783  s += v;
784  n += 1.0;
785  }
786  s += t1;
787  s += ai;
788 
789  u = a * log(x);
790  if( (a+b) < kMAXSTIR && std::abs(u) < kMAXLOG )
791  {
792  t = gamma(a+b)/(gamma(a)*gamma(b));
793  s = s * t * pow(x,a);
794  }
795  else
796  {
797  t = lgam(a+b) - lgam(a) - lgam(b) + u + std::log(s);
798  if( t < kMINLOG )
799  s = 0.0;
800  else
801  s = std::exp(t);
802  }
803  return(s);
804 }
805 
806 /*---------------------------------------------------------------------------*/
807 
808 
809 /*---------------------------------------------------------------------------*/
810 /* for evaluation of error function */
811 /*---------------------------------------------------------------------------*/
812 
813 static double erfP[] = {
814  2.46196981473530512524E-10,
815  5.64189564831068821977E-1,
816  7.46321056442269912687E0,
817  4.86371970985681366614E1,
818  1.96520832956077098242E2,
819  5.26445194995477358631E2,
820  9.34528527171957607540E2,
821  1.02755188689515710272E3,
822  5.57535335369399327526E2
823 };
824 static double erfQ[] = {
825 /* 1.00000000000000000000E0,*/
826  1.32281951154744992508E1,
827  8.67072140885989742329E1,
828  3.54937778887819891062E2,
829  9.75708501743205489753E2,
830  1.82390916687909736289E3,
831  2.24633760818710981792E3,
832  1.65666309194161350182E3,
833  5.57535340817727675546E2
834 };
835 static double erfR[] = {
836  5.64189583547755073984E-1,
837  1.27536670759978104416E0,
838  5.01905042251180477414E0,
839  6.16021097993053585195E0,
840  7.40974269950448939160E0,
841  2.97886665372100240670E0
842 };
843 static double erfS[] = {
844 /* 1.00000000000000000000E0,*/
845  2.26052863220117276590E0,
846  9.39603524938001434673E0,
847  1.20489539808096656605E1,
848  1.70814450747565897222E1,
849  9.60896809063285878198E0,
850  3.36907645100081516050E0
851 };
852 static double erfT[] = {
853  9.60497373987051638749E0,
854  9.00260197203842689217E1,
855  2.23200534594684319226E3,
856  7.00332514112805075473E3,
857  5.55923013010394962768E4
858 };
859 static double erfU[] = {
860 /* 1.00000000000000000000E0,*/
861  3.35617141647503099647E1,
862  5.21357949780152679795E2,
863  4.59432382970980127987E3,
864  2.26290000613890934246E4,
865  4.92673942608635921086E4
866 };
867 
868 /*---------------------------------------------------------------------------*/
869 /* complementary error function */
870 /* For small x, erfc(x) = 1 - erf(x); otherwise rational */
871 /* approximations are computed. */
872 
873 
874 double erfc( double a )
875 {
876  double p,q,x,y,z;
877 
878 
879  if( a < 0.0 )
880  x = -a;
881  else
882  x = a;
883 
884  if( x < 1.0 )
885  return( 1.0 - ROOT::Math::Cephes::erf(a) );
886 
887  z = -a * a;
888 
889  if( z < -kMAXLOG )
890  {
891  under:
892  if( a < 0 )
893  return( 2.0 );
894  else
895  return( 0.0 );
896  }
897 
898  z = exp(z);
899 
900  if( x < 8.0 )
901  {
902  p = Polynomialeval( x, erfP, 8 );
903  q = Polynomial1eval( x, erfQ, 8 );
904  }
905  else
906  {
907  p = Polynomialeval( x, erfR, 5 );
908  q = Polynomial1eval( x, erfS, 6 );
909  }
910  y = (z * p)/q;
911 
912  if( a < 0 )
913  y = 2.0 - y;
914 
915  if( y == 0 )
916  goto under;
917 
918  return(y);
919 }
920 
921 /*---------------------------------------------------------------------------*/
922 /* error function */
923 /* For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise */
924 /* erf(x) = 1 - erfc(x). */
925 
926 double erf( double x)
927 {
928  double y, z;
929 
930  if( std::abs(x) > 1.0 )
931  return( 1.0 - ROOT::Math::Cephes::erfc(x) );
932  z = x * x;
933  y = x * Polynomialeval( z, erfT, 4 ) / Polynomial1eval( z, erfU, 5 );
934  return( y );
935 
936 }
937 
938 } // end namespace Cephes
939 
940 
941 /*---------------------------------------------------------------------------*/
942 
943 /*---------------------------------------------------------------------------*/
944 /* Routines used within this implementation */
945 
946 
947 /*
948  * calculates a value of a polynomial of the form:
949  * a[0]x^N+a[1]x^(N-1) + ... + a[N]
950 */
951 double Polynomialeval(double x, double* a, unsigned int N)
952 {
953  if (N==0) return a[0];
954  else
955  {
956  double pom = a[0];
957  for (unsigned int i=1; i <= N; i++)
958  pom = pom *x + a[i];
959  return pom;
960  }
961 }
962 
963 /*
964  * calculates a value of a polynomial of the form:
965  * x^N+a[0]x^(N-1) + ... + a[N-1]
966 */
967 double Polynomial1eval(double x, double* a, unsigned int N)
968 {
969  if (N==0) return a[0];
970  else
971  {
972  double pom = x + a[0];
973  for (unsigned int i=1; i < N; i++)
974  pom = pom *x + a[i];
975  return pom;
976  }
977 }
978 
979 
980 } // end namespace Math
981 } // end namespace ROOT
982