-- Copyright 1993-1998, by the Cecil Project
-- Department of Computer Science and Engineering, University of Washington
-- See the LICENSE file for license information.

-- This code is an fft benchmark; adapted from the C code for the
-- Stanford integer benchmarks collected/developed originally as part of
-- the MIPS project

include "prelude.small.cecil";
include "float-vector.cecil";

object complex isa num;
  field real(@:complex):num;
  field imag(@:complex):num;

  method new_complex(r:num, i:num):complex {
      concrete object isa complex { real := r, imag := i } }

  method print_string(c@:complex):string {
      [c.real.print_string,
       if(c.imag >= 0, { "+" }, { "" }),
       c.imag.print_string,
       "i"].flatten }

  method +(c1@:complex, c2@:complex):complex {
      new_complex(c1.real + c2.real, c1.imag + c2.imag) }
  method -(c1@:complex, c2@:complex):complex {
      new_complex(c1.real - c2.real, c1.imag - c2.imag) }

  method *(c1@:complex, c2@:complex):complex {
      new_complex(c1.real * c2.real - c1.imag * c2.imag,
		  c1.real * c2.imag + c1.imag * c2.real) }
  method *(n@:num, c2@:complex):complex {
      new_complex(n * c2.real, n * c2.imag) }
  method *(c@:complex, n@:num):complex { n * c }


-- computes cosine of a float by an expansion
method cos(x@:float):float {
  let var result:float := 1.0;
  let var factor:int := 1;
  let var power:float := x;
  for(2, 10, &(i:int){
      factor := factor * i;
      power := power * x;
      if(i _bit_and 1 == 0, {
	  if(i _bit_and 3 == 0, {
	      result := result + power/factor;
	  }, {
	      result := result - power/factor;
	  })
      });
  });
  result }


-- code to generate a random sequence of float values
let var seed:int := 0;

method uniform11():float {
    seed := (4855*seed + 1731) _bit_and 8191;
    seed/8192.0 }


-- code to create the array of complexes
method exptab(n:int, n2:int):indexed[complex] {
    let theta:float := 3.1415927 (-- 3.1415925536 --);
    let var divisor:float := 4.0;

    let h:indexed[float] :=
      new_i_float_vector_init(25, &(i:int){
	  let elem := 1/(2*cos(theta/divisor));
	  divisor := divisor + divisor;
	  elem });

    let m:int := n/2;
    let var l:int := m/2;
    let var j:int := 1;

    let e:m_indexed[complex] := new_m_vector[complex](n2);
    
    e!0 := new_complex(1.0, 0.0);
    e!l := new_complex(0.0, 1.0);
    e!m := new_complex(-1.0, 0.0);

    until_false({
	let i:int := l/2;
	let var k:int := i;

	until_false({
	    e!k := (h!(j-1)) * (e!(k+i) + e!(k-i));
	    k := k + l;
	}, { k <= m });

	j := min(j+1, 25);
	l := i;
    }, { l > 1 });

    e }

-- do the fft calculation, whatever it is
let var steps:int := 0;

method fft(n:int, z:m_indexed[complex], e:indexed[complex], sqrinv:float):void {
    let m:int := n / 2;
    let var l:int := 1;

    until_false({
	let var k:int := 0;
	let var j:int := l;
	let var i:int := 0;

	let w:m_indexed[complex] := new_m_vector[complex](n);
	
	until_false({
	    until_false({
		w!(i+k) := z!i + z!m;
		w!(i+j) := e!k * (z!i - z!(i+m));
		i := i+1;
		steps := steps+1;
	    }, { i < j });
	    k := j;
	    j := k+l;
	}, { j <= m });

	-- z := w.copy, but in place
	n.do(&(index:int){
	    z!index := w!index;
	});

	l := l+l;
    }, { l <= m });

    n.do(&(i){
	z!i := fft_op(z!i, sqrinv);
    });
}

method fft_op(c:complex, sqrinv:num):complex {
    new_complex(c.real * sqrinv, c.imag * -sqrinv) }


-- the main routine, which creates the arrays and invokes the fft algorithm

let fftsize:int := 256;
let fftsize2:int := 129;

method run_fft():void {
    let e:indexed[complex] := exptab(fftsize, fftsize2);
    -- e.print_line;
    seed := 5767;
    let z:m_indexed[complex] := new_m_vector_init[complex](fftsize, &(i:int){
	new_complex(20.0 * uniform11() - 10.0, 20.0 * uniform11() - 10.0) });
    -- z.print_line;
    20.do(&(i:int){
	fft(fftsize, z, e, 0.0625);
	-- z.print_line;
    });
}

method fft() {
    "** Starting fft...".print;
    steps := 0;
    let time:int := time({
	run_fft();
    });
    ("done: " || steps.print_string || " steps; time: " ||
     time.print_string || " ms.").print_line;
}

fft();