/********************************************************************
 * Special functions. Refer to Chapter 6 of Numerical Recipe in C
 * http://www.library.cornell.edu/nr/bookcpdf.html
 ********************************************************************/

ITMAX = 100;        // maximum allowed number of iterations
EPS = 3.0e-9;       // relative accuracy
FPMIN = 1.0e-30;     // number near the smallest representable


/************************************************************
 * Returns the value ln(gamma(x)) for x > 0.
 ************************************************************/
function gammaln(x)
{
    cof = [76.18009172947146, -86.50532032941677, 24.01409824083091,
           -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];
    y = x;
    tmp = x + 5.5;
    tmp -= (x + 0.5) * log(tmp);
    ser=1.000000000190015;
    for (j = 0; j < 6; j++) {
        y++;
        ser += cof[j] / y;
    }
    return -tmp + log(2.5066282746310005 * ser / x);
}


/************************************************************
 * Returns the factorial of n, i.e., n! = n*n(-1)*...*1
 ************************************************************/
function factorial(n)
{
    if (n <= 0) error("n must be no less than 0");
    if (n > 170) error("n > 170 overflow");
    a = 1.0;
    j = 1;
    while (j < n) {
        a2 = a * j;
        a = a2;
        j++;
    }
    return a2;
}


/************************************************************
 * Returns the binomial (n, k) coefficient.
 ************************************************************/
function binomial(n, k)
{
    return floor(0.5 + factorial(n)/factorial(k)/ factorial(n-k));
}


/************************************************************
 * Returns the value of the beta function B(z,w).
 ************************************************************/
function beta(z, w)
{
    return exp(gammaln(z) + gammaln(w) - gammaln(z+w));
}

/************************************************************
 * Returns the incomplete gamma function P(a, x).
 ************************************************************/
function gammap(a, x)
{
    if (x < 0.0 || a <= 0.0) error("invalid arguments");
    
    if (x < a+1.0) return gser(a, x);

    return 1.0 - gcf(a, x);
}


/************************************************************
 * Returns the incomplete gamma function Q(a, x).
 ************************************************************/
function gammaq(a, x)
{
    if (x < 0.0 || a <= 0.0) error("invalid arguments");
    
    if (x < a+1.0) return 1.0 - gser(a, x);

    return gcf(a, x);
}


/************************************************************
 * Returns the incomplete gamma function P(a, x) evaluated
 * by its series representation as gamser.
 ************************************************************/
function gser(a, x)
{
    if (x < 0.0) error("x less than 0 in routine gser");
    
    if (x == 0.0) return 0.0;
    
    ap = a;
    sum = 1.0 / a;
    del = sum;
    
    for (n = 1; n <= ITMAX; n++) {
        ap++;
        del *= x / ap;
        sum += del;
        if (abs(del) < abs(sum)*EPS) return sum*exp(-x+a*log(x)-gammaln(a));
    }
    
    error("a too large, ITMAX too small in routine gser");
}

/************************************************************
 * Returns the incomplete gamma function Q(a, x) evaluated
 * by its continued fraction representation as gammcf.
 ************************************************************/
function gcf(a, x)
{
    b = x + 1.0 - a;
    c = 1.0 / FPMIN;
    d = 1.0 / b;
    h = d;
    
    for (i = 1; i <= ITMAX; i++) {
        an = -i * (i - a);
        b += 2.0;
        d = an * d + b;
        if (abs(d) < FPMIN) d = FPMIN;
        c = b + an / c;
        if (abs(c) < FPMIN) c = FPMIN;
        d = 1.0 / d;
        del = d * c;
        h *= del;
        if (abs(del - 1.0) < EPS) return exp(-x+a*log(x)-gammaln(a))*h;
    }
    
    error("a too large, ITMAX too small in gcf");
}

/************************************************************
 * Returns the error function erf(x).
 ************************************************************/
function erf(x)
{
    if (x < 0.0) return -gammap(0.5, x*x);
    return gammap(0.5, x*x);
}


/************************************************************
 * Returns the complementary error function erfc(x).
 ************************************************************/
function erfc(x)
{
    if (x < 0.0) return 1.0 + gammap(0.5, x*x);
    return gammaq(0.5, x*x);
}


/************************************************************
 * Returns the complementary error function erfc(x)
 * with fractional error everywhere less than 1.2 x 10e-7.
 ************************************************************/
function erfcc(x)
{
    z = abs(x);
    t = 1.0 / (1.0 + 0.5 * z);
    ans = t * exp(-z*z - 1.26551223 + 
          t * (1.00002368 + t * (0.37409196 + 
          t * (0.09678418 + t * (-0.18628806 + 
          t * (0.27886807 + t * (-1.13520398 +
          t * (1.48851587 + t * (-0.82215223 + t * 0.17087277)))))))));
   if (x >= 0.0) return ans;
   return 2.0 - ans;
}


/************************************************************
 * Returns the incomplete beta function Ix(a, b).
 ************************************************************/
function betai(a, b, x)
{
    if (x < 0.0 || x > 1.0) error("Bad x in routine betai");
    
    if (x == 0.0 || x == 1.0)
        bt = 0.0;
    else
        bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b)+a*log(x)+b*log(1.0-x));
        
    if (x < (a+1.0)/(a+b+2.0)) return bt*betacf(a,b,x)/a;
    
    return 1.0 - bt*betacf(b,a,1.0-x) / b;
}

/************************************************************
 * Evaluates continued fraction for incomplete beta function
 ************************************************************/
function betacf(a, b, x)
{
    qab = a + b;
    qap = a + 1.0;
    qam = a - 1.0;
    c = 1.0;
    d = 1.0 - qab * x / qap;
    if (abs(d) < FPMIN) d = FPMIN;
    d = 1.0 / d;
    h = d;
    for (m = 1; m <= ITMAX; m++) {
        m2 = 2 * m;
        aa = m*(b - m)*x / ((qam + m2)*(a + m2));
        d = 1.0 + aa * d;
        if (abs(d) < FPMIN) d = FPMIN;
        c = 1.0 + aa / c;
        if (abs(c) < FPMIN) c = FPMIN;
        d = 1.0 / d;
        h *= d * c;
        aa = -(a + m)*(qab + m)*x / ((a + m2)*(qap + m2));
        d = 1.0 + aa * d;
        if (abs(d) < FPMIN) d = FPMIN;
        c = 1.0 + aa/c;
        if (abs(c) < FPMIN) c = FPMIN;
        d = 1.0 / d;
        del = d * c;
        h *= del;
        if (abs(del-1.0) < EPS) return h;
    }
    error("a or b too big, or ITMAX too small in betacf");
}