(*
 %use "/homes/cbj/public_html/FISh-1.1/Samples/Matrix/fft.fsh";;
*)

(* fft.fsh - modified version of Leroy's fft.ml *)

let pi = 3.14159265358979323846
;;

let tpi = 2.0 *. pi
;;

let fft px py np =
  new i = 2
  and m = 1 
  in

  while (!i < np) do 
    i := !i + !i; 
    m := !m + 1
  done ;

  new n = !i 
  in  
    if (not (n = np)) then (* checks that np is power of 2 *)
      let f = (fun t -> 0.0)
      in 
        update f px;
        update f py
    else skip ;

    new n2 = n+n 
    in
    for (1 <= k < !m) do 
      n2 := !n2 div 2; 
      new n4 = !n2 div 4 
      in
        new e = tpi /. (int2float(!n2))
        in
	  for (1 <= j < n4 + 1) do 
	      new a = e *. (int2float(j - 1))  
	      and a3 = 3.0 *. a  
              and cc1 = cos(a) 
              and ss1 = sin(a) 
              and cc3 = cos(a3)
              and ss3 = sin(a3)
              and is = j 
              and id = 2 * !n2
              in
	        while (!is < n) do 
		new i0r = !is
                  in 
	          while (!i0r < n) do 
		    new i0 = !i0r 
                    and i1 = i0 + n4 
	            and i2 = i1 + n4 
	            and i3 = i2 + n4 
	            and r1 = px[i0] -. px[i2] 
                    in
 		      px[i0] := px[i0] +. px[i2];
	              new r2 = px[i1] -. px[i3] 
                      in
	                px[i1] := px[i1] +. px[i3];
   	                new s1 = py[i0] -. py[i2] 
                        in
	                  py[i0] := py[i0] +. py[i2];
	                    new s2 = py[i1] -. py[i3] 
                            in
	                      py[i1] := py[i1] +. py[i3];
 			      new s3 = r1 -. s2 
			      and r1 = r1 +. s2 
			      and s2 = r2 -. s1 
		              and r2 = r2 +. s1 
		  	      in
				 px[i2] := r1 *. cc1 -. s2 *. ss1; 
		  		 py[i2] := (0.0 -. s2 *. cc1) -. r1 *. ss1;
		 	         px[i3] := s3 *. cc3 +. r2*.ss3;
			         py[i3] := r2 *. cc3 -. s3*.ss3;
			         i0r := i0 + !id
	                      end (* new s3,r1,s2,r2 *)
                            end (* new s2 *) 
                          end (* new s1 *)
	                end (* new r2 *)
	              end (* new i0,i1,i2,i3,r1 *)
                  done ; (* while !i0r *)
                  is := 2 * !id - !n2 + j;
                  id := 4 * !id
                end (* new i0r *)
	      done (* while is *)
            end (* new a, a3, cc1, ss1,cc3,ss3,is,id *)
          done
          end (* new e *)
        end (* new n4 *) 
      done
      end ; (* new n2 *)

     (************************************)
     (*  Last stage, length=2 butterfly  *)
     (************************************)

      new is = 1
      and id = 4
      in
        while (!is < n) do 
          new i0r = !is 
          in
            while (!i0r <= n) do 
              new i0 = !i0r 
              and i1 = i0 + 1
              and r1 = px[i0]
	      in  
	        px[i0] := r1 +. px[i1];
                px[i1] := r1 -. px[i1];
                new r1 = py[i0]
                in
                  py[i0] := r1 +. py[i1];
                  py[i1] := r1 -. py[i1];
                  i0r := i0 + !id
	        end (* new r1 *) 
              end (* new i0,i1,r1 *)
            done ;
            is := 2 * !id - 1; 
            id := 4 * !id
          end (* new i0r *)
        done 
        end ; (* new is, id *)

	(*************************)
	(*  Bit reverse counter  *)
	(*************************)

        new j = 1 
        in
        for (1 <= i < n) do
          if (i < !j) then 
            new xt = px[!j]
            in
              px[!j] := px[i]; 
              px[i] := xt;
              new xt = py[!j] 
              in
                py[!j] := py[i];
                py[i] := xt
              end (* new xt *) 
            end (* new xt *) 
          else skip ;
          new k = n div 2 
          in
            while (!k < !j) do 
              j := !j - !k; 
              k := !k div 2
            done ;
            j := !j + !k
          end (* new k *)
        done
      end (* new j *)
    end (* new n *) 
  end (* new i,m *)  
;;

(* 
let runtest (np:size) =
  new enp = int2float(np) 
  and npm = np div 2 - 1
  and #pxr = { np + 2 : float_shape } 
  and #pxi = #pxr 
  in
       update (fun t -> 0.0) pxr;
       update (fun t -> 0.0) pxi;
       new t = pi /. enp
       in 
         pxr[1] := (enp -. 1.0) *. 0.5;
         pxi[1] := 0.0;
         new n2 = np div 2 in
           pxr[n2+1] := -.0.5;
           pxi[n2+1] := 0.0
         end ; (* new n2 *)

         for (1 <= i < npm + 1) do 
           new j = np - i 
           in
             pxr[i+1] := -.0.5;
             pxr[j+1] := -.0.5;
             new z = t *. int2float(i) 
             and y = -.0.5 *. (cos(z)/.sin(z))
             in
                 pxi[i+1] := y;
                 pxi[j+1] := 0.0 -. y
             end (* new x,y *)
           end (* new j *)
         done ;

 	 fft pxr pxi np ; 

         new zr = 0.0 
         and zi = 0.0 
         and kr = 0 
         and ki = 0 
         in
           for (0 <= i < np) do 
             new a = fabs(pxr[i+1] -. int2float(i))
             in
               if (!zr <. a) 
               then 
                 zr := a; 
                 kr := i
               else 
                 skip
             end ; (* new a *)
             new a = fabs(pxi[i+1]) 
             in
               if (!zi <. a) 
               then 
                 zi := a; 
                 ki := i
               else
                 skip 		
             end (* new a *) 
           done
         end (* new zr,zi,kr,ki *) 
       end (* new t *)	
    end (* new enp,npm,pxr,pxi *)
;;

let fftest = 
   (runtest 16; 
   runtest 32;
   runtest 64;
   runtest 128;
   runtest 256;
   runtest 512;
   runtest 1024;
   runtest 2048;
   runtest 4096;
   runtest 8192;
   runtest 16384;
   runtest 32768;
   runtest 65536)
;; 

%run -c -q fftest
;;

*)