/Users/mourrain/Devel/mmx/numerix/include/numerix/modular_int.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : modular_int.hpp
* DESCRIPTION: modulus for modular arithmetic over hardware integers
* COPYRIGHT  : (C) 2008  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__MODULAR_INT__HPP
#define __MMX__MODULAR_INT__HPP
#include <basix/operators.hpp>
#include <basix/int.hpp>
#include <numerix/integer.hpp>
#include <numerix/modular.hpp>

namespace mmx {


#define TMPL template<typename C, typename M>

/******************************************************************************
* Naive implementation relying on hardware division
******************************************************************************/
// The modulus 'p' is always taken to be positive
// The modulars are in the range [0, p-1]

template<nat size>
struct modulus_maximum_size_int {

  template<typename C, bool is_signed>
  struct _helper {
    static const nat value = (size < 8*sizeof(C) ? size : 8*sizeof(C)-1);
  };

  template<typename C>
  struct _helper<C, false> {
    static const nat value = (size <= 8*sizeof(C) ? size : 8*sizeof(C));
  };

  template<typename C>
  struct maximum_size_helper {
    static const nat value = _helper<C, is_signed_helper<C>::value>::value;
    static inline C dyn_value () { return value; }
  };

  template<typename C>
  struct maximum_value_helper {
    static const nat _size = maximum_size_helper<C>::value;
    typedef typename unsigned_of_helper<C>::type uC;
    static const C value = MMX_SAFE_LEFT_SHIFT_INT (uC, 1, _size) - 1;
    static inline C dyn_value () { return value; }
  };
};

template<typename V>
struct modulus_normalization_int_naive : public V {

  template<typename C> static inline bool
  normalize (C& p) {
    typedef typename unsigned_of_helper<C>::type uC;
    static const uC a = MMX_SAFE_LEFT_SHIFT_INT(uC, 1,
                         V::template maximum_size_helper<C>::value);
    p = abs (p);
    if ((uC) p == a || p == 0) { p = a; return true; }
    return p <= V::template maximum_value_helper<C>::dyn_value ();
  }
};

template<typename V>
struct modulus_reduction_int_naive : public V {
  template<typename C> static inline void 
  reduce_mod_core (C& dest, const C& p) {
    if (p != 0) dest %= p; }

  template<typename C> static inline void 
  reduce_mod_core (C& dest, const C& p, C& carry) {
    if (p != 0) { carry= dest / p; dest %= p; }
    else carry = 0; }

  TMPL static inline void 
  reduce_mod (C& dest, const M& m) {
    reduce_mod_core (dest, (C) m.p); }

  TMPL static inline void 
  reduce_mod (C& dest, const M& m, C& carry) {
    reduce_mod_core (dest, (C) m.p, carry); }

  TMPL static inline void 
  reduce_mod (C& dest, const C& s, const M& m) {
    dest = s;
    reduce_mod_core (dest, (C) m.p); }

  TMPL static inline void 
  reduce_mod (C& dest, const C& s, const M& m, C& carry) {
    dest = s;
    reduce_mod_core (dest, (C) m.p, carry); }
};

template<typename V>
struct modulus_add_int_naive : public V {

  TMPL static inline void
  neg_mod (C& dest, const M& m) {
    if (dest != 0) dest = m.p - dest; }

  TMPL static inline void
  neg_mod (C& dest, const M& m, C& carry) {
    VERIFY (carry == 0 || carry == 1, "unexpected large carry");
    if (dest != 0 || carry != 0) { dest= m.p - dest - carry; carry= 1; } }

  TMPL static inline void
  neg_mod (C& dest, const C& s, const M& m) {
    if (s != 0) dest = m.p - s; else dest = s; }

  TMPL static inline void
  neg_mod (C& dest, const C& s, const M& m, C& carry) {
    VERIFY (carry == 0 || carry == 1, "unexpected large carry");
    if (s != 0 || carry != 0) { dest= m.p - s - carry; carry= 1; }
    else dest= 0; }

  template<typename C> static inline void
  add_mod_without_overflow (C& dest, const C& s, const C& p) {
    dest += s;
    if (dest >= p) dest -= p; }

  template<typename C> static inline void
  add_mod_without_overflow (C& dest, const C& s, const C& p, C& carry) {
    VERIFY (carry == 0 || carry == 1, "unexpected large carry");
    dest += s + carry;
    if (dest >= p) { dest -= p; carry= 1; } else carry= 0; }

  template<typename C> static inline void
  add_mod_with_overflow (C& dest, const C& s, const C& p) {
    dest += s;
    if (dest < s) { dest -= p; return; }
    if (dest >= p) dest -= p; }

