let c1 = Complex.one in
let fft =
let (+@) = Complex.add
and (-@) = Complex.sub
and ( *@) = Complex.mul in
fun a -> let n = Array.length a in
(let j = ref 0 in for i = 0 to n-2 do
if i < !j then (let t = a.(!j) in (a.(!j) <- a.(i); a.(i) <- t));
let m = ref (n/2) in while !m <= !j do j := !j - !m; m := !m/2 done;
j := !j + !m done;
let j = ref 1 in let b = ref (Complex.neg c1) in while !j ep)
or (abs_float (y -. (sin b)) > ep) then
(Printf.fprintf stdout "er: %19.16f + %19.16fi %4d\n " (cos b) (sin b) j;
pc a.(j)) done;;