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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_crt.hpp
* DESCRIPTION: Multi-modular multiplication of matrices
* 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_CRT_HPP
#define __MMX_MATRIX_CRT_HPP
#include <numerix/modular.hpp>
#include <algebramix/matrix.hpp>
#include <algebramix/crt_blocks.hpp>
namespace mmx {

/******************************************************************************
* Special variant to avoid colisions with other global moduli
******************************************************************************/

template<typename C>
struct modular_matrix_crt {
  template<typename M>
  class modulus_storage {
    static inline M& dyn_modulus () {
      static M modulus = M ();
      return modulus; }
  public:
    static inline void set_modulus (const M& p) { dyn_modulus () = p; }
    static inline M get_modulus () { return dyn_modulus (); }
  };
};

/******************************************************************************
* Variant for matrices with CRT based operations
******************************************************************************/

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

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

/******************************************************************************
* Default non-archimedian size helper
******************************************************************************/

template<typename C>
struct matrix_crt_multiply_helper {
  static const nat dimension_threshold= 7;
  static const nat ratio_threshold= 100; // this is a percentage

  typedef crt_naive_transformer<C> crt_transformer;
  typedef moduli_helper<C,
        modulus<typename crt_transformer::modulus_base,
          typename crt_transformer::modulus_base_variant> > moduli_sequence;

  static nat size (const C* s1, nat s1_rs, nat s1_cs,
            const C* s2, nat s2_rs, nat s2_cs,
            nat r, nat l, nat c) {
    nat sz= 0;
    for (nat k= 0; k < l; k++) {
      nat sz1= 0, sz2= 0;
      const C* ss1= s1 + k * s1_cs;
      const C* ss2= s2 + k * s2_rs;
      for (nat i= 0; i < r; i++, ss1 += s1_rs) sz1= max (sz1, N (*ss1));
      for (nat j= 0; j < c; j++, ss2 += s2_cs) sz2= max (sz2, N (*ss2));
      sz= max (sz, sz1 + sz2);
    }
    return sz; }
};

/******************************************************************************
* Multi-modular multiplication of matrices
******************************************************************************/

template<typename V, typename W>
struct implementation<matrix_multiply,V,matrix_crt<W> >:
  public implementation<matrix_multiply_base,V>
{
  typedef implementation<matrix_multiply,W> Mat;

template<typename Op, typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc); }

template<typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat l, nat c) {
  Mat::template mul<mul_op> (d, s1, s2, r, r, l, l, c, c); }

template<typename C, typename I, typename MV, typename Crter>
static void
mat_direct_crt (matrix<I,MV>* dest, const C* s,
                nat s_rs, nat s_cs, nat r, nat c, Crter& crter) {
  nat n= N(crter);
  for (nat k= 0; k < n; k++)
    dest[k]= matrix<I,MV> (I (), r, c);
  I* aux= mmx_new<I> (n);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      direct_crt (aux, s[i * s_rs + j * s_cs], crter);
      for (nat k= 0; k < n ; k++) dest[k](i,j)= aux[k];
    }
  mmx_delete<I> (aux, n); }

template<typename C, typename Modulus, typename MW, typename MV, typename Crter>
static void
mat_direct_crt (matrix<modular<Modulus,MW>,MV>* dest, const C* s,
                nat s_rs, nat s_cs, nat r, nat c, Crter& crter) {
  typedef modular<Modulus,MW> Modular;
  typedef typename Modular::modulus::base I;
  nat n= N(crter);
  for (nat k= 0; k < n; k++)
    dest[k]= matrix<Modular,MV> (Modular (), r, c);
  I* aux= mmx_new<I> (n);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      direct_crt (aux, s[i * s_rs + j * s_cs], crter);
      for (nat k= 0; k < n ; k++) dest[k](i,j)= Modular (aux[k], true);
    }
  mmx_delete<I> (aux, n); }

