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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_unrolled.hpp
* DESCRIPTION: Loop unrolling for naive polynomial arithmetic
* COPYRIGHT  : (C) 2008  Gregoire Lecerf
*******************************************************************************
* 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_UNROLLED__HPP
#define __MMX__POLYNOMIAL_UNROLLED__HPP
#include <algebramix/polynomial_naive.hpp>
#include <algebramix/polynomial.hpp>
#include <algebramix/vector_unrolled.hpp>
#include <numerix/simd.hpp>

#define TMPL template<typename C>
#define TMPLK template<typename C, typename K>

namespace mmx {

/******************************************************************************
* Unrolled algorithms for polynomials
******************************************************************************/

template<typename V, nat sz>
struct polynomial_unrolled: public V {
  typedef vector_unrolled<sz> Vec;
  typedef typename V::Naive Naive;
  typedef polynomial_unrolled<typename V::Positive,sz> Positive;
  typedef polynomial_unrolled<typename V::No_simd,sz> No_simd;
  typedef polynomial_unrolled<typename V::No_thread,sz> No_thread;
  typedef polynomial_unrolled<typename V::No_scaled,sz> No_scaled;
};

template<typename F, typename V, typename W, nat sz>
struct implementation<F,V,polynomial_unrolled<W,sz> >:
  public implementation<F,V,W> {};

/******************************************************************************
* Multiplication
******************************************************************************/

template<typename V, typename C, typename K>
struct polynomial_mul_helper {
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_linear,V> Pol;

  static void
  op (C* dest, const C* s1, const K* s2, nat n1, nat n2) {
    nat l = aligned_size<K,V> (n2);
    K* rev_s2 = mmx_new<K> (l);
    Pol::reverse (rev_s2, s2, n2);
    for (nat i = 1; i <= n2; i++, dest++)
      *dest = Vec::inn_prod (s1    , rev_s2 + n2 - i, min (i, n1));
    for (nat i = 1; i <  n1; i++, dest++)
      *dest = Vec::inn_prod (s1 + i, rev_s2         , min (n2, n1 - i));
    mmx_delete (rev_s2, l);
  } 
};

template<typename V, typename C>
struct polynomial_tmul_helper {
  typedef implementation<vector_linear,V> Vec;

  static void
  op (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
    for (nat i = 0; i < n2; i++, dest++, s2++)
      *dest = Vec::inn_prod (s1, s2, n1); 
  }
};

template<typename V, typename W, nat m>
struct implementation<polynomial_multiply,V,polynomial_unrolled<W,m> >:
  public implementation<polynomial_linear,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_multiply,W> Fallback;

TMPLK static void
mul (C* dest, const C* s1, const K* s2, nat n1, nat n2) {
  Vec::clear (dest, n1+n2-1);
  polynomial_mul_helper<W,C,K>::op (dest, s1, s2, n1, n2);
}

TMPL static inline void
square (C* dest, const C* s, nat n) {
  mul (dest, s, s, n, n);
}

TMPL static void
tmul (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
  Vec::clear (dest, n2);
  polynomial_tmul_helper<W,C>::op (dest, s1, s2, n1, n2);
}

}; // implementation<polynomial_multiply,V,polynomial_unrolled<W,m> >

/******************************************************************************
* Graeffe transformations
******************************************************************************/

template<typename V, typename W, nat m>
struct implementation<polynomial_graeffe,V,polynomial_unrolled<W,m> >:
  public implementation<polynomial_multiply,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<polynomial_multiply,V> Pol;

TMPL static void
graeffe (C* dest, const C* s, nat n) {
  nat l1= (n+1) >> 1, l2= n >> 1;
  // More allocation is done to ensure memory alignment
  nat l= aligned_size<C,V> (n + 1);
  C* c1= mmx_new<C> (l);
  C* c2= mmx_new<C> (l);
  if (n>0) dest[n-1]= C(0);
  Vec::half_copy (c1, s  , l1);
  Pol::square (dest, c1, l1);
  Vec::half_copy (c1, s+1, l2);
  Pol::square (c2 , c1, l2);
  Pol::copy (c1 + 1, c2, l2 << 1);
  c1[0] = C(0);
  Pol::sub (dest, c1, l2 << 1);
  mmx_delete<C> (c1, l);
  mmx_delete<C> (c2, l);
}

}; // implementation<polynomial_graeffe,V,polynomial_unrolled<W,m> >

/******************************************************************************
* Division
******************************************************************************/

template<typename V, typename C>
struct polynomial_quo_rem_helper {
  typedef implementation<vector_linear,V> Vec;

  static void
  op (C* dest, C* s1, const C* s2, nat n1, nat n2) {
    // (s1, n1) contains the numerator on input and the remainder on output
    // (s2, n2) contains the denominator. We assume n2>0 and s2[n2-1] != 0
    // (dest, n1-n2+1) contains the quotient
    while (n1 >= n2 && n1 != 0) {
      int d= n1-n2;
      C q= quo (s1[n1-1], s2[n2-1]);
      dest[d]= q;
      Vec::template vec_binary_scalar<mul_add_op,C,C,C> (s1 + d, s2, -q, n2);
      n1--;
    }
  }
};

template<typename V, typename W, nat m>
struct implementation<polynomial_divide,V,polynomial_unrolled<W,m> >:
  public implementation<polynomial_multiply,V>
{

TMPL static void
quo_rem (C* dest, C* s1, const C* s2, nat n1, nat n2) {
  polynomial_quo_rem_helper<W,C>::op (dest, s1, s2, n1, n2);
}

}; // implementation<polynomial_divide,V,polynomial_unrolled<W,m> >

/******************************************************************************
* Efficient evaluation of polynomials
******************************************************************************/

template<typename V, typename C>
struct polynomial_evaluate_helper {
  typedef implementation<polynomial_evaluate,polynomial_naive> Pol;

  static inline C
  op (const C* p, const C& x, nat l) {
    return Pol::evaluate (p, x, l);
  }
};
#ifdef ALGEBRAMIX_ENABLE_UNSTABLE
template<typename V, typename W, nat m>
struct implementation<polynomial_evaluate,V,polynomial_unrolled<W,m> >:
  public implementation<polynomial_divide,V>
{

TMPL static void
annulator (C* d, const C* s, nat n) {
  (void) d; (void) s; (void) n;
  ERROR ("expand not yet implemented");
}

TMPL static C
evaluate (const C* p, const C& x, nat l) {
  return polynomial_evaluate_helper<W,C>::op (p, x, l);
}

TMPL static void
evaluate (C* v, const C* p, const C* x, nat l, nat n) {
  // Evaluate the polynomial (p, l) at points x[0], ..., x[n-1]
  for (nat j=0; j<n; j++)
    v[j]= evaluate (p, x[j], l);
}

TMPL static void
expand (C** v, const C* p, const C* x, const nat* nu, nat n, nat k) {
  (void) v; (void) p; (void) x; (void) nu; (void) n; (void) k;
  ERROR ("expand not yet implemented");
}

}; // implementation<polynomial_evaluate,V,polynomial_unrolled<W,m> >
#endif
} // namespace mmx

#undef TMPL
#undef TMPLK
#endif //__MMX__POLYNOMIAL_UNROLLED__HPP