double t;   t = InGamma(0.5, 0.5 * x * x);   if (x < 0.0)     return (0.5 * (1.0 - t));   else     return (0.5 * (1.0 + t));}   static double idfStandard(double u)/* ===================================  * NOTE: 0.0 < u < 1.0  * =================================== */{    double t, x = 0.0;                    /* initialize to the mean, then  */   do {                                  /* use Newton-Raphson iteration  */     t = x;     x = t + (u - cdfStandard(t)) / pdfStandard(t);   } while (fabs(x - t) >= TINY);   return (x);}   double pdfNormal(double m, double s, double x)/* ==============================================  * NOTE: x and m can be any value, but s > 0.0  * ============================================== */{    double t = (x - m) / s;   return (pdfStandard(t) / s);}   double cdfNormal(double m, double s, double x)/* ==============================================  * NOTE: x and m can be any value, but s > 0.0  * ============================================== */{    double t = (x - m) / s;   return (cdfStandard(t));}   double idfNormal(double m, double s, double u)/* =======================================================  * NOTE: m can be any value, but s > 0.0 and 0.0 < u < 1.0  * ======================================================= */{   return (m + s * idfStandard(u));}   double pdfLognormal(double a, double b, double x)/* ===================================================  * NOTE: a can have any value, but b > 0.0 and x > 0.0  * =================================================== */{    double t = (log(x) - a) / b;   return (pdfStandard(t) / (b * x));}   double cdfLognormal(double a, double b, double x)/* ===================================================  * NOTE: a can have any value, but b > 0.0 and x > 0.0  * =================================================== */{    double t = (log(x) - a) / b;   return (cdfStandard(t));}   double idfLognormal(double a, double b, double u)/* =========================================================  * NOTE: a can have any value, but b > 0.0 and 0.0 < u < 1.0  * ========================================================= */{    double t;   t = a + b * idfStandard(u);   return (exp(t));}   double pdfChisquare(long n, double x)/* =====================================  * NOTE: use n >= 1 and x > 0.0  * ===================================== */{    double t, s = n / 2.0;   t = (s - 1.0) * log(x / 2.0) - (x / 2.0) - log(2.0) - LogGamma(s);   return (exp(t));}   double cdfChisquare(long n, double x)/* =====================================  * NOTE: use n >= 1 and x > 0.0  * ===================================== */{   return (InGamma(n / 2.0, x / 2));}   double idfChisquare(long n, double u)/* =====================================  * NOTE: use n >= 1 and 0.0 < u < 1.0  * ===================================== */{    double t, x = n;                         /* initialize to the mean, then */   do {                                     /* use Newton-Raphson iteration */     t = x;     x = t + (u - cdfChisquare(n, t)) / pdfChisquare(n, t);     if (x <= 0.0)       x = 0.5 * t;   } while (fabs(x - t) >= TINY);   return (x);}   double pdfStudent(long n, double x)/* ===================================  * NOTE: use n >= 1 and x > 0.0  * =================================== */{    double s, t;   s = -0.5 * (n + 1) * log(1.0 + ((x * x) / (double) n));   t = -LogBeta(0.5, n / 2.0);   return (exp(s + t) / sqrt((double) n));}   double cdfStudent(long n, double x)/* ===================================  * NOTE: use n >= 1 and x > 0.0  * =================================== */{    double s, t;   t = (x * x) / (n + x * x);   s = InBeta(0.5, n / 2.0, t);   if (x >= 0.0)     return (0.5 * (1.0 + s));   else     return (0.5 * (1.0 - s));}   double idfStudent(long n, double u)/* ===================================  * NOTE: use n >= 1 and 0.0 < u < 1.0  * =================================== */{    double t, x = 0.0;                       /* initialize to the mean, then */   do {                                     /* use Newton-Raphson iteration */     t = x;     x = t + (u - cdfStudent(n, t)) / pdfStudent(n, t);   } while (fabs(x - t) >= TINY);   return (x);}/* =================================================================== * The six functions that follow are a 'special function' mini-library * used to support the evaluation of pdf, cdf and idf functions. * =================================================================== */   static double LogGamma(double a)/* ========================================================================  * LogGamma returns the natural log of the gamma function. * NOTE: use a > 0.0  * * The algorithm used to evaluate the natural log of the gamma function is * based on an approximation by C. Lanczos, SIAM J. Numerical Analysis, B, * vol 1, 1964.  The constants have been selected to yield a relative error * which is less than 2.0e-10 for all positive values of the parameter a.     * ========================================================================  */{    double s[6], sum, temp;   int    i;   s[0] =  76.180091729406 / a;   s[1] = -86.505320327112 / (a + 1.0);   s[2] =  24.014098222230 / (a + 2.0);   s[3] =  -1.231739516140 / (a + 3.0);   s[4] =   0.001208580030 / (a + 4.0);   s[5] =  -0.000005363820 / (a + 5.0);   sum  =   1.000000000178;   for (i = 0; i < 6; i++)      sum += s[i];   temp = (a - 0.5) * log(a + 4.5) - (a + 4.5) + log(SQRT2PI * sum);   return (temp);}   double LogFactorial(long n)/* ================================================================== * LogFactorial(n) returns the natural log of n! * NOTE: use n >= 0 * * The algorithm used to evaluate the natural log of n! is based on a * simple equation which relates the gamma and factorial functions. * ================================================================== */{   return (LogGamma(n + 1));}   static double LogBeta(double a, double b)/* ====================================================================== * LogBeta returns the natural log of the beta function. * NOTE: use a > 0.0 and b > 0.0 * * The algorithm used to evaluate the natural log of the beta function is  * based on a simple equation which relates the gamma and beta functions. * */{    return (LogGamma(a) + LogGamma(b) - LogGamma(a + b));}   double LogChoose(long n, long m)/* ======================================================================== * LogChoose returns the natural log of the binomial coefficient C(n,m). * NOTE: use 0 <= m <= n * * The algorithm used to evaluate the natural log of a binomial coefficient * is based on a simple equation which relates the beta function to a * binomial coefficient. * ======================================================================== */{   if (m > 0)     return (-LogBeta(m, n - m + 1) - log(m));   else     return (0.0);}   static double InGamma(double a, double x)/* ======================================================================== * Evaluates the incomplete gamma function. * NOTE: use a > 0.0 and x >= 0.0 * * The algorithm used to evaluate the incomplete gamma function is based on * Algorithm AS 32, J. Applied Statistics, 1970, by G. P. Bhattacharjee. * See also equations 6.5.29 and 6.5.31 in the Handbook of Mathematical * Functions, Abramowitz and Stegum (editors).  The absolute error is less  * than 1e-10 for all non-negative values of x. * ======================================================================== */{    double t, sum, term, factor, f, g, c[2], p[3], q[3];   long   n;   if (x > 0.0)     factor = exp(-x + a * log(x) - LogGamma(a));   else     factor = 0.0;   if (x < a + 1.0) {                 /* evaluate as an infinite series - */     t    = a;                        /* A & S equation 6.5.29            */     term = 1.0 / a;     sum  = term;     while (term >= TINY * sum) {     /* sum until 'term' is small */       t++;       term *= x / t;       sum  += term;     }      return (factor * sum);   }   else {                             /* evaluate as a continued fraction - */     p[0]  = 0.0;                     /* A & S eqn 6.5.31 with the extended */     q[0]  = 1.0;                     /* pattern 2-a, 2, 3-a, 3, 4-a, 4,... */     p[1]  = 1.0;                     /* - see also A & S sec 3.10, eqn (3) */     q[1]  = x;     f     = p[1] / q[1];     n     = 0;     do {                             /* recursively generate the continued */       g  = f;                        /* fraction 'f' until two consecutive */       n++;                           /* values are small                   */       if ((n % 2) > 0) {         c[0] = ((double) (n + 1) / 2) - a;         c[1] = 1.0;       }       else {         c[0] = (double) n / 2;         c[1] = x;       }       p[2] = c[1] * p[1] + c[0] * p[0];       q[2] = c[1] * q[1] + c[0] * q[0];       if (q[2] != 0.0) {             /* rescale to avoid overflow */         p[0] = p[1] / q[2];         q[0] = q[1] / q[2];         p[1] = p[2] / q[2];         q[1] = 1.0;         f    = p[1];       }     } while ((fabs(f - g) >= TINY) || (q[1] != 1.0));     return (1.0 - factor * f);   }}   static double InBeta(double a, double b, double x)/* =======================================================================  * Evaluates the incomplete beta function. * NOTE: use a > 0.0, b > 0.0 and 0.0 <= x <= 1.0 * * The algorithm used to evaluate the incomplete beta function is based on * equation 26.5.8 in the Handbook of Mathematical Functions, Abramowitz * and Stegum (editors).  The absolute error is less than 1e-10 for all x * between 0 and 1. * ======================================================================= */{    double t, factor, f, g, c, p[3], q[3];   int    swap;   long   n;   if (x > (a + 1.0) / (a + b + 1.0)) { /* to accelerate convergence   */     swap = 1;                          /* complement x and swap a & b */     x    = 1.0 - x;     t    = a;     a    = b;     b    = t;   }   else                                 /* do nothing */     swap = 0;   if (x > 0)     factor = exp(a * log(x) + b * log(1.0 - x) - LogBeta(a,b)) / a;   else     factor = 0.0;   p[0] = 0.0;   q[0] = 1.0;   p[1] = 1.0;   q[1] = 1.0;   f    = p[1] / q[1];   n    = 0;   do {                               /* recursively generate the continued */     g = f;                           /* fraction 'f' until two consecutive */     n++;                             /* values are small                   */     if ((n % 2) > 0) {       t = (double) (n - 1) / 2;       c = -(a + t) * (a + b + t) * x / ((a + n - 1.0) * (a + n));     }     else {       t = (double) n / 2;       c = t * (b - t) * x / ((a + n - 1.0) * (a + n));     }     p[2] = p[1] + c * p[0];     q[2] = q[1] + c * q[0];     if (q[2] != 0.0) {                 /* rescale to avoid overflow */       p[0] = p[1] / q[2];       q[0] = q[1] / q[2];       p[1] = p[2] / q[2];       q[1] = 1.0;       f    = p[1];     }   } while ((fabs(f - g) >= TINY) || (q[1] != 1.0));   if (swap)      return (1.0 - factor * f);   else     return (factor * f);}