#include #include #include #include "chunk_array.h" /*fast fourier transform */ /* THE COOLEY-TUKEY FAST FOURIER TRANSFORM */ /* EVALUATES COMPLEX FOURIER SERIES FOR COMPLEX OR REAL FUNCTIONS. */ /* THAT IS, IT COMPUTES */ /* FTRAN(J1,J2,...)=SUM(DATA(I1,I2,...)*W1**(I1-1)*(J1-1) */ /* *W2**(I2-1)*(J2-1)*...), */ /* WHERE W1=EXP(-2*PI*SQRT(-1)/NN(1)), W2=EXP(-2*PI*SQRT(-1)/NN(2)), */ /* ETC. AND I1 AND J1 RUN FROM 1 TO NN(1), I2 AND J2 RUN FROM 1 TO */ /* NN(2), ETC. THERE IS NO LIMIT ON THE DIMENSIONALITY (NUMBER OF */ /* SUBSCRIPTS) OF THE ARRAY OF DATA. THE PROGRAM WILL PERFORM */ /* A THREE-DIMENSIONAL FOURIER TRANSFORM AS EASILY AS A ONE-DIMEN- */ /* SIONAL ONE, THO IN A PROPORTIONATELY GREATER TIME. AN INVERSE */ /* TRANSFORM CAN BE PERFORMED, IN WHICH THE SIGN IN THE EXPONENTIALS */ /* IS +, INSTEAD OF -. IF AN INVERSE TRANSFORM IS PERFORMED UPON */ /* AN ARRAY OF TRANSFORMED DATA, THE ORIGINAL DATA WILL REAPPEAR, */ /* MULTIPLIED BY NN(1)*NN(2)*... THE ARRAY OF INPUT DATA MAY BE */ /* REAL OR COMPLEX, AT THE PROGRAMMERS OPTION, WITH A SAVING OF */ /* ABOUT THIRTY PER CENT IN RUNNING TIME FOR REAL OVER COMPLEX. */ /* (FOR FASTEST TRANSFORM OF REAL DATA, NN(1) SHOULD BE EVEN.) */ /* THE TRANSFORM VALUES ARE ALWAYS COMPLEX, AND ARE RETURNED IN THE */ /* ORIGINAL ARRAY OF DATA, REPLACING THE INPUT DATA. THE LENGTH */ /* OF EACH DIMENSION OF THE DATA ARRAY MAY BE ANY INTEGER. THE */ /* PROGRAM RUNS FASTER ON COMPOSITE INTEGERS THAN ON PRIMES, AND IS */ /* PARTICULARLY FAST ON NUMBERS RICH IN FACTORS OF TWO. */ /* TIMING IS IN FACT GIVEN BY THE FOLLOWING FORMULA. LET NTOT BE THE */ /* TOTAL NUMBER OF POINTS (REAL OR COMPLEX) IN THE DATA ARRAY, THAT */ /* IS, NTOT=NN(1)*NN(2)*... DECOMPOSE NTOT INTO ITS PRIME FACTORS, */ /* SUCH AS 2**K2 * 3**K3 * 5**K5 * ... LET SUM2 BE THE SUM OF ALL */ /* THE FACTORS OF TWO IN NTOT, THAT IS, SUM2 = 2*K2. LET SUMF BE */ /* THE SUM OF ALL OTHER FACTORS OF NTOT, THAT IS, SUMF = 3*K3+5*K5+.. */ /* THE TIME TAKEN BY A MULTIDIMENSIONAL TRANSFORM ON THESE NTOT DATA */ /* IS T = T0 + T1*NTOT + T2*NTOT*SUM2 + T3*NTOT*SUMF. FOR THE PAR- */ /* TICULAR IMPLEMENTATION FORTRAN 32 ON THE CDC 3300 (FLOATING POINT */ /* ADD TIME = SIX MICROSECONDS), */ /* T = 3000 + 600*NTOT + 50*NTOT*SUM2 + 175*NTOT*SUMF MICROSECONDS */ /* ON COMPLEX DATA. */ /* IMPLEMENTATION OF THE DEFINITION BY SUMMATION WILL RUN IN A TIME */ /* PROPORTIONAL TO NTOT**2. FOR HIGHLY COMPOSITE NTOT, THE SAVINGS */ /* OFFERED BY COOLEY-TUKEY CAN BE DRAMATIC. A MATRIX 100 BY 100 WILL */ /* BE TRANSFORMED IN TIME PROPORTIONAL TO 10000*(2+2+2+2+5+5+5+5) = */ /* 280,000 (ASSUMING T2 AND T3 TO BE ROUGHLY COMPARABLE) VERSUS */ /* 10000**2 = 100,000,000 FOR THE STRAIGHTFORWARD TECHNIQUE. */ /* THE COOLEY-TUKEY ALGORITHM PLACES TWO RESTRICTIONS UPON THE */ /* NATURE OF THE DATA BEYOND THE USUAL RESTRICTION THAT */ /* THE DATA FROM ONE CYCLE OF A PERIODIC FUNCTION. THEY ARE-- */ /* 1. THE NUMBER OF INPUT DATA AND THE NUMBER OF TRANSFORM VALUES */ /* MUST BE THE SAME. */ /* 2. CONSIDERING THE DATA TO BE IN THE TIME DOMAIN, */ /* THEY MUST BE EQUI-SPACED AT INTERVALS OF DT. FURTHER, THE TRANS- */ /* FORM VALUES, CONSIDERED TO BE IN FREQUENCY SPACE, WILL BE EQUI- */ /* SPACED FROM 0 TO 2*PI*(NN(I)-1)/(NN(I)*DT) AT INTERVALS OF */ /* 2*PI/(NN(I)*DT) FOR EACH DIMENSION OF LENGTH NN(I). OF COURSE, */ /* DT NEED NOT BE THE SAME FOR EVERY DIMENSION. */ /* THE CALLING SEQUENCE IS-- */ /* CALL FOURT(DATAR,DATAI,NN,NDIM,IFRWD,ICPLX,WORKR,WORKI) */ /* DATAR AND DATAI ARE THE ARRAYS USED TO HOLD THE REAL AND IMAGINARY */ /* PARTS OF THE INPUT DATA ON INPUT AND THE TRANSFORM VALUES ON */ /* OUTPUT. THEY ARE FLOATING POINT ARRAYS, MULTIDIMENSIONAL WITH */ /* IDENTICAL DIMENSIONALITY AND EXTENT. THE EXTENT OF EACH DIMENSION */ /* IS GIVEN IN THE INTEGER ARRAY NN, OF LENGTH NDIM. THAT IS, */ /* NDIM IS THE DIMENSIONALITY OF THE ARRAYS DATAR AND DATAI. */ /* IFRWD IS AN INTEGER USED TO INDICATE THE DIRECTION OF THE FOURIER */ /* TRANSFORM. IT IS NON-ZERO TO INDICATE A FORWARD TRANSFORM */ /* (EXPONENTIAL SIGN IS -) AND ZERO TO INDICATE AN INVERSE TRANSFORM */ /* (SIGN IS +). ICPLX IS AN INTEGER TO INDICATE WHETHER THE DATA */ /* ARE REAL OR COMPLEX. IT IS NON-ZERO FOR COMPLEX, ZERO FOR REAL. */ /* IF IT IS ZERO (REAL) THE CONTENTS OF ARRAY DATAI WILL BE ASSUMED */ /* TO BE ZERO, AND NEED NOT BE EXPLICITLY SET TO ZERO. AS EXPLAINED */ /* ABOVE, THE TRANSFORM RESULTS ARE ALWAYS COMPLEX AND ARE STORED */ /* IN DATAR AND DATAI ON RETURN. WORKR AND WORKI ARE ARRAYS USED */ /* FOR WORKING STORAGE. THEY ARE NOT NECESSARY IF ALL THE DIMENSIONS */ /* OF THE DATA ARE POWERS OF TWO. IN THIS CASE, THE ARRAYS MAY BE */ /* REPLACED BY THE NUMBER 0 IN THE CALLING SEQUENCE. THUS, USE OF */ /* POWERS OF TWO CAN FREE A GOOD DEAL OF STORAGE. IF ANY DIMENSION */ /* IS NOT A POWER OF TWO, THESE ARRAYS MUST BE SUPPLIED. THEY ARE */ /* FLOATING POINT, ONE DIMENSIONAL OF LENGTH EQUAL TO THE LARGEST */ /* ARRAY DIMENSION, THAT IS, TO THE LARGEST VALUE OF NN(I). */ /* WORKR AND WORKI, IF SUPPLIED, MUST NOT BE THE SAME ARRAYS AS DATAR */ /* OR DATAI. ALL SUBSCRIPTS OF ALL ARRAYS BEGIN AT 1. */ /* THERE ARE NO ERROR MESSAGES OR ERROR HALTS IN THIS PROGRAM. THE */ /* PROGRAM RETURNS IMMEDIATELY IF NDIM OR ANY NN(I) IS LESS THAN ONE. */ /* PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER */ /* 10-06-2000, MLR */ void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) { int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf; double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559; double valueri, valueri1, valueri3, valueii3, valuerj3, valueij3, valuerj, valueij, valueii, valuerimin, valueiimin, valuei1, valuer1; ntot = 1; for (idim = 0; idim < ndim; idim++) { ntot *= nn[idim]; } chunk_array_read(datar); chunk_array_read(datai); /*main loop for each dimension*/ np1 = 1; for (idim = 1; idim <= ndim; idim++) { n = nn[idim - 1]; np2 = np1 * n; if (n < 1) { goto L920; } else if (n == 1) { goto L900; } /*is n a power of 2 and if not, what are its factors*/ m = n; ntwo = np1; iff = 1; idiv = 2; L10: iquot = m / idiv; irem = m - idiv * iquot; if (iquot < idiv) goto L50; if (irem == 0) { ntwo *= 2; ifact[iff] = idiv; iff++; m = iquot; goto L10; } idiv = 3; inon2 = iff; L30: iquot = m / idiv; irem = m - idiv * iquot; if (iquot < idiv) goto L60; if (irem == 0) { ifact[iff] = idiv; iff++; m = iquot; goto L30; } idiv += 2; goto L30; L50: inon2 = iff; if (irem != 0) goto L60; ntwo *= 2; goto L70; L60: ifact[iff] = m; L70: non2p = np2 / ntwo; /*SEPARATE FOUR CASES-- 1. COMPLEX TRANSFORM 2. REAL TRANSFORM FOR THE 2ND, 3RD, ETC. DIMENSION. METHOD: TRANSFORM HALF THE DATA, SUPPLYING THE OTHER HALF BY CONJUGATE SYMMETRY. 3. REAL TRANSFORM FOR THE 1ST DIMENSION, N ODD. METHOD: SET THE IMAGINARY PARTS TO ZERO. 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD: TRANSFORM A COMPLEX ARRAY OF LENGTH N/2 WHOSE REAL PARTS ARE THE EVEN NUMBERED REAL VALUES AND WHOSE IMAGINARY PARTS ARE THE ODD-NUMBERED REAL VALUES. UNSCRAMBLE AND SUPPLY THE SECOND HALF BY CONJUGATE SYMMETRY. */ icase = 1; ifmin = 1; if (icplx != 0) goto L100; icase = 2; if (idim > 1) goto L100; icase = 3; if (ntwo <= np1) goto L100; icase = 4; ifmin = 2; ntwo /= 2; n /= 2; np2 /= 2; ntot /= 2; i = 1; for (j = 1; j <= ntot; j++) { chunk_array_get(datar, i, &valueri); ////printf("[1] datar[%d] = %f\n", i, valueri); chunk_array_get(datar, i+1, &valueri1); ////printf("[2] datar[%d] = %f\n", i+1, valueri1); ////printf("[48] Saving in datar the value %f in pos %d\n", valueri, j); ////printf("[49] Saving in datai the value %f in pos %d\n", valueri1, j); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueri1); i += 2; } /*shuffle data by bit reversal, since n = 2^k. As the shuffling can be done by simple interchange, no working array is needed*/ L100: if (non2p > 1) goto L200; np2hf = np2 / 2; j = 1; for (i2 = 1; i2 <= np2; i2 += np1) { if (j >= i2) goto L130; i1max = i2 + np1 - 1; for (i1 = i2; i1 <= i1max; i1++) { for (i3 = i1; i3 <= ntot; i3 += np2) { j3 = j + i3 - i2; chunk_array_get(datar, i3, &valueri3); chunk_array_get(datai, i3, &valueii3); chunk_array_get(datar, j3, &valuerj3); chunk_array_get(datai, j3, &valueij3); //printf("[3] datar[%d] = %f\n", i3, valueri3); //printf("[4] datai[%d] = %f\n", i3, valueii3); //printf("[5] datar[%d] = %f\n", j3, valuerj3); //printf("[6] datai[%d] = %f\n", j3, valueij3); tempr = valueri3; tempi = valueii3; //printf("[50] Saving in datar the value %f in pos %d\n", valuerj3, i3); //printf("[51] Saving in datai the value %f in pos %d\n", valueij3, i3); //printf("[52] Saving in datar the value %f in pos %d\n", tempr, j3); //printf("[53] Saving in datai the value %f in pos %d\n", tempi, j3); chunk_array_save(datar, i3, valuerj3); chunk_array_save(datai, i3, valueij3); chunk_array_save(datar, j3, tempr); chunk_array_save(datai, j3, tempi); } } L130: m = np2hf; L140: if (j <= m) { j += m; } else { j -= m; m /= 2; if (m >= np1) goto L140; } } goto L300; /*shuffle data by digit reversal for general n*/ L200: for (i1 = 1; i1 <= np1; i1++) { for (i3 = i1; i3 <= ntot; i3 += np2) { j = i3; for (i = 1; i <= n; i++) { chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[7] datar[%d] = %f\n", j, valuerj); //printf("[8] datai[%d] = %f\n", j, valueij); if (icase != 3) { workr[i] = valuerj; worki[i] = valueij; } else { workr[i] = valuerj; worki[i] = 0.; } ifp2 = np2; iff = ifmin; L250: ifp1 = ifp2 / ifact[iff]; j += ifp1; if (j >= i3 + ifp2) { j -= ifp2; ifp2 = ifp1; iff += 1; if (ifp2 > np1) goto L250; } } i2max = i3 + np2 - np1; i = 1; for (i2 = i3; i2 <= i2max; i2 += np1) { //printf("[54] Saving in datar the value %f in pos %d\n", workr[i], i2); //printf("[55] Saving in datai the value %f in pos %d\n", worki[i], i2); chunk_array_save(datar, i2, workr[i]); chunk_array_save(datai, i2, worki[i]); i++; } } } /*special case-- W=1*/ L300: i1rng = np1; if (icase == 2) i1rng = np0 * (1 + nprev / 2); if (ntwo <= np1) goto L600; for (i1 = 1; i1 <= i1rng; i1++) { imin = np1 + i1; istep = 2 * np1; goto L330; L310: j = i1; for (i = imin; i <= ntot; i += istep) { chunk_array_get(datar, i, &tempr); chunk_array_get(datai, i, &tempi); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[9] tempr = %f\n", i, tempr); //printf("[10] tempi = %f\n", i, tempi); //printf("[11] datar[%d] = %f\n", j, valuerj); //printf("[12] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); //printf("[56] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); //printf("[57] Saving in datai the value %f in pos %d\n", valueij - tempi, i); //printf("[58] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); //printf("[59] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } imin = 2 * imin - i1; istep *= 2; L330: if (istep <= ntwo) goto L310; /*special case-- W = -sqrt(-1)*/ imin = 3 * np1 + i1; istep = 4 * np1; goto L420; L400: j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { if (ifrwd != 0) { chunk_array_get(datai, i, &tempr); chunk_array_get(datar, i, &tempi); //printf("[13] datai[%d] = %f\n", i, tempr); //printf("[14] datar[%d] = %f\n", i, tempi); tempi = -tempi; } else { chunk_array_get(datai, i, &tempr); //printf("[15] datai[%d] = %f\n", i, tempr); tempr = -tempr; chunk_array_get(datar, i, &tempi); //printf("[16] datar[%d] = %f\n", i, tempi); } chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[17] datar[%d] = %f\n", j, valuerj); //printf("[18] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); //printf("[60] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); //printf("[61] Saving in datai the value %f in pos %d\n", valueij - tempi, i); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[19] datar[%d] = %f\n", j, valuerj); //printf("[20] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); //printf("[62] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); //printf("[63] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } imin = 2 * imin - i1; istep *= 2; L420: if (istep <= ntwo) goto L400; } /*main loop for factors of 2. W=EXP(-2*PI*SQRT(-1)*m/mmax) */ theta = -TWOPI / 8.; wstpr = 0.; wstpi = -1.; if (ifrwd == 0) { theta = -theta; wstpi = 1.; } mmax = 8 * np1; goto L540; L500: wminr = cos(theta); wmini = sin(theta); wr = wminr; wi = wmini; mmin = mmax / 2 + np1; mstep = np1 * 2; for (m = mmin; m <= mmax; m += mstep) { for (i1 = 1; i1 <= i1rng; i1++) { istep = mmax; imin = m + i1; L510: j = imin - istep / 2; for (i = imin; i <= ntot; i += istep) { chunk_array_get(datar, i, &valueri); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, i, &valueii); chunk_array_get(datai, j, &valueij); //printf("[21] datar[%d] = %f\n", i, valueri); //printf("[22] datar[%d] = %f\n", j, valuerj); //printf("[23] datai[%d] = %f\n", i, valueii); //printf("[24] datai[%d] = %f\n", j, valueij); tempr = valueri * wr - valueii * wi; tempi = valueri * wi + valueii * wr; chunk_array_save(datar, i, valuerj - tempr); chunk_array_save(datai, i, valueij - tempi); chunk_array_save(datar, j, valuerj + tempr); chunk_array_save(datai, j, valueij + tempi); //printf("[64] Saving in datar the value %f in pos %d\n", valuerj - tempr, i); //printf("[65] Saving in datai the value %f in pos %d\n", valueij - tempi, i); //printf("[66] Saving in datar the value %f in pos %d\n", valuerj + tempr, j); //printf("[67] Saving in datai the value %f in pos %d\n", valueij + tempi, j); j += istep; } imin = 2 * imin - i1; istep *= 2; if (istep <= ntwo) goto L510; } wtemp = wr * wstpi; wr = wr * wstpr - wi * wstpi; wi = wi * wstpr + wtemp; } wstpr = wminr; wstpi = wmini; theta /= 2.; mmax += mmax; L540: if (mmax <= ntwo) goto L500; /*main loop for factors not equal to 2-- W=EXP(-2*PI*SQRT(-1)*(j2-i3)/ifp2)*/ L600: if (non2p <= 1) goto L700; ifp1 = ntwo; iff = inon2; L610: ifp2 = ifact[iff] * ifp1; theta = -TWOPI / (double)ifact[iff]; if (ifrwd == 0) theta = -theta; thetm = theta / (double)(ifp1 / np1); wstpr = cos(theta); wstpi = sin(theta); wmstr = cos(thetm); wmsti = sin(thetm); wminr = 1.; wmini = 0.; for (j1 = 1; j1 <= ifp1; j1 += np1) { i1max = j1 + i1rng - 1; for (i1 = j1; i1 <= i1max; i1++) { for (i3 = i1; i3 <= ntot; i3 += np2) { i = 1; wr = wminr; wi = wmini; j2max = i3 + ifp2 - ifp1; for (j2 = i3; j2 <= j2max; j2 += ifp1) { twowr = 2. * wr; jmin = i3; j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { j = jmin + ifp2 - ifp1; chunk_array_get(datar, j, &sr); chunk_array_get(datai, j, &si); //printf("[25] datar[%d] = %f\n", j, sr); //printf("[26] datai[%d] = %f\n", j, si); oldsr = 0.; oldsi = 0.; j -= ifp1; L620: stmpr = sr; stmpi = si; chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[27] datar[%d] = %f\n", j, valuerj); //printf("[28] datai[%d] = %f\n", j, valueij); sr = twowr * sr - oldsr + valuerj; si = twowr * si - oldsi + valueij; oldsr = stmpr; oldsi = stmpi; j -= ifp1; if (j > jmin) goto L620; chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); workr[i] = wr * sr - wi * si - oldsr + valuerj; worki[i] = wi * sr + wr * si - oldsi + valueij; //printf("[85] wr = %f, sr = %f, wi = %f, si = %f, oldsr = %f, datar[j] = %f\n", wr, sr, wi, si, oldsr, valuerj); //printf("[86] wi = %f, sr = %f, wr = %f, si = %f, oldsi = %f, datai[j] = %f\n", wi, sr, wr, si, oldsi, valueij); //printf("[83] Saving in workr the value %f in pos %d\n", wr * sr - wi * si - oldsr + valuerj, i); //printf("[84] Saving in worki the value %f in pos %d\n", wi * sr + wr * si - oldsi + valueij, i); jmin += ifp2; i++; } wtemp = wr * wstpi; wr = wr * wstpr - wi * wstpi; wi = wi * wstpr + wtemp; } i = 1; for (j2 = i3; j2 <= j2max; j2 += ifp1) { j3max = j2 + np2 - ifp2; for (j3 = j2; j3 <= j3max; j3 += ifp2) { //printf("[68] Saving in datar the value %f in pos %d, i = %d\n", workr[i], j3, i); //printf("[69] Saving in datai the value %f in pos %d, i = %d\n", worki[i], j3, i); chunk_array_save(datar, j3, workr[i]); chunk_array_save(datai, j3, worki[i]); i++; } } } } wtemp = wminr * wmsti; wminr = wminr * wmstr - wmini * wmsti; wmini = wmini * wmstr + wtemp; } iff++; ifp1 = ifp2; if (ifp1 < np2) goto L610; /*complete a real transform in the 1st dimension, n even, by conjugate symmetries*/ L700: switch (icase) { case 1: goto L900; break; case 2: goto L800; break; case 3: goto L900; break; } nhalf = n; n += n; theta = -TWOPI / (double)n; if (ifrwd == 0) theta = -theta; wstpr = cos(theta); wstpi = sin(theta); wr = wstpr; wi = wstpi; imin = 2; jmin = nhalf; goto L725; L710: j = jmin; for (i = imin; i <= ntot; i += np2) { chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[29] datar[%d] = %f\n", i, valueri); //printf("[30] datai[%d] = %f\n", i, valueii); //printf("[31] datar[%d] = %f\n", j, valuerj); //printf("[32] datai[%d] = %f\n", j, valueij); sumr = (valueri + valuerj) / 2.; sumi = (valueii + valueij) / 2.; difr = (valueri - valuerj) / 2.; difi = (valueii - valueij) / 2.; tempr = wr * sumi + wi * difr; tempi = wi * sumi - wr * difr; chunk_array_save(datar, i, sumr + tempr); chunk_array_save(datai, i, difi + tempi); chunk_array_save(datar, j, sumr - tempr); chunk_array_save(datai, j, tempi - difi); //printf("[70] Saving in datar the value %f in pos %d\n", sumr + tempr, i); //printf("[71] Saving in datai the value %f in pos %d\n", difi + tempi, i); //printf("[72] Saving in datar the value %f in pos %d\n", sumr - tempr, j); //printf("[73] Saving in datai the value %f in pos %d\n", tempi - difi, j); j += np2; } imin++; jmin--; wtemp = wr * wstpi; wr = wr * wstpr - wi * wstpi; wi = wi * wstpr + wtemp; L725: if (imin < jmin) { goto L710; } else if (imin > jmin) { goto L740; } if (ifrwd == 0) goto L740; for (i = imin; i <= ntot; i += np2) { chunk_array_get(datai, i, &valueii); //printf("[33] datai[%d] = %f\n", i, valueii); chunk_array_save(datai, i, -valueii); //printf("[73] Saving in datai the value %f in pos %d\n", -valueii, i); } L740: np2 *= 2; ntot *= 2; j = ntot + 1; imax = ntot / 2 + 1; L745: imin = imax - nhalf; i = imin; goto L755; L750: ; chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); //printf("[34] datar[%d] = %f\n", i, valueri); //printf("[35] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, -valueii); //printf("[74] Saving in datar the value %f in pos %d\n", valueri, j); //printf("[75] Saving in datai the value %f in pos %d\n", -valueii, j); L755: i++; j--; if (i < imax) goto L750; chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); //printf("[36] datar[%d] = %f\n", imin, valuerimin); //printf("[37] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin - valueiimin); chunk_array_save(datai, j, 0.); //printf("[75] Saving in datar the value %f in pos %d\n", valuerimin - valueiimin, j); //printf("[76] Saving in datai the value %f in pos %d\n", 0., j); if (i >= j) { goto L780; } else { goto L770; } L765: chunk_array_get(datar, i, &valueri); chunk_array_get(datai, i, &valueii); //printf("[38] datar[%d] = %f\n", i, valueri); //printf("[39] datai[%d] = %f\n", i, valueii); chunk_array_save(datar, j, valueri); chunk_array_save(datai, j, valueii); //printf("[77] Saving in datar the value %f in pos %d\n", valueri, j); //printf("[78] Saving in datai the value %f in pos %d\n", valueii, j); L770: i--; j--; if (i > imin) goto L765; chunk_array_get(datar, imin, &valuerimin); chunk_array_get(datai, imin, &valueiimin); //printf("[40] datar[%d] = %f\n", imin, valuerimin); //printf("[41] datai[%d] = %f\n", imin, valueiimin); chunk_array_save(datar, j, valuerimin + valueiimin); chunk_array_save(datai, j, 0.); //printf("[75] Saving in datar the value %f in pos %d\n", valuerimin + valueiimin, j); //printf("[76] Saving in datai the value %f in pos %d\n", 0., j); imax = imin; goto L745; L780: ; chunk_array_get(datai, 1, &valuei1); chunk_array_get(datar, 1, &valuer1); //printf("[42] datai[1] = %f\n", valuei1); //printf("[43] datar[1] = %f\n", valuer1); chunk_array_save(datar, 1, valuei1 + valuer1); chunk_array_save(datai, 1, 0.); //printf("[77] Saving in datar the value %f in pos %d\n", valuei1 + valuer1, 1); //printf("[78] Saving in datai the value %f in pos %d\n", 0., 1); goto L900; /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/ L800: if (nprev <= 2) goto L900; for (i3 = 1; i3 <= ntot; i3 += np2) { i2max = i3 + np2 - np1; for (i2 = i3; i2 <= i2max; i2 += np1) { imax = i2 + np1 - 1; imin = i2 + i1rng; jmax = 2 * i3 + np1 - imin; if (i2 > i3) jmax += np2; if (idim > 2) { j = jmax + np0; for (i = imin; i <= imax; i++) { chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[44] datar[%d] = %f\n", j, valuerj); //printf("[45] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); //printf("[79] Saving in datar the value %f in pos %d\n", valuerj, i); //printf("[80] Saving in datai the value %f in pos %d\n", -valueij, i); j--; } } j = jmax; for (i = imin; i <= imax; i += np0) { chunk_array_get(datar, j, &valuerj); chunk_array_get(datai, j, &valueij); //printf("[46] datar[%d] = %f\n", j, valuerj); //printf("[47] datai[%d] = %f\n", j, valueij); chunk_array_save(datar, i, valuerj); chunk_array_save(datai, i, -valueij); //printf("[81] Saving in datar the value %f in pos %d\n", valuerj, i); //printf("[82] Saving in datai the value %f in pos %d\n", -valueij, i); j -= np0; } } } /*end of loop on each dimension*/ L900: np0 = np1; np1 = np2; nprev = n; } L920: return; }