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
00043
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
00055 for (unsigned long k = 0; k < N; k++)
00056 XA[k] *= XB[k] ;
00057
00058
00059 Permut(XA, N);
00060
00061
00062
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 }
00073
00074 __END_NAMESPACE_SYNAPS
00075
00076 #endif // SYNAPS_UPOL_MULT_FTT_H
00077