/*
** by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
** Rewritten into Xanadu by H. Xi, 23 May 1999
*/

/********************************************************/
/*     A Duhamel-Hollman split-radix dif fft            */
/*     Ref: Electronics Letters, Jan. 5, 1984           */
/*     Complex input and output data in arrays x and y  */
/*     Length is n.                                     */
/********************************************************/

global pi: float = 3.1415926536

{size:int | size >= 2}
int dfft(px[size+1]: float, py[size+1]: float, size: int(size))
{
  var: [i:int | i > 0] int(i) i,j,is,id,m,n;
       [i:nat] int(i) k,n2,n4;
       int i0,i1,i2,i3,tmp;
       float a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt;;

  i = 2; m = 1; 

  while (i < size) { i = i+i; m = m+1; }
  
  n = i; 
  
  assert (n == size);

  n2 = 2 * n;

  for (k = 1; k <= m-1; k = k+1) {

/*
    n2 = n2 / 2; n4 = n2 / 4; 
*/

    n4 = n2 / 8; n2 = 4 * n4;

    e  = 2.0 *. pi /. float_of_int(n2);
    a  = 0.0;

    for (j = 1; j <= n4 ; j = j+1) {

      a3 = 3.0 *. a;
      cc1 = cos(a); ss1 = sin(a);      
      cc3 = cos(a3); ss3 = sin(a3);
      a = e *. float_of_int(j); 
      is = j; id = 2 * n2;

      while (is < n) {
	  i0 = is; i1 = i0 + n4; i2 = i1 + n4; i3 = i2 + n4;

          invariant:
            [i0:int,i1:int,i2:int,i3:int | 0 < i0 <= i1 <= i2 <= i3]
	    (i0: int(i0), i1: int(i1), i2: int(i2), i3: int(i3))
	  while (i3 <= n) {
	    r1 = px[i0] -. px[i2];
	    px[i0] = px[i0] +. px[i2];
	    r2 = px[i1] -. px[i3];
	    px[i1] = px[i1] +. px[i3];
	    s1 = py[i0] -. py[i2];
	    py[i0] = py[i0] +. py[i2];
	    s2 = py[i1] -. py[i3];
	    py[i1] = py[i1] +. py[i3];
	    s3 = r1 -. s2; r1 = r1 +. s2; 
	    s2 = r2 -. s1; r2 = r2 +. s1;
	    px[i2] = r1 *. cc1 -. s2 *. ss1; 
	    py[i2] =  -.s2 *. cc1 -. r1 *. ss1;
	    px[i3] = s3 *.cc3 +. r2 *. ss3;
	    py[i3] = r2 *. cc3 -. s3 *. ss3;
	    i0 = i0 + n4; i1 = i1 + n4; i2 = i2 + n4; i3 = i3 + n4;
	  }

          tmp = 2 * id - n2 + j; 

	  assert (tmp > 0); is = tmp; 

	  id = 4 * id;
      }
    }
  }

/************************************/
/*  Last stage, length=2 butterfly  */
/************************************/
  is = 1; id = 4;
  
  while (is < n) {
    invariant: [i0:int | 0 < i0] (i0: int(i0))
    for (i0 = is; i0 < n; i0 = i0 + id) {
      i1 = i0 + 1; 
      r1 = px[i0];
      px[i0] = r1 +. px[i1];
      px[i1] = r1 -. px[i1];
      r1 = py[i0];
      py[i0] = r1 +. py[i1];
      py[i1] = r1 -. py[i1];
    }

    is = 2 * id - 1; id = 4 * id; 
  }


/*************************/
/*  Bit reverse counter  */
/*************************/
  j = 1;
  
  invariant: [j:nat | j <= size] (j: int(j))
  for (i = 1; i < n; i = i+1) {

    if (i < j) {
      xt = px[j];
      px[j] = px[i]; 
      px[i] = xt;
      xt = py[j]; 
      py[j] = py[i];
      py[i] = xt;
    }
    
    k = n / 2;
 
    invariant: [k:nat | 2 * k <= size]
               /* {metric: max(j-k)} */
               (k: int(k))
    while (k < j) {
        /* assert (k > 0); */
	j = j - k; k = k / 2;
    }

    j = j + k;
  }


/*
  for (i = 0; i < 16; i = i + 1) printf("%d  %g   %gn",i, px[i], py[i]);
*/

  return(n);
}