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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : multiplier.hpp
* DESCRIPTION: Speedups for when several products are done with a same value
* 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__MULTIPLIER__HPP
#define __MMX__MULTIPLIER__HPP
#include <basix/int.hpp>
#include <basix/wrap.hpp>
#include <numerix/modular.hpp>

namespace mmx {
#define TMPL template<typename C>
#define TMPLX template<typename C, typename X>
#define Multiplier multiplier<C>

TMPL
class multiplier {
  C rep;
public:
  inline C operator * () const { return rep; }
  inline multiplier () : rep (0) {}
  inline multiplier (const C& s) : rep (s) {}
  inline multiplier (const Multiplier& s) : rep (s.rep) {}
  inline Multiplier& operator = (const Multiplier& s) {
    rep = s.rep; return *this; }
  inline bool operator == (const Multiplier& a) const {
    return rep == a.rep; }
  inline bool operator != (const Multiplier& a) const {
    return rep != a.rep; }
  template<typename X> static inline void 
  rmul (X& dest, const Multiplier& s) { dest = dest * s.rep; }
  template<typename X> static inline void 
  lmul (X& dest, const Multiplier& s) { dest = s.rep * dest; }
};

template<typename TC, typename SC>
struct as_helper <multiplier<TC>,multiplier<SC> > {
  static inline multiplier<TC>
  cv (const multiplier<SC>& p) { return multiplier<TC> (*p); }
};

WRAP_WRAPPED_IMPL(TMPL inline,Multiplier)
WRAP_PRINT_IMPL(TMPL inline,Multiplier)

TMPLX inline X
operator * (const X& a, const Multiplier& b) {
  X c (a);
  Multiplier::rmul (c, b);
  return c; }

TMPLX inline X
operator * (const Multiplier& b, const X& a) {
  X c (a);
  Multiplier::lmul (c, b);
  return c; }

TMPLX inline X
operator *= (X& a, const Multiplier& b) {
  Multiplier::rmul (a, b);
  return a; }

#undef Multiplier
#undef TMPLX
#undef TMPL

/******************************************************************************
* Multipliers for modulars over hardware int with preinverse
******************************************************************************/

template<nat size>
struct modulus_multiplier_int_preinverse_helper {
  typedef modulus_maximum_size_int<size> size_helper;

/*
  Double size int is available
*/

  template <typename C, typename D, nat m>
  struct multiplier_helper {
    // C and D are unsigned.
    static const nat n = 8 * sizeof(C);

    static inline void
    set (C& dest, const C& src, const C& p,
         const C& q, nat r, nat s, nat t) {
      (void) t; (void) q; (void) s; (void) t;
      if (p == 0)
        dest = src;
      else
        dest = (((D) src) << r) / p;
      return;
    }

    static inline void
    mul (C& dest, const C& src, const C& p, nat r, const C& presrc) {
      C b = (C) ((((D) dest) * presrc) >> r);
      if (m + 1 <= n) {
        dest = dest * src - b * p;
        if (dest >= p) dest -= p;
        return;
      }
      D d = ((D) dest) * src - ((D) b) * p;
      if (d >= p) d -= p;
      dest = (C) d;
    }
  };

/*
  Multiplier when double size int is not available
*/

  template <typename C, nat m>
  struct multiplier_helper<C, void, m> {
    // C is unsigned.
    static const nat n = 8 * sizeof(C);
    static const nat n2 = 4 * sizeof(C);
    typedef long_int_mul_op<C>    long_mul;
    typedef long_int_lshift_op<C> long_lshift;
    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 inline void
    set (C& dest, const C& src, const C& p,
         const C& q, nat r, nat s, nat t) {
      if (r == 0) dest = 0; // means p = 1
      C d0, d1;
      long_mul::op (d1, d0, src, q);
      if (m+1 <= n) {
        long_rshift::op (d1, d0, t+s-r);
        dest = d0;
        C d = (src << r) - dest * p;
        if (d >= p) dest++;
        return;
      }
      long_rshift::op (d1, d0, r-1);
      dest = d0;
      C e0, e1;
      long_mul::op (e1, e0, dest, p);
      d0 = src;
      d1 = 0;
      long_lshift::op (d1, d0, r);
      long_sub::op (d1, d0, e1, e0);
      if (long_ge::op (d1, d0, 0, p)) {
        long_sub::op (d1, d0, 0, p);
        dest ++;
        if (long_ge::op (d1, d0, 0, p)) {
          long_sub::op (d1, d0, 0, p);
          dest ++;
        }
      }
    }

    static inline void
    mul (C& dest, const C& src, const C& p, nat r, const C& presrc) {
      C b0, b1;
      long_mul::op (b1, b0, dest, presrc);
      long_rshift::op (b1, b0, r);
      if (m + 1 <= n) {
        dest = dest * src - b0 * p;
        if (dest >= p) dest -= p;
        return;
      }
      C x0, x1, y0, y1;
      long_mul::op (y1, y0, dest, src);
      long_mul::op (x1, x0, b0, p);
      long_sub::op (y1, y0, x1, x0);
      if (long_ge::op (y1, y0, 0, p)) dest = y0-p;
      else dest = y0;
    }
  };
  
  template<typename C, typename M> static inline void 
  set (C& dest, const C& src, const M& x) {
    static const nat m = size_helper::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;
    multiplier_helper<uC,uD,m>::set (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 (C& dest, const C& src, const M& x, const C& presrc) {
    static const nat m = size_helper::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;
    multiplier_helper<uC,uD,m>::mul (tmp, (uC) src,
                                     (uC) x.p, x.r, (uC) presrc);
    dest = tmp; }
};

#define Modulus modulus<C, modulus_int_preinverse<size> >
#define Modular modular<Modulus, V>

template<typename C, typename V, nat size>
class multiplier<Modular> {
  typedef modulus_multiplier_int_preinverse_helper<size> multiplier_helper;
  C rep;
  C pre; // pre-computed shared value
public:
  inline C operator * () const { return rep; }
  inline multiplier<Modular> () : rep (0) {
    multiplier_helper::set (pre, rep, Modular::get_modulus ()); }
  inline multiplier<Modular> (const Modular& s) : rep (*s) {
    multiplier_helper::set (pre, rep, Modular::get_modulus ()); }
  inline multiplier<Modular> (const multiplier<Modular>& s)
    : rep (s.rep), pre (s.pre) {}
  inline multiplier<Modular>& operator = (const multiplier<Modular>& s) {
    rep = s.rep; pre = s.pre; return *this; }
  inline bool
  operator == (const multiplier<Modular>& a) const {
    VERIFY (pre == a.pre,"multiplier<Modular>: pre-computed value differ");
    return rep == a.rep; }
  inline bool 
  operator != (const multiplier<Modular>& a) const {
    VERIFY (pre != a.pre,"multiplier<Modular>: pre-computed value are equal");
    return rep != a.rep; }
  static inline void 
  rmul (Modular& dest, const multiplier<Modular>& s) {
    multiplier_helper::mul (dest.rep, s.rep,
                            Modular::get_modulus (), s.pre); }
  static inline void 
  lmul (Modular& dest, const multiplier<Modular>& s) {
    multiplier_helper::mul (dest.rep, s.rep,
                            Modular::get_modulus (), s.pre); }
};

#undef Modulus
#undef Modular
} // namespace mmx
#endif //__MMX__MULTIPLIER__HPP