User guide

13.6. FAST FOURIER TRANSFORM 221
LET start() = VALOF
{ writef("fft with N = %n and omega = %n modulus = %n*n*n",
N, omega, modulus)
data := getvec(upb)
UNLESS omega=1 DO // Unless doing Walsh tranform
check(omega, N) // check that omega and N are consistent
FOR i = 0 TO upb DO data!i := i
pr(data, 7)
// prints -- Original data
// 0 1 2 3 4 5 6 7
fft(data, ln, omega)
pr(data, 7)
// prints -- Transformed data
// 65017 26645 38448 37467 30114 19936 15550 42679
fft(data, ln, ovr(1,omega))
FOR i = 0 TO upb DO data!i := ovr(data!i, N)
pr(data, 7)
// prints -- Restored data
// 0 1 2 3 4 5 6 7
RESULTIS 0
}
AND fft(v, ln, w) BE // ln = log2 n w = nth root of unity
{ LET n = 1<<ln
LET vn = v+n
LET n2 = n>>1
// First do the perfect shuffle
reorder(v, n)
// Then do all the butterfly operations
FOR s = 1 TO ln DO
{ LET m = 1<<s
LET m2 = m>>1
LET wk, wkfac = 1, w
FOR i = s+1 TO ln DO wkfac := mul(wkfac, wkfac)
FOR j = 0 TO m2-1 DO
{ LET p = v+j
WHILE p<vn DO { butterfly(p, p+m2, wk); p := p+m }
wk := mul(wk, wkfac)
}
}
}
AND butterfly(p, q, wk) BE { LET a, b = !p, mul(!q, wk)
!p, !q := add(a, b), sub(a, b)
}