  template<typename C> static inline void
  add_mod_with_overflow (C& dest, const C& s, const C& p, C& carry) {
    VERIFY (carry == 0 || carry == 1, "unexpected large carry");
    dest += s + carry;
    if (dest < s || dest >= p) { dest -= p; carry= 1; }
    else carry= 0; }

  // Signed integer involve no overflow
  template<typename C, bool b> 
  struct add_mod_helper {
    static inline void
    op (C& dest, const C& s, const C& p) {
      typedef typename unsigned_of_helper<C>::type uC;
      uC t = dest, up = p, us = s;
      add_mod_without_overflow (t, us, up);
      dest = t; }
    static inline void
    op (C& dest, const C& s, const C& p, C& carry) {
      typedef typename unsigned_of_helper<C>::type uC;
      uC t = dest, up = p, us = s, ucarry= carry;
      add_mod_without_overflow (t, us, up, ucarry);
      dest = t; carry= ucarry; }
  };

  // Unsigned integer can cause overflow
  template<typename C> 
  struct add_mod_helper<C,false> {
    static inline void
    op (C& dest, const C& s, const C& p) {
      if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
        add_mod_with_overflow (dest, s, p);
      else
        add_mod_without_overflow (dest, s, p); }
    static inline void
    op (C& dest, const C& s, const C& p, C& carry) {
      if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
        add_mod_with_overflow (dest, s, p, carry);
      else
        add_mod_without_overflow (dest, s, p, carry); }
  };

  template<typename C> static inline void
  add_mod_core (C& dest, const C& s, const C& p) {
    add_mod_helper<C,is_signed_helper<C>::value>::op (dest, s, p); }

  template<typename C> static inline void
  add_mod_core (C& dest, const C& s, const C& p, C& carry) {
    add_mod_helper<C,is_signed_helper<C>::value>::op (dest, s, p, carry); }

  TMPL static inline void
  add_mod (C& dest, const C& s, const M& m) {
    add_mod_core (dest, s, m.p); }
  
  TMPL static inline void
  add_mod (C& dest, const C& s, const M& m, C& carry) {
    add_mod_core (dest, s, m.p, carry); }
  
  TMPL static inline void
  add_mod (C& dest, const C& s1, const C& s2, const M& m) {
    dest = s1;
    add_mod (dest, s2, m); }
  
  TMPL static inline void
  add_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
    dest = s1;
    add_mod (dest, s2, m, carry); }
  
  template<typename C> static inline void
  sub_mod_core (C& dest, const C& s, const C& p) {
    if (dest < s) dest += p;
    dest -= s; }

  template<typename C> static inline void
  sub_mod_core (C& dest, const C& s, const C& p, C& carry) {
    VERIFY (carry == 0 || carry == 1, "unexpected large carry");
    C t = s + carry;
    if (t == p || dest < t) { dest += p; carry= 1; } else carry= 0;
    dest -= t; }

  TMPL static inline void
  sub_mod (C& dest, const C& s, const M& m) {
    sub_mod_core (dest, s, m.p); }

  TMPL static inline void
  sub_mod (C& dest, const C& s, const M& m, C& carry) {
    sub_mod_core (dest, s, m.p, carry); }

  TMPL static inline void
  sub_mod (C& dest, const C& s1, const C& s2, const M& m) {
    dest = s1;
    sub_mod (dest, s2, m); }

  TMPL static inline void
  sub_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
    dest = s1;
    sub_mod (dest, s2, m, carry); }
};
 
template<typename V>
struct modulus_mul_int_naive : public V {

  template <typename C, typename D>
  struct mul_mod_helper {
    static inline void op (C& dest, const C& s, const C& p) {
      D a (dest);       
      a *= s;
      if (p != 0) a %= p;
      dest = (C) a;
    }
    static inline void op (C& dest, const C& s, const C& p, C& carry) {
      D a (dest);       
      a = a * s + carry;
      if (p != 0) { carry= a / p; a %= p; }
      dest = (C) a;
    }
  };

  template <typename C>
  struct mul_mod_helper<C, void> {
    static const nat half_size = 4 * sizeof (C);
    static const nat size = 2 * half_size;
    static const C lo_mask = ((C) -1) >> half_size;
    static const C hi_mask = ((C) -1) << half_size;
    // Let B = (1 << half_size)
    static inline C
    lo (const C& s) { return s & lo_mask; }
    
    static inline C
    hi (const C& s) { return (s & hi_mask) >> half_size; }

