accurate and stable variance calculation

var.c

#include <stdio.h>
#include <assert.h>
 
typedef double real;
 
static inline real
variance(const real *vec, const int len) {  
    /* iteratove variance calculation */
    register real oldm, newm, olds, news;
    assert( vec && (len > 0));
 
    if (len == 1) {
        return 0.0;
    }
 
    oldm = newm = vec[0];
    olds = 0.0;
 
    for (int i = 1; i < len; ++i) {
  	real t = vec[i];
        newm = oldm + (t - oldm) / (i + 1);
        news = olds + (t - oldm) * (t - newm);
        oldm = newm;
        olds = news;
    }
    
    return news / (len - 1);
}
 
static inline real
textbook_variance(const real *vec, const int len) {
    /* a normal way in text book, un-stable when large scale but min var */
    real ssq = .0, sqs = .0;
 
    assert( vec && (len > 0));
 
    if (len == 1) {
        return 0.0;
    }
 
    for (int i = 0; i < len; ++i) {
        ssq += vec[i] * vec[i];
        sqs += vec[i];
    }
    return ssq / (len - 1) - sqs * sqs / (len * (len - 1));
}
 
static void
tcase() {
    real vec[] = {100000, 20000000000, 3000000000, 4000000,600000,70000000,800000,90000000};
    printf("%lf\n", variance(vec, 8));
    printf("%lf\n", textbook_variance(vec, 8));
}
 
int
main(int argc, char *argv[]) {
    tcase();
    getchar();
    return 0;
}