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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : modular.hpp
* DESCRIPTION: modular arithmetic
* 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__HPP
#define __MMX__MODULAR__HPP
#include <basix/generic.hpp>
#include <basix/wrap.hpp>
#include <basix/operators.hpp>
#include <basix/int.hpp>
#include <numerix/integer.hpp>
#include <numerix/modulus.hpp>

namespace mmx {
#define Modular_variant(M) modular_variant_helper<M>::MV
#define TMPL_DEF \
  template<typename M, typename V= typename Modular_variant(M)>
#define TMPL template<typename M, typename V>
#define Modulus modulus<C,V>
#define Modular modular<M,V>

template<typename M>
struct modular_variant_helper;
TMPL class modular;

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

// One modulus is attached to each modular
struct modular_local {
  template<typename M>
  class modulus_storage {
    static inline M& dyn_modulus () {
      static M modulus = M ();
      return modulus; }
  public:
    static inline void set_modulus (const M& p) {
      if (! is_exact_zero (*p)) dyn_modulus () = p; }
    static inline M get_modulus () { return dyn_modulus (); }
  };
};

// Global storage of the modulus
struct modular_global {
  template<typename M>
  class modulus_storage {
    static inline M& dyn_modulus () {
      static M modulus = M ();
      return modulus; }
  public:
    static inline void set_modulus (const M& p) { dyn_modulus () = p; }
    static inline M get_modulus () { return dyn_modulus (); }
  };
};

// Fixed modulus
template<typename C, C p>
struct modular_fixed {
  template<typename M>
  class modulus_storage {
    static inline M& dyn_modulus () {
      static fixed_value <C, p> _p;
      static M modulus = M (_p);
      return modulus; }
  public:
    static inline void set_modulus (const M& q) {
      ASSERT (q == dyn_modulus (),
              "modular_fixed: can not change the modulus");
    }
    static inline M get_modulus () { return dyn_modulus (); }
  };
};

/******************************************************************************
* Variant helper, default storage of the modulus is "global"
******************************************************************************/

template<typename M>
struct modular_variant_helper {
  typedef modular_global MV; };

/******************************************************************************
* General modulars
******************************************************************************/

TMPL_DEF
class modular {
MMX_ALLOCATORS
public:
  typedef M modulus;
  typedef typename M::base C;
  typedef typename V::template modulus_storage<M> S;

  C rep;
  static inline M get_modulus () { return S::get_modulus (); }
  static inline void set_modulus (const M& p) { S::set_modulus (p); }

  inline C operator * () const {
    C dest;
    decode_mod (dest, rep, get_modulus ());
    return dest; }

  inline modular () { encode_mod (rep, C (), get_modulus ()); }

  inline modular (const Modular& s) : rep (s.rep) {}

  inline modular (const C& a, bool reduced) {
    if (reduced) rep = a;
    else encode_mod (rep, a, get_modulus ()); }

  inline modular (const C& s, const M& p, bool reduced=false) {
    VERIFY (p == get_modulus (), "wrong modulus");
    if (reduced) rep= s; 
    else encode_mod (rep, s, p); }

  template<typename O>
  inline modular (const O& a) {
    encode_mod (rep, as<C> (a), get_modulus ()); }

  inline Modular& operator = (const Modular& a) {
    rep = a.rep;
    return *this; }

  inline bool operator == (const Modular& a) const {
    return rep == a.rep; }