    static inline void
    half_lshift_mod(C& dest, const C& p, const C& lo_p,const C& hi_p) {
      // p = hi_p * B + lo_p, return (dest * B) % p
      // It is required that hi_p != 0
      C q = dest / hi_p;
      C r = dest - q * hi_p;
      q *= lo_p;
      q %= p;
      dest = r << half_size;
      V::sub_mod_core (dest, q, p);
    }

    static inline void
    half_lshift_mod(C& dest, const C& p, const C& lo_p,const C& hi_p,
                    C& carry) {
      // p = hi_p * B + lo_p, return (dest * B) % p and (dest * B) / p
      // It is required that hi_p != 0
      C q = dest / hi_p, c = 0; carry = q;
      C r = dest - q * hi_p;
      q *= lo_p; carry -= (q / p);
      q %= p;
      dest = r << half_size;
      V::sub_mod_core (dest, q, p, c);
      if (c != 0) carry--;
    }

    static inline void op (C& dest, const C& s, const C& p) {
      C hi_p = hi (p);
      if (hi_p) {
        C lo_p = lo (p);
        C lo_dest = lo (dest);
        C lo_s = lo (s);
        C hi_dest = hi (dest);
        C hi_s = hi (s);
        dest = lo_dest * lo_s;
        reduce_mod_core (dest, p);
        C t = lo_dest * hi_s;
        half_lshift_mod (t, p, lo_p, hi_p);
        V::add_mod_core (dest, t, p);
        t = hi_dest * lo_s;
        half_lshift_mod (t, p, lo_p, hi_p);
        V::add_mod_core(dest, t, p);
        t = hi_dest * hi_s;
        half_lshift_mod (t, p, lo_p, hi_p);
        half_lshift_mod (t, p, lo_p, hi_p);
        V::add_mod_core(dest, t, p);
      }
      else {
        dest *= s;
        if (p != 0) dest %= p;
      }
    }

    static inline void op (C& dest, const C& s, const C& p, C& carry) {
      C hi_p = hi (p);
      if (hi_p) {
        C lo_p = lo (p);
        C lo_dest = lo (dest);
        C lo_s = lo (s);
        C hi_dest = hi (dest);
        C hi_s = hi (s);
        C c= 0, cc= 0;
        dest = lo_dest * lo_s;
        reduce_mod_core (dest, p, cc);
        C t = lo_dest * hi_s;
        half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
        V::add_mod_core (dest, t, p, c);       cc += c; c= 0;
        t = hi_dest * lo_s;
        half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
        V::add_mod_core(dest, t, p, c);        cc += c; c= 0;
        t = hi_dest * hi_s;
        half_lshift_mod (t, p, lo_p, hi_p, c); cc += c << half_size; c= 0;
        half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
        V::add_mod_core(dest, t, p, c);        cc += c; c= 0;
        V::add_mod_core(dest, carry, p, c);    carry = cc + c;
      }
      else {
        dest = dest * s + carry;
        if (p != 0) { carry= dest / p; dest %= p; }
      }
    }
  };

  TMPL static inline void
  mul_mod (C& dest, const C& s, const M& m) {
    if (is_signed_helper<C>::value) {
      typedef typename unsigned_of_helper<C>::type uC;
      uC t = dest;
      mul_mod (t, (uC) s, m);
      dest = t;
    }
    else {
      typedef typename unsigned_int_with_double_size_helper<C>::type D;
      mul_mod_helper<C,D>::op (dest, s, m.p);
    }
  }

  TMPL static inline void
  mul_mod (C& dest, const C& s, const M& m, C& carry) {
    if (is_signed_helper<C>::value) {
      typedef typename unsigned_of_helper<C>::type uC;
      uC t = dest, ucarry= carry;
      mul_mod (t, (uC) s, m, ucarry);
      dest = t; carry= ucarry;
    }
    else {
      typedef typename unsigned_int_with_double_size_helper<C>::type D;
      mul_mod_helper<C,D>::op (dest, s, m.p, carry);
    }
  }

  TMPL static inline void
  mul_mod (C& dest, const C& s1, const C& s2, const M& m) {
    dest = s1;
    mul_mod (dest, s2, m); }

  TMPL static inline void
  mul_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
    dest = s1;
    mul_mod (dest, s2, m, carry); }

};

template<typename V>
struct modulus_inv_int_naive : public V {

