#include <math.h>

double mean( /* Mean */
    int n,
    double x[]
) {
    double s;
    int i;

    for (s = 0.0, i = 0; i < n; i++)
        s += x[i];
    return s/(double)n;
}

double cov( /* Covariance */
    int n,
    double x[],
    double y[]
) {
    double sxy, sx, sy;
    int i;

    for (sxy = sx = sy = 0.0, i = 0; i < n; i++) {
        sxy += x[i]*y[i];
        sx += x[i];
        sy += y[i];
    }
    return sxy/((double)(n-1))-sx*sy/((double)(n-1))/((double)n);
}

double var( /* Variance */
    int n,
    double x[]
) {
    return cov(n,x,x); /* Not exactly optimized... */
}

double corr( /* Correlation */
    int n,
    double x[],
    double y[]
) {
    double vv;

    vv = var(n,x)*var(n,y);
    return vv > 0.0? cov(n,x,y)/sqrt(vv): 0.0; /* Paranoia! Ok? */
}

double sdev( /* Standard Deviation */
    int n,
    double x[]
) {
    double v;

    v = var(n,x);
    return v > 0.0? sqrt(v): 0.0; /* Rounding error paranoia! */
}

void linreg( /* Linear Regression  y = a+b*x+c*e  where e is N(0,1) */
    int n,
    double x[],
    double y[],
    double *a,
    double *b,
    double *c
) {
    double _b, vx, v;

    vx = var(n,x);
    _b = vx > 0.0? cov(n,x,y)/vx: 0.0;

    if (a)
        *a = mean(n,y)-_b*mean(n,x);
    if (b)
        *b = _b;
    if (c)
        if (n > 2) {
            v = var(n,y)-_b*cov(n,x,y); /* The usual paranoia ... */
            *c = v > 0.0 ? sqrt(v*(double)(n-1)/(double)(n-2)): 0.0;
        } else
            *c = vx > 0.0? 0.0: sdev(n,y); /* Otherwise what? */
}

#define _DBL_RAND_MAX 32767.5

double normsim(void) { /* Very Approximate Random Normal Distribution */
    static int f;

    if (!f) {
        srand((int)time((void *)0)); /* You might have to rewrite this! */
        f = 1;
    }
    return ((double)rand()+(double)rand()+(double)rand()
           +(double)rand()+(double)rand()+(double)rand()
           -(double)rand()-(double)rand()-(double)rand()
           -(double)rand()-(double)rand()-(double)rand())
           /_DBL_RAND_MAX;
    /*
        NOTE: You may have to change the
        constant _DBL_RAND_MAX in stat.h
        to make this function work the way it's
        supposed to!
    */
}