template<typename C, typename I, typename MV, typename Crter>
static void
mat_inverse_crt (C* d, nat d_rs, nat d_cs, nat r, nat c,
                 const matrix<I,MV>* s, Crter& crter) {
  nat n= N(crter);
  I* aux= mmx_new<I> (n);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      for (nat k= 0; k < n ; k++) aux[k]= s[k](i,j);
      inverse_crt (d[i * d_rs + j * d_cs], aux, crter);
    }
  mmx_delete<I> (aux, n); }

template<typename C, typename Modulus, typename MW, typename MV, typename Crter>
static void
mat_inverse_crt (C* d, nat d_rs, nat d_cs, nat r, nat c,
                 const matrix<modular<Modulus,MW>,MV>* s, Crter& crter) {
  typedef modular<Modulus,MW> Modular;
  typedef typename Modular::modulus::base I;
  nat n= N(crter);
  I* aux= mmx_new<I> (n);
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < c; j++) {
      for (nat k= 0; k < n ; k++) aux[k]= * s[k](i,j);
      inverse_crt (d[i * d_rs + j * d_cs], aux, crter);
    }
  mmx_delete<I> (aux, n); }

template<typename C, typename Crter> static void
mul (C* d, const C* s1, const C* s2, nat r, nat l, nat c,
     Crter& crter) {
  typedef typename Crter::modulus_base I;
  typedef modulus<I,typename Crter::modulus_base_variant> Modulus;
  typedef modular<Modulus,modular_matrix_crt<C> > Modular;
  typedef matrix<Modular> Matrix_modular;
  nat n= N(crter);
  Matrix_modular* mm1= mmx_new<Matrix_modular> (n);
  Matrix_modular* mm2= mmx_new<Matrix_modular> (n);
  Matrix_modular* mmd= mmx_new<Matrix_modular> (n);
  mat_direct_crt (mm1, s1, Mat::index (1, 0, r, l), Mat::index (0, 1, r, l),
                  r, l, crter);
  mat_direct_crt (mm2, s2, Mat::index (1, 0, l, c), Mat::index (0, 1, l, c),
                  l, c, crter);
  for (nat k= 0; k < n; k++) {
    Modular::set_modulus (crter[k]);
    mmd[k]= mm1[k] * mm2[k];
  }
  mat_inverse_crt (d, Mat::index (1, 0, r, c), Mat::index (0, 1, r, c),
                   r, c, mmd, crter);
  mmx_delete<Matrix_modular> (mm1, n);
  mmx_delete<Matrix_modular> (mm2, n);
  mmx_delete<Matrix_modular> (mmd, n); }

template<typename C, typename S, typename CV> static void
mul (C* d, const C* s1, const C* s2, nat r, nat l, nat c,
     crt_naive_transformer<C,S,CV>& crter) {
  typedef implementation<crt_project,CV> Crt;
  typedef crt_naive_transformer<C,S,CV> Crter;
  typedef typename Crter::modulus_base I;
  typedef modulus<I,typename Crter::modulus_base_variant> Modulus;
  typedef modular<Modulus,modular_matrix_crt<C> > Modular;
  typedef typename Matrix_variant(Modular) MV;
  typedef implementation<matrix_multiply,MV> Mat_mod;

  nat spc1= aligned_size<Modular,V> (r * l);
  nat spc2= aligned_size<Modular,V> (l * c);
  nat spcd= aligned_size<Modular,V> (r * c);
  nat spc= spc1 + spc2 + spcd;
  Modular* x1= mmx_new<Modular> (spc);
  Modular* x2= x1 + spc1;
  Modular* xd= x2 + spc2;

  if (N(crter) == 1) {
    Modulus p (crter[0]); Modular::set_modulus (p);
    for (nat i= 0; i < r * l; i++)
      x1[i]= Modular (Crt::encode (s1[i], p), true);
    for (nat i= 0; i < l * c; i++)
      x2[i]= Modular (Crt::encode (s2[i], p), true);
    Mat_mod::mul (xd, x1, x2, r, l, c);
    for (nat i= 0; i < r * c; i++)
      d[i]= Crt::decode (C(* xd[i]), crter.P, crter.H);
  }
  else {
    for (nat k= 0; k < N(crter); k++) {
      Modulus p (crter[k]); Modular::set_modulus (p);
      for (nat i= 0; i < r * l; i++)
        x1[i]= Modular (Crt::mod (Crt::encode (s1[i], crter.P), p));
      for (nat i= 0; i < l * c; i++)
        x2[i]= Modular (Crt::mod (Crt::encode (s2[i], crter.P), p));
      Mat_mod::mul (xd, x1, x2, r, l, c);
      I m (crter.m[k]), t; C q (crter.q[k]);
      for (nat i= 0; i < r * c; i++) {
        mul_mod (t, m, * xd[i], p);
        if (k == 0) d[i]= t * q; else mul_add (d[i], t, q);
      }
    }
    for (nat i= 0; i < r * c; i++)
      d[i]= Crt::decode (Crt::mod (d[i], crter.P), crter.P, crter.H);
  }
  mmx_delete<Modular> (x1, spc); }