  inline bool operator != (const Modular& a) const {
    return rep != a.rep; }
};

STYPE_TO_TYPE(TMPL,lift_type,Modular,typename Modular::C);

TMPL inline typename Modular::C
lift (const Modular& x) { return *x; }

TMPL
struct project_helper<Modular > {
  template<typename C> static inline Modular
  op (const C& x) { return Modular (x); }
};

TMPL inline M get_modulus (const Modular& x) {
  (void) x; return Modular::get_modulus (); }
TMPL inline void set_modulus (Modular& x, const M& p) {
  (void) x; Modular::set_modulus (p); }

template<typename TM, typename TV, typename SM, typename SV>
struct as_helper<modular<TM,TV>,modular<SM,SV> > {
  static inline modular<TM,TV>
  cv (const modular<SM,SV>& p) { return modular<TM,TV> (*p); }
};

template<typename C, typename V, typename W>
struct as_helper<C,modular<modulus<C,V>,W> > {
  static inline C
  cv (const modular<modulus<C,V>,W>& p) { return *p; }
};

template<typename C, typename V, typename W>
struct as_helper<generic, modular<modulus<C,V>,W> > {
  static inline generic cv (const modular<modulus<C,V>,W> & x) {
    return new generic_concrete_rep<modular<modulus<C,V>,W> > (x); }
};

WRAP_WRAPPED_IMPL(TMPL inline,Modular)

/******************************************************************************
* Arithmetic
******************************************************************************/

TMPL inline Modular&
operator += (Modular& dest, const Modular& s) {
  add_mod (dest.rep, s.rep, Modular::get_modulus ());
  return dest; }

TMPL inline Modular&
operator -= (Modular& dest, const Modular& s) {
  sub_mod (dest.rep, s.rep, Modular::get_modulus ());
  return dest; }

TMPL inline Modular&
operator *= (Modular& dest, const Modular& s) {
  mul_mod (dest.rep, s.rep, Modular::get_modulus ());
  return dest; }

TMPL inline Modular&
operator /= (Modular& dest, const Modular& s) {
  div_mod (dest.rep, s.rep, Modular::get_modulus ());
  return dest; }

TMPL inline Modular
operator + (const Modular& a, const Modular& b) {
  Modular c;
  add_mod (c.rep, a.rep, b.rep, Modular::get_modulus ());
  return c; }

TMPL inline Modular
operator + (const Modular& a, const typename M::base& b) {
  return a + Modular (b); }

TMPL inline Modular
operator + (const typename M::base& a, const Modular& b) {
  return Modular (a) + b; }

TMPL inline Modular
operator - (const Modular& a) {
  Modular c;
  neg_mod (c.rep, a.rep, Modular::get_modulus ());
  return c; }

TMPL inline Modular
operator - (const Modular& a, const Modular& b) {
  Modular c;
  sub_mod (c.rep, a.rep, b.rep, Modular::get_modulus ());
  return c; }

TMPL inline Modular
operator - (const Modular& a, const typename M::base& b) {
  return a - Modular (b); }

TMPL inline Modular
operator - (const typename M::base& a, const Modular& b) {
  return Modular (a) - b; }

TMPL inline Modular
operator * (const Modular& a, const Modular& b) {
  Modular c;
  mul_mod (c.rep, a.rep, b.rep, Modular::get_modulus ());
  return c; }

TMPL inline Modular
operator * (const Modular& a, const typename M::base& b) {
  return a * Modular (b); }

TMPL inline Modular
operator * (const typename M::base& a, const Modular& b) {
  return Modular (a) * b; }

TMPL inline Modular
operator / (const Modular& a, const Modular& b) {
  Modular c;
  div_mod (c.rep, a.rep, b.rep, Modular::get_modulus ());
  return c; }

TMPL inline Modular
operator / (const Modular& a, const typename M::base& b) {
  return a / Modular (b); }

TMPL inline Modular
operator / (const typename M::base& a, const Modular& b) {
  return Modular (a) / b; }

TMPL inline Modular
invert (const Modular& a) {
  Modular c;
  inv_mod (c.rep, a.rep, Modular::get_modulus ());
  return c; }

/******************************************************************************
* Gcd & lcm routine
******************************************************************************/

TMPL inline Modular
gcd (const Modular& a, const Modular& b) {
  if (a == 0) return b;
  return a; }

TMPL inline Modular
lcm (const Modular& a, const Modular& b) {
  return a * b; }

/******************************************************************************
* Output routine
******************************************************************************/

WRAP_PRINT_IMPL(TMPL inline,Modular)

/******************************************************************************
* Syntactic sugar
******************************************************************************/

ARITH_INT_SUGAR(TMPL,Modular)
POOR_MAN_SQRT_SUGAR(TMPL,Modular)
POOR_MAN_ELEMENTARY_SUGAR(TMPL,Modular)

#undef TMPL_DEF
#undef TMPL
#undef Modulus
#undef Modular

/******************************************************************************
* Modulars with local moduli
******************************************************************************/

#define TMPL template<typename M>
#define Modulus modulus<C,V>
#define Modular modular<M,modular_local>

TMPL
class modular<M, modular_local> {
MMX_ALLOCATORS
public:
  typedef M modulus;
  typedef modular_local V;
  typedef typename M::base C;
  typedef typename V::template modulus_storage<M> S;

