00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__MULTIPLIER__HPP
00014 #define __MMX__MULTIPLIER__HPP
00015 #include <basix/int.hpp>
00016 #include <basix/wrap.hpp>
00017 #include <numerix/modular.hpp>
00018
00019 namespace mmx {
00020 #define TMPL template<typename C>
00021 #define TMPLX template<typename C, typename X>
00022 #define Multiplier multiplier<C>
00023
00025 TMPL
00026 class multiplier {
00027 C rep;
00028 public:
00029 inline C operator * () const { return rep; }
00030 inline multiplier () : rep (0) {}
00031 inline multiplier (const C& s) : rep (s) {}
00032 inline multiplier (const Multiplier& s) : rep (s.rep) {}
00033 inline Multiplier& operator = (const Multiplier& s) {
00034 rep = s.rep; return *this; }
00035 inline bool operator == (const Multiplier& a) const {
00036 return rep == a.rep; }
00037 inline bool operator != (const Multiplier& a) const {
00038 return rep != a.rep; }
00039 template<typename X> static inline void
00040 rmul (X& dest, const Multiplier& s) { dest = dest * s.rep; }
00041 template<typename X> static inline void
00042 lmul (X& dest, const Multiplier& s) { dest = s.rep * dest; }
00043 };
00044
00045 template<typename TC, typename SC>
00046 struct as_helper <multiplier<TC>,multiplier<SC> > {
00047 static inline multiplier<TC>
00048 cv (const multiplier<SC>& p) { return multiplier<TC> (*p); }
00049 };
00050
00051 WRAP_WRAPPED_IMPL(TMPL inline,Multiplier)
00052 WRAP_PRINT_IMPL(TMPL inline,Multiplier)
00053
00054 TMPLX inline X
00055 operator * (const X& a, const Multiplier& b) {
00056 X c (a);
00057 Multiplier::rmul (c, b);
00058 return c; }
00059
00060 TMPLX inline X
00061 operator * (const Multiplier& b, const X& a) {
00062 X c (a);
00063 Multiplier::lmul (c, b);
00064 return c; }
00065
00066 TMPLX inline X
00067 operator *= (X& a, const Multiplier& b) {
00068 Multiplier::rmul (a, b);
00069 return a; }
00070
00071 #undef Multiplier
00072 #undef TMPLX
00073 #undef TMPL
00074
00075
00076
00077
00078
00079 template<nat size>
00080 struct modulus_multiplier_int_preinverse_helper {
00081 typedef modulus_maximum_size_int<size> size_helper;
00082
00083
00084
00085
00086
00087 template <typename C, typename D, nat m>
00088 struct multiplier_helper {
00089
00090 static const nat n = 8 * sizeof(C);
00091
00092 static inline void
00093 set (C& dest, const C& src, const C& p,
00094 const C& q, nat r, nat s, nat t) {
00095 (void) t; (void) q; (void) s; (void) t;
00096 if (p == 0)
00097 dest = src;
00098 else
00099 dest = (((D) src) << r) / p;
00100 return;
00101 }
00102
00103 static inline void
00104 mul (C& dest, const C& src, const C& p, nat r, const C& presrc) {
00105 C b = (C) ((((D) dest) * presrc) >> r);
00106 if (m + 1 <= n) {
00107 dest = dest * src - b * p;
00108 if (dest >= p) dest -= p;
00109 return;
00110 }
00111 D d = ((D) dest) * src - ((D) b) * p;
00112 if (d >= p) d -= p;
00113 dest = (C) d;
00114 }
00115 };
00116
00117
00118
00119
00120
00121 template <typename C, nat m>
00122 struct multiplier_helper<C, void, m> {
00123
00124 static const nat n = 8 * sizeof(C);
00125 static const nat n2 = 4 * sizeof(C);
00126 typedef long_int_mul_op<C> long_mul;
00127 typedef long_int_lshift_op<C> long_lshift;
00128 typedef long_int_rshift_op<C> long_rshift;
00129 typedef long_int_sub_op<C> long_sub;
00130 typedef long_int_ge_op<C> long_ge;
00131
00132 static inline void
00133 set (C& dest, const C& src, const C& p,
00134 const C& q, nat r, nat s, nat t) {
00135 if (r == 0) dest = 0;
00136 C d0, d1;
00137 long_mul::op (d1, d0, src, q);
00138 if (m+1 <= n) {
00139 long_rshift::op (d1, d0, t+s-r);
00140 dest = d0;
00141 C d = (src << r) - dest * p;
00142 if (d >= p) dest++;
00143 return;
00144 }
00145 long_rshift::op (d1, d0, r-1);
00146 dest = d0;
00147 C e0, e1;
00148 long_mul::op (e1, e0, dest, p);
00149 d0 = src;
00150 d1 = 0;
00151 long_lshift::op (d1, d0, r);
00152 long_sub::op (d1, d0, e1, e0);
00153 if (long_ge::op (d1, d0, 0, p)) {
00154 long_sub::op (d1, d0, 0, p);
00155 dest ++;
00156 if (long_ge::op (d1, d0, 0, p)) {
00157 long_sub::op (d1, d0, 0, p);
00158 dest ++;
00159 }
00160 }
00161 }
00162
00163 static inline void
00164 mul (C& dest, const C& src, const C& p, nat r, const C& presrc) {
00165 C b0, b1;
00166 long_mul::op (b1, b0, dest, presrc);
00167 long_rshift::op (b1, b0, r);
00168 if (m + 1 <= n) {
00169 dest = dest * src - b0 * p;
00170 if (dest >= p) dest -= p;
00171 return;
00172 }
00173 C x0, x1, y0, y1;
00174 long_mul::op (y1, y0, dest, src);
00175 long_mul::op (x1, x0, b0, p);
00176 long_sub::op (y1, y0, x1, x0);
00177 if (long_ge::op (y1, y0, 0, p)) dest = y0-p;
00178 else dest = y0;
00179 }
00180 };
00181
00182 template<typename C, typename M> static inline void
00183 set (C& dest, const C& src, const M& x) {
00184 static const nat m = size_helper::template maximum_size_helper<C>::value;
00185 typedef typename unsigned_of_helper<C>::type uC;
00186 typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
00187 uC tmp = dest;
00188 multiplier_helper<uC,uD,m>::set (tmp, (uC) src,
00189 (uC) x.p, (uC) x.q,
00190 x.r, x.s, x.t);
00191 dest = tmp; }
00192
00193 template<typename C, typename M> static inline void
00194 mul (C& dest, const C& src, const M& x, const C& presrc) {
00195 static const nat m = size_helper::template maximum_size_helper<C>::value;
00196 typedef typename unsigned_of_helper<C>::type uC;
00197 typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
00198 uC tmp = dest;
00199 multiplier_helper<uC,uD,m>::mul (tmp, (uC) src,
00200 (uC) x.p, x.r, (uC) presrc);
00201 dest = tmp; }
00202 };
00203
00204 #define Modulus modulus<C, modulus_int_preinverse<size> >
00205 #define Modular modular<Modulus, V>
00206
00207 template<typename C, typename V, nat size>
00208 class multiplier<Modular> {
00209 typedef modulus_multiplier_int_preinverse_helper<size> multiplier_helper;
00210 C rep;
00211 C pre;
00212 public:
00213 inline C operator * () const { return rep; }
00214 inline multiplier<Modular> () : rep (0) {
00215 multiplier_helper::set (pre, rep, Modular::get_modulus ()); }
00216 inline multiplier<Modular> (const Modular& s) : rep (*s) {
00217 multiplier_helper::set (pre, rep, Modular::get_modulus ()); }
00218 inline multiplier<Modular> (const multiplier<Modular>& s)
00219 : rep (s.rep), pre (s.pre) {}
00220 inline multiplier<Modular>& operator = (const multiplier<Modular>& s) {
00221 rep = s.rep; pre = s.pre; return *this; }
00222 inline bool
00223 operator == (const multiplier<Modular>& a) const {
00224 VERIFY (pre == a.pre,"multiplier<Modular>: pre-computed value differ");
00225 return rep == a.rep; }
00226 inline bool
00227 operator != (const multiplier<Modular>& a) const {
00228 VERIFY (pre != a.pre,"multiplier<Modular>: pre-computed value are equal");
00229 return rep != a.rep; }
00230 static inline void
00231 rmul (Modular& dest, const multiplier<Modular>& s) {
00232 multiplier_helper::mul (dest.rep, s.rep,
00233 Modular::get_modulus (), s.pre); }
00234 static inline void
00235 lmul (Modular& dest, const multiplier<Modular>& s) {
00236 multiplier_helper::mul (dest.rep, s.rep,
00237 Modular::get_modulus (), s.pre); }
00238 };
00239
00240 #undef Modulus
00241 #undef Modular
00242 }
00243 #endif //__MMX__MULTIPLIER__HPP
