MODE C = LONG COMPLEX;
PROC fft = (REF[]C a)VOID: (INT n = UPB a;
  (INT j := 1; FOR i TO n-1 DO (i<j|C t = a[j]; a[j] := a[i]; a[i] := t);
    (INT m := n OVER 2; WHILE m<j DO j -:= m; m := m OVER 2 OD; j +:= m) OD);
  (INT j := 1; C b := -1; WHILE j<n DO C w := 1;
    FOR m FROM 0 TO j-1 DO FOR i FROM m BY j+j TO n-1 DO C t = w*a[i+j+1];
      a[i+j+1] := a[i+1] - t; a[i+1] +:= t OD; w *:= b OD;
    j +:= j; b := long complex sqrt(b) OD));

PROC test = (INT n, INT m, REAL ep)VOID:
# n should be a power of 2 and 0<=m<n #
([n] C a;
FOR i TO n DO a[i] := (i=m+1|1|0) OD; fft(a);
FOR i TO n DO IF ABS (a[i] - long complex exp(2 * longpi *(0 I 1)*m*(i-1)/n)) > ep THEN
  print(("Bad", (i-1), newline)) FI OD;
print(("good", newline)));
test(1024, 51, 1.0e-6); test(8192, 4096, 1.0e-24); test(8192, 7013, 1.0e-24)