(*FAST FOURIER TRANSFORM ALGORITHM 8.3
*
*  To compute the coefficients in the discrete approximation 
*   for the data (x(J),y(J)), 0<=J<=2m-1 where m=2*p and
*   x(J)=-pi+J*pi/m for 0<=J<=2m-1.
*
*  INPUT: m; y(0),y(1), ...y(2m-1).
*
* OUTPUT: complex numbers c(0), ..., c(2m-1);
*	  real numbers a(0), ..., a(m); b(1), ..., b(m-1).
*)
Input["This is the Fast Fourier Transform.\n
 The user must make provisions if the\n
 interval is not [-pi,pi].\n
 \n
 Press 1 [enter] to continue\n"];
OK = 0;
While[OK == 0,
   FLAG = Input["Choice of input method\n
   \n
   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n
   3. Generate data using a function\n"];
   If[FLAG==1 || FLAG==2 || FLAG==3,
      OK = 1;
   ];
];
If[FLAG == 1,
   OK = 0;
   While[OK == 0,
      M = Input["Input M\n"];
      If[M > 0, 
	 OK = 1;
	 n = 2*M;
	 For[JJ = 1,JJ <= n,JJ++,
	    J = JJ-1;
	    Y[JJ-1] = Input["Input y("<>ToString[J]<>")\n"];
	 ],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
];
If[FLAG == 2,
   A = InputString["Has a text file been created with the\n
   entries y(0), ..., y(2m-1) separated by a blank?\n
   \n
   Enter 'yes' or 'no'\n"];
   If[A=="y" || A=="Y" || A=="yes" || A=="YES",
      NAME = InputString["Input the file name in the form -\n
      drive:\name.ext\n
      For example A:\\DATA.DTA\n"];
      INP = OpenRead[NAME];
      OK = 0;
      While[OK == 0,
	 M = Input["Input number m\n"];
         n = 2*M;
         If[n > 0, 
	    For[JJ = 1,JJ <= n,JJ++,
	       Y[JJ-1] = Read[INP,Number];
	    ];
	    Close[INP];
	    OK = 1,
	    Input["Number must be a positive integer\n
	    \n
	    Press 1 [enter] to continue\n"];
         ];
      ];
      ,
      Input["This file will end so the input file can be created\n
      \n
      Press 1 [enter] to continue\n"];
      OK = 0;
   ];
];
If[FLAG == 3,
   TEMP = Input["Input the function F(x) in terms of x\n
   \n
   For example: Cos[x]\n"];
   F[x_] := Evaluate[TEMP];
   OK = 0;
   While[OK == 0,
      M = Input["Input the number m.\n"];
      n = 2*M;
      If[n > 0,
	 For[JJ = 1,JJ <= n,JJ++,
	    Z = -Pi+(JJ-1)*Pi/M;
 	    Y[JJ-1] = F[Z];
	 ];
	 OK = 1,
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
        For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout"
   ];
   Write[OUP,"FAST FOURIER TRANSFORM\n"];
   Write[OUP,"\n"];
   TW = Log[2.0];
   (* step 1 *)
   (* use N2 for m, NG for p, NU1 for q, WW for zeta *)
   N2 = M;
   (* step 2 *)
   For[JJ = 1,JJ <= n,JJ++,
      c[JJ-1] = Y[JJ-1];
   ];
   Z = n;
   NG = N[Round[Log[n]/TW]];
   NU1 = NG-1;
   YY = 2*Pi*I/n;
   WW = Exp[YY];
   (* step 3 *)
   (*Write[OUP,TW,N2,NG,NU1];*)
   For[JJ = 1,JJ <= N2,JJ++,
      X = 1;
      YY = 1;
      For[J = 1,J <= JJ,J++,
	 YY = X*WW;
	 X = YY;
      ];
      W[JJ-1] = X;
      W[N2+JJ-1] = -X;
   ];
   (* step 4 *)
   K = 0;
   W[-1] = 1;
   (* step 5 *)
   For[L = 1,L <= NG,L++,
   (* step 6 *)
      While[K < n-1,
      (* step 7 *)
	 For[JJ = 1,JJ <= N2,JJ++,
     	    (* step 8 *)
	    Z = Exp[NU1*TW];
	    M1 = Round[Z];
	    M1 = Floor[K/M1];
	    (* IBR does the bit reversal *)
	    R1 = M1;
            S = 0;
	    For[RR = 1,RR <= NG,RR++,
	       R2 = Floor[R1/2];
 	       S = 2*S+R1-2*R2;
	       R1 = R2;
            ];
            NP = S;
	    (* T1 is eta *)
	    T1 = c[K+N2];
	    (* step 9 *)
	    If[NP != 0,
	       X = T1;
	       T1 = X*W[NP-1];
	    ];
	    c[K+N2] = c[K]-T1;
            c[K] = c[K]+T1;
	    (* step 10 *)
	    K = K+1;
         ];
	 (*  step 12 *)
	 K = K+N2;
      ];
      K = 0;
      N2 = Floor[N2/2];
      NU1 = NU1-1;
   ];
   (* step 13 *)
   While[K < n-1,
   (* step 14 *)
      R1 = K;
      JJ = 0;
      For[RR = 1,RR <= NG,RR++,
	 R2 = Floor[R1/2];
	 JJ = 2*JJ+R1-2*R2;
	 R1 = R2;
      ];
      (* step 15 *)
      If[JJ > K,
	 T3 = c[K];
	 c[K] = c[JJ];
	 c[JJ] = T3;
      ];
      (* step 16 *)
      K = K+1;
   ];
   (* steps 17 and 18 *)
   Write[OUP,"\n"];
   Write[OUP,"Coefficients c(0), ... , c(2m-1)\n"]; 
   Write[OUP,"\n"];
   For[JJ = 1,JJ <= n,JJ++,
      YY = -(JJ-1)*Pi*I;
      X = Exp[YY];
      YY = X*c[JJ-1];
      c[JJ-1] = N[YY/(0.5*n)];
      K = JJ-1;
      Write[OUP,K,"      ",SetPrecision[c[JJ-1],9]];
   ];  
   Write[OUP,"\n"];
   Write[OUP,"Coefficients a(0), ... , a(m)\n"];
   Write[OUP,"\n"];
   For[JJ = 1,JJ <= M+1,JJ++,
      Write[OUP,SetPrecision[Re[c[JJ-1]],9]];
   ];
   Write[OUP,"\n"];
   Write[OUP,"Coefficients b(1), ... , b(m-1)\n"];
   Write[OUP,"\n"];
   For[JJ = 2,JJ <= M,JJ++,
      Write[OUP,SetPrecision[Im[c[JJ-1]],9]];
   ];
   Write[OUP,"\n"];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
]; 		
Quit[];