eulsum.c

Numerical Recipes 是国际公认的高水平的、关于数值计算的书

#include <math.h>

void eulsum(sum,term,jterm,wksp)
float *sum,term,wksp[];
int jterm;
{
	int j;
	static int nterm;
	float tmp,dum;

	if (jterm == 1) {
		nterm=1;
		*sum=0.5*(wksp[1]=term);
	} else {
		tmp=wksp[1];
		wksp[1]=term;
		for (j=1;j<=nterm-1;j++) {
			dum=wksp[j+1];
			wksp[j+1]=0.5*(wksp[j]+tmp);
			tmp=dum;
		}
		wksp[nterm+1]=0.5*(wksp[nterm]+tmp);
		if (fabs(wksp[nterm+1]) <= fabs(wksp[nterm]))
			*sum += (0.5*wksp[++nterm]);
		else
			*sum += wksp[nterm+1];
	}
}