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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_modular_int.hpp
* DESCRIPTION: modular matrices over hardware integers
* COPYRIGHT  : (C) 2009  Joris van der Hoeven and 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__MATRIX_MODULAR_INT__HPP
#define __MMX__MATRIX_MODULAR_INT__HPP
#include <numerix/modular_int.hpp>
#include <algebramix/matrix_int.hpp>
#include <algebramix/matrix_double.hpp>
#include <algebramix/matrix_modular.hpp>

namespace mmx {
#define TMPL template<typename C, typename V1, typename V2>
#define Modulus modulus<C,V1>
#define Modular modular<Modulus,V2>

/******************************************************************************
* Variant
******************************************************************************/

typedef matrix_unrolled<4,matrix_unrolled<4,matrix_naive> >
        matrix_unrolled_4_4;

DEFINE_VARIANT (matrix_modular_int,
                matrix_strassen<
                  matrix_threads<
                    matrix_unrolled_4_4> >)

#define DECLARE_HELPER(I,V1)                                            \
  template<typename V2, nat m>                                          \
  struct matrix_variant_helper<modular<modulus<I,V1>,V2> > {            \
    typedef matrix_naive Naive;                                         \
    typedef matrix_unrolled<4,matrix_unrolled<2,Naive> > Unrolled;      \
    typedef typename Matrix_simd_variant(I) Simd;                       \
    typedef matrix_modular_int MV;                                      \
  };

DECLARE_HELPER(signed char, modulus_int_naive<m>)
DECLARE_HELPER(signed char, modulus_int_preinverse<m>)
DECLARE_HELPER(unsigned char, modulus_int_naive<m>)
DECLARE_HELPER(unsigned char, modulus_int_preinverse<m>)
DECLARE_HELPER(short int, modulus_int_naive<m>)
DECLARE_HELPER(short int, modulus_int_preinverse<m>)
DECLARE_HELPER(unsigned short int, modulus_int_naive<m>)
DECLARE_HELPER(unsigned short int, modulus_int_preinverse<m>)
DECLARE_HELPER(int, modulus_int_naive<m>)
DECLARE_HELPER(int, modulus_int_preinverse<m>)
DECLARE_HELPER(unsigned int, modulus_int_naive<m>)
DECLARE_HELPER(unsigned int, modulus_int_preinverse<m>)
DECLARE_HELPER(long int, modulus_int_naive<m>)
DECLARE_HELPER(long int, modulus_int_preinverse<m>)
DECLARE_HELPER(unsigned long int, modulus_int_naive<m>)
DECLARE_HELPER(unsigned long int, modulus_int_preinverse<m>)
DECLARE_HELPER(long long int, modulus_int_naive<m>)
DECLARE_HELPER(long long int, modulus_int_preinverse<m>)
DECLARE_HELPER(unsigned long long int, modulus_int_naive<m>)
DECLARE_HELPER(unsigned long long int, modulus_int_preinverse<m>)
#undef DECLARE_HELPER

/******************************************************************************
* Encoding and decoding from/to lifting type
******************************************************************************/

template<typename D>
struct free_bits_helper {
  static const nat value= 8 * sizeof (D);
};

template<typename V, typename D, typename C, typename V1, typename V2> void
encode_modular_int (D* dest, const Modular* src,
                    nat r, nat rr, nat c, nat cc)
{
  typedef implementation<matrix_linear,V> Mat;
  nat dest_rs= Mat::index (1, 0, r, c);
  nat dest_cs= Mat::index (0, 1, r, c);
  nat src_rs = Mat::index (1, 0, rr, cc);
  nat src_cs = Mat::index (0, 1, rr, cc);
  D* dest_row= dest;
  const Modular* src_row= src;
  for (nat i=0; i<r; i++, dest_row += dest_rs, src_row += src_rs) {
    D* dest_col= dest_row;
    const Modular* src_col= src_row;
    for (nat j=0; j<c; j++, dest_col += dest_cs, src_col += src_cs)
      *dest_col= (D) (*(*src_col));
  }
}

template<typename Op, typename V, typename C, typename D,
          typename V1, typename V2> void
decode_modular_int (Modular* dest, const D* src,
                    nat r, nat rr, nat c, nat cc)
{
  typedef implementation<matrix_linear,V> Mat;
  typedef typename Op::nomul_op Set;
  C p = *Modular::get_modulus ();
  nat dest_rs= Mat::index (1, 0, rr, cc);
  nat dest_cs= Mat::index (0, 1, rr, cc);
  nat src_rs = Mat::index (1, 0, r, c);
  nat src_cs = Mat::index (0, 1, r, c);
  Modular* dest_row= dest;
  const D* src_row= src;
  for (nat i=0; i<r; i++, dest_row += dest_rs, src_row += src_rs) {
    Modular* dest_col= dest_row;
    const D* src_col= src_row;
    for (nat j=0; j<c; j++, dest_col += dest_cs, src_col += src_cs)
      Set::set_op (*dest_col, Modular ((*src_col) % p, true));
  }
}

#if defined(NUMERIX_ENABLE_SIMD)\
 && defined(ALGEBRAMIX_ENABLE_SIMD)\
 && defined(__SSE2__)
// double are faster than int's on usual platforms supporting SSE2
template<typename Op, typename V,
         typename C, typename V1, typename V2> void
decode_modular_int (Modular* dest, const double* src,
                    nat r, nat rr, nat c, nat cc)
{
  typedef implementation<matrix_linear,V> Mat;
  typedef typename Op::nomul_op Set;
  uint32_t p = (uint32_t) *Modular::get_modulus ();
  nat dest_rs= Mat::index (1, 0, rr, cc);
  nat dest_cs= Mat::index (0, 1, rr, cc);
  nat src_rs = Mat::index (1, 0, r, c);
  nat src_cs = Mat::index (0, 1, r, c);
  Modular* dest_row= dest;
  const double* src_row= src;
  for (nat i=0; i<r; i++, dest_row += dest_rs, src_row += src_rs) {
    Modular* dest_col= dest_row;
    const double* src_col= src_row;
    for (nat j=0; j<c; j++, dest_col += dest_cs, src_col += src_cs)
      Set::set_op (*dest_col, Modular (((uint64_t) (*src_col)) % p, true));
  }
}
STMPL
struct free_bits_helper<double> {
  static const nat value= 53;
};
#endif // __SSE2__


template<typename V>
struct implementation<matrix_multiply_base,V,matrix_modular_int>:
  public implementation<matrix_linear,V>
{
  typedef matrix_strassen<matrix_threads<matrix_unrolled_4_4> > W;
  typedef implementation<matrix_multiply,W> Fallback;
  typedef implementation<matrix_linear,W> Mat;
private:

#if defined(NUMERIX_ENABLE_SIMD)\
 && defined(ALGEBRAMIX_ENABLE_SIMD)\
 && defined(__SSE2__)
template<nat s>
struct lifting_int_with_size_at_least_helper {
  template<bool b, typename T, typename E>
  struct if_helper { typedef T ret; };
  template<typename T, typename E>
  struct if_helper<false,T,E> { typedef E ret; };
  typedef
    typename if_helper<8*sizeof(short unsigned int) >= s, short unsigned int,
    typename if_helper<53 >= s, double,
      long long unsigned int>::ret >::ret type;
};

#else
template<nat s>
struct lifting_int_with_size_at_least_helper {
  template<bool b, typename T, typename E>
  struct if_helper { typedef T ret; };
  template<typename T, typename E>
  struct if_helper<false,T,E> { typedef E ret; };
  typedef
    typename if_helper<8*sizeof(unsigned char)>=s, unsigned char,
    typename if_helper<8*sizeof(short unsigned int)>=s, short unsigned int,
    typename if_helper<8*sizeof(unsigned int)>=s, unsigned int,
    typename if_helper<8*sizeof(long unsigned int)>=s, long unsigned int,
    long long unsigned int>::ret >::ret >::ret >::ret type ;
};
#endif
  
template<typename Op, typename C, typename V1, typename V2>
static inline void
_mul_int (Modular* d, const Modular* s1, const Modular* s2,
         nat r, nat rr, nat l, nat ll, nat c, nat cc)
  {
    static const nat modulus_size= V1::template maximum_size_helper<C>::value;
    // 2 free bits are required
    static const nat required_size= (modulus_size << 1) + 2;
    typedef typename
      lifting_int_with_size_at_least_helper<required_size>::type D;
    // free bits in type D
    static const nat possible_free_bits=
      free_bits_helper<D>::value - (modulus_size << 1);
    static const nat free_bits= possible_free_bits >= 16
                                  ? 16 : possible_free_bits;
    const nat ul= MMX_SAFE_LEFT_SHIFT_INT(nat,1,free_bits);
    typedef typename Matrix_variant(D) W;
    typedef implementation<matrix_multiply,W> Mat_D;

    if (free_bits == 0) {
      Fallback::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc);
    }
    else if (l <= ul) {
      nat len_1= aligned_size<D,W> (r * l);
      nat len_2= aligned_size<D,W> (l * c);
      nat len_d= aligned_size<D,W> (r * c);
      nat spc = len_1 + len_2 + len_d;
      D* x1= mmx_new<D> (spc);
      D* x2= x1 + len_1;
      D* xd= x2 + len_2;
      encode_modular_int<W> (x1, s1, r, rr, l, ll);
      encode_modular_int<W> (x2, s2, l, ll, c, cc);
      Mat_D::template mul<Op> (xd, x1, x2, r, r, l, l, c, c);
      decode_modular_int<Op,W> (d, xd, r, rr, c, cc);
      mmx_delete<D> (x1, spc);
    }
    else {
      typedef typename Op::acc_op Acc;
      const nat ur= 32;
      const nat uc= 32;

      nat nr= (r + ur - 1) / ur;
      nat nl= (l + ul - 1) / ul;
      nat nc= (c + uc - 1) / uc;
      for (nat ir=0; ir<nr; ir++)
        for (nat ic=0; ic<nc; ic++) {
          nat il=0;
          for (; il<1; il++) {
            mul<Op > (d  + Mat::index (ir*ur, ic*uc, rr, cc),
                      s1 + Mat::index (ir*ur, il*ul, rr, ll),
                      s2 + Mat::index (il*ul, ic*uc, ll, cc),
                      min (ur, r - ir*ur), rr,
                      min (ul, l - il*ul), ll,
                      min (uc, c - ic*uc), cc);
          }
          for (; il<nl; il++) {
            mul<Acc> (d  + Mat::index (ir*ur, ic*uc, rr, cc),
                      s1 + Mat::index (ir*ur, il*ul, rr, ll),
                      s2 + Mat::index (il*ul, ic*uc, ll, cc),
                      min (ur, r - ir*ur), rr,
                      min (ul, l - il*ul), ll,
                      min (uc, c - ic*uc), cc);
          }
        }
    }
  }

public:

template<typename Op, typename C, typename V1, typename V2>
static inline void
mul (Modular* d, const Modular* s1, const Modular* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc)
  {
    _mul_int<Op> (d, s1, s2, r, rr, l, ll, c, cc);
  }

}; // implementation<matrix_multiply_base,V,matrix_modular_int<MV> >

#undef TMPL
#undef Modulus
#undef Modular
} // namespace mmx

#endif //__MMX__MATRIX_MODULAR_INT__HPP