EMAN2
native_fft.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
00007  * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  */
00035 
00036 #ifdef NATIVE_FFT
00037 
00038 #include <cstring>
00039 #include <cmath>
00040 #include "native_fft.h" 
00041 #include "log.h"
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <math.h>
00046 
00047 
00048 #include "util.h"
00049 
00050 using namespace EMAN;
00051 
00052 int Nativefft::fftmcf_(float *a, float *b, integer *ntot, 
00053             integer *n, integer *nspan, integer *isn)
00054 {
00055     /* System generated locals */
00056     integer i__1;
00057     float r__1, r__2;
00058     int   status = -1;
00059 
00060     /* Builtin functions */
00061     // double atan(), cos(), sin(), sqrt();
00062 
00063     /* Local variables */
00064     static integer nfac[11];
00065     static float radf;
00066     static integer maxf, maxp, i__, j, k, m, kspan;
00067     static float c1, c2, c3;
00068     static integer kspnn, k1, k2, k3, k4;
00069     static float s1, s2, s3, aa, bb, cd, aj, c72;
00070     static integer jc;
00071     static float ck[23], ak;
00072     static integer jf;
00073     static float bk, bj;
00074     static integer jj;
00075     static float at[23], bt[23], sd;
00076     static integer kk;
00077     static float s72;
00078     static integer nn, np[209];
00079     static float sk[23];
00080     static integer ks, kt, nt;
00081     static float s120, rad, ajm, akm;
00082     static integer inc;
00083     static float ajp, akp, bkp, bkm, bjp, bjm;
00084 
00085 
00086 /* ------- ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23 */
00087 
00088 /*      EQUIVALENCE (I,II) */
00089 
00090 /* ----- THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH */
00091 /* ----- THE ARRAY DIMENSIONS. */
00092 
00093     /* Parameter adjustments */
00094     --b;
00095     --a;
00096 
00097     /* Function Body */
00098     maxf = 23;
00099     maxp = 209;
00100     if (*n < 2) {
00101         return 0;
00102     }
00103     inc = *isn;
00104     rad = atan((float)1.) * (float)8.;
00105     s72 = rad / (float)5.;
00106     c72 = cos(s72);
00107     s72 = sin(s72);
00108     s120 = sqrt((float).75);
00109     if (*isn >= 0) {
00110         goto L10;
00111     }
00112     s72 = -s72;
00113     s120 = -s120;
00114     rad = -rad;
00115     inc = -inc;
00116 L10:
00117     nt = inc * *ntot;
00118     ks = inc * *nspan;
00119     kspan = ks;
00120     nn = nt - inc;
00121     jc = ks / *n;
00122     radf = rad * (float) jc * (float).5;
00123     i__ = 0;
00124     jf = 0;
00125 
00126 /* -------DETERMINE THE FACTORS OF N------------------------- */
00127 
00128     m = 0;
00129     k = *n;
00130     goto L20;
00131 L15:
00132     ++m;
00133     nfac[m - 1] = 4;
00134     k /= 16;
00135 L20:
00136     if (k - (k / 16 << 4) == 0) {
00137         goto L15;
00138     }
00139     j = 3;
00140     jj = 9;
00141     goto L30;
00142 L25:
00143     ++m;
00144     nfac[m - 1] = j;
00145     k /= jj;
00146 L30:
00147     if (k % jj == 0) {
00148         goto L25;
00149     }
00150     j += 2;
00151 /* Computing 2nd power */
00152     i__1 = j;
00153     jj = i__1 * i__1;
00154     if (jj <= k) {
00155         goto L30;
00156     }
00157     if (k > 4) {
00158         goto L40;
00159     }
00160     kt = m;
00161     nfac[m] = k;
00162     if (k != 1) {
00163         ++m;
00164     }
00165     goto L80;
00166 L40:
00167     if (k - (k / 4 << 2) != 0) {
00168         goto L50;
00169     }
00170     ++m;
00171     nfac[m - 1] = 2;
00172     k /= 4;
00173 L50:
00174     kt = m;
00175     j = 2;
00176 L60:
00177     if (k % j != 0) {
00178         goto L70;
00179     }
00180     ++m;
00181     nfac[m - 1] = j;
00182     k /= j;
00183 L70:
00184     j = ((j + 1) / 2 << 1) + 1;
00185     if (j <= k) {
00186         goto L60;
00187     }
00188 L80:
00189     if (kt == 0) {
00190         goto L100;
00191     }
00192     j = kt;
00193 L90:
00194     ++m;
00195     nfac[m - 1] = nfac[j - 1];
00196     --j;
00197     if (j != 0) {
00198         goto L90;
00199     }
00200 
00201 /* -------COMPUTE FOURIER TRANSFORM----------------------------- */
00202 
00203 L100:
00204     sd = radf / (float) kspan;
00205 /* Computing 2nd power */
00206     r__1 = sin(sd);
00207     cd = r__1 * r__1 * (float)2.;
00208     sd = sin(sd + sd);
00209     kk = 1;
00210     ++i__;
00211     if (nfac[i__ - 1] != 2) {
00212         goto L400;
00213     }
00214 
00215 /* -------TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)- */
00216 
00217     kspan /= 2;
00218     k1 = kspan + 2;
00219 L210:
00220     k2 = kk + kspan;
00221     ak = a[k2];
00222     bk = b[k2];
00223     a[k2] = a[kk] - ak;
00224     b[k2] = b[kk] - bk;
00225     a[kk] += ak;
00226     b[kk] += bk;
00227     kk = k2 + kspan;
00228     if (kk <= nn) {
00229         goto L210;
00230     }
00231     kk -= nn;
00232     if (kk <= jc) {
00233         goto L210;
00234     }
00235     if (kk > kspan) {
00236         goto L800;
00237     }
00238 L220:
00239     c1 = (float)1. - cd;
00240     s1 = sd;
00241 L230:
00242     k2 = kk + kspan;
00243     ak = a[kk] - a[k2];
00244     bk = b[kk] - b[k2];
00245     a[kk] += a[k2];
00246     b[kk] += b[k2];
00247     a[k2] = c1 * ak - s1 * bk;
00248     b[k2] = s1 * ak + c1 * bk;
00249     kk = k2 + kspan;
00250     if (kk < nt) {
00251         goto L230;
00252     }
00253     k2 = kk - nt;
00254     c1 = -c1;
00255     kk = k1 - k2;
00256     if (kk > k2) {
00257         goto L230;
00258     }
00259     ak = c1 - (cd * c1 + sd * s1);
00260     s1 = sd * c1 - cd * s1 + s1;
00261 
00262 /* ------THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION-- */
00263 /*      ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE */
00264 /*      C1 = AK */
00265 
00266 /* Computing 2nd power */
00267     r__1 = ak;
00268 /* Computing 2nd power */
00269     r__2 = s1;
00270     c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00271     s1 = c1 * s1;
00272     c1 *= ak;
00273     kk += jc;
00274     if (kk < k2) {
00275         goto L230;
00276     }
00277     k1 = k1 + inc + inc;
00278     kk = (k1 - kspan) / 2 + jc;
00279     if (kk <= jc + jc) {
00280         goto L220;
00281     }
00282     goto L100;
00283 
00284 /* -------TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)---------- */
00285 
00286 L320:
00287     k1 = kk + kspan;
00288     k2 = k1 + kspan;
00289     ak = a[kk];
00290     bk = b[kk];
00291     aj = a[k1] + a[k2];
00292     bj = b[k1] + b[k2];
00293     a[kk] = ak + aj;
00294     b[kk] = bk + bj;
00295     ak = aj * (float)-.5 + ak;
00296     bk = bj * (float)-.5 + bk;
00297     aj = (a[k1] - a[k2]) * s120;
00298     bj = (b[k1] - b[k2]) * s120;
00299     a[k1] = ak - bj;
00300     b[k1] = bk + aj;
00301     a[k2] = ak + bj;
00302     b[k2] = bk - aj;
00303     kk = k2 + kspan;
00304     if (kk < nn) {
00305         goto L320;
00306     }
00307     kk -= nn;
00308     if (kk <= kspan) {
00309         goto L320;
00310     }
00311     goto L700;
00312 
00313 /* -------TRANSFORM FOR FACTOR OF 4--------------- */
00314 
00315 L400:
00316     if (nfac[i__ - 1] != 4) {
00317         goto L600;
00318     }
00319     kspnn = kspan;
00320     kspan /= 4;
00321 L410:
00322     c1 = (float)1.;
00323     s1 = (float)0.;
00324 L420:
00325     k1 = kk + kspan;
00326     k2 = k1 + kspan;
00327     k3 = k2 + kspan;
00328     akp = a[kk] + a[k2];
00329     akm = a[kk] - a[k2];
00330     ajp = a[k1] + a[k3];
00331     ajm = a[k1] - a[k3];
00332     a[kk] = akp + ajp;
00333     ajp = akp - ajp;
00334     bkp = b[kk] + b[k2];
00335     bkm = b[kk] - b[k2];
00336     bjp = b[k1] + b[k3];
00337     bjm = b[k1] - b[k3];
00338     b[kk] = bkp + bjp;
00339     bjp = bkp - bjp;
00340     if (*isn < 0) {
00341         goto L450;
00342     }
00343     akp = akm - bjm;
00344     akm += bjm;
00345     bkp = bkm + ajm;
00346     bkm -= ajm;
00347     if (s1 == (float)0.) {
00348         goto L460;
00349     }
00350 L430:
00351     a[k1] = akp * c1 - bkp * s1;
00352     b[k1] = akp * s1 + bkp * c1;
00353     a[k2] = ajp * c2 - bjp * s2;
00354     b[k2] = ajp * s2 + bjp * c2;
00355     a[k3] = akm * c3 - bkm * s3;
00356     b[k3] = akm * s3 + bkm * c3;
00357     kk = k3 + kspan;
00358     if (kk <= nt) {
00359         goto L420;
00360     }
00361 L440:
00362     c2 = c1 - (cd * c1 + sd * s1);
00363     s1 = sd * c1 - cd * s1 + s1;
00364 
00365 /* -------THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION */
00366 /*       ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE */
00367 /*       C1 = C2 */
00368 
00369 /* Computing 2nd power */
00370     r__1 = c2;
00371 /* Computing 2nd power */
00372     r__2 = s1;
00373     c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00374     s1 = c1 * s1;
00375     c1 *= c2;
00376 /* Computing 2nd power */
00377     r__1 = c1;
00378 /* Computing 2nd power */
00379     r__2 = s1;
00380     c2 = r__1 * r__1 - r__2 * r__2;
00381     s2 = c1 * (float)2. * s1;
00382     c3 = c2 * c1 - s2 * s1;
00383     s3 = c2 * s1 + s2 * c1;
00384     kk = kk - nt + jc;
00385     if (kk <= kspan) {
00386         goto L420;
00387     }
00388     kk = kk - kspan + inc;
00389     if (kk <= jc) {
00390         goto L410;
00391     }
00392     if (kspan == jc) {
00393         goto L800;
00394     }
00395     goto L100;
00396 L450:
00397     akp = akm + bjm;
00398     akm -= bjm;
00399     bkp = bkm - ajm;
00400     bkm += ajm;
00401     if (s1 != (float)0.) {
00402         goto L430;
00403     }
00404 L460:
00405     a[k1] = akp;
00406     b[k1] = bkp;
00407     a[k2] = ajp;
00408     b[k2] = bjp;
00409     a[k3] = akm;
00410     b[k3] = bkm;
00411     kk = k3 + kspan;
00412     if (kk <= nt) {
00413         goto L420;
00414     }
00415     goto L440;
00416 
00417 /* -------TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)-------- */
00418 
00419 L510:
00420 /* Computing 2nd power */
00421     r__1 = c72;
00422 /* Computing 2nd power */
00423     r__2 = s72;
00424     c2 = r__1 * r__1 - r__2 * r__2;
00425     s2 = c72 * (float)2. * s72;
00426 L520:
00427     k1 = kk + kspan;
00428     k2 = k1 + kspan;
00429     k3 = k2 + kspan;
00430     k4 = k3 + kspan;
00431     akp = a[k1] + a[k4];
00432     akm = a[k1] - a[k4];
00433     bkp = b[k1] + b[k4];
00434     bkm = b[k1] - b[k4];
00435     ajp = a[k2] + a[k3];
00436     ajm = a[k2] - a[k3];
00437     bjp = b[k2] + b[k3];
00438     bjm = b[k2] - b[k3];
00439     aa = a[kk];
00440     bb = b[kk];
00441     a[kk] = aa + akp + ajp;
00442     b[kk] = bb + bkp + bjp;
00443     ak = akp * c72 + ajp * c2 + aa;
00444     bk = bkp * c72 + bjp * c2 + bb;
00445     aj = akm * s72 + ajm * s2;
00446     bj = bkm * s72 + bjm * s2;
00447     a[k1] = ak - bj;
00448     a[k4] = ak + bj;
00449     b[k1] = bk + aj;
00450     b[k4] = bk - aj;
00451     ak = akp * c2 + ajp * c72 + aa;
00452     bk = bkp * c2 + bjp * c72 + bb;
00453     aj = akm * s2 - ajm * s72;
00454     bj = bkm * s2 - bjm * s72;
00455     a[k2] = ak - bj;
00456     a[k3] = ak + bj;
00457     b[k2] = bk + aj;
00458     b[k3] = bk - aj;
00459     kk = k4 + kspan;
00460     if (kk < nn) {
00461         goto L520;
00462     }
00463     kk -= nn;
00464     if (kk <= kspan) {
00465         goto L520;
00466     }
00467     goto L700;
00468 
00469 /* -------TRANSFORM FOR ODD FACTORS-------------------- */
00470 
00471 L600:
00472     k = nfac[i__ - 1];
00473     kspnn = kspan;
00474     kspan /= k;
00475     if (k == 3) {
00476         goto L320;
00477     }
00478     if (k == 5) {
00479         goto L510;
00480     }
00481     if (k == jf) {
00482         goto L640;
00483     }
00484     jf = k;
00485     s1 = rad / (float) k;
00486     c1 = cos(s1);
00487     s1 = sin(s1);
00488     if (jf > maxf) {
00489         goto L998;
00490     }
00491     ck[jf - 1] = (float)1.;
00492     sk[jf - 1] = (float)0.;
00493     j = 1;
00494 L630:
00495     ck[j - 1] = ck[k - 1] * c1 + sk[k - 1] * s1;
00496     sk[j - 1] = ck[k - 1] * s1 - sk[k - 1] * c1;
00497     --k;
00498     ck[k - 1] = ck[j - 1];
00499     sk[k - 1] = -sk[j - 1];
00500     ++j;
00501     if (j < k) {
00502         goto L630;
00503     }
00504 L640:
00505     k1 = kk;
00506     k2 = kk + kspnn;
00507     aa = a[kk];
00508     bb = b[kk];
00509     ak = aa;
00510     bk = bb;
00511     j = 1;
00512     k1 += kspan;
00513 L650:
00514     k2 -= kspan;
00515     ++j;
00516     at[j - 1] = a[k1] + a[k2];
00517     ak = at[j - 1] + ak;
00518     bt[j - 1] = b[k1] + b[k2];
00519     bk = bt[j - 1] + bk;
00520     ++j;
00521     at[j - 1] = a[k1] - a[k2];
00522     bt[j - 1] = b[k1] - b[k2];
00523     k1 += kspan;
00524     if (k1 < k2) {
00525         goto L650;
00526     }
00527     a[kk] = ak;
00528     b[kk] = bk;
00529     k1 = kk;
00530     k2 = kk + kspnn;
00531     j = 1;
00532 L660:
00533     k1 += kspan;
00534     k2 -= kspan;
00535     jj = j;
00536     ak = aa;
00537     bk = bb;
00538     aj = (float)0.;
00539     bj = (float)0.;
00540     k = 1;
00541 L670:
00542     ++k;
00543     ak = at[k - 1] * ck[jj - 1] + ak;
00544     bk = bt[k - 1] * ck[jj - 1] + bk;
00545     ++k;
00546     aj = at[k - 1] * sk[jj - 1] + aj;
00547     bj = bt[k - 1] * sk[jj - 1] + bj;
00548     jj += j;
00549     if (jj > jf) {
00550         jj -= jf;
00551     }
00552     if (k < jf) {
00553         goto L670;
00554     }
00555     k = jf - j;
00556     a[k1] = ak - bj;
00557     b[k1] = bk + aj;
00558     a[k2] = ak + bj;
00559     b[k2] = bk - aj;
00560     ++j;
00561     if (j < k) {
00562         goto L660;
00563     }
00564     kk += kspnn;
00565     if (kk <= nn) {
00566         goto L640;
00567     }
00568     kk -= nn;
00569     if (kk <= kspan) {
00570         goto L640;
00571     }
00572 
00573 /* ------- MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4) */
00574 
00575 L700:
00576     if (i__ == m) {
00577         goto L800;
00578     }
00579     kk = jc + 1;
00580 L710:
00581     c2 = (float)1. - cd;
00582     s1 = sd;
00583 L720:
00584     c1 = c2;
00585     s2 = s1;
00586     kk += kspan;
00587 L730:
00588     ak = a[kk];
00589     a[kk] = c2 * ak - s2 * b[kk];
00590     b[kk] = s2 * ak + c2 * b[kk];
00591     kk += kspnn;
00592     if (kk <= nt) {
00593         goto L730;
00594     }
00595     ak = s1 * s2;
00596     s2 = s1 * c2 + c1 * s2;
00597     c2 = c1 * c2 - ak;
00598     kk = kk - nt + kspan;
00599     if (kk <= kspnn) {
00600         goto L730;
00601     }
00602     c2 = c1 - (cd * c1 + sd * s1);
00603     s1 += sd * c1 - cd * s1;
00604 
00605 /* -------THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION */
00606 /*       ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY */
00607 /*       BE DELETED. */
00608 
00609 /* Computing 2nd power */
00610     r__1 = c2;
00611 /* Computing 2nd power */
00612     r__2 = s1;
00613     c1 = (float).5 / (r__1 * r__1 + r__2 * r__2) + (float).5;
00614     s1 = c1 * s1;
00615     c2 = c1 * c2;
00616     kk = kk - kspnn + jc;
00617     if (kk <= kspan) {
00618         goto L720;
00619     }
00620     kk = kk - kspan + jc + inc;
00621     if (kk <= jc + jc) {
00622         goto L710;
00623     }
00624     goto L100;
00625 
00626 /* -------PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES- */
00627 /*       PERMUTATION FOR SQUARE FACTORS OF N */
00628 
00629 L800:
00630     np[0] = ks;
00631     if (kt == 0) {
00632         goto L890;
00633     }
00634     k = kt + kt + 1;
00635     if (m < k) {
00636         --k;
00637     }
00638     j = 1;
00639     np[k] = jc;
00640 L810:
00641     np[j] = np[j - 1] / nfac[j - 1];
00642     np[k - 1] = np[k] * nfac[j - 1];
00643     ++j;
00644     --k;
00645     if (j < k) {
00646         goto L810;
00647     }
00648     k3 = np[k];
00649     kspan = np[1];
00650     kk = jc + 1;
00651     j = 1;
00652     k2 = kspan + 1;
00653     if (*n != *ntot) {
00654         goto L850;
00655     }
00656 
00657 /* -------PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE) */
00658 
00659 L820:
00660     ak = a[kk];
00661     a[kk] = a[k2];
00662     a[k2] = ak;
00663     bk = b[kk];
00664     b[kk] = b[k2];
00665     b[k2] = bk;
00666     kk += inc;
00667     k2 = kspan + k2;
00668     if (k2 < ks) {
00669         goto L820;
00670     }
00671 L830:
00672     k2 -= np[j - 1];
00673     ++j;
00674     k2 = np[j] + k2;
00675     if (k2 > np[j - 1]) {
00676         goto L830;
00677     }
00678     j = 1;
00679 L840:
00680     if (kk < k2) {
00681         goto L820;
00682     }
00683     kk += inc;
00684     k2 = kspan + k2;
00685     if (k2 < ks) {
00686         goto L840;
00687     }
00688     if (kk < ks) {
00689         goto L830;
00690     }
00691     jc = k3;
00692     goto L890;
00693 
00694 /* -------PERMUTATION FOR MULTIVARIATE TRANSFORM------ */
00695 
00696 L850:
00697     k = kk + jc;
00698 L860:
00699     ak = a[kk];
00700     a[kk] = a[k2];
00701     a[k2] = ak;
00702     bk = b[kk];
00703     b[kk] = b[k2];
00704     b[k2] = bk;
00705     kk += inc;
00706     k2 += inc;
00707     if (kk < k) {
00708         goto L860;
00709     }
00710     kk = kk + ks - jc;
00711     k2 = k2 + ks - jc;
00712     if (kk < nt) {
00713         goto L850;
00714     }
00715     k2 = k2 - nt + kspan;
00716     kk = kk - nt + jc;
00717     if (k2 < ks) {
00718         goto L850;
00719     }
00720 L870:
00721     k2 -= np[j - 1];
00722     ++j;
00723     k2 = np[j] + k2;
00724     if (k2 > np[j - 1]) {
00725         goto L870;
00726     }
00727     j = 1;
00728 L880:
00729     if (kk < k2) {
00730         goto L850;
00731     }
00732     kk += jc;
00733     k2 = kspan + k2;
00734     if (k2 < ks) {
00735         goto L880;
00736     }
00737     if (kk < ks) {
00738         goto L870;
00739     }
00740     jc = k3;
00741 L890:
00742     if ((kt << 1) + 1 >= m) {
00743         return 0;
00744     }
00745     kspnn = np[kt];
00746 
00747 /* -------PERMUTATION FOR SQUARE-FREE FACTORS OF N----- */
00748 
00749     j = m - kt;
00750     nfac[j] = 1;
00751 L900:
00752     nfac[j - 1] *= nfac[j];
00753     --j;
00754     if (j != kt) {
00755         goto L900;
00756     }
00757     ++kt;
00758     nn = nfac[kt - 1] - 1;
00759     if (nn > maxp) {
00760         goto L998;
00761     }
00762     jj = 0;
00763     j = 0;
00764     goto L906;
00765 L902:
00766     jj -= k2;
00767     k2 = kk;
00768     ++k;
00769     kk = nfac[k - 1];
00770 L904:
00771     jj = kk + jj;
00772     if (jj >= k2) {
00773         goto L902;
00774     }
00775     np[j - 1] = jj;
00776 L906:
00777     k2 = nfac[kt - 1];
00778     k = kt + 1;
00779     kk = nfac[k - 1];
00780     ++j;
00781     if (j <= nn) {
00782         goto L904;
00783     }
00784 
00785 /* -------DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1 */
00786 
00787     j = 0;
00788     goto L914;
00789 L910:
00790     k = kk;
00791     kk = np[k - 1];
00792     np[k - 1] = -kk;
00793     if (kk != j) {
00794         goto L910;
00795     }
00796     k3 = kk;
00797 L914:
00798     ++j;
00799     kk = np[j - 1];
00800     if (kk < 0) {
00801         goto L914;
00802     }
00803     if (kk != j) {
00804         goto L910;
00805     }
00806     np[j - 1] = -j;
00807     if (j != nn) {
00808         goto L914;
00809     }
00810     maxf = inc * maxf;
00811 
00812 /* -------REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES- */
00813 
00814     goto L950;
00815 L924:
00816     --j;
00817     if (np[j - 1] < 0) {
00818         goto L924;
00819     }
00820     jj = jc;
00821 L926:
00822     kspan = jj;
00823     if (jj > maxf) {
00824         kspan = maxf;
00825     }
00826     jj -= kspan;
00827     k = np[j - 1];
00828     kk = jc * k + i__ + jj;
00829     k1 = kk + kspan;
00830     k2 = 0;
00831 L928:
00832     ++k2;
00833     at[k2 - 1] = a[k1];
00834     bt[k2 - 1] = b[k1];
00835     k1 -= inc;
00836     if (k1 != kk) {
00837         goto L928;
00838     }
00839 L932:
00840     k1 = kk + kspan;
00841     k2 = k1 - jc * (k + np[k - 1]);
00842     k = -np[k - 1];
00843 L936:
00844     a[k1] = a[k2];
00845     b[k1] = b[k2];
00846     k1 -= inc;
00847     k2 -= inc;
00848     if (k1 != kk) {
00849         goto L936;
00850     }
00851     kk = k2;
00852     if (k != j) {
00853         goto L932;
00854     }
00855     k1 = kk + kspan;
00856     k2 = 0;
00857 L940:
00858     ++k2;
00859     a[k1] = at[k2 - 1];
00860     b[k1] = bt[k2 - 1];
00861     k1 -= inc;
00862     if (k1 != kk) {
00863         goto L940;
00864     }
00865     if (jj != 0) {
00866         goto L926;
00867     }
00868     if (j != 1) {
00869         goto L924;
00870     }
00871 L950:
00872     j = k3 + 1;
00873     nt -= kspnn;
00874     i__ = nt - inc + 1;
00875     if (nt >= 0) {
00876         goto L924;
00877     }
00878     return 0;
00879 
00880 /* --------ERROR FINISH, INSUFFICIENT ARRAY STORAGE-- */
00881 
00882 L998:
00883     *isn = 0;
00884     return status;
00885 } /* fftmcf_ */
00886 
00887 #define work(i)  work[(i)-1]
00888 #define x(i)     x[(i)-1]
00889 #define y(i)     y[(i)-1]
00890 
00891 // 1d inplace FFT
00892 int Nativefft::fmrs_1rf(float *x, float *work, int nsam)
00893 {
00894         // input : x(nsam+2-nsam%2), work(nsam+2-nsam%2)
00895         // output: x(nsam+2-nsam%2), overwritten by Fourier coefficients
00896 
00897         int i, n, status;
00898         int inv=1;
00899 
00900         n = nsam;
00901 
00902         for (i=1; i<=n; i++)  work(i) = 0.0;
00903         status = fftmcf_(x,work,&n,&n,&n,&inv);
00904         // should put in appropriate exception handling here
00905         if (status == -1) {
00906                 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
00907                 fprintf(stderr, "the dimensions set to n = %d\n", nsam);  
00908                 exit(1); 
00909         }
00910 
00911         // check even/odd
00912         if (n%2 == 0) {
00913                 for (i=n+1;i>=3;i=i-2) {
00914                         x(i)   = x((i+1)/2);
00915                         x(i+1) = work((i+1)/2);
00916                 }
00917                 x(2)   = 0.0;
00918                 x(n+2) = 0.0;
00919         }
00920         else {
00921                 for (i=n;i>=3;i=i-2) {
00922                         x(i)   = x(i/2+1);
00923                         x(i+1) = work(i/2+1);
00924                 }
00925                 x(2)=0.0;
00926         }
00927         return status;
00928 }
00929 
00930 // 1d inplace IFFT
00931 int Nativefft::fmrs_1rb(float *x, float *work, int nsam)
00932 {
00933         // input:  x(nsam+2-nsam%2), work(nsam+2-nsam%2)
00934         // output: x(nsam+2-nsam%2), overwritten with Fourier coefficients
00935 
00936         int i, n, status;
00937         int inv=-1;
00938 
00939         n = nsam; 
00940 
00941         for (i=2;i<=n/2+1;i++) {
00942                 work(i)     = x(2*i)/n;
00943                 work(n-i+2) = -work(i);
00944         }
00945         work(1) = 0.0;
00946 
00947         for (i=1;i<=n/2+1;i++)  x(i) = x(2*i-1)/n;
00948         for (i=n;i>=n/2+2;i--)  x(i) = x(n-i+2);
00949 
00950         status = fftmcf_(x,work,&n,&n,&n,&inv);
00951         // should put in appropriate exception handling here
00952         if (status == -1) {
00953                 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
00954                 fprintf(stderr, "the dimensions set to n = %d\n", nsam);  
00955                 exit(1); 
00956         }
00957 
00958         return status;
00959 }
00960 
00961 #undef x
00962 #undef y
00963 #undef work
00964 
00965 
00966 //----------2D FFT INPLACE--------------------
00967 
00968 #define x(i,j) x[((j)-1)*lda + (i)-1]
00969 #define y(i,j) y[((j)-1)*lda + (i)-1]
00970 
00971 // 2d inplace fft
00972 int Nativefft::fmrs_2rf(float *x, float *work, int lda, int nsam, int nrow)
00973 {
00974         // input:  x(lda,nrow), work(lda)
00975         // output: x(lda,nrow) overwritten with Fourier coefficients
00976 
00977         int   ins, status=0, i, j;
00978 
00979         ins = lda;
00980         /*
00981         int  l;
00982         l=(int)(log2(nsam));
00983         for (j=1;j<=nrow;j++) {
00984                 Util::fftr_q(&x(1,j),l);
00985                 status = fmrs_1rf(&x(1,j), work, nsam);
00986         }
00987         */
00988 
00989         for (i=1;i<=lda;i=i+2){
00990                 status = fftmcf_(&x(i,1),&x(i+1,1),&nrow,&nrow,&nrow,&ins);
00991                 if (status == -1) {
00992                         fprintf(stderr, "Native FFT cannot be performed on images with one of ");
00993                         fprintf(stderr, "the dimensions set to n = %d\n", nrow);
00994                         exit(1); 
00995                 }
00996         }
00997 
00998         return status;
00999 }
01000 
01001 
01002 // 2d inplace IFFT
01003 int Nativefft::fmrs_2rb(float *y, float *work, int lda, int nsam, int nrow)
01004 {
01005         // input:  y(lda,nrow), work(lda)
01006         // output: y(lda,nrow), overwritten with Fourier coefficients
01007 
01008         int   ins, status=0, i, j;
01009         float q;
01010 
01011         ins=-lda;
01012 
01013         for (i=1; i<=lda; i=i+2) {
01014                 fftmcf_(&y(i,1),&y(i+1,1),&nrow,&nrow,&nrow,&ins);
01015                 if (status == -1) {
01016                         fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01017                         fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01018                         exit(1); 
01019                 }
01020         }
01021          
01022         // normalize for inverse
01023         q = 1.0/(float)(nrow);
01024         for (j=1; j<=nrow; j++)  for (i=1; i<=lda; i++) y(i,j)*=q;
01025 
01026         // need an inplace 1d ifft 
01027         for (j=1; j<=nrow; j++) status = fmrs_1rb(&y(1,j), work, nsam);
01028 
01029         return status;
01030 }
01031 #undef x
01032 #undef y
01033 
01034 //----------2D FFT out of place--------------------
01035 
01036 #define complexd(i,j) complexd[(j-1)*lda + i-1]
01037 #define reald(i,j)    reald[(j-1)*nsam + i-1]
01038 
01039 // 2D out of place fft
01040 int Nativefft::ftp_2rf(float *reald, float *complexd, int lda, int nsam, int nrow)
01041 {
01042         // input:  reald(nsam,nrow), 
01043         //  work(2*nsam+15)
01044         // output: complexd(lda,nrow) Fourier coefficients
01045 
01046         int   ins, status=0, i, j;
01047 
01048         float * work = (float*) malloc((2*nsam+15)*sizeof(float));
01049         if (!work) {
01050                  fprintf(stderr,"real_to_complex_nd(2df): failed to allocate work\n");
01051                  LOGERR("real_to_complex_nd(2df): failed to allocate work\n");
01052         }
01053         //  Do columns with fftpack     
01054         rffti(nsam, work);
01055 
01056         for (j=1; j<=nrow; j++)  {
01057                 memcpy(&complexd(2,j), &reald(1,j), nsam * sizeof(float));
01058                 rfftf(nsam, &complexd(2,j), work);
01059                 complexd(1,j) = complexd(2,j) ;
01060                 complexd(2,j) = 0.0f ;
01061                 if (nsam%2 == 0)  complexd(nsam+2,j) = 0.0f ;
01062         }
01063 
01064         free(work);
01065 
01066         ins = lda;
01067 
01068         for (i=1; i<=lda; i=i+2) {
01069            status = fftmcf_(&complexd(i,1),&complexd(i+1,1),&nrow,&nrow,&nrow,&ins);
01070            if (status == -1) {
01071                   fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01072                   fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01073                   exit(1); 
01074            }
01075         }
01076 
01077    return status;
01078 }
01079 
01080 
01081 int Nativefft::ftp_2rb(float *complexd, int lda, int nsam, int nrow)
01082 {
01083         // input:  complexd(lsd,nrow), 
01084         //  work(2*nsam+15)
01085         // output: complexd(lda,nrow) Fourier coefficients
01086 
01087         int   ins, status=0, i, j;
01088 
01089         ins = -lda;
01090 
01091         for (i=1; i<=lda; i=i+2) {
01092                 status = fftmcf_(&complexd(i,1),&complexd(i+1,1),&nrow,&nrow,&nrow,&ins);
01093                 if (status == -1) {
01094                         fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01095                         fprintf(stderr, "the dimensions set to n = %d\n", nrow);
01096                         exit(1); 
01097                 }
01098         }
01099 
01100 
01101         float * work = (float*) malloc((2*nsam+15)*sizeof(float));
01102         if (!work) {
01103                 fprintf(stderr,"real_to_complex_nd(2df): failed to allocate work\n");
01104                 LOGERR("real_to_complex_nd(2df): failed to allocate work\n");
01105         }
01106         //  Do inv columns with fftpack 
01107         rffti(nsam, work);
01108 
01109         for (j=1; j<=nrow; j++) {
01110                 memmove(&complexd(2,j), &complexd(3,j), (nsam-1) * sizeof(float));
01111                 rfftb(nsam, &complexd(1,j), work);
01112         }
01113         free(work);
01114         //  Normalize
01115         float nrm = 1.0f/float(nsam*nrow);
01116         for (int j = 1; j<=nrow; j++) for (int i = 1; i<=nsam; i++) complexd(i,j) *= nrm;
01117 
01118    return status;
01119 }
01120 #undef complexd
01121 #undef reald
01122 
01123 #define  reald(i,j,k)     reald[i-1 + (j-1+(k-1)*nrow)*nsam]
01124 #define  complexd(i,j,k)  complexd[i-1 + (j-1+(k-1)*nrow)*lda]
01125 // 3D out of place FFT
01126 int Nativefft::ftp_3rf(float *reald, float *complexd, int lda, int nsam, int nrow, int nslice)
01127 {
01128         // input :  reald(nsam,nrow,nslice)
01129         // output:  complexd(lda,nrow,nslice), overwritten with Fourier coefficients
01130 
01131         int ndr, i, j, k, status=0;
01132 
01133         for (k=1; k<=nslice;k ++) status = ftp_2rf(&reald(1,1,k), &complexd(1,1,k), lda, nsam, nrow);
01134 
01135         ndr=lda*nrow;
01136 
01137         for (j=1; j<=nrow; j++) {
01138                 for (i=1; i<=lda-1; i=i+2) {
01139                         status = fftmcf_(&complexd(i,j,1),&complexd(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01140                         if (status == -1) {
01141                                 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01142                                 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01143                                 exit(1); 
01144                         }
01145                 }
01146         }
01147 
01148         return status;
01149 }
01150 #undef reald
01151 
01152 // 3d out of place IFFT
01153 int Nativefft::ftp_3rb(float *complexd, int lda, int nsam, int nrow, int nslice)
01154 {
01155         // input:  complexd(lda,nrow,nslice)
01156         // output: complexd(lda,nrow,nslice), with Fourier coefficients 
01157 
01158         int ndr, i, j, k, status=0;
01159         float q;
01160 
01161         ndr=-lda*nrow;
01162 
01163         for (j=1; j<=nrow; j++) {
01164                 for (i=1; i<=lda-1; i=i+2) {
01165                         status = fftmcf_(&complexd(i,j,1),&complexd(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01166                         if (status == -1) {
01167                                 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01168                                 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01169                                 exit(1); 
01170                         }
01171                 }
01172         }
01173 
01174         // normalize for inverse
01175         q=1.0/(float)(nslice);
01176         for (k=1; k<=nslice; k++) for (j=1;j<=nrow;j++) for (i=1;i<=lda;i++) complexd(i,j,k) *= q;
01177 
01178         for (k=1; k<=nslice; k++) status = ftp_2rb(&complexd(1,1,k), lda, nsam, nrow);
01179 
01180         return status;
01181 }
01182 #undef complexd
01183 
01184 //---------------3D INPLACE FFT----------------------
01185 #define b(i,j,k) b[((k)-1)*lda*nrow + ((j)-1)*lda + (i) - 1]
01186 
01187 // 3d inplace FFT
01188 int Nativefft::fmrs_3rf(float *b, float *work, int lda, int nsam, int nrow, int nslice)
01189 {
01190         // input :  b(lda,nrow,nslice), work(lda)
01191         // output:  b(lda,nrow,nslice), overwritten with Fourier coefficients
01192 
01193         int ndr, i, j, k, status=0;
01194 
01195         ndr=lda*nrow;
01196 
01197         for (k=1; k<=nslice; k++) status = fmrs_2rf(&b(1,1,k), work, lda, nsam, nrow);
01198 
01199         for (j=1; j<=nrow; j++) {
01200                 for (i=1; i<=lda-1; i=i+2) {
01201                         status = fftmcf_(&b(i,j,1),&b(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01202                         if (status == -1) {
01203                                 fprintf(stderr, "Native FFT cannot be performed on images with one of ");
01204                                 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01205                                 exit(1); 
01206                         }
01207                 }
01208         }
01209         // should check for error here by looking at ndrt
01210 
01211         return status;
01212 }
01213 
01214 // 3d inplace IFFT
01215 int Nativefft::fmrs_3rb(float *b, float *work, int lda, int nsam, int nrow, int nslice)
01216 {
01217         // input:  b(lda,nrow,nslice), work(lda)
01218         // output: b(lda,nrow,nslice), overwritten with Fourier coefficients 
01219 
01220         int ndr, i, j, k, status=0;
01221         float q;
01222 
01223         ndr=-lda*nrow;
01224 
01225         for (j=1;j<=nrow;j++) {
01226                 for (i=1;i<=lda-1;i=i+2) {
01227                         status = fftmcf_(&b(i,j,1),&b(i+1,j,1),&nslice,&nslice,&nslice,&ndr);
01228                         if (status == -1) {
01229                                 fprintf(stderr, "Native IFT cannot be performed on images with one of ");
01230                                 fprintf(stderr, "the dimensions set to n = %d\n", nslice);
01231                                 exit(1); 
01232                         }
01233                 }
01234         }
01235 
01236         // should check for error here by looking at ndrt
01237 
01238         // normalize for inverse
01239         q=1.0/(float)(nslice);
01240         for (k=1;k<=nslice;k++) for (j=1;j<=nrow;j++) for (i=1;i<=lda;i++) b(i,j,k)*=q;
01241 
01242         for (k=1; k<=nslice; k++) status = fmrs_2rb(&b(1,1,k), work, lda, nsam, nrow);
01243         return status;
01244 }
01245 #undef b
01246 
01247 
01248 
01249 /*
01250 fftpack.c : A set of FFT routines in C.
01251 Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
01252 
01253 */
01254 
01255 /* isign is +1 for backward and -1 for forward transforms */
01256 
01257 #include <math.h>
01258 #include <stdio.h>
01259 //#define DOUBLE
01260 
01261 #ifdef DOUBLE
01262 #define Treal double
01263 #else
01264 #define Treal float
01265 #endif
01266 
01267 
01268 #define ref(u,a) u[a]
01269 
01270 #define MAXFAC 13    /* maximum number of factors in factorization of n */
01271 #define NSPECIAL 4   /* number of factors for which we have special-case routines */
01272 
01273 //#define __cplusplus
01274 //#ifdef __cplusplus
01275 //extern "C" {
01276 //#endif
01277 
01278 
01279 /* ----------------------------------------------------------------------
01280    passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
01281 ---------------------------------------------------------------------- */
01282 
01283 static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
01284   /* isign==+1 for backward transform */
01285   {
01286         int i, k, ah, ac;
01287         Treal ti2, tr2;
01288         if (ido <= 2) {
01289                 for (k=0; k<l1; k++) {
01290                         ah = k*ido;
01291                         ac = 2*k*ido;
01292                         ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);
01293                         ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);
01294                         ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);
01295                         ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
01296                 }
01297         } else {
01298                 for (k=0; k<l1; k++) {
01299                         for (i=0; i<ido-1; i+=2) {
01300                                 ah = i + k*ido;
01301                                 ac = i + 2*k*ido;
01302                                 ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);
01303                                 tr2      = ref(cc,ac) - ref(cc,ac + ido);
01304                                 ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
01305                                 ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
01306                                 ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
01307                                 ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
01308                         }
01309                 }
01310         }
01311   } /* passf2 */
01312 
01313 
01314 static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
01315       const Treal wa1[], const Treal wa2[], int isign)
01316   /* isign==+1 for backward transform */
01317   {
01318     static const Treal taur = -0.5;
01319     static const Treal taui = 0.866025403784439;
01320     int i, k, ac, ah;
01321     Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
01322     if (ido == 2) {
01323       for (k=1; k<=l1; k++) {
01324         ac = (3*k - 2)*ido;
01325         tr2 = ref(cc,ac) + ref(cc,ac + ido);
01326         cr2 = ref(cc,ac - ido) + taur*tr2;
01327         ah = (k - 1)*ido;
01328         ch[ah] = ref(cc,ac - ido) + tr2;
01329 
01330         ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
01331         ci2 = ref(cc,ac - ido + 1) + taur*ti2;
01332         ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
01333 
01334         cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
01335         ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
01336         ch[ah + l1*ido] = cr2 - ci3;
01337         ch[ah + 2*l1*ido] = cr2 + ci3;
01338         ch[ah + l1*ido + 1] = ci2 + cr3;
01339         ch[ah + 2*l1*ido + 1] = ci2 - cr3;
01340       }
01341     } else {
01342       for (k=1; k<=l1; k++) {
01343         for (i=0; i<ido-1; i+=2) {
01344           ac = i + (3*k - 2)*ido;
01345           tr2 = ref(cc,ac) + ref(cc,ac + ido);
01346           cr2 = ref(cc,ac - ido) + taur*tr2;
01347           ah = i + (k-1)*ido;
01348           ch[ah] = ref(cc,ac - ido) + tr2;
01349           ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
01350           ci2 = ref(cc,ac - ido + 1) + taur*ti2;
01351           ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
01352           cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
01353           ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
01354           dr2 = cr2 - ci3;
01355           dr3 = cr2 + ci3;
01356           di2 = ci2 + cr3;
01357           di3 = ci2 - cr3;
01358           ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
01359           ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
01360           ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
01361           ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
01362         }
01363       }
01364     }
01365   } /* passf3 */
01366 
01367 
01368 static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
01369       const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
01370   /* isign == -1 for forward transform and +1 for backward transform */
01371   {
01372     int i, k, ac, ah;
01373     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01374     if (ido == 2) {
01375       for (k=0; k<l1; k++) {
01376         ac = 4*k*ido + 1;
01377         ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
01378         ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
01379         tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
01380         ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
01381         tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
01382         tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
01383         ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
01384         tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
01385         ah = k*ido;
01386         ch[ah] = tr2 + tr3;
01387         ch[ah + 2*l1*ido] = tr2 - tr3;
01388         ch[ah + 1] = ti2 + ti3;
01389         ch[ah + 2*l1*ido + 1] = ti2 - ti3;
01390         ch[ah + l1*ido] = tr1 + isign*tr4;
01391         ch[ah + 3*l1*ido] = tr1 - isign*tr4;
01392         ch[ah + l1*ido + 1] = ti1 + isign*ti4;
01393         ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
01394       }
01395     } else {
01396       for (k=0; k<l1; k++) {
01397         for (i=0; i<ido-1; i+=2) {
01398           ac = i + 1 + 4*k*ido;
01399           ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
01400           ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
01401           ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
01402           tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
01403           tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
01404           tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
01405           ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
01406           tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
01407           ah = i + k*ido;
01408           ch[ah] = tr2 + tr3;
01409           cr3 = tr2 - tr3;
01410           ch[ah + 1] = ti2 + ti3;
01411           ci3 = ti2 - ti3;
01412           cr2 = tr1 + isign*tr4;
01413           cr4 = tr1 - isign*tr4;
01414           ci2 = ti1 + isign*ti4;
01415           ci4 = ti1 - isign*ti4;
01416           ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
01417           ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
01418           ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
01419           ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
01420           ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
01421           ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
01422         }
01423       }
01424     }
01425   } /* passf4 */
01426 
01427 
01428 static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
01429       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
01430   /* isign == -1 for forward transform and +1 for backward transform */
01431   {
01432     static const Treal tr11 = 0.309016994374947;
01433     static const Treal ti11 = 0.951056516295154;
01434     static const Treal tr12 = -0.809016994374947;
01435     static const Treal ti12 = 0.587785252292473;
01436     int i, k, ac, ah;
01437     Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
01438         ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01439     if (ido == 2) {
01440       for (k = 1; k <= l1; ++k) {
01441         ac = (5*k - 4)*ido + 1;
01442         ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
01443         ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
01444         ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
01445         ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
01446         tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
01447         tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
01448         tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
01449         tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
01450         ah = (k - 1)*ido;
01451         ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
01452         ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
01453         cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
01454         ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
01455         cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
01456         ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
01457         cr5 = isign*(ti11*tr5 + ti12*tr4);
01458         ci5 = isign*(ti11*ti5 + ti12*ti4);
01459         cr4 = isign*(ti12*tr5 - ti11*tr4);
01460         ci4 = isign*(ti12*ti5 - ti11*ti4);
01461         ch[ah + l1*ido] = cr2 - ci5;
01462         ch[ah + 4*l1*ido] = cr2 + ci5;
01463         ch[ah + l1*ido + 1] = ci2 + cr5;
01464         ch[ah + 2*l1*ido + 1] = ci3 + cr4;
01465         ch[ah + 2*l1*ido] = cr3 - ci4;
01466         ch[ah + 3*l1*ido] = cr3 + ci4;
01467         ch[ah + 3*l1*ido + 1] = ci3 - cr4;
01468         ch[ah + 4*l1*ido + 1] = ci2 - cr5;
01469       }
01470     } else {
01471       for (k=1; k<=l1; k++) {
01472         for (i=0; i<ido-1; i+=2) {
01473           ac = i + 1 + (k*5 - 4)*ido;
01474           ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
01475           ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
01476           ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
01477           ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
01478           tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
01479           tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
01480           tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
01481           tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
01482           ah = i + (k - 1)*ido;
01483           ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
01484           ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
01485           cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
01486 
01487           ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
01488           cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
01489 
01490           ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
01491           cr5 = isign*(ti11*tr5 + ti12*tr4);
01492           ci5 = isign*(ti11*ti5 + ti12*ti4);
01493           cr4 = isign*(ti12*tr5 - ti11*tr4);
01494           ci4 = isign*(ti12*ti5 - ti11*ti4);
01495           dr3 = cr3 - ci4;
01496           dr4 = cr3 + ci4;
01497           di3 = ci3 + cr4;
01498           di4 = ci3 - cr4;
01499           dr5 = cr2 + ci5;
01500           dr2 = cr2 - ci5;
01501           di5 = ci2 - cr5;
01502           di2 = ci2 + cr5;
01503           ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
01504           ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
01505           ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
01506           ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
01507           ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
01508           ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
01509           ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
01510           ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
01511         }
01512       }
01513     }
01514   } /* passf5 */
01515 
01516 
01517 static void passf(int *nac, int ido, int ip, int l1, int idl1,
01518       Treal cc[], Treal ch[],
01519       const Treal wa[], int isign)
01520   /* isign is -1 for forward transform and +1 for backward transform */
01521   {
01522     int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
01523     Treal wai, war;
01524 
01525     idot = ido / 2;
01526     nt = ip*idl1;
01527     ipph = (ip + 1) / 2;
01528     idp = ip*ido;
01529     if (ido >= l1) {
01530       for (j=1; j<ipph; j++) {
01531         jc = ip - j;
01532         for (k=0; k<l1; k++) {
01533           for (i=0; i<ido; i++) {
01534             ch[i + (k + j*l1)*ido] =
01535                 ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
01536             ch[i + (k + jc*l1)*ido] =
01537                 ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
01538           }
01539         }
01540       }
01541       for (k=0; k<l1; k++)
01542         for (i=0; i<ido; i++)
01543           ch[i + k*ido] = ref(cc,i + k*ip*ido);
01544     } else {
01545       for (j=1; j<ipph; j++) {
01546         jc = ip - j;
01547         for (i=0; i<ido; i++) {
01548           for (k=0; k<l1; k++) {
01549             ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
01550                 ip)*ido);
01551             ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
01552                 ip)*ido);
01553           }
01554         }
01555       }
01556       for (i=0; i<ido; i++)
01557         for (k=0; k<l1; k++)
01558           ch[i + k*ido] = ref(cc,i + k*ip*ido);
01559     }
01560 
01561     idl = 2 - ido;
01562     inc = 0;
01563     for (l=1; l<ipph; l++) {
01564       lc = ip - l;
01565       idl += ido;
01566       for (ik=0; ik<idl1; ik++) {
01567         cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
01568         cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
01569       }
01570       idlj = idl;
01571       inc += ido;
01572       for (j=2; j<ipph; j++) {
01573         jc = ip - j;
01574         idlj += inc;
01575         if (idlj > idp) idlj -= idp;
01576         war = wa[idlj - 2];
01577         wai = wa[idlj-1];
01578         for (ik=0; ik<idl1; ik++) {
01579           cc[ik + l*idl1] += war*ch[ik + j*idl1];
01580           cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
01581         }
01582       }
01583     }
01584     for (j=1; j<ipph; j++)
01585       for (ik=0; ik<idl1; ik++)
01586         ch[ik] += ch[ik + j*idl1];
01587     for (j=1; j<ipph; j++) {
01588       jc = ip - j;
01589       for (ik=1; ik<idl1; ik+=2) {
01590         ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
01591         ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
01592         ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
01593         ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
01594       }
01595     }
01596     *nac = 1;
01597     if (ido == 2) return;
01598     *nac = 0;
01599     for (ik=0; ik<idl1; ik++)
01600       cc[ik] = ch[ik];
01601     for (j=1; j<ip; j++) {
01602       for (k=0; k<l1; k++) {
01603         cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
01604         cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
01605       }
01606     }
01607     if (idot <= l1) {
01608       idij = 0;
01609       for (j=1; j<ip; j++) {
01610         idij += 2;
01611         for (i=3; i<ido; i+=2) {
01612           idij += 2;
01613           for (k=0; k<l1; k++) {
01614             cc[i - 1 + (k + j*l1)*ido] =
01615                 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
01616                 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
01617             cc[i + (k + j*l1)*ido] =
01618                 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
01619                 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
01620           }
01621         }
01622       }
01623     } else {
01624       idj = 2 - ido;
01625       for (j=1; j<ip; j++) {
01626         idj += ido;
01627         for (k = 0; k < l1; k++) {
01628           idij = idj;
01629           for (i=3; i<ido; i+=2) {
01630             idij += 2;
01631             cc[i - 1 + (k + j*l1)*ido] =
01632                 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
01633                 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
01634             cc[i + (k + j*l1)*ido] =
01635                 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
01636                 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
01637           }
01638         }
01639       }
01640     }
01641   } /* passf */
01642 
01643 
01644   /* ----------------------------------------------------------------------
01645 radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
01646 Treal FFT passes fwd and bwd.
01647 ---------------------------------------------------------------------- */
01648 
01649 static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
01650   {
01651     int i, k, ic;
01652     Treal ti2, tr2;
01653     for (k=0; k<l1; k++) {
01654       ch[2*k*ido] =
01655           ref(cc,k*ido) + ref(cc,(k + l1)*ido);
01656       ch[(2*k+1)*ido + ido-1] =
01657           ref(cc,k*ido) - ref(cc,(k + l1)*ido);
01658     }
01659     if (ido < 2) return;
01660     if (ido != 2) {
01661       for (k=0; k<l1; k++) {
01662         for (i=2; i<ido; i+=2) {
01663           ic = ido - i;
01664           tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
01665           ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
01666           ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
01667           ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
01668           ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
01669           ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
01670         }
01671       }
01672       if (ido % 2 == 1) return;
01673     }
01674     for (k=0; k<l1; k++) {
01675       ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
01676       ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
01677     }
01678   } /* radf2 */
01679 
01680 
01681 static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
01682   {
01683     int i, k, ic;
01684     Treal ti2, tr2;
01685     for (k=0; k<l1; k++) {
01686       ch[k*ido] =
01687           ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
01688       ch[(k + l1)*ido] =
01689           ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
01690     }
01691     if (ido < 2) return;
01692     if (ido != 2) {
01693       for (k = 0; k < l1; ++k) {
01694         for (i = 2; i < ido; i += 2) {
01695           ic = ido - i;
01696           ch[i-1 + k*ido] =
01697               ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
01698           tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
01699           ch[i + k*ido] =
01700               ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
01701           ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
01702           ch[i-1 + (k + l1)*ido] =
01703               wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
01704           ch[i + (k + l1)*ido] =
01705               wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
01706         }
01707       }
01708       if (ido % 2 == 1) return;
01709     }
01710     for (k = 0; k < l1; k++) {
01711       ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
01712       ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
01713     }
01714   } /* radb2 */
01715 
01716 
01717 static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
01718       const Treal wa1[], const Treal wa2[])
01719   {
01720     static const Treal taur = -0.5;
01721     static const Treal taui = 0.866025403784439;
01722     int i, k, ic;
01723     Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
01724     for (k=0; k<l1; k++) {
01725       cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
01726       ch[3*k*ido] = ref(cc,k*ido) + cr2;
01727       ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
01728       ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
01729     }
01730     if (ido == 1) return;
01731     for (k=0; k<l1; k++) {
01732       for (i=2; i<ido; i+=2) {
01733         ic = ido - i;
01734         dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
01735             wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01736         di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01737         dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
01738         di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
01739         cr2 = dr2 + dr3;
01740         ci2 = di2 + di3;
01741         ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
01742         ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
01743         tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
01744         ti2 = ref(cc,i + k*ido) + taur*ci2;
01745         tr3 = taui*(di2 - di3);
01746         ti3 = taui*(dr3 - dr2);
01747         ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
01748         ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
01749         ch[i + (3*k + 2)*ido] = ti2 + ti3;
01750         ch[ic + (3*k + 1)*ido] = ti3 - ti2;
01751       }
01752     }
01753   } /* radf3 */
01754 
01755 
01756 static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
01757       const Treal wa1[], const Treal wa2[])
01758   {
01759     static const Treal taur = -0.5;
01760     static const Treal taui = 0.866025403784439;
01761     int i, k, ic;
01762     Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
01763     for (k=0; k<l1; k++) {
01764       tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
01765       cr2 = ref(cc,3*k*ido) + taur*tr2;
01766       ch[k*ido] = ref(cc,3*k*ido) + tr2;
01767       ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
01768       ch[(k + l1)*ido] = cr2 - ci3;
01769       ch[(k + 2*l1)*ido] = cr2 + ci3;
01770     }
01771     if (ido == 1) return;
01772     for (k=0; k<l1; k++) {
01773       for (i=2; i<ido; i+=2) {
01774         ic = ido - i;
01775         tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
01776         cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
01777         ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
01778         ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
01779         ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
01780         ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
01781         cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
01782         ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
01783         dr2 = cr2 - ci3;
01784         dr3 = cr2 + ci3;
01785         di2 = ci2 + cr3;
01786         di3 = ci2 - cr3;
01787         ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
01788         ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
01789         ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
01790         ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
01791       }
01792     }
01793   } /* radb3 */
01794 
01795 
01796 static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
01797       const Treal wa1[], const Treal wa2[], const Treal wa3[])
01798   {
01799     static const Treal hsqt2 = 0.7071067811865475;
01800     int i, k, ic;
01801     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01802     for (k=0; k<l1; k++) {
01803       tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
01804       tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
01805       ch[4*k*ido] = tr1 + tr2;
01806       ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
01807       ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
01808       ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
01809     }
01810     if (ido < 2) return;
01811     if (ido != 2) {
01812       for (k=0; k<l1; k++) {
01813         for (i=2; i<ido; i += 2) {
01814           ic = ido - i;
01815           cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01816           ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01817           cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
01818               ido);
01819           ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
01820               ido);
01821           cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
01822               ido);
01823           ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
01824               ido);
01825           tr1 = cr2 + cr4;
01826           tr4 = cr4 - cr2;
01827           ti1 = ci2 + ci4;
01828           ti4 = ci2 - ci4;
01829           ti2 = ref(cc,i + k*ido) + ci3;
01830           ti3 = ref(cc,i + k*ido) - ci3;
01831           tr2 = ref(cc,i - 1 + k*ido) + cr3;
01832           tr3 = ref(cc,i - 1 + k*ido) - cr3;
01833           ch[i - 1 + 4*k*ido] = tr1 + tr2;
01834           ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
01835           ch[i + 4*k*ido] = ti1 + ti2;
01836           ch[ic + (4*k + 3)*ido] = ti1 - ti2;
01837           ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
01838           ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
01839           ch[i + (4*k + 2)*ido] = tr4 + ti3;
01840           ch[ic + (4*k + 1)*ido] = tr4 - ti3;
01841         }
01842       }
01843       if (ido % 2 == 1) return;
01844     }
01845     for (k=0; k<l1; k++) {
01846       ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
01847       tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
01848       ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
01849       ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
01850       ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
01851       ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
01852     }
01853   } /* radf4 */
01854 
01855 
01856 static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
01857       const Treal wa1[], const Treal wa2[], const Treal wa3[])
01858   {
01859     static const Treal sqrt2 = 1.414213562373095;
01860     int i, k, ic;
01861     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
01862     for (k = 0; k < l1; k++) {
01863       tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
01864       tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
01865       tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
01866       tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
01867       ch[k*ido] = tr2 + tr3;
01868       ch[(k + l1)*ido] = tr1 - tr4;
01869       ch[(k + 2*l1)*ido] = tr2 - tr3;
01870       ch[(k + 3*l1)*ido] = tr1 + tr4;
01871     }
01872     if (ido < 2) return;
01873     if (ido != 2) {
01874       for (k = 0; k < l1; ++k) {
01875         for (i = 2; i < ido; i += 2) {
01876           ic = ido - i;
01877           ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
01878           ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
01879           ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
01880           tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
01881           tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
01882           tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
01883           ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
01884           tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
01885           ch[i - 1 + k*ido] = tr2 + tr3;
01886           cr3 = tr2 - tr3;
01887           ch[i + k*ido] = ti2 + ti3;
01888           ci3 = ti2 - ti3;
01889           cr2 = tr1 - tr4;
01890           cr4 = tr1 + tr4;
01891           ci2 = ti1 + ti4;
01892           ci4 = ti1 - ti4;
01893           ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
01894           ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
01895           ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
01896           ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
01897           ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
01898           ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
01899         }
01900       }
01901       if (ido % 2 == 1) return;
01902     }
01903     for (k = 0; k < l1; k++) {
01904       ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
01905       ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
01906       tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
01907       tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
01908       ch[ido-1 + k*ido] = tr2 + tr2;
01909       ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
01910       ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
01911       ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
01912     }
01913   } /* radb4 */
01914 
01915 
01916 static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
01917       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
01918   {
01919     static const Treal tr11 = 0.309016994374947;
01920     static const Treal ti11 = 0.951056516295154;
01921     static const Treal tr12 = -0.809016994374947;
01922     static const Treal ti12 = 0.587785252292473;
01923     int i, k, ic;
01924     Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
01925         cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
01926     for (k = 0; k < l1; k++) {
01927       cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
01928       ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
01929       cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
01930       ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
01931       ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
01932       ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
01933       ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
01934       ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
01935       ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
01936     }
01937     if (ido == 1) return;
01938     for (k = 0; k < l1; ++k) {
01939       for (i = 2; i < ido; i += 2) {
01940         ic = ido - i;
01941         dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
01942         di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
01943         dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
01944         di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
01945         dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
01946         di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
01947         dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
01948         di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
01949         cr2 = dr2 + dr5;
01950         ci5 = dr5 - dr2;
01951         cr5 = di2 - di5;
01952         ci2 = di2 + di5;
01953         cr3 = dr3 + dr4;
01954         ci4 = dr4 - dr3;
01955         cr4 = di3 - di4;
01956         ci3 = di3 + di4;
01957         ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
01958         ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
01959         tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
01960         ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
01961         tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
01962         ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
01963         tr5 = ti11*cr5 + ti12*cr4;
01964         ti5 = ti11*ci5 + ti12*ci4;
01965         tr4 = ti12*cr5 - ti11*cr4;
01966         ti4 = ti12*ci5 - ti11*ci4;
01967         ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
01968         ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
01969         ch[i + (5*k + 2)*ido] = ti2 + ti5;
01970         ch[ic + (5*k + 1)*ido] = ti5 - ti2;
01971         ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
01972         ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
01973         ch[i + (5*k + 4)*ido] = ti3 + ti4;
01974         ch[ic + (5*k + 3)*ido] = ti4 - ti3;
01975       }
01976     }
01977   } /* radf5 */
01978 
01979 
01980 static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
01981       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
01982   {
01983     static const Treal tr11 = 0.309016994374947;
01984     static const Treal ti11 = 0.951056516295154;
01985     static const Treal tr12 = -0.809016994374947;
01986     static const Treal ti12 = 0.587785252292473;
01987     int i, k, ic;
01988     Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
01989         ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01990     for (k = 0; k < l1; k++) {
01991       ti5 = 2*ref(cc,(5*k + 2)*ido);
01992       ti4 = 2*ref(cc,(5*k + 4)*ido);
01993       tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
01994       tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
01995       ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
01996       cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
01997       cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
01998       ci5 = ti11*ti5 + ti12*ti4;
01999       ci4 = ti12*ti5 - ti11*ti4;
02000       ch[(k + l1)*ido] = cr2 - ci5;
02001       ch[(k + 2*l1)*ido] = cr3 - ci4;
02002       ch[(k + 3*l1)*ido] = cr3 + ci4;
02003       ch[(k + 4*l1)*ido] = cr2 + ci5;
02004     }
02005     if (ido == 1) return;
02006     for (k = 0; k < l1; ++k) {
02007       for (i = 2; i < ido; i += 2) {
02008         ic = ido - i;
02009         ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
02010         ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
02011         ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
02012         ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
02013         tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
02014         tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
02015         tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
02016         tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
02017         ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
02018         ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
02019         cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;
02020 
02021         ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
02022         cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;
02023 
02024         ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
02025         cr5 = ti11*tr5 + ti12*tr4;
02026         ci5 = ti11*ti5 + ti12*ti4;
02027         cr4 = ti12*tr5 - ti11*tr4;
02028         ci4 = ti12*ti5 - ti11*ti4;
02029         dr3 = cr3 - ci4;
02030         dr4 = cr3 + ci4;
02031         di3 = ci3 + cr4;
02032         di4 = ci3 - cr4;
02033         dr5 = cr2 + ci5;
02034         dr2 = cr2 - ci5;
02035         di5 = ci2 - cr5;
02036         di2 = ci2 + cr5;
02037         ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
02038         ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
02039         ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
02040         ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
02041         ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
02042         ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
02043         ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
02044         ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
02045       }
02046     }
02047   } /* radb5 */
02048 
02049 
02050 static void radfg(int ido, int ip, int l1, int idl1,
02051       Treal cc[], Treal ch[], const Treal wa[])
02052   {
02053     static const Treal twopi = 6.28318530717959;
02054     int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
02055     Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
02056     arg = twopi / ip;
02057     dcp = cos(arg);
02058     dsp = sin(arg);
02059     ipph = (ip + 1) / 2;
02060     nbd = (ido - 1) / 2;
02061     if (ido != 1) {
02062       for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
02063       for (j=1; j<ip; j++)
02064         for (k=0; k<l1; k++)
02065           ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
02066       if (nbd <= l1) {
02067         is = -ido;
02068         for (j=1; j<ip; j++) {
02069           is += ido;
02070           idij = is-1;
02071           for (i=2; i<ido; i+=2) {
02072             idij += 2;
02073             for (k=0; k<l1; k++) {
02074               ch[i - 1 + (k + j*l1)*ido] =
02075                   wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
02076               ch[i + (k + j*l1)*ido] =
02077                   wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
02078             }
02079           }
02080         }
02081       } else {
02082         is = -ido;
02083         for (j=1; j<ip; j++) {
02084           is += ido;
02085           for (k=0; k<l1; k++) {
02086             idij = is-1;
02087             for (i=2; i<ido; i+=2) {
02088               idij += 2;
02089               ch[i - 1 + (k + j*l1)*ido] =
02090                   wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
02091               ch[i + (k + j*l1)*ido] =
02092                   wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
02093             }
02094           }
02095         }
02096       }
02097       if (nbd >= l1) {
02098         for (j=1; j<ipph; j++) {
02099           jc = ip - j;
02100           for (k=0; k<l1; k++) {
02101             for (i=2; i<ido; i+=2) {
02102               cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02103               cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
02104               cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02105               cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
02106             }
02107           }
02108         }
02109       } else {
02110         for (j=1; j<ipph; j++) {
02111           jc = ip - j;
02112           for (i=2; i<ido; i+=2) {
02113             for (k=0; k<l1; k++) {
02114               cc[i - 1 + (k + j*l1)*ido] =
02115                   ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02116               cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
02117               cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02118               cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
02119             }
02120           }
02121         }
02122       }
02123     } else {  /* now ido == 1 */
02124       for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
02125     }
02126     for (j=1; j<ipph; j++) {
02127       jc = ip - j;
02128       for (k=0; k<l1; k++) {
02129         cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
02130         cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
02131       }
02132     }
02133 
02134     ar1 = 1;
02135     ai1 = 0;
02136     for (l=1; l<ipph; l++) {
02137       lc = ip - l;
02138       ar1h = dcp*ar1 - dsp*ai1;
02139       ai1 = dcp*ai1 + dsp*ar1;
02140       ar1 = ar1h;
02141       for (ik=0; ik<idl1; ik++) {
02142         ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
02143         ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
02144       }
02145       dc2 = ar1;
02146       ds2 = ai1;
02147       ar2 = ar1;
02148       ai2 = ai1;
02149       for (j=2; j<ipph; j++) {
02150         jc = ip - j;
02151         ar2h = dc2*ar2 - ds2*ai2;
02152         ai2 = dc2*ai2 + ds2*ar2;
02153         ar2 = ar2h;
02154         for (ik=0; ik<idl1; ik++) {
02155           ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
02156           ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
02157         }
02158       }
02159     }
02160     for (j=1; j<ipph; j++)
02161       for (ik=0; ik<idl1; ik++)
02162         ch[ik] += cc[ik + j*idl1];
02163 
02164     if (ido >= l1) {
02165       for (k=0; k<l1; k++) {
02166         for (i=0; i<ido; i++) {
02167           ref(cc,i + k*ip*ido) = ch[i + k*ido];
02168         }
02169       }
02170     } else {
02171       for (i=0; i<ido; i++) {
02172         for (k=0; k<l1; k++) {
02173           ref(cc,i + k*ip*ido) = ch[i + k*ido];
02174         }
02175       }
02176     }
02177     for (j=1; j<ipph; j++) {
02178       jc = ip - j;
02179       j2 = 2*j;
02180       for (k=0; k<l1; k++) {
02181         ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
02182             ch[(k + j*l1)*ido];
02183         ref(cc,(j2 + k*ip)*ido) =
02184             ch[(k + jc*l1)*ido];
02185       }
02186     }
02187     if (ido == 1) return;
02188     if (nbd >= l1) {
02189       for (j=1; j<ipph; j++) {
02190         jc = ip - j;
02191         j2 = 2*j;
02192         for (k=0; k<l1; k++) {
02193           for (i=2; i<ido; i+=2) {
02194             ic = ido - i;
02195             ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02196             ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
02197             ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02198             ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
02199           }
02200         }
02201       }
02202     } else {
02203       for (j=1; j<ipph; j++) {
02204         jc = ip - j;
02205         j2 = 2*j;
02206         for (i=2; i<ido; i+=2) {
02207           ic = ido - i;
02208           for (k=0; k<l1; k++) {
02209             ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
02210             ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
02211             ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
02212             ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
02213           }
02214         }
02215       }
02216     }
02217   } /* radfg */
02218 
02219 
02220 static void radbg(int ido, int ip, int l1, int idl1,
02221       Treal cc[], Treal ch[], const Treal wa[])
02222   {
02223     static const Treal twopi = 6.28318530717959;
02224     int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
02225     Treal dc2, ai1, ai2, ar1, ar2, ds2;
02226     int nbd;
02227     Treal dcp, arg, dsp, ar1h, ar2h;
02228     arg = twopi / ip;
02229     dcp = cos(arg);
02230     dsp = sin(arg);
02231     nbd = (ido - 1) / 2;
02232     ipph = (ip + 1) / 2;
02233     if (ido >= l1) {
02234       for (k=0; k<l1; k++) {
02235         for (i=0; i<ido; i++) {
02236           ch[i + k*ido] = ref(cc,i + k*ip*ido);
02237         }
02238       }
02239     } else {
02240       for (i=0; i<ido; i++) {
02241         for (k=0; k<l1; k++) {
02242           ch[i + k*ido] = ref(cc,i + k*ip*ido);
02243         }
02244       }
02245     }
02246     for (j=1; j<ipph; j++) {
02247       jc = ip - j;
02248       j2 = 2*j;
02249       for (k=0; k<l1; k++) {
02250         ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
02251             ido);
02252         ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
02253       }
02254     }
02255 
02256     if (ido != 1) {
02257       if (nbd >= l1) {
02258         for (j=1; j<ipph; j++) {
02259           jc = ip - j;
02260           for (k=0; k<l1; k++) {
02261             for (i=2; i<ido; i+=2) {
02262               ic = ido - i;
02263               ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
02264                   ic - 1 + (2*j - 1 + k*ip)*ido);
02265               ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
02266                   ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
02267               ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
02268                   + (2*j - 1 + k*ip)*ido);
02269               ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
02270                   + (2*j - 1 + k*ip)*ido);
02271             }
02272           }
02273         }
02274       } else {
02275         for (j=1; j<ipph; j++) {
02276           jc = ip - j;
02277           for (i=2; i<ido; i+=2) {
02278             ic = ido - i;
02279             for (k=0; k<l1; k++) {
02280               ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
02281                   ic - 1 + (2*j - 1 + k*ip)*ido);
02282               ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
02283                   ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
02284               ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
02285                   + (2*j - 1 + k*ip)*ido);
02286               ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
02287                   + (2*j - 1 + k*ip)*ido);
02288             }
02289           }
02290         }
02291       }
02292     }
02293 
02294     ar1 = 1;
02295     ai1 = 0;
02296     for (l=1; l<ipph; l++) {
02297       lc = ip - l;
02298       ar1h = dcp*ar1 - dsp*ai1;
02299       ai1 = dcp*ai1 + dsp*ar1;
02300       ar1 = ar1h;
02301       for (ik=0; ik<idl1; ik++) {
02302         cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
02303         cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
02304       }
02305       dc2 = ar1;
02306       ds2 = ai1;
02307       ar2 = ar1;
02308       ai2 = ai1;
02309       for (j=2; j<ipph; j++) {
02310         jc = ip - j;
02311         ar2h = dc2*ar2 - ds2*ai2;
02312         ai2 = dc2*ai2 + ds2*ar2;
02313         ar2 = ar2h;
02314         for (ik=0; ik<idl1; ik++) {
02315           cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
02316           cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
02317         }
02318       }
02319     }
02320     for (j=1; j<ipph; j++) {
02321       for (ik=0; ik<idl1; ik++) {
02322         ch[ik] += ch[ik + j*idl1];
02323       }
02324     }
02325     for (j=1; j<ipph; j++) {
02326       jc = ip - j;
02327       for (k=0; k<l1; k++) {
02328         ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
02329         ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
02330       }
02331     }
02332 
02333     if (ido == 1) return;
02334     if (nbd >= l1) {
02335       for (j=1; j<ipph; j++) {
02336         jc = ip - j;
02337         for (k=0; k<l1; k++) {
02338           for (i=2; i<ido; i+=2) {
02339             ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
02340             ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
02341             ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
02342             ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
02343           }
02344         }
02345       }
02346     } else {
02347       for (j=1; j<ipph; j++) {
02348         jc = ip - j;
02349         for (i=2; i<ido; i+=2) {
02350           for (k=0; k<l1; k++) {
02351             ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
02352             ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
02353             ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
02354             ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
02355           }
02356         }
02357       }
02358     }
02359     for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
02360     for (j=1; j<ip; j++)
02361       for (k=0; k<l1; k++)
02362         cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
02363     if (nbd <= l1) {
02364       is = -ido;
02365       for (j=1; j<ip; j++) {
02366         is += ido;
02367         idij = is-1;
02368         for (i=2; i<ido; i+=2) {
02369           idij += 2;
02370           for (k=0; k<l1; k++) {
02371             cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
02372                 ch[i + (k + j*l1)*ido];
02373             cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
02374           }
02375         }
02376       }
02377     } else {
02378       is = -ido;
02379       for (j=1; j<ip; j++) {
02380         is += ido;
02381         for (k=0; k<l1; k++) {
02382           idij = is - 1;
02383           for (i=2; i<ido; i+=2) {
02384             idij += 2;
02385             cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
02386                 ch[i + (k + j*l1)*ido];
02387             cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
02388           }
02389         }
02390       }
02391     }
02392   } /* radbg */
02393 
02394   /* ----------------------------------------------------------------------
02395 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
02396 ---------------------------------------------------------------------- */
02397 
02398 static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
02399   {
02400     int idot, i;
02401     int k1, l1, l2;
02402     int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
02403     Treal *cinput, *coutput;
02404     nf = ifac[1];
02405     na = 0;
02406     l1 = 1;
02407     iw = 0;
02408     for (k1=2; k1<=nf+1; k1++) {
02409       ip = ifac[k1];
02410       l2 = ip*l1;
02411       ido = n / l2;
02412       idot = ido + ido;
02413       idl1 = idot*l1;
02414       if (na) {
02415         cinput = ch;
02416         coutput = c;
02417       } else {
02418         cinput = c;
02419         coutput = ch;
02420       }
02421       switch (ip) {
02422       case 4:
02423         ix2 = iw + idot;
02424         ix3 = ix2 + idot;
02425         passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
02426         na = !na;
02427         break;
02428       case 2:
02429         passf2(idot, l1, cinput, coutput, &wa[iw], isign);
02430         na = !na;
02431         break;
02432       case 3:
02433         ix2 = iw + idot;
02434         passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
02435         na = !na;
02436         break;
02437       case 5:
02438         ix2 = iw + idot;
02439         ix3 = ix2 + idot;
02440         ix4 = ix3 + idot;
02441         passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
02442         na = !na;
02443         break;
02444       default:
02445         passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
02446         if (nac != 0) na = !na;
02447       }
02448       l1 = l2;
02449       iw += (ip - 1)*idot;
02450     }
02451     if (na == 0) return;
02452     for (i=0; i<2*n; i++) c[i] = ch[i];
02453   } /* cfftf1 */
02454 
02455 
02456 void cfftf(int n, Treal c[], Treal wsave[])
02457   {
02458     int iw1, iw2;
02459     if (n == 1) return;
02460     iw1 = 2*n;
02461     iw2 = iw1 + 2*n;
02462     cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
02463   } /* cfftf */
02464 
02465 
02466 void cfftb(int n, Treal c[], Treal wsave[])
02467   {
02468     int iw1, iw2;
02469     if (n == 1) return;
02470     iw1 = 2*n;
02471     iw2 = iw1 + 2*n;
02472     cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
02473   } /* cfftb */
02474 
02475 
02476 static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
02477   /* Factorize n in factors in ntryh and rest. On exit,
02478 ifac[0] contains n and ifac[1] contains number of factors,
02479 the factors start from ifac[2]. */
02480   {
02481     int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
02482 startloop:
02483     if (j < NSPECIAL)
02484       ntry = ntryh[j];
02485     else
02486       ntry+= 2;
02487     j++;
02488     do {
02489       nq = nl / ntry;
02490       nr = nl - ntry*nq;
02491       if (nr != 0) goto startloop;
02492       nf++;
02493       ifac[nf + 1] = ntry;
02494       nl = nq;
02495       if (ntry == 2 && nf != 1) {
02496         for (i=2; i<=nf; i++) {
02497           ib = nf - i + 2;
02498           ifac[ib + 1] = ifac[ib];
02499         }
02500         ifac[2] = 2;
02501       }
02502     } while (nl != 1);
02503     ifac[0] = n;
02504     ifac[1] = nf;
02505   }
02506 
02507 
02508 static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
02509   {
02510     static const Treal twopi = 6.28318530717959;
02511     Treal arg, argh, argld, fi;
02512     int idot, i, j;
02513     int i1, k1, l1, l2;
02514     int ld, ii, nf, ip;
02515     int ido, ipm;
02516 
02517     static const int ntryh[NSPECIAL] = {
02518       3,4,2,5    }; /* Do not change the order of these. */
02519 
02520     factorize(n,ifac,ntryh);
02521     nf = ifac[1];
02522     argh = twopi/(Treal)n;
02523     i = 1;
02524     l1 = 1;
02525     for (k1=1; k1<=nf; k1++) {
02526       ip = ifac[k1+1];
02527       ld = 0;
02528       l2 = l1*ip;
02529       ido = n / l2;
02530       idot = ido + ido + 2;
02531       ipm = ip - 1;
02532       for (j=1; j<=ipm; j++) {
02533         i1 = i;
02534         wa[i-1] = 1;
02535         wa[i] = 0;
02536         ld += l1;
02537         fi = 0;
02538         argld = ld*argh;
02539         for (ii=4; ii<=idot; ii+=2) {
02540           i+= 2;
02541           fi+= 1;
02542           arg = fi*argld;
02543           wa[i-1] = cos(arg);
02544           wa[i] = sin(arg);
02545         }
02546         if (ip > 5) {
02547           wa[i1-1] = wa[i-1];
02548           wa[i1] = wa[i];
02549         }
02550       }
02551       l1 = l2;
02552     }
02553   } /* cffti1 */
02554 
02555 
02556 void cffti(int n, Treal wsave[])
02557  {
02558     int iw1, iw2;
02559     if (n == 1) return;
02560     iw1 = 2*n;
02561     iw2 = iw1 + 2*n;
02562     cffti1(n, wsave+iw1, (int*)(wsave+iw2));
02563   } /* cffti */
02564 
02565   /* ----------------------------------------------------------------------
02566 rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
02567 ---------------------------------------------------------------------- */
02568 
02569 static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
02570   {
02571     int i;
02572     int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
02573     Treal *cinput, *coutput;
02574     nf = ifac[1];
02575     na = 1;
02576     l2 = n;
02577     iw = n-1;
02578     for (k1 = 1; k1 <= nf; ++k1) {
02579       kh = nf - k1;
02580       ip = ifac[kh + 2];
02581       l1 = l2 / ip;
02582       ido = n / l2;
02583       idl1 = ido*l1;
02584       iw -= (ip - 1)*ido;
02585       na = !na;
02586       if (na) {
02587         cinput = ch;
02588         coutput = c;
02589       } else {
02590         cinput = c;
02591         coutput = ch;
02592       }
02593       switch (ip) {
02594       case 4:
02595         ix2 = iw + ido;
02596         ix3 = ix2 + ido;
02597         radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
02598         break;
02599       case 2:
02600         radf2(ido, l1, cinput, coutput, &wa[iw]);
02601         break;
02602       case 3:
02603         ix2 = iw + ido;
02604         radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
02605         break;
02606       case 5:
02607         ix2 = iw + ido;
02608         ix3 = ix2 + ido;
02609         ix4 = ix3 + ido;
02610         radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
02611         break;
02612       default:
02613         if (ido == 1)
02614           na = !na;
02615         if (na == 0) {
02616           radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
02617           na = 1;
02618         } else {
02619           radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
02620           na = 0;
02621         }
02622       }
02623       l2 = l1;
02624     }
02625     if (na == 1) return;
02626     for (i = 0; i < n; i++) c[i] = ch[i];
02627   } /* rfftf1 */
02628 
02629 
02630 void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
02631   {
02632     int i;
02633     int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
02634     Treal *cinput, *coutput;
02635     nf = ifac[1];
02636     na = 0;
02637     l1 = 1;
02638     iw = 0;
02639     for (k1=1; k1<=nf; k1++) {
02640       ip = ifac[k1 + 1];
02641       l2 = ip*l1;
02642       ido = n / l2;
02643       idl1 = ido*l1;
02644       if (na) {
02645         cinput = ch;
02646         coutput = c;
02647       } else {
02648         cinput = c;
02649         coutput = ch;
02650       }
02651       switch (ip) {
02652       case 4:
02653         ix2 = iw + ido;
02654         ix3 = ix2 + ido;
02655         radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
02656         na = !na;
02657         break;
02658       case 2:
02659         radb2(ido, l1, cinput, coutput, &wa[iw]);
02660         na = !na;
02661         break;
02662       case 3:
02663         ix2 = iw + ido;
02664         radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
02665         na = !na;
02666         break;
02667       case 5:
02668         ix2 = iw + ido;
02669         ix3 = ix2 + ido;
02670         ix4 = ix3 + ido;
02671         radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
02672         na = !na;
02673         break;
02674       default:
02675         radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
02676         if (ido == 1) na = !na;
02677       }
02678       l1 = l2;
02679       iw += (ip - 1)*ido;
02680     }
02681     if (na == 0) return;
02682     for (i=0; i<n; i++) c[i] = ch[i];
02683   } /* rfftb1 */
02684 
02685 
02686 void EMAN::rfftf(int n, Treal r[], Treal wsave[])
02687   {
02688     if (n == 1) return;
02689     rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
02690   } /* rfftf */
02691 
02692 
02693 void EMAN::rfftb(int n, Treal r[], Treal wsave[])
02694   {
02695     if (n == 1) return;
02696     rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
02697   } /* rfftb */
02698 
02699 
02700 static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
02701   {
02702     static const Treal twopi = 6.28318530717959;
02703     Treal arg, argh, argld, fi;
02704     int i, j;
02705     int k1, l1, l2;
02706     int ld, ii, nf, ip, is;
02707     int ido, ipm, nfm1;
02708     static const int ntryh[NSPECIAL] = {
02709       4,2,3,5    }; /* Do not change the order of these. */
02710     factorize(n,ifac,ntryh);
02711     nf = ifac[1];
02712     argh = twopi / n;
02713     is = 0;
02714     nfm1 = nf - 1;
02715     l1 = 1;
02716     if (nfm1 == 0) return;
02717     for (k1 = 1; k1 <= nfm1; k1++) {
02718       ip = ifac[k1 + 1];
02719       ld = 0;
02720       l2 = l1*ip;
02721       ido = n / l2;
02722       ipm = ip - 1;
02723       for (j = 1; j <= ipm; ++j) {
02724         ld += l1;
02725         i = is;
02726         argld = (Treal) ld*argh;
02727         fi = 0;
02728         for (ii = 3; ii <= ido; ii += 2) {
02729           i += 2;
02730           fi += 1;
02731           arg = fi*argld;
02732           wa[i - 2] = cos(arg);
02733           wa[i - 1] = sin(arg);
02734         }
02735         is += ido;
02736       }
02737       l1 = l2;
02738     }
02739   } /* rffti1 */
02740 
02741 
02742 void EMAN::rffti(int n, Treal wsave[])
02743   {
02744         if (n == 1) return;
02745         rffti1(n, wsave+n, (int*)(wsave+2*n));
02746   } /* rffti */
02747 
02748 //#ifdef __cplusplus
02749 //}
02750 //#endif
02751 
02752 #endif  //NATIVE_FFT