Given two polynomial A(x) and B(x), find the product C(x) = A(x)*B(x). There is already an O(
A coefficient representation of a polynomial
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 point-value form and then compute the product. A point-value representation of a polynomial A(x) of degree-bound n is a set of n point-value pairs is{ (x0, y0), (x1, y1), …, (xn-1, yn-1)} such that all of the xi are distinct and yi = A(xi) for i = 0, 1, …, n-1.
Example
xi -- 0, 1, 2, 3 A(xi) -- 1, 0, 5, 22
Point-value representation of above polynomial is { (0, 1), (1, 0), (2, 5), (3, 22) }. Using Horner’s method, (discussed here), n-point evaluation takes time O(
So far we discussed,
.
This idea still solves the problem in O(
DFT
DFT is evaluating values of polynomial at n complex nth roots of unity
The product of two polynomials of degree-bound n is a polynomial of degree-bound 2n. Before evaluating the input polynomials A and B, therefore, we first double their degree-bounds to 2n by adding n high-order 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 high-order zero coefficients.
FFT
Here is the Divide-and-conquer strategy to solve this problem.
Define two new polynomials of degree-bound n/2, using even-index and odd-index coefficients of A(x) separately
The problem of evaluating A(x) at
Now combining the results by
The list
Algorithm 1. Add n higher-order zero coefficients to A(x) and B(x) 2. Evaluate A(x) and B(x) using FFT for 2n points 3. Pointwise multiplication of point-value 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, ..., an-2 ) // odd indexed coefficients A1 = (a1, a3, ..., an-1 ) y0 = Recursive_FFT(A0) // local array y1 = Recursive-FFT(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 Recursive-FFT 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, ..., yn-1.
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;
} |
from math import sin,cos,pi
# Recursive function of FFT def fft(a):
n = len (a)
# if input contains just one element
if n = = 1 :
return [a[ 0 ]]
# For storing n complex nth roots of unity
theta = - 2 * pi / n
w = list ( complex (cos(theta * i), sin(theta * i)) for i in range (n) )
# Separe coefficients
Aeven = a[ 0 :: 2 ]
Aodd = a[ 1 :: 2 ]
# Recursive call for even indexed coefficients
Yeven = fft(Aeven)
# Recursive call for odd indexed coefficients
Yodd = fft(Aodd)
# for storing values of y0, y1, y2, ..., yn-1.
Y = [ 0 ] * n
middle = n / / 2
for k in range (n / / 2 ):
w_yodd_k = w[k] * Yodd[k]
yeven_k = Yeven[k]
Y[k] = yeven_k + w_yodd_k
Y[k + middle] = yeven_k - w_yodd_k
return Y
# Driver code if __name__ = = '__main__' :
a = [ 1 , 2 , 3 , 4 ]
b = fft(a)
for B in b:
print (B)
|
// JavaScript program to implement // the approach // For storing complex values of nth roots // of unity we use complex class complex { constructor(a, b = 0)
{
this .x = a;
this .y = b;
}
} function product(c1, c2)
{ let c = new complex(0, 0);
c.x = c1.x * c2.x - c1.y * c2.y
c.y = c1.x * c2.y + c2.x * c1.y
return c
} function sum(c1, c2)
{ let c = new complex(0, 0);
c.x = c1.x + c2.x
c.y = c1.y + c2.y
return c
} function difference(c1, c2)
{ let c = new complex(0, 0);
c.x = c1.x - c2.x
c.y = c1.y - c2.y
return c
} // Recursive function of FFT function fft(a)
{ let n = a.length;
// if input contains just one element
if (n == 1)
return [a[0]]
// For storing n complex nth roots of unity
let w = new Array(n);
let alpha = -2 * Math.PI / n;
for ( var i = 0; i < n; i++) {
w[i] = new complex(Math.cos(alpha * i), Math.sin(alpha * i));
}
let A0 = new Array(Math.floor(n / 2));
let A1 = new Array(Math.floor(n / 2));
for ( var i = 0; i < Math.floor(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
let y0 = fft(A0);
// Recursive call for odd indexed coefficients
let y1 = fft(A1);
// for storing values of y0, y1, y2, ..., yn-1.
let y = new Array(n);
for ( var k = 0; k < Math.floor(n / 2); k++) {
y[k] = sum(y0[k], product(w[k], y1[k]));
y[k + Math.floor(n / 2)] = difference(y0[k], product(w[k], y1[k]));
}
return y;
} // Driver code let a = [ new complex(1, 0), new complex(2, 0), new complex(3, 0), new complex(4, 0)];
let b = fft(a); for ( var b0 of b)
console.log( "(" , b0.x, "," , b0.y, ")" )
// This code is contributed by phasing17 |
import java.util.*;
// For storing complex values of nth roots // of unity we use complex<double> class cd {
double re;
double im;
cd( double r, double i) {
re = r;
im = i;
}
cd add(cd b) {
return new cd(re + b.re, im + b.im);
}
cd sub(cd b) {
return new cd(re - b.re, im - b.im);
}
cd mul(cd b) {
return new cd(re * b.re - im * b.im, re * b.im + im * b.re);
}
} // Recursive function of FFT class FFT {
int n;
cd[] a, w;
FFT( int n) {
this .n = n;
a = new cd[n];
w = new cd[n];
}
void fft() {
// if input contains just one element
if (n == 1 )
return ;
for ( int i = 0 ; i < n; i++) {
double alpha = - 2 * Math.PI * i / n;
w[i] = new cd(Math.cos(alpha), Math.sin(alpha));
}
cd[] A0 = new cd[n / 2 ];
cd[] A1 = new cd[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
FFT y0 = new FFT(n / 2 );
y0.a = A0;
y0.fft();
// Recursive call for odd indexed coefficients
FFT y1 = new FFT(n / 2 );
y1.a = A1;
y1.fft();
// for storing values of y0, y1, y2, ..., yn-1.
cd[] y = new cd[n];
for ( int k = 0 ; k < n / 2 ; k++) {
y[k] = y0.a[k].add(w[k].mul(y1.a[k]));
y[k + n / 2 ] = y0.a[k].sub(w[k].mul(y1.a[k]));
}
System.arraycopy(y, 0 , a, 0 , n);
}
// } // Driver code // public class Main { public static void main(String[] args) {
FFT fft = new FFT( 4 );
fft.a[ 0 ] = new cd( 1 , 0 );
fft.a[ 1 ] = new cd( 2 , 0 );
fft.a[ 2 ] = new cd( 3 , 0 );
fft.a[ 3 ] = new cd( 4 , 0 );
fft.fft();
for ( int i = 0 ; i < 4 ; i++)
System.out.println(fft.a[i].re + " + " + fft.a[i].im + "i" );
}
} |
using System;
// For storing complex values of nth roots // of unity we use a struct named Complex struct Complex
{ public double Re;
public double Im;
public Complex( double r, double i)
{
Re = r;
Im = i;
}
public Complex Add(Complex b)
{
return new Complex(Re + b.Re, Im + b.Im);
}
public Complex Sub(Complex b)
{
return new Complex(Re - b.Re, Im - b.Im);
}
public Complex Mul(Complex b)
{
return new Complex(Re * b.Re - Im * b.Im,
Re * b.Im + Im * b.Re);
}
} // Recursive function of FFT class FFT {
int n;
Complex[] a, w;
public FFT( int n)
{
this .n = n;
a = new Complex[n];
w = new Complex[n];
}
public void Fft()
{
// if input contains just one element
if (n == 1)
return ;
for ( int i = 0; i < n; i++) {
double alpha = -2 * Math.PI * i / n;
w[i] = new Complex(Math.Cos(alpha),
Math.Sin(alpha));
}
Complex[] A0 = new Complex[n / 2];
Complex[] A1 = new Complex[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
FFT y0 = new FFT(n / 2);
y0.a = A0;
y0.Fft();
// Recursive call for odd indexed coefficients
FFT y1 = new FFT(n / 2);
y1.a = A1;
y1.Fft();
// for storing values of y0, y1, y2, ..., yn-1.
Complex[] y = new Complex[n];
for ( int k = 0; k < n / 2; k++) {
y[k] = y0.a[k].Add(w[k].Mul(y1.a[k]));
y[k + n / 2] = y0.a[k].Sub(w[k].Mul(y1.a[k]));
}
Array.Copy(y, a, n);
}
public static void Main()
{
FFT fft = new FFT(4);
fft.a[0] = new Complex(1, 0);
fft.a[1] = new Complex(2, 0);
fft.a[2] = new Complex(3, 0);
fft.a[3] = new Complex(4, 0);
fft.Fft();
// css
// Copy code
for ( int i = 0; i < 4; i++)
Console.WriteLine(fft.a[i].Re + " + "
+ fft.a[i].Im + "i" );
}
} |
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).