html/doxygen/fft__double_8cpp_source.html

00001
00002 /******************************************************************************
00003 * MODULE     : fft_double.hpp
00004 * DESCRIPTION: special initialization of FFT coefficients for doubles
00005 * COPYRIGHT  : (C) 2007  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012
00013 #include <numerix/floating.hpp>
00014 #include <numerix/complex.hpp>
00015 #include <algebramix/fft_naive.hpp>
00016 #include <algebramix/polynomial_tft.hpp>
00017
00018 namespace mmx {
00019
00020 #ifdef ACCURATE_ROOTS
00021 complex<double>
00022 primitive_root_helper<complex<double> >::op (nat n, nat i) {
00023   nat old= mmx_bit_precision;
00024   mmx_bit_precision= 64;
00025   complex<floating<> > zf= primitive_root<complex<floating<> > > (n, i);
00026   complex<double> zd (as_double (Re (zf)), as_double (Im (zf)));
00027   mmx_bit_precision= old;
00028   return zd;
00029 }
00030 #endif // ACCURATE_ROOTS
00031
00032 void
00033 fft_mul (double* dest, const double* s1, const double *s2, nat n1, nat n2) {
00034   typedef polynomial_dicho<
00035             polynomial_tft<
00036               polynomial_karatsuba<
00037                 polynomial_naive> > > PV;
00038   typedef implementation<vector_linear,PV> Vec;
00039   typedef implementation<polynomial_multiply,PV> Pol;
00040   //mmout << "s1= "; Vec::print (mmout, s1, n1); mmout << "\n";
00041   //mmout << "s2= "; Vec::print (mmout, s2, n2); mmout << "\n";
00042   nat nd = n1 + n2 - 1;
00043   nat spc= n1 + n2 + nd;
00044   complex<double>* m1= mmx_new<complex<double> > (spc);
00045   complex<double>* m2= m1 + n1;
00046   complex<double>* md= m2 + n2;
00047   Vec::cast (m1, s1, n1);
00048   Vec::cast (m2, s2, n2);
00049   Pol::mul (md, m1, m2, n1, n2);
00050   Vec::vec_unary<Re_op> (dest, md, nd);
00051   mmx_delete<complex<double> > (m1, spc);
00052   //mmout << "dd= "; Vec::print (mmout, dest, nd); mmout << "\n";
00053 }
00054
00055 } // namespace mmx