// This is pseudo-code for the fast fourier transform.
  // omega(n) denotes the principle n-th root of unity.
  // omega(i,n) denotes omega(n)**i.
  // Input condition: vals of omega(i,n) have been pre-computed for
  // all n which are powers of 2 (up to some maximum n) and all 0 <= i < n.
  // Note that omega(i,n) = omega(2*i,2*n).
 
 evaluate(int n, complex(n) coef) {
  // input condition: n is a power of 2
  // given the coefs of a polynomial f of degree less than n,
  // evaluate returns the vals {f(omega(i,n))} for 0 <= i < n
    pnt_val = new complex(n);
    if ( n == 1)
      pnt_val = coef;
    else
      {
        coef_even = new complex(n/2);
        coef_odd = new complex(n/2);
        for (int i = 0 , i++, i < n/2)
         {
           coef_even(i) = coef(2*i);
           coef_odd(i) = coef(2*i+1);
         }
        pnt_val_even = new complex(n/2) = evaluate(n/2,coef_even);
        pnt_val_odd = new complex(n/2) = evaluate(n/2,coef_odd);
        for (int i = 0, i++, i < n/2)
         {
           pnt_val(i) = pnt_val_even(i) + omega(i,n)*pnt_val_odd(i);
           pnt_val(n/2 + i) = pnt_val_even(i) - omega(i,n)*pnt_val_odd(i);
            // note that omega(n/2+i,n) == - omega(i,n)
         }
      }
    return pnt_val;
  }
 
 interpolate(int n, complex(n) pnt_val) {
  // input condition: n is a power of 2
  // given the point-values f(omega(i,n)), for 0 <= i < n, of a polynomial f
  // of degree less than n, interpolate the coefficients of f.
    coef = new complex(n);
    if ( n == 1)
      coef = pnt_val;
    else
      {
        pnt_val_even = new complex(n/2);
        pnt_val_odd = new complex(n/2);
        for (int i = 0 , i++, i < n/2)
         {
           pnt_val_even(i) = (pnt_val(i) + pnt_val(i + n/2))/2;
           pnt_val_odd(i) = omega((n-i)%n,n)*(pnt_val(i) - pnt_val(i + n/2))/2;
         }
        coef_even = new complex(n/2) = interpolate(n/2,pnt_val_even);
        coef_odd = new complex(n/2) = interpolate(n/2,pnt_val_odd);
        for (int i = 0, i++, i < n/2)
         {
           coef(2*i) = coef_even(i);
           coef(2*i+1) = coef_odd(i);
         }
      }
    return coef;
  }