/Users/mourrain/Devel/mmx/algebramix/include/algebramix/polynomial_complex.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_complex.hpp
* DESCRIPTION: Complex multiplication of polynomials
* 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.
******************************************************************************/

#ifndef __MMX__POLYNOMIAL_COMPLEX__HPP
#define __MMX__POLYNOMIAL_COMPLEX__HPP
#include <algebramix/polynomial.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <numerix/complex.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Variant for complex polynomials
******************************************************************************/

template<typename V>
struct polynomial_complex {
  typedef typename V::Vec Vec;
  typedef typename V::Naive Naive;
  typedef typename V::Positive Positive;
  typedef polynomial_complex<typename V::No_simd> No_simd;
  typedef polynomial_complex<typename V::No_thread> No_thread;
  typedef polynomial_complex<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,polynomial_complex<W> >:
  public implementation<F,V,W> {};

template<typename C>
struct polynomial_variant_helper<complex<C> > {
  typedef typename Polynomial_variant(C) CPV;
  typedef polynomial_complex<CPV> PV;
};

/******************************************************************************
* Complex multiplication
******************************************************************************/

template<typename V, typename CV>
struct implementation<polynomial_multiply,V,polynomial_complex<CV> >:
  public implementation<polynomial_linear,V>
{
  typedef implementation<vector_linear,CV> Vec_C;
  typedef implementation<polynomial_multiply,CV> Pol_C;

TMPL static inline void
mul (complex<C>* dest, const complex<C>* s1, const complex<C>* s2,
     nat n1, nat n2)
{
  nat nd    = n1 + n2 - 1;
  nat spc1  = aligned_size<C,CV> (n1);
  nat spc2  = aligned_size<C,CV> (n2);
  nat spcd  = aligned_size<C,CV> (nd);
  nat spc   = 5 * (spc1 + spc2);
  C* real_x1= mmx_new<C> (spc);
  C* imag_x1= real_x1 + spc1;
  C* real_x2= imag_x1 + spc1;
  C* imag_x2= real_x2 + spc2;
  C* real_d = imag_x2 + spc2;
  C* imag_d = real_d  + spcd;
  C* temp_d = imag_d  + spcd;
  
  Vec_C::half_copy (real_x1, (C*) ((void*) s1), n1);
  Vec_C::half_copy (imag_x1, ((C*) ((void*) s1)) + 1, n1);
  Vec_C::half_copy (real_x2, (C*) ((void*) s2), n2);
  Vec_C::half_copy (imag_x2, ((C*) ((void*) s2)) + 1, n2);
  Pol_C::mul (real_d, real_x1, real_x2, n1, n2);
  Pol_C::mul (temp_d , imag_x1, imag_x2, n1, n2);
  Pol_C::add (real_x1, imag_x1, n1);
  Pol_C::add (real_x2, imag_x2, n2);
  Pol_C::mul (imag_d, real_x1, real_x2, n1, n2);
  Pol_C::sub (imag_d, real_d, nd);
  Pol_C::sub (imag_d, temp_d, nd);
  Pol_C::sub (real_d, temp_d, nd);
  Vec_C::double_copy ((C*) ((void*) dest), real_d, nd);
  Vec_C::double_copy (((C*) ((void*) dest)) + 1, imag_d, nd);

  mmx_delete<C> (real_x1, spc);
}

TMPL static inline void
square (complex<C>* dest, const complex<C>* s1, nat n1) {
  nat nd    = n1 + n1 - 1;
  nat spc1  = aligned_size<C,CV> (n1);
  nat spcd  = aligned_size<C,CV> (nd);
  nat spc   = 8 * spc1;
  C* real_x1= mmx_new<C> (spc);
  C* imag_x1= real_x1 + spc1;
  C* real_d = imag_x1 + spc1;
  C* imag_d = real_d  + spcd;
  C* temp_d = imag_d  + spcd;
  
  Vec_C::half_copy (real_x1, (C*) ((void*) s1), n1);
  Vec_C::half_copy (imag_x1, ((C*) ((void*) s1)) + 1, n1);
  Pol_C::square (real_d, real_x1, n1);
  Pol_C::square (temp_d, imag_x1, n1);
  Pol_C::add (real_x1, imag_x1, n1);
  Pol_C::square (imag_d, real_x1, n1);
  Pol_C::sub (imag_d, real_d, nd);
  Pol_C::sub (imag_d, temp_d, nd);
  Pol_C::sub (real_d, temp_d, nd);
  Vec_C::double_copy ((C*) ((void*) dest), real_d, nd);
  Vec_C::double_copy (((C*) ((void*) dest)) + 1, imag_d, nd);

  mmx_delete<C> (real_x1, spc);
}

TMPL static inline void
tmul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  (void) dest; (void) s1; (void) s2; (void) n1; (void) n2;
  ERROR ("transposed complexified multiplication not implemented");
}

}; // implementation<polynomial_multiply,V,polynomial_complex<CV> >

#undef TMPL
} // namespace mmx
#endif //__MMX__POLYNOMIAL_COMPLEX__HPP