fftpack.c

00001 /*
00002  *  This file is part of libfftpack.
00003  *
00004  *  libfftpack is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  libfftpack is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with libfftpack; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  */
00018 
00019 /*
00020  *  libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
00021  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00022  *  (DLR).
00023  */
00024 
00025 /*
00026  * This file was originally part of tela the Tensor Language.
00027  * Copyright (C) 1994-1995 Pekka Janhunen
00028  */
00029 
00030 /*
00031   fftpack.c : A set of FFT routines in C.
00032   Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
00033   (Version 4, 1985).
00034 
00035   Pekka Janhunen 23.2.1995
00036 
00037   (reformatted by joerg arndt)
00038 
00039   reformatted and slightly enhanced by Martin Reinecke (2004)
00040  */
00041 
00042 #include <math.h>
00043 #include <stdlib.h>
00044 #include <string.h>
00045 #include "fftpack.h"
00046 
00047 static void passf2(int ido, int l1, const double *cc, double *ch,
00048   const double *wa1)
00049   {
00050   int i, k, ah, ac;
00051   double ti2, tr2;
00052 
00053   if(ido<=2)
00054     {
00055     for(k=0; k<l1; k++)
00056       {
00057       ah=k*ido;
00058       ac=2*k*ido;
00059       ch[ah]=cc[ac]+cc[ac+ido];
00060       ch[ah+ido*l1]=cc[ac]-cc[ac+ido];
00061       ch[ah+1]=cc[ac+1]+cc[ac+ido+1];
00062       ch[ah+ido*l1+1]=cc[ac+1]-cc[ac+ido+1];
00063       }
00064     }
00065   else
00066     {
00067     for(k=0; k<l1; k++)
00068       {
00069       for(i=0; i<ido-1; i+=2)
00070         {
00071         ah=i+k*ido;
00072         ac=i+2*k*ido;
00073         ch[ah]=cc[ac]+cc[ac+ido];
00074         tr2=cc[ac]-cc[ac+ido];
00075         ch[ah+1]=cc[ac+1]+cc[ac+1+ido];
00076         ti2=cc[ac+1]-cc[ac+1+ido];
00077         ch[ah+l1*ido+1]=wa1[i]*ti2-wa1[i+1]*tr2;
00078         ch[ah+l1*ido]=wa1[i]*tr2+wa1[i+1]*ti2;
00079         }
00080       }
00081     }
00082   }
00083 
00084 static void passb2(int ido, int l1, const double *cc, double *ch,
00085   const double *wa1)
00086   {
00087   int i, k, ah, ac;
00088   double ti2, tr2;
00089 
00090   if(ido<=2)
00091     {
00092     for(k=0; k<l1; k++)
00093       {
00094       ah=k*ido;
00095       ac=2*k*ido;
00096       ch[ah]=cc[ac]+cc[ac+ido];
00097       ch[ah+ido*l1]=cc[ac]-cc[ac+ido];
00098       ch[ah+1]=cc[ac+1]+cc[ac+ido+1];
00099       ch[ah+ido*l1+1]=cc[ac+1]-cc[ac+ido+1];
00100       }
00101     }
00102   else
00103     {
00104     for(k=0; k<l1; k++)
00105       {
00106       for(i=0; i<ido-1; i+=2)
00107         {
00108         ah=i+k*ido;
00109         ac=i+2*k*ido;
00110         ch[ah]=cc[ac]+cc[ac+ido];
00111         tr2=cc[ac]-cc[ac+ido];
00112         ch[ah+1]=cc[ac+1]+cc[ac+1+ido];
00113         ti2=cc[ac+1]-cc[ac+1+ido];
00114         ch[ah+l1*ido+1]=wa1[i]*ti2+wa1[i+1]*tr2;
00115         ch[ah+l1*ido]=wa1[i]*tr2-wa1[i+1]*ti2;
00116         }
00117       }
00118     }
00119   }
00120 
00121 static void passf3(int ido, int l1, const double *cc, double *ch,
00122   const double *wa1, const double *wa2)
00123   {
00124   static const double taur=-0.5, taui=0.86602540378443864676;
00125   int i, k, ac, ah;
00126   double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
00127 
00128   if(ido==2)
00129     {
00130     for(k=1; k<=l1; k++)
00131       {
00132       ac=(3*k-2)*ido;
00133       tr2=cc[ac]+cc[ac+ido];
00134       cr2=cc[ac-ido]+taur*tr2;
00135       ah=(k-1)*ido;
00136       ch[ah]=cc[ac-ido]+tr2;
00137 
00138       ti2=cc[ac+1]+cc[ac+ido+1];
00139       ci2=cc[ac-ido+1]+taur*ti2;
00140       ch[ah+1]=cc[ac-ido+1]+ti2;
00141 
00142       cr3=-taui*(cc[ac]-cc[ac+ido]);
00143       ci3=-taui*(cc[ac+1]-cc[ac+ido+1]);
00144       ch[ah+l1*ido]=cr2-ci3;
00145       ch[ah+2*l1*ido]=cr2+ci3;
00146       ch[ah+l1*ido+1]=ci2+cr3;
00147       ch[ah+2*l1*ido+1]=ci2-cr3;
00148       }
00149     }
00150   else
00151     {
00152     for(k=1; k<=l1; k++)
00153       {
00154       for(i=0; i<ido-1; i+=2)
00155         {
00156         ac=i+(3*k-2)*ido;
00157         tr2=cc[ac]+cc[ac+ido];
00158         cr2=cc[ac-ido]+taur*tr2;
00159         ah=i+(k-1)*ido;
00160         ch[ah]=cc[ac-ido]+tr2;
00161         ti2=cc[ac+1]+cc[ac+ido+1];
00162         ci2=cc[ac-ido+1]+taur*ti2;
00163         ch[ah+1]=cc[ac-ido+1]+ti2;
00164         cr3=-taui*(cc[ac]-cc[ac+ido]);
00165         ci3=-taui*(cc[ac+1]-cc[ac+ido+1]);
00166         dr2=cr2-ci3;
00167         dr3=cr2+ci3;
00168         di2=ci2+cr3;
00169         di3=ci2-cr3;
00170         ch[ah+l1*ido+1]=wa1[i]*di2-wa1[i+1]*dr2;
00171         ch[ah+l1*ido]=wa1[i]*dr2+wa1[i+1]*di2;
00172         ch[ah+2*l1*ido+1]=wa2[i]*di3-wa2[i+1]*dr3;
00173         ch[ah+2*l1*ido]=wa2[i]*dr3+wa2[i+1]*di3;
00174         }
00175       }
00176     }
00177   }
00178 
00179 static void passb3(int ido, int l1, const double *cc, double *ch,
00180   const double *wa1, const double *wa2)
00181   {
00182   static const double taur=-0.5, taui=0.86602540378443864676;
00183   int i, k, ac, ah;
00184   double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
00185 
00186   if(ido==2)
00187     {
00188     for(k=1; k<=l1; k++)
00189       {
00190       ac=(3*k-2)*ido;
00191       tr2=cc[ac]+cc[ac+ido];
00192       cr2=cc[ac-ido]+taur*tr2;
00193       ah=(k-1)*ido;
00194       ch[ah]=cc[ac-ido]+tr2;
00195 
00196       ti2=cc[ac+1]+cc[ac+ido+1];
00197       ci2=cc[ac-ido+1]+taur*ti2;
00198       ch[ah+1]=cc[ac-ido+1]+ti2;
00199 
00200       cr3=taui*(cc[ac]-cc[ac+ido]);
00201       ci3=taui*(cc[ac+1]-cc[ac+ido+1]);
00202       ch[ah+l1*ido]=cr2-ci3;
00203       ch[ah+2*l1*ido]=cr2+ci3;
00204       ch[ah+l1*ido+1]=ci2+cr3;
00205       ch[ah+2*l1*ido+1]=ci2-cr3;
00206       }
00207     }
00208   else
00209     {
00210     for(k=1; k<=l1; k++)
00211       {
00212       for(i=0; i<ido-1; i+=2)
00213         {
00214         ac=i+(3*k-2)*ido;
00215         tr2=cc[ac]+cc[ac+ido];
00216         cr2=cc[ac-ido]+taur*tr2;
00217         ah=i+(k-1)*ido;
00218         ch[ah]=cc[ac-ido]+tr2;
00219         ti2=cc[ac+1]+cc[ac+ido+1];
00220         ci2=cc[ac-ido+1]+taur*ti2;
00221         ch[ah+1]=cc[ac-ido+1]+ti2;
00222         cr3=taui*(cc[ac]-cc[ac+ido]);
00223         ci3=taui*(cc[ac+1]-cc[ac+ido+1]);
00224         dr2=cr2-ci3;
00225         dr3=cr2+ci3;
00226         di2=ci2+cr3;
00227         di3=ci2-cr3;
00228         ch[ah+l1*ido+1]=wa1[i]*di2+wa1[i+1]*dr2;
00229         ch[ah+l1*ido]=wa1[i]*dr2-wa1[i+1]*di2;
00230         ch[ah+2*l1*ido+1]=wa2[i]*di3+wa2[i+1]*dr3;
00231         ch[ah+2*l1*ido]=wa2[i]*dr3-wa2[i+1]*di3;
00232         }
00233       }
00234     }
00235   }
00236 
00237 static void passf4(int ido, int l1, const double *cc, double *ch,
00238   const double *wa1, const double *wa2, const double *wa3)
00239   {
00240   int i, k, ac, ah;
00241   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
00242         
00243   if(ido==2)
00244     {
00245     for(k=0; k<l1; k++)
00246       {
00247       ac=4*k*ido+1;
00248       ti1=cc[ac]-cc[ac+2*ido];
00249       ti2=cc[ac]+cc[ac+2*ido];
00250       tr4=cc[ac+3*ido]-cc[ac+ido];
00251       ti3=cc[ac+ido]+cc[ac+3*ido];
00252       tr1=cc[ac-1]-cc[ac+2*ido-1];
00253       tr2=cc[ac-1]+cc[ac+2*ido-1];
00254       ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
00255       tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
00256       ah=k*ido;
00257       ch[ah]=tr2+tr3;
00258       ch[ah+2*l1*ido]=tr2-tr3;
00259       ch[ah+1]=ti2+ti3;
00260       ch[ah+2*l1*ido+1]=ti2-ti3;
00261       ch[ah+l1*ido]=tr1-tr4;
00262       ch[ah+3*l1*ido]=tr1+tr4;
00263       ch[ah+l1*ido+1]=ti1-ti4;
00264       ch[ah+3*l1*ido+1]=ti1+ti4;
00265       }
00266     }
00267   else
00268     {
00269     for(k=0; k<l1; k++)
00270       {
00271       for(i=0; i<ido-1; i+=2)
00272         {
00273         ac=i+1+4*k*ido;
00274         ti1=cc[ac]-cc[ac+2*ido];
00275         ti2=cc[ac]+cc[ac+2*ido];
00276         ti3=cc[ac+ido]+cc[ac+3*ido];
00277         tr4=cc[ac+3*ido]-cc[ac+ido];
00278         tr1=cc[ac-1]-cc[ac+2*ido-1];
00279         tr2=cc[ac-1]+cc[ac+2*ido-1];
00280         ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
00281         tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
00282         ah=i+k*ido;
00283         ch[ah]=tr2+tr3;
00284         cr3=tr2-tr3;
00285         ch[ah+1]=ti2+ti3;
00286         ci3=ti2-ti3;
00287         cr2=tr1-tr4;
00288         cr4=tr1+tr4;
00289         ci2=ti1-ti4;
00290         ci4=ti1+ti4;
00291         ch[ah+l1*ido]=wa1[i]*cr2+wa1[i+1]*ci2;
00292         ch[ah+l1*ido+1]=wa1[i]*ci2-wa1[i+1]*cr2;
00293         ch[ah+2*l1*ido]=wa2[i]*cr3+wa2[i+1]*ci3;
00294         ch[ah+2*l1*ido+1]=wa2[i]*ci3-wa2[i+1]*cr3;
00295         ch[ah+3*l1*ido]=wa3[i]*cr4+wa3[i+1]*ci4;
00296         ch[ah+3*l1*ido+1]=wa3[i]*ci4-wa3[i+1]*cr4;
00297         }
00298       }
00299     }
00300   }
00301 
00302 static void passb4(int ido, int l1, const double *cc, double *ch,
00303   const double *wa1, const double *wa2, const double *wa3)
00304   {
00305   int i, k, ac, ah;
00306   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
00307         
00308   if(ido==2)
00309     {
00310     for(k=0; k<l1; k++)
00311       {
00312       ac=4*k*ido+1;
00313       ti1=cc[ac]-cc[ac+2*ido];
00314       ti2=cc[ac]+cc[ac+2*ido];
00315       tr4=cc[ac+3*ido]-cc[ac+ido];
00316       ti3=cc[ac+ido]+cc[ac+3*ido];
00317       tr1=cc[ac-1]-cc[ac+2*ido-1];
00318       tr2=cc[ac-1]+cc[ac+2*ido-1];
00319       ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
00320       tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
00321       ah=k*ido;
00322       ch[ah]=tr2+tr3;
00323       ch[ah+2*l1*ido]=tr2-tr3;
00324       ch[ah+1]=ti2+ti3;
00325       ch[ah+2*l1*ido+1]=ti2-ti3;
00326       ch[ah+l1*ido]=tr1+tr4;
00327       ch[ah+3*l1*ido]=tr1-tr4;
00328       ch[ah+l1*ido+1]=ti1+ti4;
00329       ch[ah+3*l1*ido+1]=ti1-ti4;
00330       }
00331     }
00332   else
00333     {
00334     for(k=0; k<l1; k++)
00335       {
00336       for(i=0; i<ido-1; i+=2)
00337         {
00338         ac=i+1+4*k*ido;
00339         ti1=cc[ac]-cc[ac+2*ido];
00340         ti2=cc[ac]+cc[ac+2*ido];
00341         ti3=cc[ac+ido]+cc[ac+3*ido];
00342         tr4=cc[ac+3*ido]-cc[ac+ido];
00343         tr1=cc[ac-1]-cc[ac+2*ido-1];
00344         tr2=cc[ac-1]+cc[ac+2*ido-1];
00345         ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
00346         tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
00347         ah=i+k*ido;
00348         ch[ah]=tr2+tr3;
00349         cr3=tr2-tr3;
00350         ch[ah+1]=ti2+ti3;
00351         ci3=ti2-ti3;
00352         cr2=tr1+tr4;
00353         cr4=tr1-tr4;
00354         ci2=ti1+ti4;
00355         ci4=ti1-ti4;
00356         ch[ah+l1*ido]=wa1[i]*cr2-wa1[i+1]*ci2;
00357         ch[ah+l1*ido+1]=wa1[i]*ci2+wa1[i+1]*cr2;
00358         ch[ah+2*l1*ido]=wa2[i]*cr3-wa2[i+1]*ci3;
00359         ch[ah+2*l1*ido+1]=wa2[i]*ci3+wa2[i+1]*cr3;
00360         ch[ah+3*l1*ido]=wa3[i]*cr4-wa3[i+1]*ci4;
00361         ch[ah+3*l1*ido+1]=wa3[i]*ci4+wa3[i+1]*cr4;
00362         }
00363       }
00364     }
00365   }
00366 
00367 static void passf5(int ido, int l1, const double *cc, double *ch,
00368   const double *wa1, const double *wa2, const double *wa3,
00369   const double *wa4)
00370   {
00371   static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212;
00372   static const double tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
00373   int i, k, ac, ah;
00374   double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
00375          ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
00376 
00377   if(ido==2)
00378     {
00379     for(k=1; k<=l1;++k)
00380       {
00381       ac=(5*k-4)*ido+1;
00382       ti5=cc[ac]-cc[ac+3*ido];
00383       ti2=cc[ac]+cc[ac+3*ido];
00384       ti4=cc[ac+ido]-cc[ac+2*ido];
00385       ti3=cc[ac+ido]+cc[ac+2*ido];
00386       tr5=cc[ac-1]-cc[ac+3*ido-1];
00387       tr2=cc[ac-1]+cc[ac+3*ido-1];
00388       tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
00389       tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
00390       ah=(k-1)*ido;
00391       ch[ah]=cc[ac-ido-1]+tr2+tr3;
00392       ch[ah+1]=cc[ac-ido]+ti2+ti3;
00393       cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
00394       ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
00395       cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
00396       ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
00397       cr5=-(ti11*tr5+ti12*tr4);
00398       ci5=-(ti11*ti5+ti12*ti4);
00399       cr4=-(ti12*tr5-ti11*tr4);
00400       ci4=-(ti12*ti5-ti11*ti4);
00401       ch[ah+l1*ido]=cr2-ci5;
00402       ch[ah+4*l1*ido]=cr2+ci5;
00403       ch[ah+l1*ido+1]=ci2+cr5;
00404       ch[ah+2*l1*ido+1]=ci3+cr4;
00405       ch[ah+2*l1*ido]=cr3-ci4;
00406       ch[ah+3*l1*ido]=cr3+ci4;
00407       ch[ah+3*l1*ido+1]=ci3-cr4;
00408       ch[ah+4*l1*ido+1]=ci2-cr5;
00409       }
00410     }
00411   else
00412     {
00413     for(k=1; k<=l1; k++)
00414       {
00415       for(i=0; i<ido-1; i+=2)
00416         {
00417         ac=i+1+(k*5-4)*ido;
00418         ti5=cc[ac]-cc[ac+3*ido];
00419         ti2=cc[ac]+cc[ac+3*ido];
00420         ti4=cc[ac+ido]-cc[ac+2*ido];
00421         ti3=cc[ac+ido]+cc[ac+2*ido];
00422         tr5=cc[ac-1]-cc[ac+3*ido-1];
00423         tr2=cc[ac-1]+cc[ac+3*ido-1];
00424         tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
00425         tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
00426         ah=i+(k-1)*ido;
00427         ch[ah]=cc[ac-ido-1]+tr2+tr3;
00428         ch[ah+1]=cc[ac-ido]+ti2+ti3;
00429         cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
00430 
00431         ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
00432         cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
00433 
00434         ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
00435         cr5=-(ti11*tr5+ti12*tr4);
00436         ci5=-(ti11*ti5+ti12*ti4);
00437         cr4=-(ti12*tr5-ti11*tr4);
00438         ci4=-(ti12*ti5-ti11*ti4);
00439         dr3=cr3-ci4;
00440         dr4=cr3+ci4;
00441         di3=ci3+cr4;
00442         di4=ci3-cr4;
00443         dr5=cr2+ci5;
00444         dr2=cr2-ci5;
00445         di5=ci2-cr5;
00446         di2=ci2+cr5;
00447         ch[ah+l1*ido]=wa1[i]*dr2+wa1[i+1]*di2;
00448         ch[ah+l1*ido+1]=wa1[i]*di2-wa1[i+1]*dr2;
00449         ch[ah+2*l1*ido]=wa2[i]*dr3+wa2[i+1]*di3;
00450         ch[ah+2*l1*ido+1]=wa2[i]*di3-wa2[i+1]*dr3;
00451         ch[ah+3*l1*ido]=wa3[i]*dr4+wa3[i+1]*di4;
00452         ch[ah+3*l1*ido+1]=wa3[i]*di4-wa3[i+1]*dr4;
00453         ch[ah+4*l1*ido]=wa4[i]*dr5+wa4[i+1]*di5;
00454         ch[ah+4*l1*ido+1]=wa4[i]*di5-wa4[i+1]*dr5;
00455         }
00456       }
00457     }
00458   }
00459 
00460 static void passb5(int ido, int l1, const double *cc, double *ch,
00461   const double *wa1, const double *wa2, const double *wa3,
00462   const double *wa4)
00463   {
00464   static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212;
00465   static const double tr12=-0.8090169943749474241, ti12=0.58778525229247312917;
00466   int i, k, ac, ah;
00467   double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
00468          ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
00469 
00470   if(ido==2)
00471     {
00472     for(k=1; k<=l1;++k)
00473       {
00474       ac=(5*k-4)*ido+1;
00475       ti5=cc[ac]-cc[ac+3*ido];
00476       ti2=cc[ac]+cc[ac+3*ido];
00477       ti4=cc[ac+ido]-cc[ac+2*ido];
00478       ti3=cc[ac+ido]+cc[ac+2*ido];
00479       tr5=cc[ac-1]-cc[ac+3*ido-1];
00480       tr2=cc[ac-1]+cc[ac+3*ido-1];
00481       tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
00482       tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
00483       ah=(k-1)*ido;
00484       ch[ah]=cc[ac-ido-1]+tr2+tr3;
00485       ch[ah+1]=cc[ac-ido]+ti2+ti3;
00486       cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
00487       ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
00488       cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
00489       ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
00490       cr5=ti11*tr5+ti12*tr4;
00491       ci5=ti11*ti5+ti12*ti4;
00492       cr4=ti12*tr5-ti11*tr4;
00493       ci4=ti12*ti5-ti11*ti4;
00494       ch[ah+l1*ido]=cr2-ci5;
00495       ch[ah+4*l1*ido]=cr2+ci5;
00496       ch[ah+l1*ido+1]=ci2+cr5;
00497       ch[ah+2*l1*ido+1]=ci3+cr4;
00498       ch[ah+2*l1*ido]=cr3-ci4;
00499       ch[ah+3*l1*ido]=cr3+ci4;
00500       ch[ah+3*l1*ido+1]=ci3-cr4;
00501       ch[ah+4*l1*ido+1]=ci2-cr5;
00502       }
00503     }
00504   else
00505     {
00506     for(k=1; k<=l1; k++)
00507       {
00508       for(i=0; i<ido-1; i+=2)
00509         {
00510         ac=i+1+(k*5-4)*ido;
00511         ti5=cc[ac]-cc[ac+3*ido];
00512         ti2=cc[ac]+cc[ac+3*ido];
00513         ti4=cc[ac+ido]-cc[ac+2*ido];
00514         ti3=cc[ac+ido]+cc[ac+2*ido];
00515         tr5=cc[ac-1]-cc[ac+3*ido-1];
00516         tr2=cc[ac-1]+cc[ac+3*ido-1];
00517         tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
00518         tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
00519         ah=i+(k-1)*ido;
00520         ch[ah]=cc[ac-ido-1]+tr2+tr3;
00521         ch[ah+1]=cc[ac-ido]+ti2+ti3;
00522         cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
00523 
00524         ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
00525         cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
00526 
00527         ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
00528         cr5=ti11*tr5+ti12*tr4;
00529         ci5=ti11*ti5+ti12*ti4;
00530         cr4=ti12*tr5-ti11*tr4;
00531         ci4=ti12*ti5-ti11*ti4;
00532         dr3=cr3-ci4;
00533         dr4=cr3+ci4;
00534         di3=ci3+cr4;
00535         di4=ci3-cr4;
00536         dr5=cr2+ci5;
00537         dr2=cr2-ci5;
00538         di5=ci2-cr5;
00539         di2=ci2+cr5;
00540         ch[ah+l1*ido]=wa1[i]*dr2-wa1[i+1]*di2;
00541         ch[ah+l1*ido+1]=wa1[i]*di2+wa1[i+1]*dr2;
00542         ch[ah+2*l1*ido]=wa2[i]*dr3-wa2[i+1]*di3;
00543         ch[ah+2*l1*ido+1]=wa2[i]*di3+wa2[i+1]*dr3;
00544         ch[ah+3*l1*ido]=wa3[i]*dr4-wa3[i+1]*di4;
00545         ch[ah+3*l1*ido+1]=wa3[i]*di4+wa3[i+1]*dr4;
00546         ch[ah+4*l1*ido]=wa4[i]*dr5-wa4[i+1]*di5;
00547         ch[ah+4*l1*ido+1]=wa4[i]*di5+wa4[i+1]*dr5;
00548         }
00549       }
00550     }
00551   }
00552 
00553 static void passfg(int *nac, int ido, int ip, int l1, int idl1,
00554   double *cc, double *ch, const double *wa)
00555   {
00556   int idij, idlj, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp, idx;
00557   double wai, war;
00558 
00559   ipph=(ip+1)/ 2;
00560   idp=ip*ido;
00561   for(j=1; j<ipph; j++)
00562     {
00563     jc=ip-j;
00564     for(k=0; k<l1; k++)
00565       {
00566       for(i=0; i<ido; i++)
00567         {
00568         ch[i+(k+j*l1)*ido] = cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
00569         ch[i+(k+jc*l1)*ido]= cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
00570         }
00571       }
00572     }
00573   for(k=0; k<l1; k++)
00574     memcpy (ch+k*ido, cc+k*ip*ido, ido*sizeof(double));
00575 
00576   idl=2-ido;
00577   inc=0;
00578   for(l=1; l<ipph; l++)
00579     {
00580     lc=ip-l;
00581     idl+=ido;
00582     for(ik=0; ik<idl1; ik++)
00583       {
00584       cc[ik+l*idl1]=ch[ik]+wa[idl-2]*ch[ik+idl1];
00585       cc[ik+lc*idl1]=-wa[idl-1]*ch[ik+(ip-1)*idl1];
00586       }
00587     idlj=idl;
00588     inc+=ido;
00589     for(j=2; j<ipph; j++)
00590       {
00591       jc=ip-j;
00592       idlj+=inc;
00593       if(idlj>idp)
00594         idlj-=idp;
00595       war=wa[idlj-2];
00596       wai=wa[idlj-1];
00597       for(ik=0; ik<idl1; ik++)
00598         {
00599         cc[ik+l*idl1]+=war*ch[ik+j*idl1];
00600         cc[ik+lc*idl1]-=wai*ch[ik+jc*idl1];
00601         }
00602       }
00603     }
00604   for(j=1; j<ipph; j++)
00605     for(ik=0; ik<idl1; ik++)
00606       ch[ik]+=ch[ik+j*idl1];
00607   for(j=1; j<ipph; j++)
00608     {
00609     jc=ip-j;
00610     for(ik=1; ik<idl1; ik+=2)
00611       {
00612       ch[ik-1+j *idl1]=cc[ik-1+j*idl1]-cc[ik  +jc*idl1];
00613       ch[ik-1+jc*idl1]=cc[ik-1+j*idl1]+cc[ik  +jc*idl1];
00614       ch[ik  +j *idl1]=cc[ik  +j*idl1]+cc[ik-1+jc*idl1];
00615       ch[ik  +jc*idl1]=cc[ik  +j*idl1]-cc[ik-1+jc*idl1];
00616       }
00617     }
00618   *nac=1;
00619   if(ido==2)
00620     return;
00621   *nac=0;
00622   for(ik=0; ik<idl1; ik++)
00623     cc[ik]=ch[ik];
00624   for(j=1; j<ip; j++)
00625     {
00626     for(k=0; k<l1; k++)
00627       {
00628       cc[(k+j*l1)*ido  ]=ch[(k+j*l1)*ido  ];
00629       cc[(k+j*l1)*ido+1]=ch[(k+j*l1)*ido+1];
00630       }
00631     }
00632 
00633   idj=2-ido;
00634   for(j=1; j<ip; j++)
00635     {
00636     idj+=ido;
00637     for(k=0; k<l1; k++)
00638       {
00639       idij=idj;
00640       for(i=3; i<ido; i+=2)
00641         {
00642         idij+=2;
00643         idx = (k+j*l1)*ido;
00644         cc[i-1+idx] = wa[idij-2]*ch[i-1+idx]+wa[idij-1]*ch[i  +idx];
00645         cc[i  +idx] = wa[idij-2]*ch[i  +idx]-wa[idij-1]*ch[i-1+idx];
00646         }
00647       }
00648     }
00649   }
00650 
00651 static void passbg(int *nac, int ido, int ip, int l1, int idl1,
00652   double *cc, double *ch, const double *wa)
00653   {
00654   int idij, idlj, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp, idx;
00655   double wai, war;
00656 
00657   ipph=(ip+1)/ 2;
00658   idp=ip*ido;
00659   for(j=1; j<ipph; j++)
00660     {
00661     jc=ip-j;
00662     for(k=0; k<l1; k++)
00663       {
00664       for(i=0; i<ido; i++)
00665         {
00666         ch[i+(k+j*l1)*ido] = cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
00667         ch[i+(k+jc*l1)*ido]= cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
00668         }
00669       }
00670     }
00671   for(k=0; k<l1; k++)
00672     memcpy (ch+k*ido, cc+k*ip*ido, ido*sizeof(double));
00673 
00674   idl=2-ido;
00675   inc=0;
00676   for(l=1; l<ipph; l++)
00677     {
00678     lc=ip-l;
00679     idl+=ido;
00680     for(ik=0; ik<idl1; ik++)
00681       {
00682       cc[ik+l*idl1]=ch[ik]+wa[idl-2]*ch[ik+idl1];
00683       cc[ik+lc*idl1]=wa[idl-1]*ch[ik+(ip-1)*idl1];
00684       }
00685     idlj=idl;
00686     inc+=ido;
00687     for(j=2; j<ipph; j++)
00688       {
00689       jc=ip-j;
00690       idlj+=inc;
00691       if(idlj>idp)
00692         idlj-=idp;
00693       war=wa[idlj-2];
00694       wai=wa[idlj-1];
00695       for(ik=0; ik<idl1; ik++)
00696         {
00697         cc[ik+l*idl1]+=war*ch[ik+j*idl1];
00698         cc[ik+lc*idl1]+=wai*ch[ik+jc*idl1];
00699         }
00700       }
00701     }
00702   for(j=1; j<ipph; j++)
00703     for(ik=0; ik<idl1; ik++)
00704       ch[ik]+=ch[ik+j*idl1];
00705   for(j=1; j<ipph; j++)
00706     {
00707     jc=ip-j;
00708     for(ik=1; ik<idl1; ik+=2)
00709       {
00710       ch[ik-1+j *idl1]=cc[ik-1+j*idl1]-cc[ik  +jc*idl1];
00711       ch[ik-1+jc*idl1]=cc[ik-1+j*idl1]+cc[ik  +jc*idl1];
00712       ch[ik  +j *idl1]=cc[ik  +j*idl1]+cc[ik-1+jc*idl1];
00713       ch[ik  +jc*idl1]=cc[ik  +j*idl1]-cc[ik-1+jc*idl1];
00714       }
00715     }
00716   *nac=1;
00717   if(ido==2)
00718     return;
00719   *nac=0;
00720   for(ik=0; ik<idl1; ik++)
00721     cc[ik]=ch[ik];
00722   for(j=1; j<ip; j++)
00723     {
00724     for(k=0; k<l1; k++)
00725       {
00726       cc[(k+j*l1)*ido  ]=ch[(k+j*l1)*ido  ];
00727       cc[(k+j*l1)*ido+1]=ch[(k+j*l1)*ido+1];
00728       }
00729     }
00730 
00731   idj=2-ido;
00732   for(j=1; j<ip; j++)
00733     {
00734     idj+=ido;
00735     for(k=0; k<l1; k++)
00736       {
00737       idij=idj;
00738       for(i=3; i<ido; i+=2)
00739         {
00740         idij+=2;
00741         idx = (k+j*l1)*ido;
00742         cc[i-1+idx] = wa[idij-2]*ch[i-1+idx]-wa[idij-1]*ch[i  +idx];
00743         cc[i  +idx] = wa[idij-2]*ch[i  +idx]+wa[idij-1]*ch[i-1+idx];
00744         }
00745       }
00746     }
00747   }
00748 
00749 
00750 static void radf2 (int ido, int l1, const double *cc, double *ch,
00751   const double *wa1)
00752   {
00753   int i, k, ic;
00754   double ti2, tr2;
00755 
00756   for(k=0; k<l1; k++)
00757     {
00758     ch[2*k*ido] = cc[k*ido]+cc[(k+l1)*ido];
00759     ch[(2*k+1)*ido+ido-1] = cc[k*ido]-cc[(k+l1)*ido];
00760     }
00761   if(ido<2)
00762     return;
00763   if(ido !=2)
00764     {
00765     for(k=0; k<l1; k++)
00766       {
00767       for(i=2; i<ido; i+=2)
00768         {
00769         ic=ido-i;
00770         tr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
00771         ti2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
00772         ch[i+2*k*ido]=cc[i+k*ido]+ti2;
00773         ch[ic+(2*k+1)*ido]=ti2-cc[i+k*ido];
00774         ch[i-1+2*k*ido]=cc[i-1+k*ido]+tr2;
00775         ch[ic-1+(2*k+1)*ido]=cc[i-1+k*ido]-tr2;
00776         }
00777       }
00778     if(ido%2==1)
00779       return;
00780     }
00781   for(k=0; k<l1; k++)
00782     {
00783     ch[(2*k+1)*ido] = -cc[ido-1+(k+l1)*ido];
00784     ch[ido-1+2*k*ido] = cc[ido-1+k*ido];
00785     }
00786   }
00787 
00788 static void radb2(int ido, int l1, const double *cc, double *ch,
00789   const double *wa1)
00790   {
00791   int i, k, ic;
00792   double ti2, tr2;
00793 
00794   for(k=0; k<l1; k++)
00795     {
00796     ch[k*ido] = cc[2*k*ido]+cc[ido-1+(2*k+1)*ido];
00797     ch[(k+l1)*ido] = cc[2*k*ido]-cc[ido-1+(2*k+1)*ido];
00798     }
00799   if(ido<2)
00800     return;
00801   if(ido !=2)
00802     {
00803     for(k=0; k<l1;++k)
00804       {
00805       for(i=2; i<ido; i+=2)
00806         {
00807         ic=ido-i;
00808         ch[i-1+k*ido] = cc[i-1+2*k*ido]+cc[ic-1+(2*k+1)*ido];
00809         tr2 = cc[i-1+2*k*ido]-cc[ic-1+(2*k+1)*ido];
00810         ch[i+k*ido] = cc[i+2*k*ido]-cc[ic+(2*k+1)*ido];
00811         ti2 = cc[i+(2*k)*ido]+cc[ic+(2*k+1)*ido];
00812         ch[i-1+(k+l1)*ido] = wa1[i-2]*tr2-wa1[i-1]*ti2;
00813         ch[i+(k+l1)*ido] = wa1[i-2]*ti2+wa1[i-1]*tr2;
00814         }
00815       }
00816     if(ido%2==1)
00817       return;
00818     }
00819   for(k=0; k<l1; k++)
00820     {
00821     ch[ido-1+k*ido]=2*cc[ido-1+2*k*ido];
00822     ch[ido-1+(k+l1)*ido]=-2*cc[(2*k+1)*ido];
00823     }
00824   }
00825 
00826 static void radf3(int ido, int l1, const double *cc, double *ch,
00827   const double *wa1, const double *wa2)
00828   {
00829   static const double taur=-0.5, taui=0.86602540378443864676;
00830   int i, k, ic;
00831   double ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
00832 
00833   for(k=0; k<l1; k++)
00834     {
00835     cr2=cc[(k+l1)*ido]+cc[(k+2*l1)*ido];
00836     ch[3*k*ido]=cc[k*ido]+cr2;
00837     ch[(3*k+2)*ido]=taui*(cc[(k+l1*2)*ido]-cc[(k+l1)*ido]);
00838     ch[ido-1+(3*k+1)*ido]=cc[k*ido]+taur*cr2;
00839     }
00840   if(ido==1)
00841     return;
00842   for(k=0; k<l1; k++)
00843     {
00844     for(i=2; i<ido; i+=2)
00845       {
00846       ic=ido-i;
00847       dr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+
00848           wa1[i-1]*cc[i+(k+l1)*ido];
00849       di2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
00850       dr3=wa2[i-2]*cc[i-1+(k+l1*2)*ido]+wa2[i-1]*cc[i+(k+l1*2)*ido];
00851       di3=wa2[i-2]*cc[i+(k+l1*2)*ido]-wa2[i-1]*cc[i-1+(k+l1*2)*ido];
00852       cr2=dr2+dr3;
00853       ci2=di2+di3;
00854       ch[i-1+3*k*ido]=cc[i-1+k*ido]+cr2;
00855       ch[i+3*k*ido]=cc[i+k*ido]+ci2;
00856       tr2=cc[i-1+k*ido]+taur*cr2;
00857       ti2=cc[i+k*ido]+taur*ci2;
00858       tr3=taui*(di2-di3);
00859       ti3=taui*(dr3-dr2);
00860       ch[i-1+(3*k+2)*ido]=tr2+tr3;
00861       ch[ic-1+(3*k+1)*ido]=tr2-tr3;
00862       ch[i+(3*k+2)*ido]=ti2+ti3;
00863       ch[ic+(3*k+1)*ido]=ti3-ti2;
00864       }
00865     }
00866   }
00867 
00868 static void radb3(int ido, int l1, const double *cc, double *ch,
00869   const double *wa1, const double *wa2)
00870   {
00871   static const double taur=-0.5, taui=0.86602540378443864676;
00872   int i, k, ic;
00873   double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
00874 
00875   for(k=0; k<l1; k++)
00876     {
00877     tr2=2*cc[ido-1+(3*k+1)*ido];
00878     cr2=cc[3*k*ido]+taur*tr2;
00879     ch[k*ido]=cc[3*k*ido]+tr2;
00880     ci3=2*taui*cc[(3*k+2)*ido];
00881     ch[(k+l1)*ido]=cr2-ci3;
00882     ch[(k+2*l1)*ido]=cr2+ci3;
00883     }
00884   if(ido==1)
00885     return;
00886   for(k=0; k<l1; k++)
00887     {
00888     for(i=2; i<ido; i+=2)
00889       {
00890       ic=ido-i;
00891       tr2=cc[i-1+(3*k+2)*ido]+cc[ic-1+(3*k+1)*ido];
00892       cr2=cc[i-1+3*k*ido]+taur*tr2;
00893       ch[i-1+k*ido]=cc[i-1+3*k*ido]+tr2;
00894       ti2=cc[i+(3*k+2)*ido]-cc[ic+(3*k+1)*ido];
00895       ci2=cc[i+3*k*ido]+taur*ti2;
00896       ch[i+k*ido]=cc[i+3*k*ido]+ti2;
00897       cr3=taui*(cc[i-1+(3*k+2)*ido]-cc[ic-1+(3*k+1)*ido]);
00898       ci3=taui*(cc[i+(3*k+2)*ido]+cc[ic+(3*k+1)*ido]);
00899       dr2=cr2-ci3;
00900       dr3=cr2+ci3;
00901       di2=ci2+cr3;
00902       di3=ci2-cr3;
00903       ch[i-1+(k+l1)*ido]=wa1[i-2]*dr2-wa1[i-1]*di2;
00904       ch[i+(k+l1)*ido]=wa1[i-2]*di2+wa1[i-1]*dr2;
00905       ch[i-1+(k+2*l1)*ido]=wa2[i-2]*dr3-wa2[i-1]*di3;
00906       ch[i+(k+2*l1)*ido]=wa2[i-2]*di3+wa2[i-1]*dr3;
00907       }
00908     }
00909   }
00910 
00911 static void radf4(int ido, int l1, const double *cc, double *ch,
00912   const double *wa1, const double *wa2, const double *wa3)
00913   {
00914   static const double hsqt2=0.70710678118654752440;
00915   int i, k, ic;
00916   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
00917 
00918   for(k=0; k<l1; k++)
00919     {
00920     tr1=cc[(k+l1)*ido]+cc[(k+3*l1)*ido];
00921     tr2=cc[k*ido]+cc[(k+2*l1)*ido];
00922     ch[4*k*ido]=tr1+tr2;
00923     ch[ido-1+(4*k+3)*ido]=tr2-tr1;
00924     ch[ido-1+(4*k+1)*ido]=cc[k*ido]-cc[(k+2*l1)*ido];
00925     ch[(4*k+2)*ido]=cc[(k+3*l1)*ido]-cc[(k+l1)*ido];
00926     }
00927   if(ido<2)
00928     return;
00929   if(ido !=2)
00930     {
00931     for(k=0; k<l1; k++)
00932       {
00933       for(i=2; i<ido; i+=2)
00934         {
00935         ic=ido-i;
00936         cr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
00937         ci2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
00938         cr3=wa2[i-2]*cc[i-1+(k+2*l1)*ido]+wa2[i-1]*cc[i+(k+2*l1)*ido];
00939         ci3=wa2[i-2]*cc[i+(k+2*l1)*ido]-wa2[i-1]*cc[i-1+(k+2*l1)*ido];
00940         cr4=wa3[i-2]*cc[i-1+(k+3*l1)*ido]+wa3[i-1]*cc[i+(k+3*l1)*ido];
00941         ci4=wa3[i-2]*cc[i+(k+3*l1)*ido]-wa3[i-1]*cc[i-1+(k+3*l1)*ido];
00942         tr1=cr2+cr4;
00943         tr4=cr4-cr2;
00944         ti1=ci2+ci4;
00945         ti4=ci2-ci4;
00946         ti2=cc[i+k*ido]+ci3;
00947         ti3=cc[i+k*ido]-ci3;
00948         tr2=cc[i-1+k*ido]+cr3;
00949         tr3=cc[i-1+k*ido]-cr3;
00950         ch[i-1+4*k*ido]=tr1+tr2;
00951         ch[ic-1+(4*k+3)*ido]=tr2-tr1;
00952         ch[i+4*k*ido]=ti1+ti2;
00953         ch[ic+(4*k+3)*ido]=ti1-ti2;
00954         ch[i-1+(4*k+2)*ido]=ti4+tr3;
00955         ch[ic-1+(4*k+1)*ido]=tr3-ti4;
00956         ch[i+(4*k+2)*ido]=tr4+ti3;
00957         ch[ic+(4*k+1)*ido]=tr4-ti3;
00958         }
00959       }
00960     if(ido%2==1)
00961       return;
00962     }
00963   for(k=0; k<l1; k++)
00964     {
00965     ti1=-hsqt2*(cc[ido-1+(k+l1)*ido]+cc[ido-1+(k+3*l1)*ido]);
00966     tr1=hsqt2*(cc[ido-1+(k+l1)*ido]-cc[ido-1+(k+3*l1)*ido]);
00967     ch[ido-1+4*k*ido]=tr1+cc[ido-1+k*ido];
00968     ch[ido-1+(4*k+2)*ido]=cc[ido-1+k*ido]-tr1;
00969     ch[(4*k+1)*ido]=ti1-cc[ido-1+(k+2*l1)*ido];
00970     ch[(4*k+3)*ido]=ti1+cc[ido-1+(k+2*l1)*ido];
00971     }
00972   }
00973 
00974 static void radb4(int ido, int l1, const double *cc, double *ch,
00975   const double *wa1, const double *wa2, const double *wa3)
00976   {
00977   static const double sqrt2=1.41421356237309504880;
00978   int i, k, ic;
00979   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
00980 
00981   for(k=0; k<l1; k++)
00982     {
00983     tr1=cc[4*k*ido]-cc[ido-1+(4*k+3)*ido];
00984     tr2=cc[4*k*ido]+cc[ido-1+(4*k+3)*ido];
00985     tr3=cc[ido-1+(4*k+1)*ido]+cc[ido-1+(4*k+1)*ido];
00986     tr4=cc[(4*k+2)*ido]+cc[(4*k+2)*ido];
00987     ch[k*ido]=tr2+tr3;
00988     ch[(k+l1)*ido]=tr1-tr4;
00989     ch[(k+2*l1)*ido]=tr2-tr3;
00990     ch[(k+3*l1)*ido]=tr1+tr4;
00991     }
00992   if(ido<2)
00993     return;
00994   if(ido !=2)
00995     {
00996     for(k=0; k<l1;++k)
00997       {
00998       for(i=2; i<ido; i+=2)
00999         {
01000         ic=ido-i;
01001         ti1=cc[i+4*k*ido]+cc[ic+(4*k+3)*ido];
01002         ti2=cc[i+4*k*ido]-cc[ic+(4*k+3)*ido];
01003         ti3=cc[i+(4*k+2)*ido]-cc[ic+(4*k+1)*ido];
01004         tr4=cc[i+(4*k+2)*ido]+cc[ic+(4*k+1)*ido];
01005         tr1=cc[i-1+4*k*ido]-cc[ic-1+(4*k+3)*ido];
01006         tr2=cc[i-1+4*k*ido]+cc[ic-1+(4*k+3)*ido];
01007         ti4=cc[i-1+(4*k+2)*ido]-cc[ic-1+(4*k+1)*ido];
01008         tr3=cc[i-1+(4*k+2)*ido]+cc[ic-1+(4*k+1)*ido];
01009         ch[i-1+k*ido]=tr2+tr3;
01010         cr3=tr2-tr3;
01011         ch[i+k*ido]=ti2+ti3;
01012         ci3=ti2-ti3;
01013         cr2=tr1-tr4;
01014         cr4=tr1+tr4;
01015         ci2=ti1+ti4;
01016         ci4=ti1-ti4;
01017         ch[i-1+(k+l1)*ido]=wa1[i-2]*cr2-wa1[i-1]*ci2;
01018         ch[i+(k+l1)*ido]=wa1[i-2]*ci2+wa1[i-1]*cr2;
01019         ch[i-1+(k+2*l1)*ido]=wa2[i-2]*cr3-wa2[i-1]*ci3;
01020         ch[i+(k+2*l1)*ido]=wa2[i-2]*ci3+wa2[i-1]*cr3;
01021         ch[i-1+(k+3*l1)*ido]=wa3[i-2]*cr4-wa3[i-1]*ci4;
01022         ch[i+(k+3*l1)*ido]=wa3[i-2]*ci4+wa3[i-1]*cr4;
01023         }
01024       }
01025     if(ido%2==1)
01026       return;
01027     }
01028   for(k=0; k<l1; k++)
01029     {
01030     ti1=cc[(4*k+1)*ido]+cc[(4*k+3)*ido];
01031     ti2=cc[(4*k+3)*ido]-cc[(4*k+1)*ido];
01032     tr1=cc[ido-1+4*k*ido]-cc[ido-1+(4*k+2)*ido];
01033     tr2=cc[ido-1+4*k*ido]+cc[ido-1+(4*k+2)*ido];
01034     ch[ido-1+k*ido]=tr2+tr2;
01035     ch[ido-1+(k+l1)*ido]=sqrt2*(tr1-ti1);
01036     ch[ido-1+(k+2*l1)*ido]=ti2+ti2;
01037     ch[ido-1+(k+3*l1)*ido]=-sqrt2*(tr1+ti1);
01038     }
01039   }
01040 
01041 static void radf5(int ido, int l1, const double *cc, double *ch,
01042   const double *wa1, const double *wa2, const double *wa3, const double *wa4)
01043   {
01044   static const double tr11=0.3090169943749474241;
01045   static const double ti11=0.95105651629515357212;
01046   static const double tr12=-0.8090169943749474241;
01047   static const double ti12=0.58778525229247312917;
01048   int i, k, ic;
01049   double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
01050          dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
01051 
01052   for(k=0; k<l1; k++)
01053     {
01054     cr2=cc[(k+4*l1)*ido]+cc[(k+l1)*ido];
01055     ci5=cc[(k+4*l1)*ido]-cc[(k+l1)*ido];
01056     cr3=cc[(k+3*l1)*ido]+cc[(k+2*l1)*ido];
01057     ci4=cc[(k+3*l1)*ido]-cc[(k+2*l1)*ido];
01058     ch[5*k*ido]=cc[k*ido]+cr2+cr3;
01059     ch[ido-1+(5*k+1)*ido]=cc[k*ido]+tr11*cr2+tr12*cr3;
01060     ch[(5*k+2)*ido]=ti11*ci5+ti12*ci4;
01061     ch[ido-1+(5*k+3)*ido]=cc[k*ido]+tr12*cr2+tr11*cr3;
01062     ch[(5*k+4)*ido]=ti12*ci5-ti11*ci4;
01063     }
01064   if(ido==1)
01065     return;
01066   for(k=0; k<l1;++k)
01067     {
01068     for(i=2; i<ido; i+=2)
01069       {
01070       ic=ido-i;
01071       dr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
01072       di2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
01073       dr3=wa2[i-2]*cc[i-1+(k+2*l1)*ido]+wa2[i-1]*cc[i+(k+2*l1)*ido];
01074       di3=wa2[i-2]*cc[i+(k+2*l1)*ido]-wa2[i-1]*cc[i-1+(k+2*l1)*ido];
01075       dr4=wa3[i-2]*cc[i-1+(k+3*l1)*ido]+wa3[i-1]*cc[i+(k+3*l1)*ido];
01076       di4=wa3[i-2]*cc[i+(k+3*l1)*ido]-wa3[i-1]*cc[i-1+(k+3*l1)*ido];
01077       dr5=wa4[i-2]*cc[i-1+(k+4*l1)*ido]+wa4[i-1]*cc[i+(k+4*l1)*ido];
01078       di5=wa4[i-2]*cc[i+(k+4*l1)*ido]-wa4[i-1]*cc[i-1+(k+4*l1)*ido];
01079       cr2=dr2+dr5;
01080       ci5=dr5-dr2;
01081       cr5=di2-di5;
01082       ci2=di2+di5;
01083       cr3=dr3+dr4;
01084       ci4=dr4-dr3;
01085       cr4=di3-di4;
01086       ci3=di3+di4;
01087       ch[i-1+5*k*ido]=cc[i-1+k*ido]+cr2+cr3;
01088       ch[i+5*k*ido]=cc[i+k*ido]+ci2+ci3;
01089       tr2=cc[i-1+k*ido]+tr11*cr2+tr12*cr3;
01090       ti2=cc[i+k*ido]+tr11*ci2+tr12*ci3;
01091       tr3=cc[i-1+k*ido]+tr12*cr2+tr11*cr3;
01092       ti3=cc[i+k*ido]+tr12*ci2+tr11*ci3;
01093       tr5=ti11*cr5+ti12*cr4;
01094       ti5=ti11*ci5+ti12*ci4;
01095       tr4=ti12*cr5-ti11*cr4;
01096       ti4=ti12*ci5-ti11*ci4;
01097       ch[i-1+(5*k+2)*ido]=tr2+tr5;
01098       ch[ic-1+(5*k+1)*ido]=tr2-tr5;
01099       ch[i+(5*k+2)*ido]=ti2+ti5;
01100       ch[ic+(5*k+1)*ido]=ti5-ti2;
01101       ch[i-1+(5*k+4)*ido]=tr3+tr4;
01102       ch[ic-1+(5*k+3)*ido]=tr3-tr4;
01103       ch[i+(5*k+4)*ido]=ti3+ti4;
01104       ch[ic+(5*k+3)*ido]=ti4-ti3;
01105       }
01106     }
01107   }
01108 
01109 static void radb5(int ido, int l1, const double *cc, double *ch,
01110   const double *wa1, const double *wa2, const double *wa3, const double *wa4)
01111   {
01112   static const double tr11=0.3090169943749474241;
01113   static const double ti11=0.95105651629515357212;
01114   static const double tr12=-0.8090169943749474241;
01115   static const double ti12=0.58778525229247312917;
01116   int i, k, ic;
01117   double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
01118           ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01119 
01120   for(k=0; k<l1; k++)
01121     {
01122     ti5=2*cc[(5*k+2)*ido];
01123     ti4=2*cc[(5*k+4)*ido];
01124     tr2=2*cc[ido-1+(5*k+1)*ido];
01125     tr3=2*cc[ido-1+(5*k+3)*ido];
01126     ch[k*ido]=cc[5*k*ido]+tr2+tr3;
01127     cr2=cc[5*k*ido]+tr11*tr2+tr12*tr3;
01128     cr3=cc[5*k*ido]+tr12*tr2+tr11*tr3;
01129     ci5=ti11*ti5+ti12*ti4;
01130     ci4=ti12*ti5-ti11*ti4;
01131     ch[(k+l1)*ido]=cr2-ci5;
01132     ch[(k+2*l1)*ido]=cr3-ci4;
01133     ch[(k+3*l1)*ido]=cr3+ci4;
01134     ch[(k+4*l1)*ido]=cr2+ci5;
01135     }
01136   if(ido==1)
01137     return;
01138   for(k=0; k<l1;++k)
01139     {
01140     for(i=2; i<ido; i+=2)
01141       {
01142       ic=ido-i;
01143       ti5=cc[i+(5*k+2)*ido]+cc[ic+(5*k+1)*ido];
01144       ti2=cc[i+(5*k+2)*ido]-cc[ic+(5*k+1)*ido];
01145       ti4=cc[i+(5*k+4)*ido]+cc[ic+(5*k+3)*ido];
01146       ti3=cc[i+(5*k+4)*ido]-cc[ic+(5*k+3)*ido];
01147       tr5=cc[i-1+(5*k+2)*ido]-cc[ic-1+(5*k+1)*ido];
01148       tr2=cc[i-1+(5*k+2)*ido]+cc[ic-1+(5*k+1)*ido];
01149       tr4=cc[i-1+(5*k+4)*ido]-cc[ic-1+(5*k+3)*ido];
01150       tr3=cc[i-1+(5*k+4)*ido]+cc[ic-1+(5*k+3)*ido];
01151       ch[i-1+k*ido]=cc[i-1+5*k*ido]+tr2+tr3;
01152       ch[i+k*ido]=cc[i+5*k*ido]+ti2+ti3;
01153       cr2=cc[i-1+5*k*ido]+tr11*tr2+tr12*tr3;
01154 
01155       ci2=cc[i+5*k*ido]+tr11*ti2+tr12*ti3;
01156       cr3=cc[i-1+5*k*ido]+tr12*tr2+tr11*tr3;
01157 
01158       ci3=cc[i+5*k*ido]+tr12*ti2+tr11*ti3;
01159       cr5=ti11*tr5+ti12*tr4;
01160       ci5=ti11*ti5+ti12*ti4;
01161       cr4=ti12*tr5-ti11*tr4;
01162       ci4=ti12*ti5-ti11*ti4;
01163       dr3=cr3-ci4;
01164       dr4=cr3+ci4;
01165       di3=ci3+cr4;
01166       di4=ci3-cr4;
01167       dr5=cr2+ci5;
01168       dr2=cr2-ci5;
01169       di5=ci2-cr5;
01170       di2=ci2+cr5;
01171       ch[i-1+(k+l1)*ido]=wa1[i-2]*dr2-wa1[i-1]*di2;
01172       ch[i+(k+l1)*ido]=wa1[i-2]*di2+wa1[i-1]*dr2;
01173       ch[i-1+(k+2*l1)*ido]=wa2[i-2]*dr3-wa2[i-1]*di3;
01174       ch[i+(k+2*l1)*ido]=wa2[i-2]*di3+wa2[i-1]*dr3;
01175       ch[i-1+(k+3*l1)*ido]=wa3[i-2]*dr4-wa3[i-1]*di4;
01176       ch[i+(k+3*l1)*ido]=wa3[i-2]*di4+wa3[i-1]*dr4;
01177       ch[i-1+(k+4*l1)*ido]=wa4[i-2]*dr5-wa4[i-1]*di5;
01178       ch[i+(k+4*l1)*ido]=wa4[i-2]*di5+wa4[i-1]*dr5;
01179       }
01180     }
01181   }
01182 
01183 static void radfg(int ido, int ip, int l1, int idl1,
01184   double *cc, double *ch, const double *wa)
01185   {
01186   static const double twopi=6.28318530717958647692;
01187   int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
01188   double ai1, ai2, ar1, ar2, arg;
01189   double *csarr;
01190   int aidx;
01191 
01192   ipph=(ip+1)/ 2;
01193   if(ido !=1)
01194     {
01195     for(ik=0; ik<idl1; ik++)
01196       ch[ik]=cc[ik];
01197     for(j=1; j<ip; j++)
01198       for(k=0; k<l1; k++)
01199         ch[(k+j*l1)*ido]=cc[(k+j*l1)*ido];
01200 
01201     is=-ido;
01202     for(j=1; j<ip; j++)
01203       {
01204       is+=ido;
01205       for(k=0; k<l1; k++)
01206         {
01207         idij=is-1;
01208         for(i=2; i<ido; i+=2)
01209           {
01210           idij+=2;
01211           ch[i-1+(k+j*l1)*ido]=
01212             wa[idij-1]*cc[i-1+(k+j*l1)*ido]+wa[idij]*cc[i+(k+j*l1)*ido];
01213           ch[i+(k+j*l1)*ido]=
01214             wa[idij-1]*cc[i+(k+j*l1)*ido]-wa[idij]*cc[i-1+(k+j*l1)*ido];
01215           }
01216         }
01217       }
01218 
01219     for(j=1; j<ipph; j++)
01220       {
01221       jc=ip-j;
01222       for(k=0; k<l1; k++)
01223         {
01224         for(i=2; i<ido; i+=2)
01225           {
01226           cc[i-1+(k+j*l1)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
01227           cc[i-1+(k+jc*l1)*ido]=ch[i+(k+j*l1)*ido]-ch[i+(k+jc*l1)*ido];
01228           cc[i+(k+j*l1)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
01229           cc[i+(k+jc*l1)*ido]=ch[i-1+(k+jc*l1)*ido]-ch[i-1+(k+j*l1)*ido];
01230           }
01231         }
01232       }
01233     }
01234   else
01235     {                           /*now ido==1*/
01236     for(ik=0; ik<idl1; ik++)
01237       cc[ik]=ch[ik];
01238     }
01239   for(j=1; j<ipph; j++)
01240     {
01241     jc=ip-j;
01242     for(k=0; k<l1; k++)
01243       {
01244       cc[(k+j*l1)*ido]=ch[(k+j*l1)*ido]+ch[(k+jc*l1)*ido];
01245       cc[(k+jc*l1)*ido]=ch[(k+jc*l1)*ido]-ch[(k+j*l1)*ido];
01246       }
01247     }
01248 
01249   csarr=(double *)malloc(2*ip*sizeof(double));
01250   arg=twopi / ip;
01251   csarr[0]=1.;
01252   csarr[1]=0.;
01253   csarr[2]=csarr[2*ip-2]=cos(arg);
01254   csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3];
01255   for (i=2; i<=ip/2; ++i)
01256     {
01257     csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg);
01258     csarr[2*i+1]=sin(i*arg);
01259     csarr[2*ip-2*i+1]=-csarr[2*i+1];
01260     }
01261   for(l=1; l<ipph; l++)
01262     {
01263     lc=ip-l;
01264     ar1=csarr[2*l];
01265     ai1=csarr[2*l+1];
01266     for(ik=0; ik<idl1; ik++)
01267       {
01268       ch[ik+l*idl1]=cc[ik]+ar1*cc[ik+idl1];
01269       ch[ik+lc*idl1]=ai1*cc[ik+(ip-1)*idl1];
01270       }
01271     aidx=2*l;
01272     for(j=2; j<ipph; j++)
01273       {
01274       jc=ip-j;
01275       aidx+=2*l;
01276       if (aidx>=2*ip) aidx-=2*ip;
01277       ar2=csarr[aidx];
01278       ai2=csarr[aidx+1];
01279       for(ik=0; ik<idl1; ik++)
01280         {
01281         ch[ik+l*idl1]+=ar2*cc[ik+j*idl1];
01282         ch[ik+lc*idl1]+=ai2*cc[ik+jc*idl1];
01283         }
01284       }
01285     }
01286   free(csarr);
01287 
01288   for(j=1; j<ipph; j++)
01289     for(ik=0; ik<idl1; ik++)
01290       ch[ik]+=cc[ik+j*idl1];
01291 
01292   for(k=0; k<l1; k++)
01293     for(i=0; i<ido; i++)
01294       cc[i+k*ip*ido]=ch[i+k*ido];
01295   for(j=1; j<ipph; j++)
01296     {
01297     jc=ip-j;
01298     j2=2*j;
01299     for(k=0; k<l1; k++)
01300       {
01301       cc[ido-1+(j2-1+k*ip)*ido] = ch[(k+j*l1)*ido];
01302       cc[(j2+k*ip)*ido] = ch[(k+jc*l1)*ido];
01303       }
01304     }
01305   if(ido==1)
01306     return;
01307 
01308   for(j=1; j<ipph; j++)
01309     {
01310     jc=ip-j;
01311     j2=2*j;
01312     for(k=0; k<l1; k++)
01313       {
01314       for(i=2; i<ido; i+=2)
01315         {
01316         ic=ido-i;
01317         cc[i-1+(j2+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]+ch[i-1+(k+jc*l1)*ido];
01318         cc[ic-1+(j2-1+k*ip)*ido]=ch[i-1+(k+j*l1)*ido]-ch[i-1+(k+jc*l1)*ido];
01319         cc[i+(j2+k*ip)*ido]=ch[i+(k+j*l1)*ido]+ch[i+(k+jc*l1)*ido];
01320         cc[ic+(j2-1+k*ip)*ido]=ch[i+(k+jc*l1)*ido]-ch[i+(k+j*l1)*ido];
01321         }
01322       }
01323     }
01324   }
01325 
01326 static void radbg(int ido, int ip, int l1, int idl1,
01327   double *cc, double *ch, const double *wa)
01328   {
01329   static const double twopi=6.28318530717958647692;
01330   int     idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
01331   double ai1, ai2, ar1, ar2, arg;
01332   double *csarr;
01333   int aidx;
01334 
01335   ipph=(ip+1)/ 2;
01336   for(k=0; k<l1; k++)
01337     for(i=0; i<ido; i++)
01338       ch[i+k*ido]=cc[i+k*ip*ido];
01339   for(j=1; j<ipph; j++)
01340     {
01341     jc=ip-j;
01342     j2=2*j;
01343     for(k=0; k<l1; k++)
01344       {
01345       ch[(k+j*l1)*ido]=cc[ido-1+(j2-1+k*ip)*ido]+cc[ido-1+(j2-1+k*ip)*ido];
01346       ch[(k+jc*l1)*ido]=cc[(j2+k*ip)*ido]+cc[(j2+k*ip)*ido];
01347       }
01348     }
01349 
01350   if(ido !=1)
01351     {
01352     for(j=1; j<ipph; j++)
01353       {
01354       jc=ip-j;
01355       for(k=0; k<l1; k++)
01356         {
01357         for(i=2; i<ido; i+=2)
01358           {
01359           ic=ido-i;
01360           ch[i-1+(k+j*l1)*ido] =
01361             cc[i-1+(2*j+k*ip)*ido]+cc[ic-1+(2*j-1+k*ip)*ido];
01362           ch[i-1+(k+jc*l1)*ido] =
01363             cc[i-1+(2*j+k*ip)*ido]-cc[ic-1+(2*j-1+k*ip)*ido];
01364           ch[i+(k+j*l1)*ido]=cc[i+(2*j+k*ip)*ido]-cc[ic+(2*j-1+k*ip)*ido];
01365           ch[i+(k+jc*l1)*ido]=cc[i+(2*j+k*ip)*ido]+cc[ic+(2*j-1+k*ip)*ido];
01366           }
01367         }
01368       }
01369     }
01370 
01371   csarr=(double *)malloc(2*ip*sizeof(double));
01372   arg=twopi / ip;
01373   csarr[0]=1.;
01374   csarr[1]=0.;
01375   csarr[2]=csarr[2*ip-2]=cos(arg);
01376   csarr[3]=sin(arg); csarr[2*ip-1]=-csarr[3];
01377   for (i=2; i<=ip/2; ++i)
01378     {
01379     csarr[2*i]=csarr[2*ip-2*i]=cos(i*arg);
01380     csarr[2*i+1]=sin(i*arg);
01381     csarr[2*ip-2*i+1]=-csarr[2*i+1];
01382     }
01383   for(l=1; l<ipph; l++)
01384     {
01385     lc=ip-l;
01386     ar1=csarr[2*l];
01387     ai1=csarr[2*l+1];
01388     for(ik=0; ik<idl1; ik++)
01389       {
01390       cc[ik+l*idl1]=ch[ik]+ar1*ch[ik+idl1];
01391       cc[ik+lc*idl1]=ai1*ch[ik+(ip-1)*idl1];
01392       }
01393     aidx=2*l;
01394     for(j=2; j<ipph; j++)
01395       {
01396       jc=ip-j;
01397       aidx+=2*l;
01398       if (aidx>=2*ip) aidx-=2*ip;
01399       ar2=csarr[aidx];
01400       ai2=csarr[aidx+1];
01401       for(ik=0; ik<idl1; ik++)
01402         {
01403         cc[ik+l*idl1]+=ar2*ch[ik+j*idl1];
01404         cc[ik+lc*idl1]+=ai2*ch[ik+jc*idl1];
01405         }
01406       }
01407     }
01408   free(csarr);
01409 
01410   for(j=1; j<ipph; j++)
01411     for(ik=0; ik<idl1; ik++)
01412       ch[ik]+=ch[ik+j*idl1];
01413 
01414   for(j=1; j<ipph; j++)
01415     {
01416     jc=ip-j;
01417     for(k=0; k<l1; k++)
01418       {
01419       ch[(k+j*l1)*ido]=cc[(k+j*l1)*ido]-cc[(k+jc*l1)*ido];
01420       ch[(k+jc*l1)*ido]=cc[(k+j*l1)*ido]+cc[(k+jc*l1)*ido];
01421       }
01422     }
01423 
01424   if(ido==1)
01425     return;
01426   for(j=1; j<ipph; j++)
01427     {
01428     jc=ip-j;
01429     for(k=0; k<l1; k++)
01430       {
01431       for(i=2; i<ido; i+=2)
01432         {
01433         ch[i-1+(k+j*l1)*ido]=cc[i-1+(k+j*l1)*ido]-cc[i+(k+jc*l1)*ido];
01434         ch[i-1+(k+jc*l1)*ido]=cc[i-1+(k+j*l1)*ido]+cc[i+(k+jc*l1)*ido];
01435         ch[i+(k+j*l1)*ido]=cc[i+(k+j*l1)*ido]+cc[i-1+(k+jc*l1)*ido];
01436         ch[i+(k+jc*l1)*ido]=cc[i+(k+j*l1)*ido]-cc[i-1+(k+jc*l1)*ido];
01437         }
01438       }
01439     }
01440   for(ik=0; ik<idl1; ik++)
01441     cc[ik]=ch[ik];
01442   for(j=1; j<ip; j++)
01443     for(k=0; k<l1; k++)
01444       cc[(k+j*l1)*ido]=ch[(k+j*l1)*ido];
01445 
01446   is=-ido;
01447   for(j=1; j<ip; j++)
01448     {
01449     is+=ido;
01450     for(k=0; k<l1; k++)
01451       {
01452       idij=is-1;
01453       for(i=2; i<ido; i+=2)
01454         {
01455         idij+=2;
01456         cc[i-1+(k+j*l1)*ido]=
01457           wa[idij-1]*ch[i-1+(k+j*l1)*ido]-wa[idij]*ch[i+(k+j*l1)*ido];
01458         cc[i+(k+j*l1)*ido]=
01459           wa[idij-1]*ch[i+(k+j*l1)*ido]+wa[idij]*ch[i-1+(k+j*l1)*ido];
01460         }
01461       }
01462     }
01463   }
01464 
01465 
01466 /*----------------------------------------------------------------------
01467    cfftf1, cfftb1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
01468   ----------------------------------------------------------------------*/
01469 
01470 static void cfftf1(int n, double c[], double ch[], const double wa[],
01471   const int ifac[])
01472   {
01473   int idot, k1, l1, l2, na, nf, ip, iw, nac, ido, idl1;
01474   double *p1, *p2;
01475 
01476   nf=ifac[1];
01477   na=0;
01478   l1=1;
01479   iw=0;
01480   for(k1=2; k1<=nf+1; k1++)
01481     {
01482     ip=ifac[k1];
01483     l2=ip*l1;
01484     ido=n / l2;
01485     idot=ido+ido;
01486     idl1=idot*l1;
01487     p1 = (na==0) ? c : ch;
01488     p2 = (na==0) ? ch : c;
01489     if(ip==4)
01490       passf4(idot, l1, p1, p2, wa+iw, wa+iw+idot, wa+iw+2*idot);
01491     else if(ip==2)
01492       passf2(idot, l1, p1, p2, wa+iw);
01493     else if(ip==3)
01494       passf3(idot, l1, p1, p2, wa+iw, wa+iw+idot);
01495     else if(ip==5)
01496       passf5(idot, l1, p1, p2, wa+iw, wa+iw+idot, wa+iw+2*idot, wa+iw+3*idot);
01497     else
01498       {
01499       passfg(&nac, idot, ip, l1, idl1, p1, p2, &wa[iw]);
01500       if(nac==0)
01501         na=1-na;
01502       }
01503     na=1-na;
01504     l1=l2;
01505     iw+=(ip-1)*idot;
01506     }
01507   if(na!=0)
01508     memcpy (c,ch,2*n*sizeof(double));
01509   }
01510 
01511 static void cfftb1(int n, double c[], double ch[], const double wa[],
01512   const int ifac[])
01513   {
01514   int idot, k1, l1, l2, na, nf, ip, iw, nac, ido, idl1;
01515   double *p1, *p2;
01516 
01517   nf=ifac[1];
01518   na=0;
01519   l1=1;
01520   iw=0;
01521   for(k1=2; k1<=nf+1; k1++)
01522     {
01523     ip=ifac[k1];
01524     l2=ip*l1;
01525     ido=n / l2;
01526     idot=ido+ido;
01527     idl1=idot*l1;
01528     p1 = (na==0) ? c : ch;
01529     p2 = (na==0) ? ch : c;
01530     if(ip==4)
01531       passb4(idot, l1, p1, p2, wa+iw, wa+iw+idot, wa+iw+2*idot);
01532     else if(ip==2)
01533       passb2(idot, l1, p1, p2, wa+iw);
01534     else if(ip==3)
01535       passb3(idot, l1, p1, p2, wa+iw, wa+iw+idot);
01536     else if(ip==5)
01537       passb5(idot, l1, p1, p2, wa+iw, wa+iw+idot, wa+iw+2*idot, wa+iw+3*idot);
01538     else
01539       {
01540       passbg(&nac, idot, ip, l1, idl1, p1, p2, &wa[iw]);
01541       if(nac==0)
01542         na=1-na;
01543       }
01544     na=1-na;
01545     l1=l2;
01546     iw+=(ip-1)*idot;
01547     }
01548   if(na!=0)
01549     memcpy (c,ch,2*n*sizeof(double));
01550   }
01551 
01552 void cfftf(int n, double c[], double wsave[])
01553   {
01554   if(n!=1)
01555     cfftf1(n, c, wsave, wsave+2*n,(int*)(wsave+4*n));
01556   }
01557 
01558 void cfftb(int n, double c[], double wsave[])
01559   {
01560   if(n!=1)
01561     cfftb1(n, c, wsave, wsave+2*n,(int*)(wsave+4*n));
01562   }
01563 
01564 static void cffti1(int n, double wa[], int ifac[])
01565   {
01566   static const int ntryh[4]= {3, 4, 2, 5};
01567   static const double twopi=6.28318530717958647692;
01568   double argh, argld, arg, fi;
01569   int idot, ntry=0, i, j, i1, k1, l1, l2, ib;
01570   int ld, ii, nf, ip, nl, nq, nr, ido, ipm;
01571 
01572   nl=n;
01573   nf=0;
01574   j=0;
01575 startloop:
01576   j++;
01577   if(j<=4)
01578     ntry=ntryh[j-1];
01579   else
01580     ntry+=2;
01581   do
01582     {
01583     nq=nl / ntry;
01584     nr=nl-ntry*nq;
01585     if(nr !=0)
01586       goto startloop;
01587     nf++;
01588     ifac[nf+1]=ntry;
01589     nl=nq;
01590     if(ntry==2 && nf !=1)
01591       {
01592       for(i=2; i<=nf; i++)
01593         {
01594         ib=nf-i+2;
01595         ifac[ib+1]=ifac[ib];
01596         }
01597       ifac[2]=2;
01598       }
01599     }
01600   while(nl !=1);
01601   ifac[0]=n;
01602   ifac[1]=nf;
01603   argh=twopi /(double)n;
01604   i=1;
01605   l1=1;
01606   for(k1=1; k1<=nf; k1++)
01607     {
01608     ip=ifac[k1+1];
01609     ld=0;
01610     l2=l1*ip;
01611     ido=n / l2;
01612     idot=ido+ido+2;
01613     ipm=ip-1;
01614     for(j=1; j<=ipm; j++)
01615       {
01616       i1=i;
01617       wa[i-1]=1;
01618       wa[i]=0;
01619       ld+=l1;
01620       fi=0;
01621       argld=ld*argh;
01622       for(ii=4; ii<=idot; ii+=2)
01623         {
01624         i+=2;
01625         fi+=1;
01626         arg=fi*argld;
01627         wa[i-1]=cos(arg);
01628         wa[i]=sin(arg);
01629         }
01630       if(ip>5)
01631         {
01632         wa[i1-1]=wa[i-1];
01633         wa[i1]=wa[i];
01634         }
01635       }
01636     l1=l2;
01637     }
01638   }
01639 
01640 void cffti(int n, double wsave[])
01641   {
01642   if (n!=1)
01643     cffti1(n, wsave+2*n,(int*)(wsave+4*n));
01644   }
01645 
01646 
01647 /*----------------------------------------------------------------------
01648    rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Real FFTs.
01649   ----------------------------------------------------------------------*/
01650 
01651 static void rfftf1(int n, double c[], double ch[], const double wa[],
01652   const int ifac[])
01653   {
01654   int i, k1, l1, l2, na, kh, nf, ip, iw, ido, idl1;
01655   double *p1, *p2;
01656 
01657   nf=ifac[1];
01658   na=1;
01659   l2=n;
01660   iw=n-1;
01661   for(k1=1; k1<=nf;++k1)
01662     {
01663     kh=nf-k1;
01664     ip=ifac[kh+2];
01665     l1=l2 / ip;
01666     ido=n / l2;
01667     idl1=ido*l1;
01668     iw-=(ip-1)*ido;
01669     na=1-na;
01670     p1 = (na==0) ? c : ch;
01671     p2 = (na==0) ? ch : c;
01672     if(ip==4)
01673       radf4(ido, l1, p1, p2, wa+iw, wa+iw+ido, wa+iw+2*ido);
01674     else if(ip==2)
01675       radf2(ido, l1, p1, p2, wa+iw);
01676     else if(ip==3)
01677       radf3(ido, l1, p1, p2, wa+iw, wa+iw+ido);
01678     else if(ip==5)
01679       radf5(ido, l1, p1, p2, wa+iw, wa+iw+ido, wa+iw+2*ido, wa+iw+3*ido);
01680     else
01681       {
01682       if(ido==1)
01683         na=1-na;
01684       if(na==0)
01685         radfg(ido, ip, l1, idl1, c, ch, wa+iw);
01686       else
01687         radfg(ido, ip, l1, idl1, ch, c, wa+iw);
01688       na=1-na;
01689       }
01690     l2=l1;
01691     }
01692   if(na==1)
01693     return;
01694   for(i=0; i<n; i++)
01695     c[i]=ch[i];
01696 }
01697 
01698 static void rfftb1(int n, double c[], double ch[], const double wa[],
01699   const int ifac[])
01700   {
01701   int k1, l1, l2, na, nf, ip, iw, ido, idl1;
01702   double *p1, *p2;
01703 
01704   nf=ifac[1];
01705   na=0;
01706   l1=1;
01707   iw=0;
01708   for(k1=1; k1<=nf; k1++)
01709     {
01710     ip=ifac[k1+1];
01711     l2=ip*l1;
01712     ido=n / l2;
01713     idl1=ido*l1;
01714     p1 = (na==0) ? c : ch;
01715     p2 = (na==0) ? ch : c;
01716     if(ip==4)
01717       radb4(ido, l1, p1, p2, wa+iw, wa+iw+ido, wa+iw+2*ido);
01718     else if(ip==2)
01719       radb2(ido, l1, p1, p2, wa+iw);
01720     else if(ip==3)
01721       radb3(ido, l1, p1, p2, wa+iw, wa+iw+ido);
01722     else if(ip==5)
01723       radb5(ido, l1, p1, p2, wa+iw, wa+iw+ido, wa+iw+2*ido, wa+iw+3*ido);
01724     else
01725       {
01726       radbg(ido, ip, l1, idl1, p1, p2, wa+iw);
01727       if(ido!=1)
01728         na=1-na;
01729       }
01730     na=1-na;
01731     l1=l2;
01732     iw+=(ip-1)*ido;
01733     }
01734   if(na!=0)
01735     memcpy (c,ch,n*sizeof(double));
01736   }
01737 
01738 void rfftf(int n, double r[], double wsave[])
01739   {
01740   if(n!=1)
01741     rfftf1(n, r, wsave, wsave+n,(int*)(wsave+2*n));
01742   }
01743 
01744 void rfftb(int n, double r[], double wsave[])
01745   {
01746   if(n!=1)
01747     rfftb1(n, r, wsave, wsave+n,(int*)(wsave+2*n));
01748   }
01749 
01750 static void rffti1(int n, double wa[], int ifac[])
01751   {
01752   static const int ntryh[4]={4, 2, 3, 5};
01753   static const double twopi=6.28318530717958647692;
01754   double argh, argld, arg, fi;
01755   int ntry=0, i, j, k1, l1, l2, ib, ld, ii, nf, ip, nl, is, nq, nr;
01756   int ido, ipm, nfm1;
01757 
01758   nl=n;
01759   nf=0;
01760   j=0;
01761 startloop:
01762  ++j;
01763   if(j<=4)
01764     ntry=ntryh[j-1];
01765   else
01766     ntry+=2;
01767   do
01768     {
01769     nq=nl / ntry;
01770     nr=nl-ntry*nq;
01771     if(nr !=0)
01772       goto startloop;
01773     ++nf;
01774     ifac[nf+1]=ntry;
01775     nl=nq;
01776     if(ntry==2 && nf !=1)
01777       {
01778       for(i=2; i<=nf; i++)
01779         {
01780         ib=nf-i+2;
01781         ifac[ib+1]=ifac[ib];
01782         }
01783       ifac[2]=2;
01784       }
01785     }
01786   while(nl !=1);
01787   ifac[0]=n;
01788   ifac[1]=nf;
01789   argh=twopi /(double)(n);
01790   is=0;
01791   nfm1=nf-1;
01792   l1=1;
01793   if(nfm1==0)
01794     return;
01795   for(k1=1; k1<=nfm1; k1++)
01796     {
01797     ip=ifac[k1+1];
01798     ld=0;
01799     l2=l1*ip;
01800     ido=n / l2;
01801     ipm=ip-1;
01802     for(j=1; j<=ipm;++j)
01803       {
01804       ld+=l1;
01805       i=is;
01806       argld=(double)ld*argh;
01807 
01808       fi=0;
01809       for(ii=3; ii<=ido; ii+=2)
01810         {
01811         i+=2;
01812         fi+=1;
01813         arg=fi*argld;
01814         wa[i-2]=cos(arg);
01815         wa[i-1]=sin(arg);
01816         }
01817       is+=ido;
01818       }
01819     l1=l2;
01820     }
01821   }
01822 
01823 void rffti(int n, double wsave[])
01824   {
01825   if (n!=1)
01826     rffti1(n, wsave+n,(int*)(wsave+2*n));
01827   }