/*


  POLYNOMIAL ARITHMETIC USING FAST FOURIER TRANSFORM

  Submitted by
  Ashish Gupta
  98131
  Group 3
  
  This implements five functions :
  Add
  Subtract
  Multiply
  Divide
  Reciprocal


*/

#include <iostream.h>
#include <complex>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>



const double pi=3.141592654;
const int MAXSIZE=10000;   // maximum degree of resultant polynomial allowed

complex <double>  A[MAXSIZE],B[MAXSIZE],C[MAXSIZE];
int s1,s2,bound;
int lnbound;


// rounds a complex number to 6th decimal place
complex<double> round(complex<double> a)
{
  double i,j;
  i=a.real();
  j=a.imag();
  i= floor(i*100000+0.5)/100000;
  j= floor(j*100000+0.5)/100000;
  return complex<double>(i,j);
}


// rounds a double to 6th decimal place
double round(double a)
{
  a= floor(a*100000+0.5)/100000;
  return a;
}


void display(complex<double> a[],int n)
{
  cout << endl;
  for(int i=0;i<n;i++)
    {
      cout << round(a[i]) << " , " ;
    }

}

// add two complex numbers
void add(complex<double> a[],complex <double> b[],complex <double> dest[],int n)
{
  for(int i=0;i<n;i++)
    dest[i]=a[i]+b[i];  
}


// subtract two complex numbers
void sub(complex<double> a[],complex <double> b[],complex <double> dest[],int n)
{
  for(int i=0;i<n;i++)
    dest[i]=a[i]-b[i];
}


// Bit reverse an integer
int rev(int k,int b)
{
  char a[b],i;
  for(i=0;i<b;i++)
    {
      a[i]=k%2;
      k/=2;
    }
  k=0;
  int s=1;
  for(i=b-1;i>=0;i--)
    {
      k+=a[i]*s;
      s*=2;
      
    }
  return k;
}

//____________Bit reverse an array___
void bitrev(complex <double> a[],complex <double> A[],int n)
{
  int ln=(int)ceil(log(n)/log(2));
  for(int i=0;i<n;i++)
      A[rev(i,ln)]=a[i];
}


//_______Fast Fourier Transform ( also includes inverse )
// mode = -1  ==> Inverse FFT
void fft(complex <double> a[],complex <double> A[],int n,int mode=1)
{
  bitrev(a,A,n);
  int ln=(int)ceil(log(n)/log(2));
  complex<double> t,u;
  int s,j,k;

  for(s=1;s<=ln;s++)
    {
      int m=(int)pow(2,s);
      complex<double> wm(cos(2*pi/m),sin(mode*2*pi/m));
      complex<double> w(1,0);
      for(j=0;j<=m/2-1;j++)
	{

	  for(k=j;k<=n-1;k+=m)
	    {

	      t=w*A[k+m/2];
	      u=A[k];
	      A[k]=u+t;
	      A[k+m/2]=u-t;
	    }

	  w=w*wm;
	}
    }
  if(mode==-1)
    {
      for(int i=0;i<n;i++)
	A[i]/=n;
    }  
 
}
  

//______Multiply two polynomials , where n is the resultant bound___________
// Takes time : O(n*log(n))
void mult(complex<double> a[],complex<double> b[],complex<double> c[],int n)
{
  complex<double> t1[n],t2[n];

  // Take FFT of both polynomials
  fft(a,t1,n);
  fft(b,t2,n);

  // Pointwise multiplication 
  for(int i=0;i<n;i++)
      t1[i]=t1[i]*t2[i];

  // Take Inverse FFT 
  fft(t1,c,n,-1);
}

// Left Shift a polynomial by s places
void lshift(complex<double> a[],int n,int s)
{
  int i=0;
  for(i=0;i<n-s;i++)
    a[i]=a[i+s];
  for(i=n-s;i<n;i++)
    a[i]=complex<double>(0,0);
}

// Reciprocal of a polynomial
// Takes time : O(n*log(n)) using very efficient implementation involving shifting of bits
// n is not the bound but 1+degree of polynomial 


//complex<double> t1[MAXSIZE*2],t2[MAXSIZE*4];


void rec(complex<double> a[],complex<double> d[],int n)
{
  complex<double> t1[MAXSIZE*2],t2[MAXSIZE*4];
  int i,j; 
  d[n-1]=1.0/a[n-1];

  // Iterate log(n)+1 times to get the reciprocal
  for(i=0;i<=(int)ceil(log(n)/log(2));i++)
    {
      int m=(int)pow(2,i);
      mult(&d[n-m],&d[n-m],t1,m*2);
      mult(&t1[0],&a[n-m*2],&t2[n],4*m);
      for(j=0;j<n;j++)
	d[j]=complex<double>(2,0)*d[j]-t2[j+2*(n-1)-3*n+4*m+n];
    }
}


