Category : Files from Magazines
Archive   : CUJ9305.ZIP
Filename : 1105068A

 
Output of file : 1105068A contained in archive : CUJ9305.ZIP

#include
#define SWAP(a,b) tempr=a;a=b;b=tempr
void four1(float *data, int *nn, int *isign)
{ /* altered for consistency
* with original FORTRAN */
/* Press, Flannery, Teukolsky, Vettering "Numerical
* Recipes in C" tuned up ; Code works only when *nn is
* a power of 2 */
int n, mmax, m, j, i;
double wtemp, wr, wpr, wpi, wi, theta, wpin;
float tempr, tempi, datar, datai;
n = *nn * 2;
j = 0;
for (i = 0; i < n; i += 2) {
if (j > i) { /* could use j>i+1 to help
* compiler analysis */
SWAP(data[j], data[i]);
SWAP(data[j + 1], data[i + 1]);
}
m = *nn;
while (m >= 2 && j >= m) {
j -= m;
m >>= 1;
}
j += m;
}
theta = 3.141592653589795 * .5;
if (*isign < 0)
theta = -theta;
wpin = 0; /* sin(+-PI) */
for (mmax = 2; n > mmax; mmax *= 2) {
wpi = wpin;
wpin = sin(theta);
wpr = 1 - wpin * wpin - wpin * wpin; /* cos(theta*2) */
theta *= .5;
wr = 1;
wi = 0;
for (m = 0; m < mmax; m += 2) {
for (i = m; i < n; i += mmax * 2) {
j = i + mmax;
/* mixed precision not significantly more
* accurate here; if removing float casts,
* tempr and tempi should be double */
tempr = (float) wr *data[j] - (float) wi *data[j + 1];
tempi = (float) wr *data[j + 1] + (float) wi *data[j];
/* don't expect compiler to analyze j >
* i+1; original code stored data[j..]
* first and avoided using data[ri] */
data[i] = (datar = data[i]) + tempr;
data[i + 1] = (datai = data[i + 1]) + tempi;
data[j] = datar - tempr;
data[j + 1] = datai - tempi;
}
wr = (wtemp = wr) * wpr - wi * wpi;
wi = wtemp * wpi + wi * wpr;
}
}
}