fft: proc(a); dcl (a(*), t, (b, w) binary(40)) float complex;
dcl i, j, n, m, (sqrt, hbound) builtin;
 n = hbound(a, 1);
 j = 1;
twist: do i = 1 to n-1;
  if i < j then do; t =a(j); a(j) = a(i); a(i) = t; end;
  m = n/2; do while(m<j); j = j-m; m = m/2; end; j = j+m;
end twist;

j = 1; b = -1;
do while(j < n);
 w = 1;
 do m = 1 to j;
  do i = m to n by j*2;
    t = w*a(i+j); a(i+j) = a(i) - t; a(i) = a(i) + t;
   end; w = w*b;
  end; j = j*2; b = sqrt(b);
 end;
end fft;