fourfs.c

Numerical Recipes in C的源代码

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define KBF 128

void fourfs(FILE *file[5], unsigned long nn[], int ndim, int isign)
{
	void fourew(FILE *file[5], int *na, int *nb, int *nc, int *nd);
	unsigned long j,j12,jk,k,kk,n=1,mm,kc=0,kd,ks,kr,nr,ns,nv;
	int cc,na,nb,nc,nd;
	float tempr,tempi,*afa,*afb,*afc;
	double wr,wi,wpr,wpi,wtemp,theta;
	static int mate[5] = {0,2,1,4,3};

	afa=vector(1,KBF);
	afb=vector(1,KBF);
	afc=vector(1,KBF);
	for (j=1;j<=ndim;j++) {
		n *= nn[j];
		if (nn[j] <= 1) nrerror("invalid float or wrong ndim in fourfs");
	}
	nv=1;
	jk=nn[nv];
	mm=n;
	ns=n/KBF;
	nr=ns >> 1;
	kd=KBF >> 1;
	ks=n;
	fourew(file,&na,&nb,&nc,&nd);
	for (;;) {
		theta=isign*3.141592653589793/(n/mm);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		mm >>= 1;
		for (j12=1;j12<=2;j12++) {
			kr=0;
			do {
				cc=fread(&afa[1],sizeof(float),KBF,file[na]);
				if (cc != KBF) nrerror("read error in fourfs");
				cc=fread(&afb[1],sizeof(float),KBF,file[nb]);
				if (cc != KBF) nrerror("read error in fourfs");
				for (j=1;j<=KBF;j+=2) {
					tempr=((float)wr)*afb[j]-((float)wi)*afb[j+1];
					tempi=((float)wi)*afb[j]+((float)wr)*afb[j+1];
					afb[j]=afa[j]-tempr;
					afa[j] += tempr;
					afb[j+1]=afa[j+1]-tempi;
					afa[j+1] += tempi;
				}
				kc += kd;
				if (kc == mm) {
					kc=0;
					wr=(wtemp=wr)*wpr-wi*wpi+wr;
					wi=wi*wpr+wtemp*wpi+wi;
				}
				cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]);
				if (cc != KBF) nrerror("write error in fourfs");
				cc=fwrite(&afb[1],sizeof(float),KBF,file[nd]);
				if (cc != KBF) nrerror("write error in fourfs");
			} while (++kr < nr);
			if (j12 == 1 && ks != n && ks == KBF) {
				na=mate[na];
				nb=na;
			}
			if (nr == 0) break;
		}
		fourew(file,&na,&nb,&nc,&nd);
		jk >>= 1;
		while (jk == 1) {
			mm=n;
			jk=nn[++nv];
		}
		ks >>= 1;
		if (ks > KBF) {
			for (j12=1;j12<=2;j12++) {
				for (kr=1;kr<=ns;kr+=ks/KBF) {
					for (k=1;k<=ks;k+=KBF) {
						cc=fread(&afa[1],sizeof(float),KBF,file[na]);
						if (cc != KBF) nrerror("read error in fourfs");
						cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]);
						if (cc != KBF) nrerror("write error in fourfs");
					}
					nc=mate[nc];
				}
				na=mate[na];
			}
			fourew(file,&na,&nb,&nc,&nd);
		} else if (ks == KBF) nb=na;
		else break;
	}
	j=1;
	for (;;) {
		theta=isign*3.141592653589793/(n/mm);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		mm >>= 1;
		ks=kd;
		kd >>= 1;
		for (j12=1;j12<=2;j12++) {
			for (kr=1;kr<=ns;kr++) {
				cc=fread(&afc[1],sizeof(float),KBF,file[na]);
				if (cc != KBF) nrerror("read error in fourfs");
				kk=1;
				k=ks+1;
				for (;;) {
					tempr=((float)wr)*afc[kk+ks]-((float)wi)*afc[kk+ks+1];
					tempi=((float)wi)*afc[kk+ks]+((float)wr)*afc[kk+ks+1];
					afa[j]=afc[kk]+tempr;
					afb[j]=afc[kk]-tempr;
					afa[++j]=afc[++kk]+tempi;
					afb[j++]=afc[kk++]-tempi;
					if (kk < k) continue;
					kc += kd;
					if (kc == mm) {
						kc=0;
						wr=(wtemp=wr)*wpr-wi*wpi+wr;
						wi=wi*wpr+wtemp*wpi+wi;
					}
					kk += ks;
					if (kk > KBF) break;
					else k=kk+ks;
				}
				if (j > KBF) {
					cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]);
					if (cc != KBF) nrerror("write error in fourfs");
					cc=fwrite(&afb[1],sizeof(float),KBF,file[nd]);
					if (cc != KBF) nrerror("write error in fourfs");
					j=1;
				}
			}
			na=mate[na];
		}
		fourew(file,&na,&nb,&nc,&nd);
		jk >>= 1;
		if (jk > 1) continue;
		mm=n;
		do {
			if (nv < ndim) jk=nn[++nv];
			else {
				free_vector(afc,1,KBF);
				free_vector(afb,1,KBF);
				free_vector(afa,1,KBF);
				return;
			}
		} while (jk == 1);
	}
}
#undef KBF
#undef NRANSI