  C rep;
  M mod;
  static inline M get_modulus () { return S::get_modulus (); }
  static inline void set_modulus (const M& p) { S::set_modulus (p); }

  inline M get_local_modulus () const { return mod; }
  inline void set_local_modulus (const M& p) { mod= p; S::set_modulus (p); }

  inline C operator * () const {
    C dest;
    decode_mod (dest, rep, mod);
    return dest; }

  inline modular () {
    mod= get_modulus ();
    encode_mod (rep, C (), mod); }

  inline modular (const Modular& s) : rep (s.rep), mod (s.mod) {}

  inline modular (const C& s, bool reduced) {
    mod= get_modulus ();
    if (reduced) rep= s;
    else encode_mod (rep, s, mod); }

  inline modular (const C& s, const M& p, bool reduced=false) {
    mod= p;
    set_modulus (p);
    if (reduced) rep= s; 
    else encode_mod (rep, s, mod); }

  template<typename O>
  inline modular (const O& a) {
    mod= get_modulus ();
    encode_mod (rep, as<C> (a), mod); }

  inline Modular& operator = (const Modular& a) {
    rep = a.rep;
    mod = a.mod;
    return *this; }

  inline bool operator == (const Modular& a) const {
    if (mod   == a.mod) return rep == a.rep;
    if (mod   == C(0))  return Modular (* (*this), a.mod) == a;
    if (a.mod == C(0))  return Modular (* a, mod) == *this;
    return false; }

