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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_naive.hpp
* DESCRIPTION: Naive implementations for Chinese remaindering
* COPYRIGHT  : (C) 2009  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__CRT_NAIVE__HPP
#define __MMX__CRT_NAIVE__HPP
#include <basix/vector.hpp>
#include <numerix/modulus.hpp>

namespace mmx {

/******************************************************************************
* Sequence of coprime moduli
******************************************************************************/

struct coprime_moduli_sequence_naive {
  template<typename M> static bool
  extend (vector<M>& v, nat i) { return false; }
};

template<typename M>
struct coprime_moduli_helper {
  typedef coprime_moduli_sequence_naive sequence;
};

#define Coprime_moduli_variant(M) coprime_moduli_helper<M>::sequence

template<typename M, typename V= typename Coprime_moduli_variant(M)>
struct coprime_moduli_sequence {
  vector<M> v;

  coprime_moduli_sequence () : v () {}
  coprime_moduli_sequence (const coprime_moduli_sequence& x) : v (x.v) {}

  M operator [] (nat i) {
    if (i >= N(v)) {
      bool b= V::extend (v, i+1);
      if (!b) return M(0);
    }
    return v[i]; }
};

template<typename M, typename V> vector<M>
range (coprime_moduli_sequence<M, V>& seq, nat beg, nat end) {
  if (end <= beg) return vector<M> (M (), 0);
  vector<M> v (M (), end - beg);
  for (nat i= beg; i < end; i++) {
    if (seq[i] == M(0)) return vector<M> ();
    v[i-beg]= seq[i];
  }
  return v;
}

/******************************************************************************
* Covering in size s for type C with moduli of type M
******************************************************************************/

template<typename C, typename M=modulus<C>,
         typename W= typename Coprime_moduli_variant(M)>
struct moduli_helper {
  template<typename V> static bool
  covering (vector<M, V>& v, nat s) {
    static coprime_moduli_sequence<M,W> seq;
    v= vector<M,V> (M (0), s);
    for (nat i= 0; i < s; i++) {
      v[i]= seq[i];
      if (v[i] == 0) return false;
    }
    return true; }
};

/******************************************************************************
* Standard parameters
******************************************************************************/

template<typename C>
struct std_crt_naive {
  typedef C                           base;
  typedef C                           modulus_base;
  typedef typename Modulus_variant(C) modulus_base_variant;
  typedef typename Modulus_variant(C) modulus_variant;
};

/******************************************************************************
* Naive variant for Crt
******************************************************************************/

#define Crt_naive_variant(C) crt_naive_variant_helper<C>::CV

struct crt_naive {};

template<typename C>
struct crt_naive_variant_helper {
  typedef crt_naive CV;
};

/******************************************************************************
* Projection and lifting
******************************************************************************/

struct crt_project {};

template<typename V>
struct implementation<crt_project,V,crt_naive> {

  template<typename M> static inline typename M::base
  half (const M& P) {
    // return *P >> 1 for signed integer types
    return typename M::base (0); }

  template<typename C, typename M> static inline typename M::base
  encode (const C& a, const M& p) {
    (void) p;
    return as<typename M::base> (a); }

  template<typename C, typename M> static inline C
  decode (const C& a, const M& P, const C& H) {
    (void) P; (void) H;
    return a; }

  template<typename C, typename M> static inline typename M::base
  mod (const C& a, const M& p) {
    return as<typename M::base> (rem (a, *p)); }

};

/******************************************************************************
* Basic operations for signed structures
******************************************************************************/

template<typename V> struct crt_signed {};

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

template<typename V,typename W>
struct implementation<crt_project,V,crt_signed<W> > :
    public implementation<crt_project,V,W> {

  template<typename M> static inline typename M::base 
  half (const M& P) {
    return *P >> 1; }

  template<typename C, typename M> static inline typename M::base
  encode (const C& a, const M& p) {
    typename M::base b;
    if (a < 0) {
      b= as<typename M::base> (-a); neg_mod (b, p);
    }
    else 
      b= as<typename M::base> (a);
    return b; }

  template<typename C, typename M> static inline C
  decode (const C& a, const M& P, const C& H) {
    return a > H ? a - *P : a; }

};

/******************************************************************************
* Direct and inverse transforms
******************************************************************************/

struct crt_transform {};

template<typename V>
struct implementation<crt_transform,V,crt_naive> :
    public implementation<crt_project,V> {
  typedef implementation<crt_project,V> Crt;

  template<typename C, typename M, typename I> static inline void
  direct (I* c, const C& a, const M* p, nat n) {
    for (nat i= 0; i < n; i++)
      c[i]= Crt::mod (a, p[i]); }

  template<typename C, typename I> static inline void
  combine (C& a, const I* c, const C* q, nat n) {
    a= 0;
    for (nat i= 0; i < n; i++)
      mul_add (a, c[i], q[i]); }

  template<typename C, typename M, typename I, typename K> static inline void
  inverse (C& a, const I* c, const M* p, const K& P,
           const C* q, const I* m, nat n) {
    a= 0; I t;
    for (nat i= 0; i < n; i++) {
      mul_mod (t, m[i], c[i], p[i]);
      mul_add (a, t, q[i]);
    }
    a= Crt::mod (a, P); }
};

/******************************************************************************
* The naive CRT transformer class
******************************************************************************/

template<typename C, typename S=std_crt_naive<C>,
         typename V=typename Crt_naive_variant(C) >
struct crt_naive_transformer : public S {

  typedef typename S::modulus_base I;
  typedef modulus<I, typename S::modulus_base_variant> M;
  typedef modulus<C, typename S::modulus_variant> Modulus;
  typedef implementation<vector_linear, vector_naive> Vec;
  typedef implementation<crt_transform,V> Crt;

