BEGIN {
    ITMAX = 100
    EPS7 = 3.0e-7
}

function fatalError(msg, errcode) {
    print msg > "/dev/null"
    exit(errcode+0)
}

function gammaln(xx,
                 x,y,tmp,ser,cof,j) {
    cof[0] = 76.18009172947146
    cof[1] = -86.50532032941677
    cof[2] =  24.01409824083091
    cof[3] =  -1.231739572450155
    cof[4] =   0.1208650973866179e-2
    cof[5] =  -0.5395239384953e-5

    if(xx == 0.0)
        fatalError("gammaln: ERR_INVALID_PARAMETER_0, parameter was: " xx, 2)

    y = x = xx
    tmp = x + 5.5
    tmp -= (x+0.5) * log(tmp)
    ser = 1.000000000190015

    for (j=0;j<=5;j++) ser += cof[j]/++y

    return -tmp + log(2.5066282746310005*ser/x);
}

function fabs(x) {
    return x >= 0.0 ? x : -x
}

function betacf(a, b, x,
                qap,qam,qab,em,tem,d,
                bz,bm,bp,bpp,
                az,am,ap,app,aold,
                m)
{

    bm = az = am = 1.0
    qab=a+b
    qap=a+1.0
    qam=a-1.0
    bz=1.0-qab*x/qap
    for (m=1;m<=ITMAX;m++)
    {
        em=m
        tem=em+em
        d=em*(b-em)*x/((qam+tem)*(a+tem))
        ap=az+d*am
        bp=bz+d*bm
        d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
        app=ap+d*az
        bpp=bp+d*bz
        aold=az
        am=ap/bpp
        bm=bp/bpp
        az=app/bpp
        bz=1.0
        if (fabs(az-aold) < (EPS7 * fabs(az))) return az
    }
    return(paramError = 1)
}

function betai(a, b, x,
               bt) {
    if (x < 0.0 || x > 1.0)
    {
        paramError = 1  # TODO: find out what this is for.
        return(0.0)
    }

    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)
    else
        return (1.0 - bt * betacf(b,a,1.0-x) / b)
}

function probT(t, df,
               bta)
{
    bta = betai(df/2.0, 0.5, 1.0/(1.0 + t*t/df));
    if(t > 0.0) return(1.0 - 0.5 * bta);
    else if(t < 0.0) return(0.5 * bta);
    return(0.5);
}

function critical_value(p, df1, df2, max_val, type,
                        NEWTON_EPSILON, minval,
                        maxval, critval, prob)
{
    NEWTON_EPSILON = 0.000001   # Accuracy of Newton approximation
    minval = 0.0;
    maxval = max_val;
    critval;
    prob = 0.0;

    if (p <= 0.0) return 0.0;
    else if (p >= 1.0) return maxval;

    critval = (df1 + df2) / sqrt(p);    # fair first value
    while ((maxval - minval) > NEWTON_EPSILON)
    {
        prob = probT(critval, df1);

        if (prob < p) minval = critval;
        else maxval = critval;
        critval = (maxval + minval) * 0.5;
    }

    return critval;
}

function criticalX(p, df1,
                   df2, x)
{
    df2 = 0;

    x = critical_value((1.0 - p), df1, df2, 99999.0, type);
    if(x < NEWTON_EPSILON) x = 0.0;
    return(x);
}

function criticalT(p, df) {
    return criticalX(p, df)
}

function confidenceT(alpha, stdev, samplesize) {
    return sqrt(stdev/samplesize) * criticalT(alpha/2, samplesize-1)
}