  inline bool operator != (const Modular& a) const {
    if (mod   == a.mod) return rep != a.rep;
    if (mod   == C(0))  return Modular (* (*this), a.mod) != a;
    if (a.mod == C(0))  return Modular (* a, mod) != *this;
    return true; }
};

template<typename M> inline M
get_modulus (const modular<M, modular_local>& x) {
  return x.get_local_modulus (); }

template<typename M> inline void
set_modulus (Modular& x, const M& p) {
  x.set_local_modulus (p); }

template<typename TM, typename SM, typename SV>
struct as_helper<modular<TM,modular_local>,modular<SM,SV> > {
  static inline modular<TM,modular_local>
  cv (const modular<SM,SV>& p) {
    return modular<TM,modular_local> (* p, get_modulus (p)); }
};

TMPL nat 
hash (const Modular& c) {
  nat h= hash (*c);
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ hash (get_modulus (c)); }

TMPL nat
exact_hash (const Modular& c) {
  nat h= exact_hash (*c);
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ exact_hash (get_modulus (c)); }

TMPL nat hard_hash (const Modular& c) {
  nat h= hard_hash (*c);
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ hard_hash (get_modulus (c)); }

TMPL bool exact_eq (const Modular& c1, const Modular& c2) {
  return exact_eq (*c1, *c2) && 
         exact_eq (get_modulus (c1), get_modulus (c2)); }

TMPL bool exact_neq (const Modular& c1, const Modular& c2) {
  return ! exact_eq (c1, c2); }

TMPL bool hard_eq (const Modular& c1, const Modular& c2) {
  return hard_eq (*c1, *c2) && 
         hard_eq (get_modulus (c1), get_modulus (c2)); }

TMPL bool hard_neq (const Modular& c1, const Modular& c2) {
  return ! hard_eq (c1, c2); }

TMPL inline Modular&
operator += (Modular& dest, const Modular& s) {
  M md= get_modulus (dest), ms= get_modulus (s);
  M mod= md == 0 ? ms : md;
  VERIFY (md == 0 || ms == 0 || md == ms, "incompatible moduli");
  add_mod (dest.rep, s.rep, mod);
  return dest; }

TMPL inline Modular&
operator -= (Modular& dest, const Modular& s) {
  M md= get_modulus (dest), ms= get_modulus (s);
  M mod= md == 0 ? ms : md;
  VERIFY (md == 0 || ms == 0 || md == ms, "incompatible moduli");
  sub_mod (dest.rep, s.rep, mod);
  return dest; }

TMPL inline Modular&
operator *= (Modular& dest, const Modular& s) {
  M md= get_modulus (dest), ms= get_modulus (s);
  M mod= md == 0 ? ms : md;
  VERIFY (md == 0 || ms == 0 || md == ms, "incompatible moduli");
  mul_mod (dest.rep, s.rep, mod);
  return dest; }

TMPL inline Modular&
operator /= (Modular& dest, const Modular& s) {
  M md= get_modulus (dest), ms= get_modulus (s);
  M mod= md == 0 ? ms : md;
  VERIFY (md == 0 || ms == 0 || md == ms, "incompatible moduli");
  div_mod (dest.rep, s.rep, mod);
  return dest; }

TMPL inline Modular
operator + (const Modular& a, const Modular& b) {
  M ma= get_modulus (a), mb= get_modulus (b);
  M mod= ma == 0 ? mb : ma;
  VERIFY (ma == 0 || mb == 0 || ma == mb, "incompatible moduli");
  Modular c;
  add_mod (c.rep, a.rep, b.rep, mod);
  c.mod= mod;
  return c; }

TMPL inline Modular
operator - (const Modular& a) {
  M mod= get_modulus (a);
  Modular c;
  neg_mod (c.rep, a.rep, mod);
  c.mod= mod;
  return c; }

TMPL inline Modular
operator - (const Modular& a, const Modular& b) {
  M ma= get_modulus (a), mb= get_modulus (b);
  M mod= ma == 0 ? mb : ma;
  VERIFY (ma == 0 || mb == 0 || ma == mb, "incompatible moduli");
  Modular c;
  sub_mod (c.rep, a.rep, b.rep, mod);
  c.mod= mod;
  return c; }

TMPL inline Modular
operator * (const Modular& a, const Modular& b) {
  M ma= get_modulus (a), mb= get_modulus (b);
  M mod= ma == 0 ? mb : ma;
  VERIFY (ma == 0 || mb == 0 || ma == mb, "incompatible moduli");
  Modular c;
  mul_mod (c.rep, a.rep, b.rep, mod);
  c.mod= mod;
  return c; }

TMPL inline Modular
operator / (const Modular& a, const Modular& b) {
  M ma= get_modulus (a), mb= get_modulus (b);
  M mod= ma == 0 ? mb : ma;
  VERIFY (ma == 0 || mb == 0 || ma == mb, "incompatible moduli");
  Modular c;
  div_mod (c.rep, a.rep, b.rep, mod);
  c.mod= mod;
  return c; }

TMPL inline Modular
invert (const Modular& a) {
  M mod= get_modulus (a);
  Modular c;
  inv_mod (c.rep, a.rep, mod);
  c.mod= mod;
  return c; }

/******************************************************************************
* Extras for glue
******************************************************************************/

#define mmx_modular(R) modular<modulus<R >,modular_local>

#undef TMPL
#undef Modulus
#undef Modular
} // namespace mmx
#endif //__MMX__MODULAR__HPP