(* 
%use "Matrix/fourier.fsh";; 
*)

%use "Preludes/complex_prelude.fsh";;
%use "BMF/scan.fsh";;

	(* nth root of unity *)

let root (n:size) = 
  let pi = 3.1415927 in 
    from_mod_arg 
    	1.0 
	(2.0 *. pi /. (int2float n)) 
;;  


let fourier_pr (z:var b) (x: var a) = 
  let n = lendim #x 
  and omega = root n in 
  new twiddle = scan times_complex one_complex (copy {n:float_shape} omega) 
  in 
     z := copy {n:int_shape} zero_complex ;
     for i < n do 
	for j < n do
	  z[i] := plus_complex z[i] (times_complex twiddle[i*j mod n] x[j])
	done
     done
  end
;;

let fourier = proc2fun fourier_pr identity ;;