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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_blocks.hpp
* DESCRIPTION: Crt by blocks
* 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_BLOCKS__HPP
#define __MMX__CRT_BLOCKS__HPP
#include <algebramix/crt_dicho.hpp>

namespace mmx {

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

#define Crt_blocks_variant(C) crt_blocks_variant_helper<C>::CV

template<typename V>
struct crt_blocks: public V {};

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

template<typename C>
struct crt_blocks_variant_helper {
  typedef crt_blocks<typename crt_naive_variant_helper<C>::CV> CV;
};

/******************************************************************************
* Stack of two CRTs
******************************************************************************/

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

  template<typename I, typename C, typename CL, typename CH>
  static inline void
  direct (I* c, const C& a, CL** low, CH* high,
          C* aux, nat m, nat s) {
    high -> direct_transform (aux, a);
    for (nat i= 0; i < m; i++) 
      low[i] -> direct_transform (c + i * s, aux[i]); }
  
  template<typename C, typename I, typename CL, typename CH>
  static inline void
  combine (C& a, const I* c, CL* low, CH* high, C* aux, nat m, nat s) {
    for (nat i= 0; i < m; i++) 
      low[i] -> combine (aux[i], c + i * s);
    high -> combine (a, aux); }

  template<typename C, typename I, typename CL, typename CH>
  static inline void
  inverse (C& a, const I* c, CL* low, CH* high, C* aux, nat m, nat s) {
    for (nat i= 0; i < m; i++) 
      low[i] -> inverse_transform (aux[i], c + i * s);
    high -> inverse_transform (a, aux); }
};

/******************************************************************************
* Default threshold
******************************************************************************/

struct crt_blocks_threshold {};

template<typename C>
struct threshold_helper<C,crt_blocks_threshold> {
  typedef fixed_value<nat,32> impl;
};

/******************************************************************************
* Transformer
******************************************************************************/
  
template<typename WL, typename WH,
         nat s=Threshold(typename WH::base,crt_blocks_threshold),
         typename V=typename Crt_blocks_variant(typename WH::base) >
struct crt_blocks_transformer {

  typedef typename WL::base base;
  typedef typename WL::modulus_variant modulus_variant;
  typedef typename WL::modulus_base modulus_base;
  typedef typename WL::modulus_base_variant modulus_base_variant;

private:
  typedef base C;
  typedef modulus_base I;
  typedef modulus<I, modulus_base_variant> M;
  typedef modulus<C, modulus_variant> Modulus;
  typedef implementation<crt_transform,V> Crt;

public:
  
  nat n, m;
  WL** low;
  WH* high;
  C* aux;
  Modulus P;
  C H;
  bool setup;

  inline void setup_comoduli () {}

  inline void
  setup_inverse () {
    if (setup) return;
    high -> setup_inverse ();
    for (nat i= 0; i < m; i++)
      low[i] -> setup_inverse ();
    setup= true; }

  inline crt_blocks_transformer (const vector<M>& p, bool lazy= true) {
    n= N(p); m= (n + s - 1) / s; setup= !lazy;
    if (n == 0) { P= 1; H= 0; return; }
    low= mmx_new<WL*> (m);
    vector<Modulus> v (Modulus (), m);
    for (nat i= 0; i < m; i++) {
      low[i]= mmx_new<WL> (1, range (p, i * s, min (n, (i+1) * s)), lazy);
      v[i]= low[i] -> product ();
    }
    high= mmx_new<WH> (1, v, lazy);
    P= high -> product ();
    H= Crt::half (P);
    aux= mmx_new<C> (m); }

  inline ~crt_blocks_transformer () {
    if (n == 0) return;
    for (nat i= 0; i < m; i++) mmx_delete<WL> (low[i], 1);
    mmx_delete<WL*> (low, m);
    mmx_delete<WH> (high, 1);
    mmx_delete<C> (aux, m); }

  inline M operator[] (nat i) const {
    VERIFY (i < n, "index out of range");
    return low[i / s]-> operator [] (i % s); }

  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 high -> comodulus(i / s) * low[i / s] -> comodulus(i % s); }
    
  inline void direct_transform (I* c, const C& a) {
    if (n == 0) return;
    C b= Crt::encode (a, P);
    Crt::direct (c, b, low, high, aux, m, s); }

  inline void combine (C& a, const I* c) {
    if (n == 0) { a= 0; return; }
    setup_comoduli ();
    Crt::combine (a, c, low, high, aux, m, s); }

  inline void inverse_transform (C& a, const I* c) {
    if (n == 0) { a= 0; return; }
    setup_inverse ();
    Crt::inverse (a, c, low, high, aux, m, s);
    a= Crt::decode (a, P, H); }

};

template<typename WL, typename WH, nat s, typename V> inline nat
N (const crt_blocks_transformer<WL,WH,s,V>& crter) {
  return crter.size (); 
}

} // namespace mmx
#endif //__MMX__CRT_BLOCKS__HPP