//
// Assignment 2 : Program to compute the product of two polynomials in O(nlogn) time using FFT.
//
// Author : Saket Soni, M.Tech. CS, Sem II, Session 2004-06
// Roll : mtc0420
// Subject : Design and Analysis of Algorithms
// Submitted to : Shri. Krishnendu Mukhopadhyay
// Date : March 14, 2005
///////////////////////////////////////////////////////////////////////////////////////////////
// Header Files ...
///////////////////////////////////////////////////////////////////////////////////////////////
//
#include
#include
#include
#include
#include
using namespace std;
///////////////////////////////////////////////////////////////////////////////////////////////
// Constant Declarations ...
///////////////////////////////////////////////////////////////////////////////////////////////
//
const double PI = 3.1415926535;
///////////////////////////////////////////////////////////////////////////////////////////////
// Function Declarations ...
///////////////////////////////////////////////////////////////////////////////////////////////
//
void fft(int size_of_arr, int distance_between_consecutive_elements,
complex array_of_polynomial_coefficients_with_lowest_order_term_first [],
complex nth_primitive_root_of_unity,
complex output_array []);
// computes the fft of the array_of polynomial_coefficients and stores the output
// int output_array, size_of_arr is the number of the terms in the polynomial,
// distance_between_consecutive_elements is the as the name suggests,
// nth primitive root of unity is required by the function to work correctly.
void multPolyFft(int size_of_arrays,
complex a[], complex b[],
complex output_product_of_a_and_b[]);
// function to do point by point complex multiplication of complex vectors a and b
// and generate result in c
complex evalPoly(int size_of_arr,
complex array_of_polynomial_coefficients_with_lowest_order_term_first [],
complex x);
// function to evaluate a polynomial at x
///////////////////////////////////////////////////////////////////////////////////////////////
// Function Definitions ...
///////////////////////////////////////////////////////////////////////////////////////////////
//
void fft(int n, int j, complex arr[], complex w, complex output[])
{
if ( n == 1 )
output[0] = arr[0];
else
{
complex *outE, *outO, temp;
int nn = n>>1;
outE = new complex[nn]; // for storing the output
outO = new complex[nn];
temp = w*w;
fft(nn,j<<1,arr,temp,outE); // computing fft of the even terms
fft(nn,j<<1,arr+j,temp,outO); // computing fft of the odd terms
for ( int i = 0 ; i < nn ; i++ )
{
temp = pow(w,i);
temp = temp*outO[i];
output[i] = outE[i] + temp; // generating the output
output[i+nn] = outE[i] - temp;
}
delete[] outE;
delete[] outO;
}
}
void multPolyFft(int n, complex a[], complex b[], complex c[])
{
for ( int i = 0 ; i < n ; i++ ) // point by point complex vector multiplication
c[i] = a[i]*b[i];
}
complex evalPoly(int n, complex arr[], complex x)
{
complex result = 0; // polynomial evaluation according to
for ( int i = 0 ; i < n ; i++ ) // Horner's rule
result = (result*x) + arr[n-1-i];
return result;
}
///////////////////////////////////////////////////////////////////////////////////////////////
// Application Entry Point ...
///////////////////////////////////////////////////////////////////////////////////////////////
//
void main()
{
// variable declarations ...
int i, n1, n2, n, nn;
double temp;
complex *A, *B, *C, *fftA, *fftB, *fftC, w;
// setting output format flags ...
cout << setiosflags(ios::scientific) << setprecision(4);
// reading order of the polynomials ...
cout << "Enter the number of terms in the polynomial A : ";
cin >> n1;
cout << "Enter the number of terms in the polynomial B : ";
cin >> n2;
// allocating space for the polynomials ...
if ( n1 > n2 ) n = n1;
else n = n2;
i = 1;
while ( i < n ) (i = i<<1);
n = i;
nn = n << 1;
A = new complex[nn];
B = new complex[nn];
C = new complex[nn];
fftA = new complex[nn];
fftB = new complex[nn];
fftC = new complex[nn];
// reading the polynomials from the user ...
cout << "Enter the polynomial A (highest degree coeff first) :\n";
for ( i = n1-1 ; i >= 0 ; i-- ) {
cin >> temp;
A[i] = temp;
}
for ( i = n1 ; i < nn ; i++ )
A[i] = 0;
cout << "Enter the polynomial B (highest degree coeff first) :\n";
for ( i = n2-1 ; i >= 0 ; i-- ) {
cin >> temp;
B[i] = temp;
}
for ( i = n2 ; i < nn ; i++ )
B[i] = 0;
// finding the nnth complex root of the unity ...
w = complex(cos(-2.0*PI/nn),sin(-2.0*PI/nn));
// finding the ffts of A and B polynomilas ...
fft(nn,1,A,w,fftA);
fft(nn,1,B,w,fftB);
// doing point by poing multiplication to find the fft of the product polynomial ...
multPolyFft(nn,fftA,fftB,fftC);
// doing inverse fft transform to find the product polynomial ...
w = 1.0/w;
fft(nn,1,fftC,w,C);
for ( i = 0 ; i < nn ; i++ )
C[i] = C[i]/double(nn);
// displaying the product polynomial ...
cout << "The product polynomial is (highest degree coeff first) :\n";
for ( i = nn-1; i >=0 ; i-- )
cout << setw(20) << C[i].real();
// now emperically testing the computed polynomial ...
cout << "Now testing the product polynomial emperically for random values...\n";
complex x;
complex y1, y2;
srand(0);
int errorcntR=0, errorcntI=0;
double errorR, errorI;
for ( i = 0 ; i < 100000 ; i++ )
{
x = (float(rand())/RAND_MAX)/5000; // generating random value
if ( i%2 )
x = -1.0*x;
y1 = evalPoly(n,A,x)*evalPoly(n,B,x); // actual value
y2 = evalPoly(nn,C,x); // emperical value
errorR = (y1.real()-y2.real())/y1.real(); // relative error in real part
errorI = (y1.imag()-y2.imag())/y2.imag(); // relative error in imag part
if ( errorR < 0 ) // finding the absolute
errorR = errorR * -1.0;
if ( errorI < 0 ) // finding the absolute
errorI = errorI * -1.0;
if ( errorR > 0.001 ) // count the error
errorcntR++;
if ( errorI > 0.00001 ) // count the error
errorcntI++;
}
// display the result of emperical testing ...
cout << "For 100000 random values, error counts are ...\n";
cout << "Error Count in real part : " << errorcntR << endl
<< "Error Count in imaginary part : " << errorcntI << endl;
}