  template<typename C> static void
  inv_mod_signed (C& a, const C& p) {
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename signed_of_helper<C>::type sC;
    sC coa0=0, coa1=1, r0=p, r1=a, q, t;        
    if ((p == 0) && (r1 != 0)) {
      q = (((uC) r0) - ((uC) r1)) / ((uC) r1) + 1;
      t = r0 - q * r1;
      r0 = r1;
      r1 = t;
      t = coa1;
      coa1 = coa0 - q * coa1;
      coa0 = t;
    }
    while (r1 != 0) {
      q = r0 / r1;
      t = r0 - q * r1;
      r0 = r1;
      r1 = t;
      t = coa1;
      coa1 = coa0 - q * coa1;
      coa0 = t;
    }
    if (r0 != 1)
      ERROR ("inv_mod: argument is not invertible");
    a = coa0 < 0 ? p+coa0 : coa0;
  }

  template<typename C> static void
  inv_mod_unsigned (C& a, const C& p) {
    C coa0=0, coa1=1, r0=p, r1=a, q, t; 
    bool toggle = true;
    if ((p == 0) && (r1 != 0)) {
      q = (r0-r1) / r1 + 1;
      t = r0 - q * r1;
      r0 = r1;
      r1 = t;
      t = coa1;
      coa1 = coa0 - q * coa1;
      coa0 = t;
      toggle = !toggle;
    }
    while (r1 != 0) {
      q = r0 / r1;
      t = r0 - q * r1;
      r0 = r1;
      r1 = t;
      t = coa1;
      coa1 = coa0 - q * coa1;
      coa0 = t;
      toggle = !toggle;
    }
    if (r0 != 1)
      ERROR ("inv_mod: argument is not invertible");
    a = (toggle && (coa0 != 0)) ? p+coa0 : coa0;
  }

  TMPL static inline void
  inv_mod (C& a, const M& m) {
    if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
      inv_mod_unsigned(a, m.p);
    else
      inv_mod_signed(a, m.p); }

  TMPL static inline void
  inv_mod (C& dest, const C& s, const M& m) {
    dest = s;
    inv_mod (dest, m); }
};

template<typename V>
struct modulus_div_int_naive : public V {

  TMPL static inline void 
  div_mod (C& dest, const C& s, const M& m) {
    C t = s;
    V::inv_mod (t, m);
    V::mul_mod (dest, t, m); }

  TMPL static inline void 
  div_mod (C& dest, const C& s1, const C& s2, const M& m) {
    dest = s1;
    div_mod (dest, s2, m); }
};

template<typename V>
struct modulus_encoding_int_naive : public V {
  TMPL static inline void 
  encode_mod (C& dest, const C& s, const M& m) {
    typedef typename unsigned_of_helper<C>::type uC;
    if (sign (s) < 0) {
      uC tmp = - ((uC) s);
      V::reduce_mod (tmp, m);
      dest = (C) (((uC) m.p) - tmp);
    }
    else
      dest = s;
      V::reduce_mod (dest, m);
  }
  TMPL static inline void 
  decode_mod (C& dest, const C& s, const M& m) {
    (void) m;
    dest = s; }
};

template<int size>
struct modulus_int_naive:
  public modulus_encoding_int_naive<
         modulus_div_int_naive<
         modulus_inv_int_naive<
         modulus_mul_int_naive<
         modulus_add_int_naive<
         modulus_reduction_int_naive<
         modulus_normalization_int_naive<
         modulus_maximum_size_int<size> > > > > > > > {};

#undef TMPL

/******************************************************************************
* Fast implementation relying on hardware product (see doc/texmacs/modular.tm)
******************************************************************************/

#define TMPL template<typename C>

template<typename V>
struct modulus_mul_int_preinverse : public V {

/*
  Dynamic computation of the auxiliaries
*/
  template <typename C, typename D> 
  struct _dynamic_inverse_helper {
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename unsigned_of_helper<D>::type uD;
    static const nat m = V::template maximum_size_helper<C>::value;
    // p must be normalized
    static inline uC
    op (const uC& p, nat r, nat s, nat t) {
      (void) r;
      if (p == 1 && s+t == (nat) -1) return 0; // only case when s+t < 0
      uD loc_p = (p == 0) ? ((uD) 1) << m : (uD) ((uC) p);
      return (((uD) 1) << (s+t)) / loc_p; }
  };
  
  template <typename C>
  struct _dynamic_inverse_helper<C,void> {
    static const nat n2 = 4 * sizeof(C);
    static const nat n  = 8 * sizeof(C);
    static const nat m  = V::template maximum_size_helper<C>::value;
    typedef typename unsigned_of_helper<C>::type uC;
    typedef long_int_mul_op<C>    long_mul;
    typedef long_int_rshift_op<C> long_rshift;
    typedef long_int_sub_op<C>    long_sub;
    