  nat n;             // number of the moduli
  nat l_p, l_m, l_q; // allocated sizes
  M* p;              // sequence of the moduli
  Modulus P;      // product of all the p[i]
  C H;        // half of P -- used for integers
  C* q;       // sequence of P / p[i]
  I* m;       // sequence of inverse (q[i]) mod p[i]
  bool product_done, comoduli_done, inverse_done;

  inline void setup_product () {
    if (product_done) return;
    q= mmx_formatted_new<C> (l_q, format<C> (/*FIXME*/));
    for (nat i= 0; i < n; i++) q[i]= * p[i];
    P= Vec::template vec_unary_big<mul_op> (q, n);
    product_done= true; }
  
  inline void setup_comoduli () {
    if (comoduli_done) return;
    setup_product ();
    for (nat i= 0; i < n; i++) q[i]= * P / * p[i];
    comoduli_done= true; }

  inline void setup_inverse () {
    if (inverse_done) return;
    setup_comoduli ();
    H= Crt::half (P);
    m= mmx_formatted_new<I> (l_m, format<I> (/*FIXME*/));
    for (nat i= 0; i < n; i++) {
      I a (1), b;
      for (nat j= 0; j < n; j++) {
        reduce_mod (b, * p[j], p[i]);
        if (j != i) mul_mod (a, b, p[i]);
      }
      inv_mod (m[i], a, p[i]);
      ASSERT (m[i] != 0, "moduli must be pairwise coprime"); }
    inverse_done= true; }

  template<typename NV,typename VV> 
  inline crt_naive_transformer (const vector<NV,VV>& v, bool lazy= true) 
    : n(N(v)) {
    l_p= default_aligned_size<M> (n);
    l_m= default_aligned_size<I> (n);
    l_q= default_aligned_size<C> (n);
    p= mmx_formatted_new<M> (l_p, format<M> (/*FIXME*/));
    for (nat i= 0; i < n; i++) p[i]= as<M> (v[i]);
    q= NULL; m= NULL; P= C(1);
    product_done= false; comoduli_done= false; inverse_done= false;
    if (! lazy) setup_inverse (); }

  inline ~crt_naive_transformer () {
    mmx_delete<M> (p, l_p);
    if (q != NULL) mmx_delete<C> (q, l_q);
    if (m != NULL) mmx_delete<I> (m, l_m); }

  inline M operator[] (nat i) const {
    VERIFY (i < n, "index out of range");
    return p[i]; }

  inline nat size () const {
    return n; }

  inline Modulus product () {
    setup_product (); return P; }

  inline C comodulus (nat i) {
    VERIFY (i < n, "index out of range");
    setup_comoduli (); return q[i]; }
    
  inline void direct_transform (I* c, const C& a) const {
    C b= Crt::encode (a, P);
    Crt::direct (c, b, p, n); }

  inline void combine (C& a, const I* c) {
    setup_comoduli ();
    Crt::combine (a, c, q, n); }

  inline void inverse_transform (C& a, const I* c) {
    setup_inverse ();
    Crt::inverse (a, c, p, P, q, m, n);
    a= Crt::decode (a, P, H); }
};

template<typename C, typename M, typename V> inline nat
N (const crt_naive_transformer<C,M,V>& crter) {
  return crter.size (); 
}

/******************************************************************************
* Higher interface level -- common to all transformers
******************************************************************************/

#define C typename Crter::base
#define I typename Crter::modulus_base
#define M modulus<typename Crter::modulus_base,\
                  typename Crter::modulus_base_variant>
#define Modulus modulus<C, typename Crter::modulus_variant>

template<typename Crter> inline Modulus
moduli_product (Crter& crter) {
  return crter.product (); }

template<typename Crter> inline typename Crter::base
comodulus (const Crter& crter, nat i) {
  return crter.comodulus (i);
}

template<typename Crter> inline void
direct_crt (I* dest, const C& s, Crter& crter) {
  crter.direct_transform (dest, s);
}

template<typename Crter, typename W> inline void
direct_crt (vector<I,W>& dest, const C& s, Crter& crter) {
  nat l= aligned_size<I,W> (N(crter));
  I* tmp= mmx_formatted_new<I> (l, CF(dest));
  crter.direct_transform (tmp, s);
  dest= vector<I,W> (tmp, N(crter), l, CF(dest));
}

template<typename Crter> inline vector<I>
direct_crt (const C& s, Crter& crter) {
  vector<I> dest;
  direct_crt (dest, s, crter);
  return dest;
}

template<typename Crter> inline void
combine_crt (C& d, const I* src, Crter& crter) {
  crter.combine (d, src);
}

template<typename Crter, typename W> inline void
combine_crt (C& d, const vector<I,W>& src, Crter& crter) {
  crter.combine (d, seg (src));
}

template<typename Crter, typename W> inline C
combine_crt (const vector<I,W>& src, Crter& crter) {
  C d; combine_crt (d, seg (src), crter);
  return d;
}

template<typename Crter> inline void
inverse_crt (C& d, const I* src, Crter& crter) {
  crter.inverse_transform (d, src);
}

template<typename Crter, typename W> inline void
inverse_crt (C& d, const vector<I,W>& src, Crter& crter) {
  crter.inverse_transform (d, seg (src));
}

template<typename Crter, typename W> inline C
inverse_crt (const vector<I,W>& src, Crter& crter) {
  C d; inverse_crt (d, seg (src), crter);
  return d;
}

#undef C
#undef I
#undef M
#undef Modulus
} // namespace mmx
#endif //__MMX__CRT_NAIVE__HPP