// Divides polynomials
// Time taken  : O(n*log(n))
void div(complex<double> a[],complex<double> b[],complex<double> c[],int bound,int s)
{
  complex<double> t1[MAXSIZE];
  rec(b,t1,s);
  mult(t1,a,c,bound);
  lshift(c,bound,2*(s-1));
}

// For generating random polynomials of size n
void rpoly(complex<double> a[],int n)
{
  int i;
  for(i=0;i<n;i++)
    a[i]=complex<double>(random()%1000,0);
}

//________ Main functions
void main()
{
  int i,j,k;
  double c;
  char op[15];
  
  cout << "\n\n#\tPOLYNOMIAL ARITHMETIC USING FAST FOURIER TRANSFORM\n\n";

  /*
    
    // This commented out part can be used to run the operations on random polynomials
    // of increasing size for generating graphs etc.
 
  for(j=2;j<2000;j+=32)
    {
      s1=s2=j;
      s2/=2;
      rpoly(A,s1);
      rpoly(B,s2);

  
      int maxs=s1<s2 ? s2 : s1; 
      bound=2*(int)pow(2,ceil(log(maxs)/log(2)));
      
      //      cout << "Bound : " << bound;
      for(i=s1;i<bound;i++)
	A[i]=complex<double>(0,0);
      for(i=s2;i<bound;i++)
	B[i]=complex<double>(0,0);
      
      for(i=0;i<bound;i++)
	C[i]=complex<double>(0,0);
      
      timeval t1,t2;
      gettimeofday(&t1,NULL);
  

      //cout << "mult " << bound << flush;
      mult(A,B,C,bound);
      gettimeofday(&t2,NULL);
  
      long timediff=(t2.tv_sec-t1.tv_sec)*1000000 + (t2.tv_usec-t1.tv_usec);

      cout << "\n" <<  j << "\t" << timediff;
  
    }

  exit(0);
  
  */  
  
  
  while(1)
    {     
      cout << "\n------------------------------------------------------------";      
      cout << "\n\nEnter Operation ( add , sub , mult , div , rec  or end ): ";
      cin >> op;
      cout << op;
      if(strcmp(op,"end")==0)
	  break;
      
      cout << "\nPolynomial 1 : \n";
      cout << "Enter Size : ";
      cin >> s1;
      cout << "Enter coefficients seperated by space in decreasing order of degree :\n";
      for(i=s1-1;i>=0;i--)
	{
	  cin >> c;
	  A[i]=complex <double> (c,0);
	  cout << c << "  ";
	}       
 

      if(strcmp(op,"rec"))
	{
	  cout << "\n\nPolynomial 2 : \n";
	  cout << "Enter Size : ";
	  cin >> s2;
	  cout << "Enter coefficients seperated by space in decreasing order of degree :\n";
	  
	  for(i=s2-1;i>=0;i--)
	    {
	      cin >> c;
	      B[i]=complex <double> (c,0);   
	      cout << c << "  ";
	    }       
	}
      else
	s2=0;
 
      int maxs=s1<s2 ? s2 : s1; 
      bound=2*(int)pow(2,ceil(log(maxs)/log(2)));
 
      //      cout << "Bound : " << bound;
      for(i=s1;i<bound;i++)
	A[i]=complex<double>(0,0);
      for(i=s2;i<bound;i++)
	B[i]=complex<double>(0,0);

     for(i=0;i<bound;i++)
	C[i]=complex<double>(0,0);

      lnbound=1+(int)ceil(log(maxs)/log(2));
      if(strcmp(op,"add")==0)
	add(A,B,C,bound);
      if(strcmp(op,"sub")==0)
	sub(A,B,C,bound);
      if(strcmp(op,"mult")==0)
	mult(A,B,C,bound);
      if(strcmp(op,"div")==0)
	div(A,B,C,bound,s2);
      if(strcmp(op,"rec")==0)
	rec(A,C,s1);

      cout << endl;
      cout << "\nAnswer is : \n";
      for(i=0;i<bound;i++)
	  cout << round(C[bound-i-1].real()) << "  ";
    }
  
  cout << endl  << "\n---------- Finished --------------\n";  
  
}