    static uC
    _q0 (const uC& p, nat r, nat s, nat t, nat s0, nat s1) {
      // p != 1 or 2
      uC q0;
      if (r <= s0)
        q0 = (((uC) 1) << (s0+r-1)) / p;
      else {
        uC p1 = (p == 0) ? ((uC) 1) << s0 : p >> (r-s0);
        uC p0 = p & (((uC) -1) >> (n+s0-r));
        uC d = ((uC) 1) << s1;
        q0 = d / p1;
        uC b = (d % p1) << (r-s0);
        uC c = q0 * p0;
        if (c >= p) { q0--; c -= p; }
        if (b < c) q0--;
      }
      return q0;
    }
    
    static uC
    op (const uC& p, nat r, nat s, nat t) {
      // If p is 1 or 2 then s0 == 0.
      // p must be normalized.
      if (p == 1) {
        if (m+2 <= n) return 2;
        if (m+1 <= n) return 1;
        return 0;
      }
      if (p == 2) {
        if (m+2 <= n) return 4;
        if (m+1 <= n) return 2;
        return 1;
      }
      nat s0 = (s+t-(r-1)) >> 1;
      nat s1 = 2*s0 - 1;
      uC q0 = _q0 (p, r, s, t, s0, s1);
      uC a0 = 0, a1 = 0;
      long_mul::op (a1, a0, q0 * q0, p);
      if (long_rshift::op_b (a1, a0, r)) a0++;
      uC c = (q0 << s0) - a0;
      uC d;
      if (m+1 <= n || s1+r <= n) {
        if (s1+r-1 >= n) d = 0; else d = (((uC) 1) << (s1+r-1));
        d -= c*p;
        if (d >= p) { c++; d-=p; }
      }
      else {
        uC d1, e0, e1;
        long_mul::op (e1, e0, c, p);
        d = 0;
        d1 = ((uC) 1) << (s1+r-1-n);
        long_sub::op (d1, d, e1, e0);
        if ((d1 != 0) || (d >= p)) { c++; d-=p; }
      }
      uC d2 = d << 1;
      c <<= 1;
      if ((d2 >= p) || (d2 < d)) { c++; d = d2-p; } else d = d2;
      if (s1+1 == s+t-(r-1)) return c;
      d2 = d << 1;
      c <<= 1;
      if ((d2 >= p) || (d2 < d)) c++;
      return c;
    }
  };
  
  TMPL static inline nat
  dyn_r (const C& p) {
    typedef typename unsigned_of_helper<C>::type uC;
    static const nat m = V::template maximum_size_helper<C>::value;
    C abs_p = abs (p);
    nat r = bit_size (abs_p);
    return (r == 0) ? m : ( (uC) abs_p == (((uC) 1) << (r-1)) ? r-1 : r );
  }

  TMPL static inline nat
  dyn_s (const C& p, nat r) {
    static const nat m = V::template maximum_size_helper<C>::value;
    return m+2 <= 8 * sizeof(C) ? r-2 : r-1;
  }

  TMPL static inline nat
  dyn_t (const C& p, nat r) {
    static const nat n = 8 * sizeof(C);
    static const nat m = V::template maximum_size_helper<C>::value;
    return (m+2 <= n) ? r+3 : ((m+1 <= n) ? r+1 : r);
  }

  TMPL static inline C
  dyn_q (const C& p, nat r, nat s, nat t) {
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename unsigned_int_with_double_size_helper<C>::type D;
    return _dynamic_inverse_helper<C,D>::op ((uC) p, r, s, t);
  }

/*
  Static computation of the auxiliaries
*/
  template<typename C, C p>
  struct auxiliaries_helper {
    static const nat m = V::template maximum_size_helper<C>::value;
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename unsigned_int_with_double_size_helper<uC>::type uD;

    static const nat n2= 4*sizeof(C);
    static const nat n = 8*sizeof(C);
    static const uC abs_p = (uC) ((p < 0) ? -p : p);
    static const uC a  = MMX_SAFE_LEFT_SHIFT_INT(uC,1,m);
    static const uC up = (abs_p == 0) ? a : abs_p;
    static const nat bitsize = int_bitsize_helper<uC, up>::value;
    static const nat _b = (bitsize == 0) ? n : bitsize-1;
    static const nat r  = (bitsize == 0) ? n :
      ( (up == MMX_SAFE_LEFT_SHIFT_INT(uC,1,_b)) ?
        bitsize-1 : bitsize );
    static const nat s = (m+2 <= n) ? r-2 : r-1;
    static const nat t = (m+2 <= n) ? r+3 : ((m+1 <= n) ? r+1 : r);

