If you need numerical stability, I'd use Gary Perlmann's [|stat](http://oldwww.acm.org/perlman/stat/history.html). It's older and somewhat harder to get a copy of, but it's as reliably correct as a piece of software can be...
I've just had a look at the |STAT source, and it computes the variance with
double M = Sum/N; /* mean */
double var = (s2 - M*Sum)/(N-1); /* variance */
where s2 is the sum of squares. In most reasonable situations, this approach will work fine. It just doesn't take as wide a variety of inputs as is easily possible to achieve with normal floating point doubles. In fairness, the |STAT terms and conditions state:
|STAT PROGRAMS HAVE NOT BEEN VALIDATED FOR LARGE DATASETS, HIGHLY VARIABLE DATA, NOR VERY LARGE NUMBERS.
|STAT PROGRAMS HAVE NOT BEEN VALIDATED FOR LARGE DATASETS, HIGHLY VARIABLE DATA, NOR VERY LARGE NUMBERS.