76static double s2pi = 2.50662827463100050242E0;
78static double P0[5] = {
79-5.99633501014107895267E1,
80 9.80010754185999661536E1,
81-5.66762857469070293439E1,
82 1.39312609387279679503E1,
83-1.23916583867381258016E0,
85static double Q0[8] = {
86 1.95448858338141759834E0,
87 4.67627912898881538453E0,
88 8.63602421390890590575E1,
89-2.25462687854119370527E2,
90 2.00260212380060660359E2,
91-8.20372256168333339912E1,
92 1.59056225126211695515E1,
93-1.18331621121330003142E0,
95static double P1[9] = {
96 4.05544892305962419923E0,
97 3.15251094599893866154E1,
98 5.71628192246421288162E1,
99 4.40805073893200834700E1,
100 1.46849561928858024014E1,
101 2.18663306850790267539E0,
102-1.40256079171354495875E-1,
103-3.50424626827848203418E-2,
104-8.57456785154685413611E-4,
106static double Q1[8] = {
107 1.57799883256466749731E1,
108 4.53907635128879210584E1,
109 4.13172038254672030440E1,
110 1.50425385692907503408E1,
111 2.50464946208309415979E0,
112-1.42182922854787788574E-1,
113-3.80806407691578277194E-2,
114-9.33259480895457427372E-4,
116static double P2[9] = {
117 3.23774891776946035970E0,
118 6.91522889068984211695E0,
119 3.93881025292474443415E0,
120 1.33303460815807542389E0,
121 2.01485389549179081538E-1,
122 1.23716634817820021358E-2,
123 3.01581553508235416007E-4,
124 2.65806974686737550832E-6,
125 6.23974539184983293730E-9,
127static double Q2[8] = {
128 6.02427039364742014255E0,
129 3.67983563856160859403E0,
130 1.37702099489081330271E0,
131 2.16236993594496635890E-1,
132 1.34204006088543189037E-2,
133 3.28014464682127739104E-4,
134 2.89247864745380683936E-6,
135 6.79019408009981274425E-9,
139 double x,
y, z,
y2, x0,
x1;
142 return( - std::numeric_limits<double>::infinity() );
144 return( + std::numeric_limits<double>::infinity() );
147 if(
y > (1.0 - 0.13533528323661269189) )
152 if(
y > 0.13533528323661269189 )
160 x = std::sqrt( -2.0 * std::log(
y) );
161 x0 =
x - std::log(
x)/
x;
227 double x0,
x1,
x, yl, yh,
y,
d, lgm, dithresh;
232 MATH_ERROR_MSG(
"Cephes::igami",
"Wrong domain for parameter a (must be > 0)");
237 return std::numeric_limits<double>::infinity();
246 static double kMAXNUM = std::numeric_limits<double>::max();
255 y = ( 1.0 -
d -
ndtri(y0) * std::sqrt(
d) );
260 for( i=0; i<10; i++ )
262 if(
x > x0 ||
x <
x1 )
265 if( y < yl || y > yh )
278 d = (
a - 1.0) * std::log(
x) -
x - lgm;
297 while( x0 == kMAXNUM )
313 for( i=0; i<400; i++ )
317 lgm = (x0 -
x1)/(
x1 + x0);
318 if( std::fabs(lgm) < dithresh )
321 if( std::fabs(lgm) < dithresh )
337 d = (y0 - yl)/(yh - yl);
352 d = (y0 - yl)/(yh - yl);
411double incbi(
double aa,
double bb,
double yy0 )
413 double a,
b, y0,
d,
y,
x, x0,
x1, lgm, yp, di, dithresh, yl, yh, xt;
414 int i, rflg, dir, nflg;
418 MATH_ERROR_MSG(
"Cephes::incbi",
"Wrong domain for parameter a (must be > 0)");
422 MATH_ERROR_MSG(
"Cephes::incbi",
"Wrong domain for parameter b (must be > 0)");
438 if( aa <= 1.0 || bb <= 1.0 )
473 lgm = (yp * yp - 3.0)/6.0;
474 x = 2.0/( 1.0/(2.0*
a-1.0) + 1.0/(2.0*
b-1.0) );
475 d = yp * std::sqrt(
x + lgm ) /
x
476 - ( 1.0/(2.0*
b-1.0) - 1.0/(2.0*
a-1.0) )
477 * (lgm + 5.0/6.0 - 2.0/(3.0*
x));
484 x =
a/(
a +
b * std::exp(
d) );
487 if( std::fabs(yp) < 0.2 )
495 for( i=0; i<100; i++ )
499 x = x0 + di * (
x1 - x0);
505 x = x0 + di * (
x1 - x0);
510 yp = (
x1 - x0)/(
x1 + x0);
511 if( std::fabs(yp) < dithresh )
514 if( std::fabs(yp) < dithresh )
527 di = 1.0 - (1.0 - di) * (1.0 - di);
531 di = (y0 -
y)/(yh - yl);
577 di = (
y - y0)/(yh - yl);
627 if(
x == 1.0 ||
x == 0.0 )
630 d = (
a - 1.0) * std::log(
x) + (
b - 1.0) * std::log(1.0-
x) + lgm;
641 y = (
x - x0) / (
x1 - x0);
642 xt = x0 + 0.5 *
y * (
x - x0);
649 xt =
x1 - 0.5 *
y * (
x1 -
x);
#define MATH_ERROR_MSG(loc, str)
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Namespace for new Math classes and functions.
double incbet(double aa, double bb, double xx)
DESCRIPTION:
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
double igami(double a, double y)
double incbi(double a, double b, double y)
double Polynomial1eval(double x, double *a, unsigned int N)
double Polynomialeval(double x, double *a, unsigned int N)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...