    template <typename D, typename Dummy> 
    struct _inverse_helper {
      static const uD _up = (up == 0) ? ((uD) 1) << m : (uD) up;
      static const nat _s = (s+t >= 2*n) ? 2*n : s+t;
      static const uD _t  = (uD) (MMX_SAFE_LEFT_SHIFT_INT(uD,1,_s));
      static const uC value = (up == 1 && m >= n) ?
                              (uC) 0 : (uC) (_t / (uD) _up);
      static inline C dyn_value () { return value; }
    };

    template <typename Dummy> 
    struct _inverse_helper <void, Dummy> {
      static inline C dyn_value () { return dyn_q ((C) up, r, s, t); }
    };

    static inline C q () {
      return _inverse_helper<uD, void>::dyn_value (); }
  };

/*
  Modular product when int of double size are available
*/

  template <typename C, typename D, nat m>
  struct mul_mod_helper {
    static const nat n = 8 * sizeof(C);
    // C and D are unsigned.
    static inline void
    op (C& dest, const C& src, const C& p,
        const C& q, nat r, nat s, nat t) {
      D a (dest);
      a *= src;
      if (m + 2 <= n) {
        C hi_a = a >> s;
        C b = (((D) hi_a) * q) >> t;
        dest = ((C) a) - b * p;
        if (dest >= p) dest -=p;
        return;
      }
      D hi_a = a >> s;
      C b = (hi_a * q) >> t;
      D c = a - ((D) b) * p;
      if (m + 1 <= n) {
        if (c >= p) {
          dest = c - p;
          if (dest >= p)
            dest -= p;
        }
        else
          dest = c;
        return;
      }
      if (c >= p) {
        c -= p;
        if (c >= p) {
          c -= p;
          if (c >= p) 
            c -= p;
        }
      }
      dest = c;
    }

    static inline void
    op (C& dest, const C& src, const C& p, C& carry,
        const C& q, nat r, nat s, nat t) {
      D a (dest);
      a *= src;
      if (m + 2 <= n) {
        C hi_a = a >> s;
        C b = (((D) hi_a) * q) >> t, c = 0;;
        dest = ((C) a) - b * p; 
        if (dest >= p) { dest -= p; b++; }
        V::add_mod_core (dest, carry, p, c); carry = b + c;
        return;
      }
      D hi_a = a >> s;
      C b = (hi_a * q) >> t, z= 0;
      D c = a - ((D) b) * p;
      if (m + 1 <= n) {
        if (c >= p) {
          dest = c - p; b++;
          if (dest >= p) {
            dest -= p; b++; }
        }
        else
          dest = c;
        V::add_mod_core (dest, carry, p, z); carry = b + z;
        return;
      }
      if (c >= p) {
        c -= p; b++;
        if (c >= p) {
          c -= p; b++;
          if (c >= p) {
            c -= p; b++; }
        }
      }
      dest = c;
      V::add_mod_core (dest, carry, p, z); carry = b + z;
    }
  };

/*
  Modular product when int of double size are not available
*/

  template <typename C, nat m>
  struct mul_mod_helper<C, void, m> {
    // C is unsigned.
    static const nat n2 = 4 * sizeof(C);
    static const nat n  = 8 * sizeof(C);
    typedef long_int_mul_op<C>    long_mul;
    typedef long_int_rshift_op<C> long_rshift;
    typedef long_int_sub_op<C>    long_sub;
    typedef long_int_ge_op<C>     long_ge;

    static void
    op (C& dest, const C& src, const C& p,
        const C& q, nat r, nat s, nat t) {

      C a1, a0, hi_a1, hi_a0;
      C tmp = dest;
      
      long_mul::op (a1, a0, tmp, src);
      hi_a1 = a1;
      hi_a0 = a0;
      long_rshift::op (hi_a1, hi_a0, s);
      C b1, b0;
      long_mul::op (b1, b0, hi_a0, q);
      if (m >= n) {if (hi_a1 != 0) b1 += q;}
      if (m+1 <= n) {
        if (t >= n)
          if (n >= t) b0 = b1 << (n-t); else  b0 = b1 >> (t-n);
        else
          b0 = (b0 >> t) | (b1 << (n-t));
      }
      else
        long_rshift::op (b1, b0, t);
      if (m+2 <= n) {
        dest = a0 - b0*p;
        dest = min (dest, dest - p);
        return;
      }
      C d1, d0;
      long_mul::op (d1, d0, b0, p);
      long_sub::op (a1, a0, d1, d0);
      if (m+1 <= n) {
        if (long_ge::op (a1, a0, 0, p)) {
          dest = a0 - p;
          if (dest >= p)
            dest -=p; }
        else dest = a0;
        return;
      }
      if (long_ge::op (a1, a0, 0, p)) {
        long_sub::op (a1, a0, 0, p);
        if (long_ge::op (a1, a0, 0, p)) {
          long_sub::op (a1, a0, 0, p);
          if (long_ge::op (a1, a0, 0, p)) 
            a0 -= p;
        }
      }
      dest = a0;
    }
 
