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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_integer.hpp
* DESCRIPTION: Subroutines for efficient Chinese remaindering
* COPYRIGHT  : (C) 2008  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_CRT_INTEGER_HPP
#define __MMX_CRT_INTEGER_HPP
#include <numerix/modular_integer.hpp>
#include <algebramix/crt_naive.hpp>
#include <algebramix/crt_dicho.hpp>
#include <algebramix/crt_blocks.hpp>
#include <algebramix/crt_int.hpp>

namespace mmx {

/******************************************************************************
* Naive variant
******************************************************************************/

DEFINE_VARIANT(crt_naive_integer,crt_signed<crt_naive>)

#if defined(__GNU_MP__)
template<typename V>
struct implementation<crt_project,V,crt_naive_integer> :
    public implementation<crt_project,crt_signed<crt_naive> >
{
  template<typename C, typename I, typename W> static inline I
  mod (const C& a, const modulus<I,W>& p) {
    static integer r;
    return (I) mpz_fdiv_r_ui (*r, *a, *p); }

  template<typename C, typename W> static inline C
  mod (const C& a, const modulus<integer,W>& p) {
    return rem (a, *p); }
};
#endif

STMPL
struct crt_naive_variant_helper<integer> {
  typedef crt_naive_integer CV; };

DEFINE_VARIANT(crt_dicho_integer,crt_dicho<crt_naive_integer>)

STMPL
struct crt_dicho_variant_helper<integer> {
  typedef crt_dicho_integer CV; };

/******************************************************************************
* Default parameters
******************************************************************************/

STMPL
struct std_crt_naive<integer> {
  typedef integer base;
  typedef int modulus_base;
  typedef modulus_int_preinverse<8*sizeof(int)-9> modulus_base_variant;
  typedef Modulus_variant(integer) modulus_variant;
};

STMPL
struct std_crt_dicho<integer> {
  typedef integer base;
  typedef integer modulus_base;
  typedef Modulus_variant(integer) modulus_base_variant;
  typedef Modulus_variant(integer) modulus_variant;
};

/******************************************************************************
* Covering
******************************************************************************/

template<typename M, typename W>
struct moduli_helper<integer,M,W> :
    moduli_signed_integer_helper<integer,M,W> {};

template <typename M, nat t>
struct moduli_helper<integer, M, fft_prime_sequence_int<t> > {
  template<typename V> static bool
  covering (vector<M, V>& v, nat s) {
    static coprime_moduli_sequence<M, fft_prime_sequence_int<t> > seq;
    if (t < 2) {
      v= vector<M, V> ();
      return s == 0;
    }
    nat n= (s + t) / (t - 1);
    v= range (seq, 0, n);
    return N(v) == n; }
};

/******************************************************************************
* Special variant for GMP
******************************************************************************/

#if defined(__GNU_MP__)
#define C integer

vector<mp_limb_t> mpz_setup_crt (const vector<mp_limb_t>& mods);
// Build auxiliary conversion information for the moduli p[i]

integer mpz_moduli_product (const vector<mp_limb_t>& mods, 
                            const vector<mp_limb_t>& cv);
// Product of all the moduli

void mpz_encode_crt (mp_limb_t* dest, const integer& src,
                     const vector<mp_limb_t>& mods,
                     const vector<mp_limb_t>& cv);
// Compute src mod p[i] for each of the moduli p[i]

integer mpz_decode_crt (const mp_limb_t* src,
                        const vector<mp_limb_t>& mods,
                        const vector<mp_limb_t>& cv);
// Reconstruct the integer from its remainders w.r.t. p[i] for each i

template<typename I>
struct crt_gmp_transformer {
  typedef C base;
  typedef int modulus_base;
  typedef modulus_int_preinverse<23> modulus_base_variant;
  typedef typename Modulus_variant(C) modulus_variant;
private:
  typedef modulus<modulus_base,modulus_base_variant> M;

  nat n;
  vector<mp_limb_t> mods;
  vector<mp_limb_t> cv;
  modulus<C> P;
  mp_limb_t* aux;  

public:
  inline void setup_inverse () {}
  
  inline crt_gmp_transformer (const vector<M>& p, bool lazy=true) {
    (void) lazy;
    ASSERT (sizeof (I) <= sizeof (mp_limb_t), "moduli too large");
    n= N(p); mods= vector<mp_limb_t> ((mp_limb_t) 0, n);
    for (nat i= 0; i < n; i++) mods[i]= (mp_limb_t) * p[i];
    cv= mpz_setup_crt (mods);
    P= mpz_moduli_product (mods, cv);
    if (sizeof (I) != sizeof (mp_limb_t))
      aux= mmx_new<mp_limb_t> (n); }

  inline ~crt_gmp_transformer () {
    if (sizeof (I) != sizeof (mp_limb_t))
      mmx_delete<mp_limb_t> (aux, n); }

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

  inline nat size () const {
    return n; }

  inline modulus<C> product () const {
    return P; }

  inline C comodulus (nat i) {
    VERIFY (i < n, "index out of range");
    return (* P) / mods[i]; }

  inline void direct_transform (I* c, const C& a) const {
    mp_limb_t* tmp= sizeof (I) == sizeof (mp_limb_t) ?
      (mp_limb_t*) (void*) c : aux;
    mpz_encode_crt (tmp, a, mods, cv);
    if (sizeof (I) != sizeof (mp_limb_t))
      for (nat i= 0; i < n; i++) c[i]= (I) tmp[i]; }
  
  inline void inverse_transform (C& a, const I* c) {
    mp_limb_t* tmp;
    if (sizeof (I) == sizeof (mp_limb_t))
      tmp= (mp_limb_t*) (void*) c;
    else {
      for (nat i= 0; i < n; i++) aux[i]= (mp_limb_t) c[i];
      tmp= aux;
    }
    a= mpz_decode_crt (tmp, mods, cv); }
};

template<typename I> inline nat
N (const crt_gmp_transformer<I>& crter) {
  return crter.size (); 
}

#undef C
#endif // __GNU_MP__

} // namespace mmx
#endif // __MMX_CRT_INTEGER_HPP