Fast Fourier Transformation for poynomial multiplication
Given two polynomial A(x) and B(x), find the product C(x) = A(x)*B(x). There is already an O() naive approach to solve this problem. here. This approach uses the coefficient form of the polynomial to calculate the product.
A coefficient representation of a polynomial is a = a0, a1, …, an1.
• Example
Coefficient representation of A(x) = (9, 10, 7, 6)
Coefficient representation of B(x) = (5, 4, 0, 2)
Input : A[] = {9, 10, 7, 6} B[] = {5, 4, 0, 2} Output :
We can do better, if we represent the polynomial in another form.
Idea is to represent polynomial in pointvalue form and then compute the product. A pointvalue representation of a polynomial A(x) of degreebound n is a set of n pointvalue pairs is{ (x0, y0), (x1, y1), …, (xn1, yn1)} such that all of the xi are distinct and yi = A(xi) for i = 0, 1, …, n1.
Example
xi  0, 1, 2, 3 A(xi)  1, 0, 5, 22
Pointvalue representation of above polynomial is { (0, 1), (1, 0), (2, 5), (3, 22) }. Using Horner’s method, (discussed here), npoint evaluation takes time O(). It’s just calculation of values of A(x) at some x for n different points, so time complexity is O(). Now that the polynomial is converted into point value, it can be easily calculated C(x) = A(x)*B(x) again using horner’s method. This takes O(n) time. An important point here is C(x) has degree bound 2n, then n points will give only n points of C(x), so for that case we need 2n different values of x to calculate 2n different values of y. Now that the product is calculated, the answer can to be converted back into coefficient vector form. To get back to coefficient vector form we use inverse of this evaluation. The inverse of evaluation is called interpolation. Interpolation using Lagrange’s formula gives point valueform to coefficient vector form of the polynomial.Lagrange’s formula is –
So far we discussed,
.
This idea still solves the problem in O() time complexity. We can use any points we want as evaluation points, but by choosing the evaluation points carefully, we can convert between representations in only O(n log n) time. If we choose “complex roots of unity” as the evaluation points, we can produce a pointvalue representation by taking the discrete Fourier transform (DFT) of a coefficient vector. We can perform the inverse operation, interpolation, by taking the “inverse DFT” of pointvalue pairs, yielding a coefficient vector. Fast Fourier Transform (FFT) can perform DFT and inverse DFT in time O(nlogn).
DFT
DFT is evaluating values of polynomial at n complex nth roots of unity . So, for k = 0, 1, 2, …, n1, y = (y0, y1, y2, …, yn1) is Discrete fourier Transformation (DFT) of given polynomial.
The product of two polynomials of degreebound n is a polynomial of degreebound 2n. Before evaluating the input polynomials A and B, therefore, we first double their degreebounds to 2n by adding n highorder coefficients of 0. Because the vectors have 2n elements, we use “complex 2nth roots of unity, ” which are denoted by the W2n (omega 2n). We assume that n is a power of 2; we can always meet this requirement by adding highorder zero coefficients.
FFT
Here is the Divideandconquer strategy to solve this problem.

Define two new polynomials of degreebound n/2, using evenindex and oddindex coefficients of A(x) separately

The problem of evaluating A(x) at reduces to evaluating the degreebound n/2 polynomials A0(x) and A1(x) at the points
 Iterative Fast Fourier Transformation for polynomial multiplication
 Exponential Squaring (Fast Modulo Multiplication)
 Karatsuba algorithm for fast multiplication using Divide and Conquer algorithm
 Window to Viewport Transformation in Computer Graphics with Implementation
 Fast inverse square root
 How does Floyd's slow and fast pointers approach work?
 Fast method to calculate inverse square root of a floating point number in IEEE 754 format
 Booth’s Multiplication Algorithm
 Sum of multiplication of triplet of divisors of a number
 Right most nonzero digit in multiplication of array elements
 How to avoid overflow in modular multiplication?
 Count divisors of array multiplication
 Check for integer overflow on multiplication
 Multiplication of two complex numbers given as strings
 Program to print multiplication table of a number
 Now combining the results by
The list does not contain n distinct values, but n/2 complex n/2th roots of unity. Polynomials A0 and A1 are recursively evaluated at the n complex nth roots of unity. Subproblems have exactly the same form as the original problem, but are half the size. So recurrence formed is T(n) = 2T(n/2) + O(n), i.e complexity O(nlogn).
Algorithm 1. Add n higherorder zero coefficients to A(x) and B(x) 2. Evaluate A(x) and B(x) using FFT for 2n points 3. Pointwise multiplication of pointvalue forms 4. Interpolate C(x) using FFT to compute inverse DFT
Pseudo code of recursive FFT
Recursive_FFT(a){ n = length(a) // a is the input coefficient vector if n = 1 then return a // wn is principle complex nth root of unity. wn = e^(2*pi*i/n) w = 1 // even indexed coefficients A0 = (a0, a2, ..., an2 ) // odd indexed coefficients A1 = (a1, a3, ..., an1 ) y0 = Recursive_FFT(A0) // local array y1 = RecursiveFFT(A1) // local array for k = 0 to n/2  1 // y array stores values of the DFT // of given polynomial. do y[k] = y0[k] + w*y1[k] y[k+(n/2)] = y0[k]  w*y1[k] w = w*wn return y } Recursion Tree of Above Execution
Why does this work ?
since,
Thus, the vector y returned by RecursiveFFT is indeed the DFT of the input
vector a.
#include <bits/stdc++.h> using namespace std; // For storing complex values of nth roots // of unity we use complex<double> typedef complex< double > cd; // Recursive function of FFT vector<cd> fft(vector<cd>& a) { int n = a.size(); // if input contains just one element if (n == 1) return vector<cd>(1, a[0]); // For storing n complex nth roots of unity vector<cd> w(n); for ( int i = 0; i < n; i++) { double alpha = 2 * M_PI * i / n; w[i] = cd( cos (alpha), sin (alpha)); } vector<cd> A0(n / 2), A1(n / 2); for ( int i = 0; i < n / 2; i++) { // even indexed coefficients A0[i] = a[i * 2]; // odd indexed coefficients A1[i] = a[i * 2 + 1]; } // Recursive call for even indexed coefficients vector<cd> y0 = fft(A0); // Recursive call for odd indexed coefficients vector<cd> y1 = fft(A1); // for storing values of y0, y1, y2, ..., yn1. vector<cd> y(n); for ( int k = 0; k < n / 2; k++) { y[k] = y0[k] + w[k] * y1[k]; y[k + n / 2] = y0[k]  w[k] * y1[k]; } return y; } // Driver code int main() { vector<cd> a{1, 2, 3, 4}; vector<cd> b = fft(a); for ( int i = 0; i < 4; i++) cout << b[i] << endl; return 0; } 
Input: 1 2 3 4 Output: (10, 0) (2, 2) (2, 0) (2, 2)
Interpolation
Switch the roles of a and y.
Replace wn by wn^1.
Divide each element of the result by n.
Time Complexity : O(nlogn).
Recommended Posts:
If you like GeeksforGeeks and would like to contribute, you can also write an article using contribute.geeksforgeeks.org or mail your article to contribute@geeksforgeeks.org. See your article appearing on the GeeksforGeeks main page and help other Geeks.
Please Improve this article if you find anything incorrect by clicking on the "Improve Article" button below.