#include <math.h>
#include <stdio.h>

typedef double gnm_float;
#define	GNM_FLOAT_PRINTF	"f"

/* Arithmetic sum.  */
/* We are using the double pass algorithm 2.3.1 described in */
/* Thisted: Elements of Statistical Computing                */
/* Note: mathematically sum_2 should be zero, except that it */
/* contains some of the numerical errors from the first pass */
static int
range_sum (const gnm_float *xs, int n, gnm_float *res)
{
	gnm_float sum_1 = 0.;
	gnm_float sum_2 = 0.;
	gnm_float xbar = 0.;
	int i;

	if (n == 0) {
		*res = 0.0;
		return 0;
	}

	for (i = 0; i < n; i++)
		sum_1 += xs[i];
	xbar = sum_1 / n;
	for (i = 0; i < n; i++)
		sum_2 += (xs[i] - xbar);
	printf ("sum_2 %.20" GNM_FLOAT_PRINTF "\n", sum_2);
	*res = sum_1 + sum_2;

	return 0;
}

main()
{
    static gnm_float data[] =
    {
	1., 1., 1., 0., 0.,
	0., 0., 0., 0., 0.,
	0., 0., 0., 0., 0.,
	0., 0., 0., 0., 0.,
	0., 0., 0., 0., 0.
    };
    gnm_float res;

    range_sum (data, sizeof (data)/ sizeof (gnm_float), &res);
}
