src/fft_double.cpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_double.hpp
* DESCRIPTION: special initialization of FFT coefficients for doubles
* COPYRIGHT  : (C) 2007  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <numerix/floating.hpp>
#include <numerix/complex.hpp>
#include <algebramix/fft_naive.hpp>
#include <algebramix/polynomial_tft.hpp>

namespace mmx {

#ifdef ACCURATE_ROOTS
complex<double>
primitive_root_helper<complex<double> >::op (nat n, nat i, const format<complex<double> >&) {
  nat old= mmx_bit_precision;
  mmx_bit_precision= 64;
  format<complex<floating<> > > ffm;
  complex<floating<> > zf= primitive_root<complex<floating<> > > (n, i, ffm);
  complex<double> zd (as_double (Re (zf)), as_double (Im (zf)));
  mmx_bit_precision= old;
  return zd;
}
#endif // ACCURATE_ROOTS

void
fft_mul (double* dest, const double* s1, const double *s2, nat n1, nat n2) {
  typedef polynomial_dicho<
            polynomial_tft<
              polynomial_karatsuba<
                polynomial_naive> > > PV;
  typedef implementation<vector_linear,PV> Vec;
  typedef implementation<polynomial_multiply,PV> Pol;
  //mmout << "s1= "; Vec::print (mmout, s1, n1); mmout << "\n";
  //mmout << "s2= "; Vec::print (mmout, s2, n2); mmout << "\n";
  nat nd = n1 + n2 - 1;
  nat spc= n1 + n2 + nd;
  complex<double>* m1= mmx_new<complex<double> > (spc);
  complex<double>* m2= m1 + n1;
  complex<double>* md= m2 + n2;
  Vec::cast (m1, s1, n1);
  Vec::cast (m2, s2, n2);
  Pol::mul (md, m1, m2, n1, n2);
  Vec::vec_unary<Re_op> (dest, md, nd);
  mmx_delete<complex<double> > (m1, spc);
  //mmout << "dd= "; Vec::print (mmout, dest, nd); mmout << "\n";
}

} // namespace mmx

---