Logo ROOT   6.18/05
Reference Guide
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
20namespace ROOT {
21namespace Math {
22
23namespace Cephes {
24
25
26static double kBig = 4.503599627370496e15;
27static double kBiginv = 2.22044604925031308085e-16;
28
29/* log( sqrt( 2*pi ) ) */
30static 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
51double 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
127double 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
170static double A[] = {
171 8.11614167470508450300E-4,
172 -5.95061904284301438324E-4,
173 7.93650340457716943945E-4,
174 -2.77777777730099687205E-3,
175 8.33333333333331927722E-2
176};
177
178static double B[] = {
179 -1.37825152569120859100E3,
180 -3.88016315134637840924E4,
181 -3.31612992738871184744E5,
182 -1.16237097492762307383E6,
183 -1.72173700820839662146E6,
184 -8.53555664245765465627E5
185};
186
187static 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
197double 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/*---------------------------------------------------------------------------*/
285static 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};
294static 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 */
306static 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 */
316static 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
339double 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
415small:
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 */
428double beta(double z, double w)
429{
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
484double 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;
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
562done:
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
581double 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
663cdone:
664 return(ans);
665}
666
667
668/*---------------------------------------------------------------------------*/
669
670/* Continued fraction expansion #2
671 * for incomplete beta integral
672 */
673
674double 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 );
756cdone:
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
766double 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
813static 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};
824static 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};
835static double erfR[] = {
836 5.64189583547755073984E-1,
837 1.27536670759978104416E0,
838 5.01905042251180477414E0,
839 6.16021097993053585195E0,
840 7.40974269950448939160E0,
841 2.97886665372100240670E0
842};
843static double erfS[] = {
844/* 1.00000000000000000000E0,*/
845 2.26052863220117276590E0,
846 9.39603524938001434673E0,
847 1.20489539808096656605E1,
848 1.70814450747565897222E1,
849 9.60896809063285878198E0,
850 3.36907645100081516050E0
851};
852static double erfT[] = {
853 9.60497373987051638749E0,
854 9.00260197203842689217E1,
855 2.23200534594684319226E3,
856 7.00332514112805075473E3,
857 5.55923013010394962768E4
858};
859static 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
874double 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
926double 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*/
951double 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*/
967double 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
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define SQTPI
#define kMACHEP
#define kMAXLOG
#define kMAXLGM
#define kMAXSTIR
#define kMINLOG
#define N
float * q
Definition: THbookFile.cxx:87
double pow(double, double)
double floor(double)
double sin(double)
double exp(double)
double log(double)
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
static double kBiginv
double incbd(double a, double b, double x)
static double B[]
static double STIR[5]
double erfc(double a)
static double erfU[]
static double erfR[]
static double kBig
static double erfS[]
double erf(double x)
double incbcf(double a, double b, double x)
static double Q[]
double incbet(double aa, double bb, double xx)
DESCRIPTION:
double pseries(double a, double b, double x)
double igam(double a, double x)
static double P[]
double lgam(double x)
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
static double erfT[]
static double A[]
double beta(double z, double w)
static double LS2PI
static double erfP[]
static double erfQ[]
static double C[]
static double stirf(double x)
double gamma(double x)
double Polynomial1eval(double x, double *a, unsigned int N)
double Pi()
Mathematical constants.
Definition: Math.h:88
double Polynomialeval(double x, double *a, unsigned int N)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static constexpr double s
auto * a
Definition: textangle.C:12
auto * t1
Definition: textangle.C:20