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

Generated on Thu May 3 10:06:25 2012 for EMAN2 by  doxygen 1.4.7