synaps/upol/mult_fft.h

00001 #ifndef SYNAPS_UPOL_MULT_FFT_H
00002 #define SYNAPS_UPOL_MULT_FFT_H
00003 //======================================================================
00004 #include <synaps/init.h>
00005 #include <synaps/arithm/RealOf.h>
00006 #include <synaps/linalg/FFT.h>
00007 //======================================================================
00008 
00009 __BEGIN_NAMESPACE_SYNAPS
00010 
00011 
00012 namespace UPOL
00013 {
00021 template<class R, class TA, class TB>
00022 void mul_fft(R & r, const TA & a, const TB & b)
00023 {
00024   #ifdef VERBATIM
00025         std::cout <<"Using fft for the multiplication of polynomials"<<std::endl;
00026   #endif
00027 
00028   typedef typename TA::value_type C;
00029   typedef std::complex<typename RealOf<C>::TYPE> Cpl;
00030   unsigned long degree = degree(a) + degree(b);
00031   unsigned long N = (int) pow(2, 1 + (int) (log(degree)/log(2)));
00032   Cpl  w ;
00033 
00034   Root_Unity(w, N);
00035 
00036   Cpl * XA = new Cpl[N];
00037   Cpl * XB = new Cpl[N];
00038 
00039   for (long int i = 0; i <= degree(a) ; i++)
00040     XA[i] = a[i];
00041 
00042   // Bit reversal permutation, to allow Butterfly algorithm to calculate FFT
00043   // on the two polynomials coefficients.
00044 
00045   Permut(XA, N);
00046   FFT(N, w, XA,0);
00047 
00048   for (long int i = 0; i <= degree(b) ; i++)
00049     XB[i] = b[i];
00050 
00051   Permut(XB, N);
00052   FFT(N, w, XB, 0);
00053 
00054   // multiplication of FFT(polynomial a) by FFT(polynomial b)
00055   for (unsigned long k = 0; k < N; k++)
00056     XA[k] *= XB[k] ;
00057 
00058   // result of this multiplication is also permuted for the same reason
00059   Permut(XA, N);
00060 
00061   // Backward FFT , to retrieve the result of multiplication
00062   // of the two polynomials
00063   FFT(N, pow(w,-1), XA, 1);
00064 
00065   for (unsigned long k = 0; k <= degree; k++)
00066     RealOf<C>::assign(r[k], (XA[k]/N));
00067 
00068   delete [] XA;
00069   delete [] XB;
00070 }
00071 //----------------------------------------------------------------------
00072 }//UPOL
00073 
00074 __END_NAMESPACE_SYNAPS
00075 
00076 #endif // SYNAPS_UPOL_MULT_FTT_H
00077 

SYNAPS DOCUMENTATION
logo