17 void nrerror(
const char * msg) {
27 static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
31 tmp -= (x+0.5)*log(tmp);
32 ser=1.000000000190015;
33 for (j=0;j<=5;j++) ser += cof[j]/++y;
34 return -tmp+log(2.5066282746310005*ser/x);
39 double gammp(
double a,
double x)
41 void gcf(
double *gammcf,
double a,
double x,
double *gln);
42 void gser(
double *gamser,
double a,
double x,
double *gln);
44 double gamser,gammcf,gln;
45 if (x < 0.0 || a <= 0.0) nrerror(
"Invalid arguments in routine gammp");
47 gser(&gamser,a,x,&gln);
50 gcf(&gammcf,a,x,&gln);
56 double gammq(
double a,
double x)
58 void gcf(
double *gammcf,
double a,
double x,
double *gln);
59 void gser(
double *gamser,
double a,
double x,
double *gln);
61 double gamser,gammcf,gln;
62 if (x < 0.0 || a <= 0.0) nrerror(
"Invalid arguments in routine gammq");
64 gser(&gamser,a,x,&gln);
67 gcf(&gammcf,a,x,&gln);
74 void gser(
double *gamser,
double a,
double x,
double *gln)
82 if (x < 0.0) nrerror(
"x less than 0 in routine gser");
88 for (n=1;n<=
ITMAX;n++) {
92 if (fabs(del) < fabs(sum)*
EPS) {
93 *gamser=sum*
exp(-x+a*log(x)-(*gln));
97 nrerror(
"a too large, ITMAX too small in routine gser");
103 void gcf(
double *gammcf,
double a,
double x,
double *gln)
119 for (i=1;i<=
ITMAX;i++) {
129 if (fabs(del-1.0) <
EPS)
break;
131 if (i >
ITMAX) nrerror(
"a too large, ITMAX too small in gcf");
132 *gammcf=
exp(-x+a*log(x)-(*gln))*h;
139 double gammp(
double a,
double x);
140 return x < 0.0 ? -
gammp(0.5,x*x) :
gammp(0.5,x*x);
#define EPS
Relative accuracy.
#define FPMIN
Number near the smallest representable doubleing-point number.
double erf(double x)
Returns the error function erf(x).
double gammln(double xx)
Returns the value ln[?(xx)] for xx > 0.
static const double exp
exponential
bool str2bool(string str)
double gammq(double a, double x)
Returns the incomplete gamma function Q(a, x) ? 1 ? P(a, x).
#define ITMAX
Maximum allowed number of iterations.
void gcf(double *gammcf, double a, double x, double *gln)
Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as g...
void gser(double *gamser, double a, double x, double *gln)
double gammp(double a, double x)
Returns the incomplete gamma function P(a, x).