template<typename C, typename Low, typename High, nat s, typename CV>
static void
mul (C* d, const C* s1, const C* s2, nat r, nat l, nat c,
     crt_blocks_transformer<Low,High,s,CV>& crter) {
  typedef matrix<C,matrix_crt<W> > Matrix;
  nat n= crter.high -> size ();
  if (n == 1) {
    mul (d, s1, s2, r, l, c, * crter.low[0]);
    return;
  }
  Matrix* mm1= mmx_new<Matrix> (n);
  Matrix* mm2= mmx_new<Matrix> (n);
  Matrix* mmd= mmx_new<Matrix> (n);
  mat_direct_crt (mm1, s1, Mat::index (1, 0, r, l), Mat::index (0, 1, r, l),
                  r, l, * crter.high);
  mat_direct_crt (mm2, s2, Mat::index (1, 0, l, c), Mat::index (0, 1, l, c),
                  l, c, * crter.high);
  for (nat k= 0; k < n; k++)
    mul (tab(mmd[k]), tab(mm1[k]), tab(mm2[k]), r, l, c, * crter.low[k]);
  mat_inverse_crt (d, Mat::index (1, 0, r, c), Mat::index (0, 1, r, c),
                   r, c, mmd, * crter.high);
  mmx_delete<Matrix> (mm1, n);
  mmx_delete<Matrix> (mm2, n);
  mmx_delete<Matrix> (mmd, n); }

template<typename C> static void
mul (C* d, const C* s1, const C* s2, nat r, nat l, nat c) {
  typedef matrix_crt_multiply_helper<C> Matrix_crt;
  typedef typename Matrix_crt::crt_transformer Crter;
  typedef typename Matrix_crt::moduli_sequence Sequence;
  typedef typename Crter::modulus_base I;
  typedef modulus<I,typename Crter::modulus_base_variant> Modulus;
  static const nat dim_thr= Matrix_crt::dimension_threshold;
  static const nat ratio_thr= Matrix_crt::ratio_threshold;;
  
  if (r <= dim_thr || l <= dim_thr || c <= dim_thr) {
    Mat::mul (d, s1, s2, r, l, c);
    return;
  }
  nat sz= matrix_crt_multiply_helper<C>
    ::size (s1, Mat::index (1, 0, r, l), Mat::index (0, 1, r, l),
            s2, Mat::index (1, 0, l, c), Mat::index (0, 1, l, c),
            r, l, c);
  nat wd= (ratio_thr * sz) / 100; 

  if (r < wd || l < wd || c < wd) {
    Mat::mul (d, s1, s2, r, l, c);
    return;
  }
  vector<Modulus> mods;
  if (! Sequence::covering (mods, sz))
    // Not enough moduli can be generated
    Mat::mul (d, s1, s2, r, l, c);
  else {
    Crter crter (mods, false);
    mul (d, s1, s2, r, l, c, crter); } }

}; // implementation<matrix_multiply,V,matrix_crt<W> >

} // namespace mmx
#endif // __MMX_MATRIX_CRT_HPP