The following implementation uses a nonrecursive Fast Fourier Transformation algorithm to multiply polynomials. If we view the elements of a vector (a0,...,an-1) to be transformed as coefficients of a polynomial a(x), then the Fourier transform is the same as computing a(wj) for 0 <= j < N. This is also equivalent to computing the remainder by making use of a principle called balancing. In this case we compute the remainders with the help of a process which is structured like a binary tree. Since the complex description of this algorithm is beyond the scope of this document, refer to the bibliography for further information.
1 procedure fft(arr: in out coeff_array; m: Integer) 2 is 3 i,j,k,l,n,ndiv2,pow2,pow2m1,index : Integer; 4 r,s,t : complex; 5 begin 6 n := 2**m; 7 ndiv2 := n/2; 8 j := 1; 9 for i in 1..n-1 loop 11 if i<j then 12 t := arr(j); 13 arr(j) := arr(i); 14 arr(i) := t; 15 end if; 16 k := ndiv2; 17 discrete a:=k in reverse j..k 18 new a:=a/2 19 loop 20 j := j-k; 21 a := a/2; 22 k := a; 23 end loop; 24 j := j+k; 25 end loop; 26 for l in 1..m loop 27 pow2 := 2**l; 28 pow2m1 := pow2/2; 29 r := (1.0,0.0); 30 s := (cos(pi/double(pow2m1)),sin(pi/double(pow2m1))); 31 for j in 1..pow2m1 loop 32 discrete i:=j in j..n 33 new i:=i+pow2 34 loop 35 index := i+pow2m1; 36 t := arr(index)*r; 37 arr(index) := arr(i)-t; 38 arr(i) := arr(i)+t; 39 i := i+pow2; 40 end loop; 41 r := r*s; 42 end loop; 43 end loop; 44 end fft;
[ Sourcefile ]
In our example above we assume complex arithmetic and w=e2*PI*i/j is expressed in terms of sines and cosines. 'arr' contains the input coefficients and the elements of the transform are computed in-place.
The complexity of the non-recursive Fast Fourier Transform is
N log N.
[Sav76, p. 342]
[ Samples Overview | Bibliography | WOOP Project ]