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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_dicho.hpp
* DESCRIPTION: Multimod and Chinese remaindering with dichotomic algorithms
* 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_DICHO__HPP
#define __MMX__CRT_DICHO__HPP
#include <algebramix/crt_naive.hpp>

namespace mmx {

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

template<typename C>
struct std_crt_dicho : std_crt_naive<C> {};

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

#define Crt_dicho_variant(C) crt_dicho_variant_helper<C>::CV

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

template<typename C>
struct crt_dicho_variant_helper {
  typedef crt_dicho<typename crt_naive_variant_helper<C>::CV> CV;
};

/******************************************************************************
* Dichotomic CRT implementations
******************************************************************************/

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

  template<typename C, typename Modulus> static inline void 
  down (C* c, const C& a, Modulus** tree, nat* k, nat depth, nat n) {
    if (n == 0) return;
    c[0]= Crt::mod (a, tree[depth-1][0]);
    if (depth == 1) return;
    for (nat i= depth-1; i > 0; i--) {
      if (k[i-1] & 1) c[k[i-1] - 1]= c[k[i] - 1];
      for (nat j= (k[i-1] >> 1) - 1; j + 1 > 0; j--) {
        c[(j << 1) + 1]= Crt::mod (c[j], tree[i-1][(j << 1) + 1]);
        c[(j << 1)    ]= Crt::mod (c[j], tree[i-1][(j << 1)    ]); } } }

  template<typename C, typename Modulus> static inline void
  up (C& a, C* c, Modulus** tree, nat* k, nat depth, nat n) {
    if (n == 0) return;
    for (nat i= 0; i < depth - 1; i++) {
      for (nat j= 0; j < (k[i] >> 1); j++)
        c[j]= c[j << 1]       * * tree[i][(j << 1) + 1] +
              c[(j << 1) + 1] * * tree[i][j << 1];
      if (k[i] & 1) c[k[i+1] - 1]= c[k[i] - 1]; }
    a= c[0]; }

  template<typename C, typename J, typename Modulus> static inline void
  direct (J* c, const C& a, 
          Modulus** tree, nat* k, nat depth, nat n, C* aux) {
    down (aux, a, tree, k, depth, n);
    for (nat i= 0; i < n; i++) c[i]= as<J> (aux[i]); }
  
  template<typename C, typename Modulus> static inline void
  direct (C* c, const C& a,
          Modulus** tree, nat* k, nat depth, nat n, C* aux) {
    (void) aux; down (c, a, tree, k, depth, n); }
  
  template<typename C, typename I, typename Modulus> static inline void
  combine (C& a, const I* c, Modulus** tree, nat* k, nat depth, nat n,
           C* aux) {
    if (n == 0) { a= C(0); return; }
    for (nat i= 0; i < n; i++) aux[i]= C (c[i]);
    up (a, aux, tree, k, depth, n); }

  template<typename C, typename I, typename Modulus> static inline void
  inverse (C& a, const I* c,
           Modulus** tree, nat* k, nat depth,
           C* m, nat n, const Modulus& P, C* aux) {
    if (n == 0) { a= C(0); return; }
    for (nat i= 0; i < n; i++)
      aux[i]= Crt::mod (C (c[i]) * m[i], tree[0][i]);
    up (a, aux, tree, k, depth, n);
    a= Crt::mod (a, P); }

};

/******************************************************************************
* Transformer
******************************************************************************/
    
template<typename C, typename S=std_crt_dicho<C>,
         typename V=typename Crt_dicho_variant(C)>
struct crt_dicho_transformer : public S {
private:
  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;
  M* p;     // moduli
  modulus<C>** tree; // subproduct tree of depth 'depth' with 'n' leafs
  nat depth;// depth of the tree
  nat* k;   // size of each level of the tree
  Modulus P;
  C H;
  C* m;     // sequence of inverse P/p[i] mod p[i]
  C* aux;

  inline void setup_direct (const M* _p, nat nn) {
    m= NULL; // setup only if necessary
    n= nn;
    if (n == 0) { depth= 0; P= 1; H= 0; return; }
    p= mmx_new<M> (n);
    for (nat i= 0; i < n; i++) p[i]= _p[i];
    depth= log_2 (next_power_of_two (n)) + 1;
    tree= mmx_new<Modulus*> (depth);
    k= mmx_new<nat> (depth);
    k[0]= n;
    tree[0]= mmx_new<Modulus> (k[0]);
    for (nat j= 0; j < k[0]; j++) 
      tree[0][j]= C(* p[j]);
    for (nat i= 1; i < depth; i++) {
      k[i]= (k[i-1] >> 1) + (k[i-1] & 1);
      tree[i]= mmx_new<Modulus> (k[i]);
      for (nat j= 0; j < (k[i-1] >> 1); j++)
        tree[i][j]= * tree[i-1][j << 1] * * tree[i-1][(j << 1) + 1];
      if (k[i-1] & 1) tree[i][k[i] - 1]= * tree[i-1][k[i-1] - 1];
    }
    P= tree[depth-1][0];
    H= Crt::half (P); }

public:

  inline void setup_inverse () {
    typedef implementation<vector_linear,vector_naive> Vec;
    if (n == 0 || m != NULL) return;
    m= mmx_new<C> (n);
    C c;
    Vec::set (m, C(1), n);
    Crt::up (c, m, tree, k, depth, n);
    Crt::down (m, c, tree, k, depth, n);
    for (nat i= 0; i < n; i++)
      inv_mod (m[i], tree[0][i]); }

  inline crt_dicho_transformer (const vector<M>& p, bool lazy= true)
    : tree(NULL), k(NULL) {
    setup_direct (seg (p), N(p));
    aux= mmx_new<C> (n);
    if (! lazy) setup_inverse (); }

  inline ~crt_dicho_transformer () {
    if (n == 0) return;
    mmx_delete<M> (p, n);
    mmx_delete<C> (aux, n);
    for (nat i= 0; i < depth; i++)
      mmx_delete<Modulus> (tree[i], k[i]);
    mmx_delete<Modulus*> (tree, depth);
    mmx_delete<nat> (k, depth);
    if (m != NULL) mmx_delete<C> (m, n); }

  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 () const {
    return P; }

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

  inline void combine (C& a, const I* c) {
    Crt::combine (a, c, tree, k, depth, n, aux); }

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

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

} // namespace mmx
#endif //__MMX__CRT_DICHO__HPP