www.mapleapps.com/categories/mathematics/algebra/html/FFT-Polynomial.html
An Analogy *
A Series of Needed Procedures *
Initial Data *
The Traditional Multiplication of Two Polynomials *
Explanation of the DFT/FFT Method *
List Padding as Preparation for DFT/FFT * 8 The DFT Approach *
The FFT Approach * 10 Time Comparisons * 11 References POLYNOMIALS MULTIPLICATION USING THE FAST FOURIER TRANSFORM Bruno Guerrieri Florida A&M University We will be multiplying polynomials using two different ways, (traditional and FFT), and see whether one method is consistently faster than the other one. We will let the polynomials A(x) and B(x) be of the same degree for simplicity and we will generate the polynomials coefficients randomly. For a more detailed explanation of what follows, please consult
in the Reference section. An Analogy Consider the following table: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8182 16384 32768 You have, of course, noticed that entry i in row 1 corresponds to entry 2^i in row 2. Or, to put it another way, the top row lists the logarithms in base 2 of the bottom list. What if we wanted to calculate (64)(512)? One way would be to follow the traditional algorithm that we all learned in grammar school and perform the standard multiplication. The other way would be to locate the corresponding logarithms, namely 6 and 9, ADD these two logarithms to get 15, and then determine the antilog, our answer, namely 32768. In what follows, the strategy is the same. The quantities "worked on" and the operations are different. Instead of multiplying numbers, we will be multiplying polynomials. The equivalent of taking logarithms and antilogarithms will be polynomial evaluations and interpolations and the equivalent of our ADD will be POINTWISE multiplication. RETURN(c); The procedure "padding" "right-pads" a list of coefficients with "deg" zeros. RETURN(d); Vector "c" is obtained by "convolving" the input coefficients vectors "a" and "b". RETURN(c); The procedure "dft" (Discrete Fourier Transform) is present here since we wanted to, in fact compare the three processes for multiplication of two polynomials, namely the traditional, DFT, and FFT (Fast Fourier Transform) processes. One has to get into high degrees to see the FFT overtake the traditional method. The DFT is present here to make us appreciate the speed improvement that the FFT brings to the situation. The procedure "dft" which follows takes two arguments. The first one is the list to be transformed. The second argument is set so that a = 1 for DFT and a = -1 for IDFT (Inverse Discrete Fourier Transform). The last part of the code takes very small values (usually on the inverse transformation when dealing with original real data) located in the imaginary part and zero them out. So first the list is "split" into real and imaginary parts, Fourier transformation occurs, zeroing out is done where needed, and then the two lists are "sewn" back together. LS1 := nops(h_in); LS1 do > sum1 := 0; LS1 do > sum1 := sum1 + h_in
*exp(a*I*2*Pi*(j-1)*(k-1)/LS1); RETURN(h_out); The FFT/IFFT procedure call is to be found below. The only reason the procedure shown below looks lengthy is because real part and imaginary part of lists going through the FFT or the IFFT have to be kept separated. LS1,nombre,re_list,im_list,re_array,im_array,L_out; LS1 := nops(L_in); FFT(nombre,re_array,im_array) else iFFT(nombre,re_array,im_array); L_out := zip((a,b)->simplify((a+b*I)),re_list,im_list); RETURN(L_out); The "power_of_two" procedure determines the nearest power greater than deg(A(x)) + deg(B(x) + 1. This is needed, for padding purposes, since the FFT/IFFT needs "power of two" list sizes. RETURN(2^i); The " zeros_plot" procedure "reconstructs" a polynomial from its coefficients so that its zeros can be determined and plotted if necessary. N := nops(coef_list); N do > y := y + coef_list
*x^(i-1); R := fsolve(y,x,complex) ; Re(z),Im(z) ,R); RETURN(p); Initial Data Now that all the procedures we will need have been written, let us randomly generate the "deg+1" coefficients of two polynomials A(x) and B(x) of degree "deg". For simplicity, we will take the two polynomials to be of the same degree. Otherwise, we would have to pad with zeros. We are saving these two lists, in c_list and d_list, so that we can use the same polynomials with the DFT and then the FFT approaches. The Traditional Multiplication of Two Polynomials In the straightforward multiplication (convolution) method that follows, list sizes are handled automatically. We will notice, when we use the DFT approach that the padding will always be up to a power of 2. The final answer provides a list of the coefficients of the resulting polynomials C(x). Explanation of the DFT/FFT Method Again, for a more detailed explanation we recommend that you refer to
. We would like to multiply the same two polynomials seen above using the DFT. Very quickly, we will come to realize that this process requires a lot of computing time. In fact the DFT is no match for the traditional method. We wanted to investigate the DFT approach for the simple reason that it is a routine that you can write on your own by using the definition for the DFT. This way, anybody unfamiliar with the intricacies of the FFT, can feel somewhat at ease when told that the FFT is nothing but a "souped up" DFT - (Apologies to Cooley and Tukey
). Therefore we recommend that you bypass the DFT process and go directly to the FFT one whenever your polynomial degrees start to get above 50. Note that the polynomial representation we used above, where polynomials are identified by their coefficient lists, is called the "coefficient representation (CR)". One can also use the "point-value representation (PVR)". While multiplying polynomials in CR form takes O( n^2 ), the alternative multiplicative option of using the three steps of a) conversion to PVR, b) pointwise multiplication, c) reconversion to CR only takes O( n*log(n) ). The evaluation of a polynomial at a particular value of x is usually carried out using Horner's method. PVR form by relying on Lagrange polynomials. Thus n -points evaluation and interpolation are therefore well-defined operations that allow one to move from CR to PVR and back again. Let us assume in the rest of the discussion that n = deg( A(x) = deg( B(x) ). Since the resulting polynomial is of degree 2*n , we must start with "extended" 2*n+1 PVR form for both polynomials A and B. The quantity (+ 1 ) that appears at times has to do with the fact that a polynomial of degree n has n+1 coefficients. We must also, since the FFT is at its best when processing lists whose lengths are powers of 2, let 2 n be a power of 2. Be aware that there are two sorts of padding going on. One was seen at work in the traditional multiplication process. But the one we are referring to now has to do with bringing all lists to lengths that are powers of 2. One can be clever in the choice of the x
's (there will be 2 n of them) at which polynomials will be evaluated so as to determine the corresponding PVR form. In fact, if we choose "complex roots of unity" as evaluation points, we can produce a PVR form of a given polynomial in CR form by taking the DFT of the coefficient vector. IDFT takes care of interpolation so as to bring us back to a CR form, from a PVR form. Here is the algorithm we will be following : Start with A(x) and B(x) in coefficient (CR) form. Create coefficient representation of A(x) and B(x) as polynomials of degree 2 n . Time = O( n ) . Evaluate: Compute point-value representation of A(x) and B(x) of length 2*n through applications of the DFT (later the FFT) of order 2*n to each coefficient vector. These representations contains the values of the two polynomials at the ( 2*n )th roots of unity. Time = O( nlog(n) ) . The representation obtained contains the value of C(x) at each ( 2*n )th roots of unity. Time = O( n ) . Interpolate: Create the coefficient representation of the polynomial C(x) through a single application of a IDFT (later IFFT) on 2*n point-value pairs. Time = O( nlog(n) ) . What has been done is as follows: {a} convolved with {b} = IDFT DFT(a)*DFT(b) , where vectors {a} and {b} have been padded wi...
|
http://csua.com/feed/