    static void
    op (C& dest, const C& src, const C& p, C& carry,
        const C& q, nat r, nat s, nat t) {

      C a1, a0, hi_a1, hi_a0;
      C tmp = dest;
      
      long_mul::op (a1, a0, tmp, src);
      hi_a1 = a1;
      hi_a0 = a0;
      long_rshift::op (hi_a1, hi_a0, s);
      C b1, b0, c = 0;
      long_mul::op (b1, b0, hi_a0, q);
      if (m >= n) {if (hi_a1 != 0) b1 += q;}
      if (m+1 <= n) {
        if (t >= n)
          if (n >= t) b0 = b1 << (n-t); else  b0 = b1 >> (t-n);
        else
          b0 = (b0 >> t) | (b1 << (n-t));
      }
      else
        long_rshift::op (b1, b0, t);
      if (m+2 <= n) {
        dest = a0 - b0 * p;
        if (dest >= p) { dest -= p; b0++; }
        V::add_mod_core (dest, carry, p, c); carry = b0 + c;
        return;
      }
      C d1, d0;
      long_mul::op (d1, d0, b0, p);
      long_sub::op (a1, a0, d1, d0);
      if (m+1 <= n) {
        if (long_ge::op (a1, a0, 0, p)) {
          dest = a0 - p; b0++;
          if (dest >= p) {
            dest -=p; b0++; } }
        else dest = a0;
        V::add_mod_core (dest, carry, p, c); carry = b0 + c;
        return;
      }
      if (long_ge::op (a1, a0, 0, p)) {
        long_sub::op (a1, a0, 0, p); b0++;
        if (long_ge::op (a1, a0, 0, p)) {
          long_sub::op (a1, a0, 0, p); b0++;
          if (long_ge::op (a1, a0, 0, p)) {
            a0 -= p; b0++; }
        }
      }
      dest = a0;
      V::add_mod_core (dest, carry, p, c); carry = b0 + c;
    }
  };
  
  template<typename C, typename M> static inline void
  mul_mod (C& dest, const C& src, const M& x) {
    static const nat m = V::template maximum_size_helper<C>::value;
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
    uC tmp = dest;
    mul_mod_helper<uC,uD,m>::op (tmp, (uC) src, (uC) x.p,
                                 (uC) x.q, x.r, x.s, x.t);
    dest = tmp; }

  template<typename C, typename M> static inline void
  mul_mod (C& dest, const C& src, const M& x, C& carry) {
    static const nat m = V::template maximum_size_helper<C>::value;
    typedef typename unsigned_of_helper<C>::type uC;
    typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
    uC tmp = dest, ucarry = carry;
    mul_mod_helper<uC,uD,m>::op (tmp, (uC) src, (uC) x.p, ucarry,
                                 (uC) x.q, x.r, x.s, x.t);
    dest = tmp; carry = ucarry; }

  template<typename C, typename M> static inline void
  mul_mod (C& dest, const C& s1, const C& s2, const M& x) {
    dest = s1;
    mul_mod (dest, s2, x); }

  template<typename C, typename M> static inline void
  mul_mod (C& dest, const C& s1, const C& s2, const M& x, C& carry) {
    dest = s1;
    mul_mod (dest, s2, x, carry); }
};

template<nat m>
struct modulus_int_preinverse:
  public modulus_encoding_int_naive<
         modulus_div_int_naive<
         modulus_inv_int_naive<
         modulus_mul_int_preinverse<
         modulus_add_int_naive<
         modulus_reduction_int_naive<
         modulus_normalization_int_naive<
         modulus_maximum_size_int<m> > > > > > > > {};

/******************************************************************************
* Specialized modulus
******************************************************************************/

template<typename C, nat m>
class modulus <C, modulus_int_preinverse<m> > {
MMX_ALLOCATORS
public:
  C p;
  C q;
  nat r;
  nat s;
  nat t;
  C h; // p >> 1
  typedef C base;
  typedef modulus_int_preinverse<m> variant;
  typedef modulus <C, modulus_int_preinverse<m> > self_type;

