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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : polynomial_simd.hpp
* DESCRIPTION: Simd support for polynomials
* 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_SIMD__HPP
#define __MMX__POLYNOMIAL_SIMD__HPP
#include <numerix/simd.hpp>
#ifdef NUMERIX_ENABLE_SIMD
#include <algebramix/polynomial_unrolled.hpp>

namespace mmx {

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

#ifdef __SSE2__

template<typename V, typename C, typename K>
struct polynomial_mul_simd_helper {
  typedef implementation<polynomial_linear,V> Pol;
  static const nat m = Simd_size (C);

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

#define POLYNOMIAL_MUL_IMPL(C)                                          \
  template<typename V, typename K>                                      \
  struct polynomial_mul_helper<V,C,K> {                                 \
    static inline void                                                  \
    op (C* dest, const C* s1, const C* s2, nat n1, nat n2) {            \
      polynomial_mul_simd_helper<V,C,K>::op (dest, s1, s2, n1, n2); }   \
  };

template<typename V, typename C>
struct polynomial_tmul_simd_helper {
  typedef implementation<polynomial_linear,V> Pol;
  static const nat m = Simd_size (C);
  static void
  op (C* dest, const C* s1, const C* s2, nat n1, nat n2) {
    nat l = aligned_size<C,V> (n1+m);
    C* aligned_s1 = mmx_new<C> (l);
    for (nat j = 0; j < m; j++) {
      C* d = dest + j;
      const C* t2 = s2;
      Pol::copy (aligned_s1 + j, s1, n1);
      if (j > 0) aligned_s1[j-1] = (C) 0;
      for (nat i = j; i < n2; i+=m, d+=m, t2+=m)
        *d = Pol::Vec::inn_prod (aligned_s1, t2, n1+j); 
    }
    mmx_delete<C> (aligned_s1, l);
  } 
};

#define POLYNOMIAL_TMUL_IMPL(C)                                         \
  template<typename V>                                                  \
  struct polynomial_tmul_helper<V,C> {                                  \
    static inline void                                                  \
    op (C* dest, const C* s1, const C* s2, nat n1, nat n2) {            \
      polynomial_tmul_simd_helper<V,C>::op (dest, s1, s2, n1, n2); }    \
  };

POLYNOMIAL_MUL_IMPL (double);
POLYNOMIAL_TMUL_IMPL (double);

#undef POLYNOMIAL_MUL_IMPL
#undef POLYNOMIAL_TMUL_IMPL

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

template<typename V, typename C>
struct polynomial_quo_rem_simd_helper {
  typedef implementation<polynomial_linear,V> Pol;
  static const nat m = Simd_size (C);
  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
    C* aligned_s2 [m];
    for (nat j = 1; j < m; j++) {
      aligned_s2[j] = mmx_new<C> (aligned_size<C,V> (n2+j));
      for (nat k = 0; k < j; k++)
        aligned_s2[j][k] = C (0);
      Pol::copy (aligned_s2[j] + j, s2, n2);
    }
    while (n1 >= n2 && n1 != 0) {
      int d= n1-n2;
      nat j= d % m; 
      C q= quo (s1[n1-1], s2[n2-1]);
      dest[d]= q;
      if (j == 0)
        Pol::Vec::template vec_binary_scalar<mul_add_op,C,C,C>
          (s1 + d, s2, -q, n2);
      else
        Pol::Vec::template vec_binary_scalar<mul_add_op,C,C,C>
          (s1 + d - j, aligned_s2 [j], -q, n2 + j);
      n1--;
    }
    for (nat j = 1; j < m; j++)
      mmx_delete<C> (aligned_s2[j], aligned_size<C,V> (n2+j));
  }
};

#define POLYNOMIAL_QUO_REM_IMPL(C)                                      \
  template<typename V>                                                  \
  struct polynomial_quo_rem_helper<V,C> {                               \
    static inline void                                                  \
    op (C* dest, C* s1, const C* s2, nat n1, nat n2) {                  \
      polynomial_quo_rem_simd_helper<V,C>::op (dest, s1, s2, n1, n2); } \
  };

POLYNOMIAL_QUO_REM_IMPL (double);

#undef POLYNOMIAL_QUO_REM_IMPL

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

template<typename V, typename C>
struct polynomial_evaluate_simd_helper {
  static const nat m = Simd_size (C);
  typedef typename Simd_type (C) sC;
  static inline C
  op (const C* p, const C& x, nat l) {
    if (l == 0) return 0;
    nat k = (l+m-1) / m;
    const C* q = p + (k-1) * m;
    sC sr = simd_load_aligned (q);
    nat i = 1; q -= m;
    sC sxm= simd_set_duplicate (binpow (x, m));
    for (; i + 4 < k; i += 4, q -= 4*m) {
      sr = simd_load_aligned (q)     + sxm * sr;
      sr = simd_load_aligned (q-m)   + sxm * sr;
      sr = simd_load_aligned (q-2*m) + sxm * sr;
      sr = simd_load_aligned (q-3*m) + sxm * sr;
    }
    for (; i < k; i++, q -= m)
      sr = simd_load_aligned (q)     + sxm * sr;
    C r = ((C*) &sr) [m-1];
    for (i = m-2; i != ((nat) -1); i--)
      r = r * x + ((C*) &sr) [i];
    return r;
  }
};
  
#define POLYNOMIAL_EVALUATE_IMPL(C)                                     \
  template<typename V>                                                  \
  struct polynomial_evaluate_helper<V,C> {                              \
    static inline C                                                     \
    op (const C* p, const C& x, nat l) {                                \
      return polynomial_evaluate_simd_helper<V,C>::op (p, x, l); }      \
  };

POLYNOMIAL_EVALUATE_IMPL (double);

#undef POLYNOMIAL_EVALUATE_IMPL

#endif // __SSE2__

} // namespace mmx

#endif // NUMERIX_ENABLE_SIMD
#endif //__MMX__POLYNOMIAL_SIMD__HPP