12.6 External Storage or Memory-Local FFTs 533 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[i]; if (nn[j]<1)nrerror("invalid float or wrong ndim in fourfs"); nv=1; http://www.nr. readable files Permission is jk=nn[nv]; mm=n: ns=n/KBF nr=ns >1: granted fori kd=KBF >>1; ks=n; .com or call 1-800-872- (including this one) fourew (file,&na,&nb,&nc,&nd); internet The first phase of the transform starts here. for(;;)[ Start of the computing pass. theta=1s1gn*3.141592653589793/(n/mm); wtemp=sin(0.5*theta); wpr =-2.0*wtemp*wtemp; wpi=sin(theta); r=1.0 -7423(North America w1=0.0; to any server computer,is strictly prohibited. t users to make one paper copy for their 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN THE mm>>=1; for(j12=1;j12<=2;j12++) kr=0; only),or 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[i+1]: tempi=((float)wi)*afb[j]+((float)wr)*afb[j+1]; afb[j]-afa[j]-tempr; afa[j]+tempr; afb[i+1]safa[i+1]-tempi; afa[j+1]+tempi; kc +kd; rsend email to directcustserv@cambridge.org ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) if (kc ==mm){ kc=0; personal use.Further reproduction,or 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") (outside North America) 之 Software. cc=fwrite(&afb[1],sizeof(float),KBF,file [nd]); if (cc !KBF)nrerror("write error in fourfs"); while (++kr nr); if (i12==1 ks !n&&ks =KBF){ ying of machine na=mate [na]; nb=na; if (nr ==0)break; fourew(file,&na,&nb,&nc,&nd); Start of the permutation pass jk>>=1; while (jk ==1){ mm=ni12.6 External Storage or Memory-Local FFTs 533 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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); The first phase of the transform starts here. for (;;) { Start of the computing pass. 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); Start of the permutation pass. jk >>= 1; while (jk == 1) { mm=n;