  inline C operator * () const { return p; }

  inline modulus () : p (0) {
    r = variant::dyn_r (p);
    s = variant::dyn_s (p, r);
    t = variant::dyn_t (p, r);
    q = variant::dyn_q (p, r, s, t);
    h = p >> 1;
  }

  template<typename X>
  inline modulus (const X& x) {
    p = (C) x;
    ASSERT (variant::normalize (p), "modulus: bad value");
    r = variant::dyn_r (p);
    s = variant::dyn_s (p, r);
    t = variant::dyn_t (p, r);
    q = variant::dyn_q (p, r, s, t);
    h = p >> 1;
  }

  template<typename X, X _p>
  inline modulus (const fixed_value<X, _p>& x) {
    typedef typename variant::template auxiliaries_helper<C, _p> A;
    C loc_p = _p;
    ASSERT (variant::normalize (loc_p), "modulus: bad value");
    p = A::up;
    q = A::q ();
    r = A::r;
    s = A::s;
    t = A::t;
    h = p >> 1; }
    
  inline modulus (const self_type& x) : p (x.p), q (x.q),
                                        r (x.r), s (x.s),
                                        t (x.t), h (x.h) {}

  inline self_type&
  operator = (const self_type& x) {
    p = x.p;
    q = x.q;
    r = x.r;
    s = x.s;
    t = x.t;
    h = x.h;
    return *this; }

  inline bool operator == (const self_type& a) const {
    return p == a.p; }

  inline bool operator != (const self_type& a) const {
    return p != a.p; }
};

#undef TMPL

/******************************************************************************
* Helpers specializations for hardware integers
******************************************************************************/

// Default variant 
#define Variant modulus_int_preinverse

template<>
struct modulus_variant_helper<unsigned char> {
  typedef Variant<8*sizeof(unsigned char)> MV; };

template<>
struct modulus_variant_helper<short unsigned int> {
  typedef Variant<8*sizeof(short unsigned int)> MV; };

template<>
struct modulus_variant_helper<unsigned int> {
  typedef Variant<8*sizeof(unsigned int)> MV; };

template<>
struct modulus_variant_helper<long unsigned int> {
  typedef Variant<8*sizeof(long unsigned int)> MV; };

template<>
struct modulus_variant_helper<long long unsigned int> {
  typedef Variant<8*sizeof(long long unsigned int)> MV; };

template<>
struct modulus_variant_helper<char> {
  typedef Variant<8*sizeof(char)-1> MV; };

template<>
struct modulus_variant_helper<signed char> {
  typedef Variant<8*sizeof(char)-1> MV; };

template<>
struct modulus_variant_helper<short int> {
  typedef Variant<8*sizeof(short int)-1> MV; };

template<>
struct modulus_variant_helper<int> {
  typedef Variant<8*sizeof(int)-1> MV; };

template<>
struct modulus_variant_helper<long int> {
  typedef Variant<8*sizeof(long int)-1> MV; };

template<>
struct modulus_variant_helper<long long int> {
  typedef Variant<8*sizeof(long long int)-1> MV; };

#undef Variant

/******************************************************************************
* Lift and project 
******************************************************************************/

#define TMPLVW template<typename V, typename W>
#define Modular(I) modular<modulus<I,V>,W>
#define DECLARE_HELPER(I) \
STYPE_TO_TYPE(TMPLVW,lift_type,Modular(I),integer);     \
TMPLVW inline integer \
lift (const Modular(I)& x) { return integer (*x); } \
TMPLVW \
struct project_helper<Modular(I)> { \
  template<typename C> static inline Modular(I) \
  op (const C& x) { return Modular(I) (x); } \
  static inline Modular(I) \
  op (const integer& x) { \
    integer p= *Modular(I)::get_modulus (); \
    return (x >= 0) ? Modular(I) (as<I> (rem (x, p))) \
      : - Modular(I) (as<I> (rem (-x, p))); } \
};

DECLARE_HELPER(unsigned char)
DECLARE_HELPER(signed char)
DECLARE_HELPER(unsigned short int)
DECLARE_HELPER(signed short int)
DECLARE_HELPER(unsigned int)
DECLARE_HELPER(int)
DECLARE_HELPER(unsigned long int)
DECLARE_HELPER(long int)
DECLARE_HELPER(unsigned long long int)
DECLARE_HELPER(long long int)
#undef DECLARE_HELPER
#undef Modular
#undef TMPLVW

} // namespace mmx
#endif //__MMX__MODULAR_INT__HPP