00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__MODULAR_INT__HPP
00014 #define __MMX__MODULAR_INT__HPP
00015 #include <basix/operators.hpp>
00016 #include <basix/int.hpp>
00017 #include <numerix/rational.hpp>
00018 #include <numerix/modular.hpp>
00019
00020 namespace mmx {
00021
00025
00026 #define TMPL template<typename C, typename M>
00027
00028
00029
00030
00031
00032
00033
00034 template<nat size>
00035 struct modulus_maximum_size_int {
00036
00037 template<typename C, bool is_signed>
00038 struct _helper {
00039 static const nat value = (size < 8*sizeof(C) ? size : 8*sizeof(C)-1);
00040 };
00041
00042 template<typename C>
00043 struct _helper<C, false> {
00044 static const nat value = (size <= 8*sizeof(C) ? size : 8*sizeof(C));
00045 };
00046
00047 template<typename C>
00048 struct maximum_size_helper {
00049 static const nat value = _helper<C, is_signed_helper<C>::value>::value;
00050 static inline C dyn_value () { return value; }
00051 };
00052
00053 template<typename C>
00054 struct maximum_value_helper {
00055 static const nat _size = maximum_size_helper<C>::value;
00056 typedef typename unsigned_of_helper<C>::type uC;
00057 static const C value = MMX_SAFE_LEFT_SHIFT_INT (uC, 1, _size) - 1;
00058 static inline C dyn_value () { return value; }
00059 };
00060 };
00061
00062 template<typename V>
00063 struct modulus_normalization_int_naive : public V {
00064
00065 template<typename C> static inline bool
00066 normalize (C& p) {
00067 typedef typename unsigned_of_helper<C>::type uC;
00068 static const uC a = MMX_SAFE_LEFT_SHIFT_INT(uC, 1,
00069 V::template maximum_size_helper<C>::value);
00070 p = abs (p);
00071 if ((uC) p == a || p == 0) { p = a; return true; }
00072 return p <= V::template maximum_value_helper<C>::dyn_value ();
00073 }
00074 };
00075
00076 template<typename V>
00077 struct modulus_reduction_int_naive : public V {
00078 template<typename C> static inline void
00079 reduce_mod_core (C& dest, const C& p) {
00080 if (p != 0) dest %= p; }
00081
00082 template<typename C> static inline void
00083 reduce_mod_core (C& dest, const C& p, C& carry) {
00084 if (p != 0) { carry= dest / p; dest %= p; }
00085 else carry = 0; }
00086
00087 TMPL static inline void
00088 reduce_mod (C& dest, const M& m) {
00089 reduce_mod_core (dest, (C) m.p); }
00090
00091 TMPL static inline void
00092 reduce_mod (C& dest, const M& m, C& carry) {
00093 reduce_mod_core (dest, (C) m.p, carry); }
00094
00095 TMPL static inline void
00096 reduce_mod (C& dest, const C& s, const M& m) {
00097 dest = s;
00098 reduce_mod_core (dest, (C) m.p); }
00099
00100 TMPL static inline void
00101 reduce_mod (C& dest, const C& s, const M& m, C& carry) {
00102 dest = s;
00103 reduce_mod_core (dest, (C) m.p, carry); }
00104 };
00105
00106 template<typename V>
00107 struct modulus_add_int_naive : public V {
00108
00109 TMPL static inline void
00110 neg_mod (C& dest, const M& m) {
00111 if (dest != 0) dest = m.p - dest; }
00112
00113 TMPL static inline void
00114 neg_mod (C& dest, const M& m, C& carry) {
00115 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00116 if (dest != 0 || carry != 0) { dest= m.p - dest - carry; carry= 1; } }
00117
00118 TMPL static inline void
00119 neg_mod (C& dest, const C& s, const M& m) {
00120 if (s != 0) dest = m.p - s; else dest = s; }
00121
00122 TMPL static inline void
00123 neg_mod (C& dest, const C& s, const M& m, C& carry) {
00124 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00125 if (s != 0 || carry != 0) { dest= m.p - s - carry; carry= 1; }
00126 else dest= 0; }
00127
00128 template<typename C> static inline void
00129 add_mod_without_overflow (C& dest, const C& s, const C& p) {
00130 dest += s;
00131 if (dest >= p) dest -= p; }
00132
00133 template<typename C> static inline void
00134 add_mod_without_overflow (C& dest, const C& s, const C& p, C& carry) {
00135 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00136 dest += s + carry;
00137 if (dest >= p) { dest -= p; carry= 1; } else carry= 0; }
00138
00139 template<typename C> static inline void
00140 add_mod_with_overflow (C& dest, const C& s, const C& p) {
00141 dest += s;
00142 if (dest < s) { dest -= p; return; }
00143 if (dest >= p) dest -= p; }
00144
00145 template<typename C> static inline void
00146 add_mod_with_overflow (C& dest, const C& s, const C& p, C& carry) {
00147 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00148 dest += s + carry;
00149 if (dest < s || dest >= p) { dest -= p; carry= 1; }
00150 else carry= 0; }
00151
00152
00153 template<typename C, bool b>
00154 struct add_mod_helper {
00155 static inline void
00156 op (C& dest, const C& s, const C& p) {
00157 typedef typename unsigned_of_helper<C>::type uC;
00158 uC t = dest, up = p, us = s;
00159 add_mod_without_overflow (t, us, up);
00160 dest = t; }
00161 static inline void
00162 op (C& dest, const C& s, const C& p, C& carry) {
00163 typedef typename unsigned_of_helper<C>::type uC;
00164 uC t = dest, up = p, us = s, ucarry= carry;
00165 add_mod_without_overflow (t, us, up, ucarry);
00166 dest = t; carry= ucarry; }
00167 };
00168
00169
00170 template<typename C>
00171 struct add_mod_helper<C,false> {
00172 static inline void
00173 op (C& dest, const C& s, const C& p) {
00174 if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
00175 add_mod_with_overflow (dest, s, p);
00176 else
00177 add_mod_without_overflow (dest, s, p); }
00178 static inline void
00179 op (C& dest, const C& s, const C& p, C& carry) {
00180 if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
00181 add_mod_with_overflow (dest, s, p, carry);
00182 else
00183 add_mod_without_overflow (dest, s, p, carry); }
00184 };
00185
00186 template<typename C> static inline void
00187 add_mod_core (C& dest, const C& s, const C& p) {
00188 add_mod_helper<C,is_signed_helper<C>::value>::op (dest, s, p); }
00189
00190 template<typename C> static inline void
00191 add_mod_core (C& dest, const C& s, const C& p, C& carry) {
00192 add_mod_helper<C,is_signed_helper<C>::value>::op (dest, s, p, carry); }
00193
00194 TMPL static inline void
00195 add_mod (C& dest, const C& s, const M& m) {
00196 add_mod_core (dest, s, m.p); }
00197
00198 TMPL static inline void
00199 add_mod (C& dest, const C& s, const M& m, C& carry) {
00200 add_mod_core (dest, s, m.p, carry); }
00201
00202 TMPL static inline void
00203 add_mod (C& dest, const C& s1, const C& s2, const M& m) {
00204 dest = s1;
00205 add_mod (dest, s2, m); }
00206
00207 TMPL static inline void
00208 add_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
00209 dest = s1;
00210 add_mod (dest, s2, m, carry); }
00211
00212 template<typename C> static inline void
00213 sub_mod_core (C& dest, const C& s, const C& p) {
00214 if (dest < s) dest += p;
00215 dest -= s; }
00216
00217 template<typename C> static inline void
00218 sub_mod_core (C& dest, const C& s, const C& p, C& carry) {
00219 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00220 C t = s + carry;
00221 if (t == p || dest < t) { dest += p; carry= 1; } else carry= 0;
00222 dest -= t; }
00223
00224 TMPL static inline void
00225 sub_mod (C& dest, const C& s, const M& m) {
00226 sub_mod_core (dest, s, m.p); }
00227
00228 TMPL static inline void
00229 sub_mod (C& dest, const C& s, const M& m, C& carry) {
00230 sub_mod_core (dest, s, m.p, carry); }
00231
00232 TMPL static inline void
00233 sub_mod (C& dest, const C& s1, const C& s2, const M& m) {
00234 dest = s1;
00235 sub_mod (dest, s2, m); }
00236
00237 TMPL static inline void
00238 sub_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
00239 dest = s1;
00240 sub_mod (dest, s2, m, carry); }
00241 };
00242
00243 template<typename V>
00244 struct modulus_mul_int_naive : public V {
00245
00246 template <typename C, typename D>
00247 struct mul_mod_helper {
00248 static inline void op (C& dest, const C& s, const C& p) {
00249 D a (dest);
00250 a *= s;
00251 if (p != 0) a %= p;
00252 dest = (C) a;
00253 }
00254 static inline void op (C& dest, const C& s, const C& p, C& carry) {
00255 D a (dest);
00256 a = a * s + carry;
00257 if (p != 0) { carry= a / p; a %= p; }
00258 dest = (C) a;
00259 }
00260 };
00261
00262 template <typename C>
00263 struct mul_mod_helper<C, void> {
00264 static const nat half_size = 4 * sizeof (C);
00265 static const nat size = 2 * half_size;
00266 static const C lo_mask = ((C) -1) >> half_size;
00267 static const C hi_mask = ((C) -1) << half_size;
00268
00269 static inline C
00270 lo (const C& s) { return s & lo_mask; }
00271
00272 static inline C
00273 hi (const C& s) { return (s & hi_mask) >> half_size; }
00274
00275 static inline void
00276 half_lshift_mod(C& dest, const C& p, const C& lo_p,const C& hi_p) {
00277
00278
00279 C q = dest / hi_p;
00280 C r = dest - q * hi_p;
00281 q *= lo_p;
00282 q %= p;
00283 dest = r << half_size;
00284 V::sub_mod_core (dest, q, p);
00285 }
00286
00287 static inline void
00288 half_lshift_mod(C& dest, const C& p, const C& lo_p,const C& hi_p,
00289 C& carry) {
00290
00291
00292 C q = dest / hi_p, c = 0; carry = q;
00293 C r = dest - q * hi_p;
00294 q *= lo_p; carry -= (q / p);
00295 q %= p;
00296 dest = r << half_size;
00297 V::sub_mod_core (dest, q, p, c);
00298 if (c != 0) carry--;
00299 }
00300
00301 static inline void op (C& dest, const C& s, const C& p) {
00302 C hi_p = hi (p);
00303 if (hi_p) {
00304 C lo_p = lo (p);
00305 C lo_dest = lo (dest);
00306 C lo_s = lo (s);
00307 C hi_dest = hi (dest);
00308 C hi_s = hi (s);
00309 dest = lo_dest * lo_s;
00310 V::reduce_mod_core (dest, p);
00311 C t = lo_dest * hi_s;
00312 half_lshift_mod (t, p, lo_p, hi_p);
00313 V::add_mod_core (dest, t, p);
00314 t = hi_dest * lo_s;
00315 half_lshift_mod (t, p, lo_p, hi_p);
00316 V::add_mod_core(dest, t, p);
00317 t = hi_dest * hi_s;
00318 half_lshift_mod (t, p, lo_p, hi_p);
00319 half_lshift_mod (t, p, lo_p, hi_p);
00320 V::add_mod_core(dest, t, p);
00321 }
00322 else {
00323 dest *= s;
00324 if (p != 0) dest %= p;
00325 }
00326 }
00327
00328 static inline void op (C& dest, const C& s, const C& p, C& carry) {
00329 C hi_p = hi (p);
00330 if (hi_p) {
00331 C lo_p = lo (p);
00332 C lo_dest = lo (dest);
00333 C lo_s = lo (s);
00334 C hi_dest = hi (dest);
00335 C hi_s = hi (s);
00336 C c= 0, cc= 0;
00337 dest = lo_dest * lo_s;
00338 reduce_mod_core (dest, p, cc);
00339 C t = lo_dest * hi_s;
00340 half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
00341 V::add_mod_core (dest, t, p, c); cc += c; c= 0;
00342 t = hi_dest * lo_s;
00343 half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
00344 V::add_mod_core(dest, t, p, c); cc += c; c= 0;
00345 t = hi_dest * hi_s;
00346 half_lshift_mod (t, p, lo_p, hi_p, c); cc += c << half_size; c= 0;
00347 half_lshift_mod (t, p, lo_p, hi_p, c); cc += c; c= 0;
00348 V::add_mod_core(dest, t, p, c); cc += c; c= 0;
00349 V::add_mod_core(dest, carry, p, c); carry = cc + c;
00350 }
00351 else {
00352 dest = dest * s + carry;
00353 if (p != 0) { carry= dest / p; dest %= p; }
00354 }
00355 }
00356 };
00357
00358 TMPL static inline void
00359 mul_mod (C& dest, const C& s, const M& m) {
00360 if (is_signed_helper<C>::value) {
00361 typedef typename unsigned_of_helper<C>::type uC;
00362 uC t = dest;
00363 mul_mod (t, (uC) s, m);
00364 dest = t;
00365 }
00366 else {
00367 typedef typename unsigned_int_with_double_size_helper<C>::type D;
00368 mul_mod_helper<C,D>::op (dest, s, m.p);
00369 }
00370 }
00371
00372 TMPL static inline void
00373 mul_mod (C& dest, const C& s, const M& m, C& carry) {
00374 if (is_signed_helper<C>::value) {
00375 typedef typename unsigned_of_helper<C>::type uC;
00376 uC t = dest, ucarry= carry;
00377 mul_mod (t, (uC) s, m, ucarry);
00378 dest = t; carry= ucarry;
00379 }
00380 else {
00381 typedef typename unsigned_int_with_double_size_helper<C>::type D;
00382 mul_mod_helper<C,D>::op (dest, s, m.p, carry);
00383 }
00384 }
00385
00386 TMPL static inline void
00387 mul_mod (C& dest, const C& s1, const C& s2, const M& m) {
00388 dest = s1;
00389 mul_mod (dest, s2, m); }
00390
00391 TMPL static inline void
00392 mul_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
00393 dest = s1;
00394 mul_mod (dest, s2, m, carry); }
00395
00396 };
00397
00398 template<typename V>
00399 struct modulus_inv_int_naive : public V {
00400
00401 template<typename C> static void
00402 inv_mod_signed (C& a, const C& p) {
00403 typedef typename unsigned_of_helper<C>::type uC;
00404 typedef typename signed_of_helper<C>::type sC;
00405 sC coa0=0, coa1=1, r0=p, r1=a, q, t;
00406 if ((p == 0) && (r1 != 0)) {
00407 q = (((uC) r0) - ((uC) r1)) / ((uC) r1) + 1;
00408 t = r0 - q * r1;
00409 r0 = r1;
00410 r1 = t;
00411 t = coa1;
00412 coa1 = coa0 - q * coa1;
00413 coa0 = t;
00414 }
00415 while (r1 != 0) {
00416 q = r0 / r1;
00417 t = r0 - q * r1;
00418 r0 = r1;
00419 r1 = t;
00420 t = coa1;
00421 coa1 = coa0 - q * coa1;
00422 coa0 = t;
00423 }
00424 if (r0 != 1)
00425 ERROR ("inv_mod: argument is not invertible");
00426 a = coa0 < 0 ? p+coa0 : coa0;
00427 }
00428
00429 template<typename C> static void
00430 inv_mod_unsigned (C& a, const C& p) {
00431 C coa0=0, coa1=1, r0=p, r1=a, q, t;
00432 bool toggle = true;
00433 if ((p == 0) && (r1 != 0)) {
00434 q = (r0-r1) / r1 + 1;
00435 t = r0 - q * r1;
00436 r0 = r1;
00437 r1 = t;
00438 t = coa1;
00439 coa1 = coa0 - q * coa1;
00440 coa0 = t;
00441 toggle = !toggle;
00442 }
00443 while (r1 != 0) {
00444 q = r0 / r1;
00445 t = r0 - q * r1;
00446 r0 = r1;
00447 r1 = t;
00448 t = coa1;
00449 coa1 = coa0 - q * coa1;
00450 coa0 = t;
00451 toggle = !toggle;
00452 }
00453 if (r0 != 1)
00454 ERROR ("inv_mod: argument is not invertible");
00455 a = (toggle && (coa0 != 0)) ? p+coa0 : coa0;
00456 }
00457
00458 TMPL static inline bool
00459 is_invertible_mod (const C& s, const M& m) {
00460 return abs (gcd (s, m.p)) == 1; }
00461
00462 TMPL static inline void
00463 inv_mod (C& a, const M& m) {
00464 if (V::template maximum_size_helper<C>::value >= 8*sizeof(C))
00465 inv_mod_unsigned(a, m.p);
00466 else
00467 inv_mod_signed(a, m.p); }
00468
00469 TMPL static inline void
00470 inv_mod (C& dest, const C& s, const M& m) {
00471 dest = s;
00472 inv_mod (dest, m); }
00473 };
00474
00475 template<typename V>
00476 struct modulus_div_int_naive : public V {
00477
00478 TMPL static inline void
00479 div_mod (C& dest, const C& s, const M& m) {
00480 C t = s;
00481 V::inv_mod (t, m);
00482 V::mul_mod (dest, t, m); }
00483
00484 TMPL static inline void
00485 div_mod (C& dest, const C& s1, const C& s2, const M& m) {
00486 dest = s1;
00487 div_mod (dest, s2, m); }
00488 };
00489
00490 template<typename V>
00491 struct modulus_encoding_int_naive : public V {
00492 TMPL static inline void
00493 encode_mod (C& dest, const C& s, const M& m) {
00494 typedef typename unsigned_of_helper<C>::type uC;
00495 if (sign (s) < 0) {
00496 uC tmp = - ((uC) s);
00497 V::reduce_mod (tmp, m);
00498 dest = (C) (((uC) m.p) - tmp);
00499 }
00500 else
00501 dest = s;
00502 V::reduce_mod (dest, m);
00503 }
00504 TMPL static inline void
00505 decode_mod (C& dest, const C& s, const M& m) {
00506 (void) m;
00507 dest = s; }
00508 };
00509
00510 template<int size>
00511 struct modulus_int_naive:
00512 public modulus_encoding_int_naive<
00513 modulus_div_int_naive<
00514 modulus_inv_int_naive<
00515 modulus_mul_int_naive<
00516 modulus_add_int_naive<
00517 modulus_reduction_int_naive<
00518 modulus_normalization_int_naive<
00519 modulus_maximum_size_int<size> > > > > > > > {};
00520
00521 #undef TMPL
00522
00523
00524
00525
00526
00527 #define TMPL template<typename C>
00528
00529 template<typename V>
00530 struct modulus_mul_int_preinverse : public V {
00531
00532
00533
00534
00535 template <typename C, typename D>
00536 struct _dynamic_inverse_helper {
00537 typedef typename unsigned_of_helper<C>::type uC;
00538 typedef typename unsigned_of_helper<D>::type uD;
00539 static const nat m = V::template maximum_size_helper<C>::value;
00540
00541 static inline uC
00542 op (const uC& p, nat r, nat s, nat t) {
00543 (void) r;
00544 if (p == 1 && s+t == (nat) -1) return 0;
00545 uD loc_p = (p == 0) ? ((uD) 1) << m : (uD) ((uC) p);
00546 return (((uD) 1) << (s+t)) / loc_p; }
00547 };
00548
00549 template <typename C>
00550 struct _dynamic_inverse_helper<C,void> {
00551 static const nat n2 = 4 * sizeof(C);
00552 static const nat n = 8 * sizeof(C);
00553 static const nat m = V::template maximum_size_helper<C>::value;
00554 typedef typename unsigned_of_helper<C>::type uC;
00555 typedef long_int_mul_op<C> long_mul;
00556 typedef long_int_rshift_op<C> long_rshift;
00557 typedef long_int_sub_op<C> long_sub;
00558
00559 static uC
00560 _q0 (const uC& p, nat r, nat s, nat t, nat s0, nat s1) {
00561
00562 uC q0;
00563 if (r <= s0)
00564 q0 = (((uC) 1) << (s0+r-1)) / p;
00565 else {
00566 uC p1 = (p == 0) ? ((uC) 1) << s0 : p >> (r-s0);
00567 uC p0 = p & (((uC) -1) >> (n+s0-r));
00568 uC d = ((uC) 1) << s1;
00569 q0 = d / p1;
00570 uC b = (d % p1) << (r-s0);
00571 uC c = q0 * p0;
00572 if (c >= p) { q0--; c -= p; }
00573 if (b < c) q0--;
00574 }
00575 return q0;
00576 }
00577
00578 static uC
00579 op (const uC& p, nat r, nat s, nat t) {
00580
00581
00582 if (p == 1) {
00583 if (m+2 <= n) return 2;
00584 if (m+1 <= n) return 1;
00585 return 0;
00586 }
00587 if (p == 2) {
00588 if (m+2 <= n) return 4;
00589 if (m+1 <= n) return 2;
00590 return 1;
00591 }
00592 nat s0 = (s+t-(r-1)) >> 1;
00593 nat s1 = 2*s0 - 1;
00594 uC q0 = _q0 (p, r, s, t, s0, s1);
00595 uC a0 = 0, a1 = 0;
00596 long_mul::op (a1, a0, q0 * q0, p);
00597 if (long_rshift::op_b (a1, a0, r)) a0++;
00598 uC c = (q0 << s0) - a0;
00599 uC d;
00600 if (m+1 <= n || s1+r <= n) {
00601 if (s1+r-1 >= n) d = 0; else d = (((uC) 1) << (s1+r-1));
00602 d -= c*p;
00603 if (d >= p) { c++; d-=p; }
00604 }
00605 else {
00606 uC d1, e0, e1;
00607 long_mul::op (e1, e0, c, p);
00608 d = 0;
00609 d1 = ((uC) 1) << (s1+r-1-n);
00610 long_sub::op (d1, d, e1, e0);
00611 if ((d1 != 0) || (d >= p)) { c++; d-=p; }
00612 }
00613 uC d2 = d << 1;
00614 c <<= 1;
00615 if ((d2 >= p) || (d2 < d)) { c++; d = d2-p; } else d = d2;
00616 if (s1+1 == s+t-(r-1)) return c;
00617 d2 = d << 1;
00618 c <<= 1;
00619 if ((d2 >= p) || (d2 < d)) c++;
00620 return c;
00621 }
00622 };
00623
00624 TMPL static inline nat
00625 dyn_r (const C& p) {
00626 typedef typename unsigned_of_helper<C>::type uC;
00627 static const nat m = V::template maximum_size_helper<C>::value;
00628 C abs_p = abs (p);
00629 nat r = bit_size (abs_p);
00630 return (r == 0) ? m : ( (uC) abs_p == (((uC) 1) << (r-1)) ? r-1 : r );
00631 }
00632
00633 TMPL static inline nat
00634 dyn_s (const C& p, nat r) {
00635 static const nat m = V::template maximum_size_helper<C>::value;
00636 return m+2 <= 8 * sizeof(C) ? r-2 : r-1;
00637 }
00638
00639 TMPL static inline nat
00640 dyn_t (const C& p, nat r) {
00641 static const nat n = 8 * sizeof(C);
00642 static const nat m = V::template maximum_size_helper<C>::value;
00643 return (m+2 <= n) ? r+3 : ((m+1 <= n) ? r+1 : r);
00644 }
00645
00646 TMPL static inline C
00647 dyn_q (const C& p, nat r, nat s, nat t) {
00648 typedef typename unsigned_of_helper<C>::type uC;
00649 typedef typename unsigned_int_with_double_size_helper<C>::type D;
00650 return _dynamic_inverse_helper<C,D>::op ((uC) p, r, s, t);
00651 }
00652
00653
00654
00655
00657 template<typename C, C p>
00658 struct auxiliaries_helper {
00659 static const nat m = V::template maximum_size_helper<C>::value;
00660 typedef typename unsigned_of_helper<C>::type uC;
00661 typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
00662
00663 static const nat n2= 4*sizeof(C);
00664 static const nat n = 8*sizeof(C);
00665 static const uC abs_p = (uC) ((p < 0) ? -p : p);
00666 static const uC a = MMX_SAFE_LEFT_SHIFT_INT(uC,1,m);
00667 static const uC up = (abs_p == 0) ? a : abs_p;
00668 static const nat bitsize = int_bitsize_helper<uC, up>::value;
00669 static const nat _b = (bitsize == 0) ? n : bitsize-1;
00670 static const nat r = (bitsize == 0) ? n :
00671 ( (up == MMX_SAFE_LEFT_SHIFT_INT(uC,1,_b)) ?
00672 bitsize-1 : bitsize );
00673 static const nat s = (m+2 <= n) ? r-2 : r-1;
00674 static const nat t = (m+2 <= n) ? r+3 : ((m+1 <= n) ? r+1 : r);
00675
00676 template <typename D, typename Dummy>
00677 struct _inverse_helper {
00678 static const uD _up = (up == 0) ? ((uD) 1) << m : (uD) up;
00679 static const nat _s = (s+t >= 2*n) ? 2*n : s+t;
00680 static const uD _t = (uD) (MMX_SAFE_LEFT_SHIFT_INT(uD,1,_s));
00681 static const uC value = (up == 1 && m >= n) ?
00682 (uC) 0 : (uC) (_t / (uD) _up);
00683 static inline C dyn_value () { return value; }
00684 };
00685
00686 template <typename Dummy>
00687 struct _inverse_helper <void, Dummy> {
00688 static inline C dyn_value () { return dyn_q ((C) up, r, s, t); }
00689 };
00690
00691 static inline C q () {
00692 return _inverse_helper<uD, void>::dyn_value (); }
00693 };
00694
00695
00696
00697
00698
00699 template <typename C, typename D, nat m>
00700 struct mul_mod_helper {
00701 static const nat n = 8 * sizeof(C);
00702
00703 static inline void
00704 op (C& dest, const C& src, const C& p,
00705 const C& q, nat r, nat s, nat t) {
00706 D a (dest);
00707 a *= src;
00708 if (m + 2 <= n) {
00709 C hi_a = a >> s;
00710 C b = (((D) hi_a) * q) >> t;
00711 dest = ((C) a) - b * p;
00712 if (dest >= p) dest -=p;
00713 return;
00714 }
00715 D hi_a = a >> s;
00716 C b = (hi_a * q) >> t;
00717 D c = a - ((D) b) * p;
00718 if (m + 1 <= n) {
00719 if (c >= p) {
00720 dest = c - p;
00721 if (dest >= p)
00722 dest -= p;
00723 }
00724 else
00725 dest = c;
00726 return;
00727 }
00728 if (c >= p) {
00729 c -= p;
00730 if (c >= p) {
00731 c -= p;
00732 if (c >= p)
00733 c -= p;
00734 }
00735 }
00736 dest = c;
00737 }
00738
00739 static inline void
00740 op (C& dest, const C& src, const C& p, C& carry,
00741 const C& q, nat r, nat s, nat t) {
00742 D a (dest);
00743 a *= src;
00744 if (m + 2 <= n) {
00745 C hi_a = a >> s;
00746 C b = (((D) hi_a) * q) >> t, c = 0;;
00747 dest = ((C) a) - b * p;
00748 if (dest >= p) { dest -= p; b++; }
00749 V::add_mod_core (dest, carry, p, c); carry = b + c;
00750 return;
00751 }
00752 D hi_a = a >> s;
00753 C b = (hi_a * q) >> t, z= 0;
00754 D c = a - ((D) b) * p;
00755 if (m + 1 <= n) {
00756 if (c >= p) {
00757 dest = c - p; b++;
00758 if (dest >= p) {
00759 dest -= p; b++; }
00760 }
00761 else
00762 dest = c;
00763 V::add_mod_core (dest, carry, p, z); carry = b + z;
00764 return;
00765 }
00766 if (c >= p) {
00767 c -= p; b++;
00768 if (c >= p) {
00769 c -= p; b++;
00770 if (c >= p) {
00771 c -= p; b++; }
00772 }
00773 }
00774 dest = c;
00775 V::add_mod_core (dest, carry, p, z); carry = b + z;
00776 }
00777 };
00778
00779
00780
00781
00782
00783 template <typename C, nat m>
00784 struct mul_mod_helper<C, void, m> {
00785
00786 static const nat n2 = 4 * sizeof(C);
00787 static const nat n = 8 * sizeof(C);
00788 typedef long_int_mul_op<C> long_mul;
00789 typedef long_int_rshift_op<C> long_rshift;
00790 typedef long_int_sub_op<C> long_sub;
00791 typedef long_int_ge_op<C> long_ge;
00792
00793 static void
00794 op (C& dest, const C& src, const C& p,
00795 const C& q, nat r, nat s, nat t) {
00796
00797 C a1, a0, hi_a1, hi_a0;
00798 C tmp = dest;
00799
00800 long_mul::op (a1, a0, tmp, src);
00801 hi_a1 = a1;
00802 hi_a0 = a0;
00803 long_rshift::op (hi_a1, hi_a0, s);
00804 C b1, b0;
00805 long_mul::op (b1, b0, hi_a0, q);
00806 if (m >= n) {if (hi_a1 != 0) b1 += q;}
00807 if (m+1 <= n) {
00808 if (t >= n)
00809 if (n >= t) b0 = b1 << (n-t); else b0 = b1 >> (t-n);
00810 else
00811 b0 = (b0 >> t) | (b1 << (n-t));
00812 }
00813 else
00814 long_rshift::op (b1, b0, t);
00815 if (m+2 <= n) {
00816 dest = a0 - b0*p;
00817 dest = min (dest, dest - p);
00818 return;
00819 }
00820 C d1, d0;
00821 long_mul::op (d1, d0, b0, p);
00822 long_sub::op (a1, a0, d1, d0);
00823 if (m+1 <= n) {
00824 if (long_ge::op (a1, a0, 0, p)) {
00825 dest = a0 - p;
00826 if (dest >= p)
00827 dest -=p; }
00828 else dest = a0;
00829 return;
00830 }
00831 if (long_ge::op (a1, a0, 0, p)) {
00832 long_sub::op (a1, a0, 0, p);
00833 if (long_ge::op (a1, a0, 0, p)) {
00834 long_sub::op (a1, a0, 0, p);
00835 if (long_ge::op (a1, a0, 0, p))
00836 a0 -= p;
00837 }
00838 }
00839 dest = a0;
00840 }
00841
00842 static void
00843 op (C& dest, const C& src, const C& p, C& carry,
00844 const C& q, nat r, nat s, nat t) {
00845
00846 C a1, a0, hi_a1, hi_a0;
00847 C tmp = dest;
00848
00849 long_mul::op (a1, a0, tmp, src);
00850 hi_a1 = a1;
00851 hi_a0 = a0;
00852 long_rshift::op (hi_a1, hi_a0, s);
00853 C b1, b0, c = 0;
00854 long_mul::op (b1, b0, hi_a0, q);
00855 if (m >= n) {if (hi_a1 != 0) b1 += q;}
00856 if (m+1 <= n) {
00857 if (t >= n)
00858 if (n >= t) b0 = b1 << (n-t); else b0 = b1 >> (t-n);
00859 else
00860 b0 = (b0 >> t) | (b1 << (n-t));
00861 }
00862 else
00863 long_rshift::op (b1, b0, t);
00864 if (m+2 <= n) {
00865 dest = a0 - b0 * p;
00866 if (dest >= p) { dest -= p; b0++; }
00867 V::add_mod_core (dest, carry, p, c); carry = b0 + c;
00868 return;
00869 }
00870 C d1, d0;
00871 long_mul::op (d1, d0, b0, p);
00872 long_sub::op (a1, a0, d1, d0);
00873 if (m+1 <= n) {
00874 if (long_ge::op (a1, a0, 0, p)) {
00875 dest = a0 - p; b0++;
00876 if (dest >= p) {
00877 dest -=p; b0++; } }
00878 else dest = a0;
00879 V::add_mod_core (dest, carry, p, c); carry = b0 + c;
00880 return;
00881 }
00882 if (long_ge::op (a1, a0, 0, p)) {
00883 long_sub::op (a1, a0, 0, p); b0++;
00884 if (long_ge::op (a1, a0, 0, p)) {
00885 long_sub::op (a1, a0, 0, p); b0++;
00886 if (long_ge::op (a1, a0, 0, p)) {
00887 a0 -= p; b0++; }
00888 }
00889 }
00890 dest = a0;
00891 V::add_mod_core (dest, carry, p, c); carry = b0 + c;
00892 }
00893 };
00894
00895 template<typename C, typename M> static inline void
00896 mul_mod (C& dest, const C& src, const M& x) {
00897 static const nat m = V::template maximum_size_helper<C>::value;
00898 typedef typename unsigned_of_helper<C>::type uC;
00899 typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
00900 uC tmp = dest;
00901 mul_mod_helper<uC,uD,m>::op (tmp, (uC) src, (uC) x.p,
00902 (uC) x.q, x.r, x.s, x.t);
00903 dest = tmp; }
00904
00905 template<typename C, typename M> static inline void
00906 mul_mod (C& dest, const C& src, const M& x, C& carry) {
00907 static const nat m = V::template maximum_size_helper<C>::value;
00908 typedef typename unsigned_of_helper<C>::type uC;
00909 typedef typename unsigned_int_with_double_size_helper<uC>::type uD;
00910 uC tmp = dest, ucarry = carry;
00911 mul_mod_helper<uC,uD,m>::op (tmp, (uC) src, (uC) x.p, ucarry,
00912 (uC) x.q, x.r, x.s, x.t);
00913 dest = tmp; carry = ucarry; }
00914
00915 template<typename C, typename M> static inline void
00916 mul_mod (C& dest, const C& s1, const C& s2, const M& x) {
00917 dest = s1;
00918 mul_mod (dest, s2, x); }
00919
00920 template<typename C, typename M> static inline void
00921 mul_mod (C& dest, const C& s1, const C& s2, const M& x, C& carry) {
00922 dest = s1;
00923 mul_mod (dest, s2, x, carry); }
00924 };
00925
00926 template<nat m>
00927 struct modulus_int_preinverse:
00928 public modulus_encoding_int_naive<
00929 modulus_div_int_naive<
00930 modulus_inv_int_naive<
00931 modulus_mul_int_preinverse<
00932 modulus_add_int_naive<
00933 modulus_reduction_int_naive<
00934 modulus_normalization_int_naive<
00935 modulus_maximum_size_int<m> > > > > > > > {};
00936
00937
00938
00939
00940
00941 template<typename C, nat m>
00942 class modulus <C, modulus_int_preinverse<m> > {
00943 MMX_ALLOCATORS
00944 public:
00945 C p;
00946 C q;
00947 nat r;
00948 nat s;
00949 nat t;
00950 C h;
00951 typedef C base;
00952 typedef modulus_int_preinverse<m> variant;
00953 typedef modulus <C, modulus_int_preinverse<m> > self_type;
00954
00955 inline C operator * () const { return p; }
00956
00957 inline modulus () : p (0) {
00958 r = variant::dyn_r (p);
00959 s = variant::dyn_s (p, r);
00960 t = variant::dyn_t (p, r);
00961 q = variant::dyn_q (p, r, s, t);
00962 h = p >> 1;
00963 }
00964
00965 template<typename X>
00966 inline modulus (const X& x) {
00967 p = (C) x;
00968 ASSERT (variant::normalize (p), "modulus: bad value");
00969 r = variant::dyn_r (p);
00970 s = variant::dyn_s (p, r);
00971 t = variant::dyn_t (p, r);
00972 q = variant::dyn_q (p, r, s, t);
00973 h = p >> 1;
00974 }
00975
00976 template<typename X, X _p>
00977 inline modulus (const fixed_value<X, _p>& x) {
00978 typedef typename variant::template auxiliaries_helper<C, _p> A;
00979 C loc_p = _p;
00980 ASSERT (variant::normalize (loc_p), "modulus: bad value");
00981 p = A::up;
00982 q = A::q ();
00983 r = A::r;
00984 s = A::s;
00985 t = A::t;
00986 h = p >> 1; }
00987
00988 inline modulus (const self_type& x) : p (x.p), q (x.q),
00989 r (x.r), s (x.s),
00990 t (x.t), h (x.h) {}
00991
00992 inline self_type&
00993 operator = (const self_type& x) {
00994 p = x.p;
00995 q = x.q;
00996 r = x.r;
00997 s = x.s;
00998 t = x.t;
00999 h = x.h;
01000 return *this; }
01001
01002 inline bool operator == (const self_type& a) const {
01003 return p == a.p; }
01004
01005 inline bool operator != (const self_type& a) const {
01006 return p != a.p; }
01007 };
01008
01009 #undef TMPL
01010
01011
01012
01013
01014
01015
01016 #define Variant modulus_int_preinverse
01017
01018 template<>
01019 struct modulus_variant_helper<unsigned char> {
01020 typedef Variant<8*sizeof(unsigned char)> MV; };
01021
01022 template<>
01023 struct modulus_variant_helper<short unsigned int> {
01024 typedef Variant<8*sizeof(short unsigned int)> MV; };
01025
01026 template<>
01027 struct modulus_variant_helper<unsigned int> {
01028 typedef Variant<8*sizeof(unsigned int)> MV; };
01029
01030 template<>
01031 struct modulus_variant_helper<long unsigned int> {
01032 typedef Variant<8*sizeof(long unsigned int)> MV; };
01033
01034 template<>
01035 struct modulus_variant_helper<long long unsigned int> {
01036 typedef Variant<8*sizeof(long long unsigned int)> MV; };
01037
01038 template<>
01039 struct modulus_variant_helper<char> {
01040 typedef Variant<8*sizeof(char)-1> MV; };
01041
01042 template<>
01043 struct modulus_variant_helper<signed char> {
01044 typedef Variant<8*sizeof(char)-1> MV; };
01045
01046 template<>
01047 struct modulus_variant_helper<short int> {
01048 typedef Variant<8*sizeof(short int)-1> MV; };
01049
01050 template<>
01051 struct modulus_variant_helper<int> {
01052 typedef Variant<8*sizeof(int)-1> MV; };
01053
01054 template<>
01055 struct modulus_variant_helper<long int> {
01056 typedef Variant<8*sizeof(long int)-1> MV; };
01057
01058 template<>
01059 struct modulus_variant_helper<long long int> {
01060 typedef Variant<8*sizeof(long long int)-1> MV; };
01061
01062 #undef Variant
01063
01064
01065
01066
01067
01068 #define TMPLVW template<typename V, typename W>
01069 #define Modular(I) modular<modulus<I,V>,W>
01070 #define DECLARE_HELPER(I) \
01071 UNARY_RETURN_TYPE(TMPLVW,lift_op,Modular(I),integer); \
01072 UNARY_RETURN_TYPE(STMPL,project_op,I,modular<modulus<I> >); \
01073 TMPLVW \
01074 struct lift_helper<Modular(I)> { \
01075 template<typename R> static inline void \
01076 set_op (R&y, const Modular(I)& x) { \
01077 set_as (y, *x); } \
01078 static inline integer \
01079 op (const Modular(I)& x) { \
01080 return as<integer> (*x); } \
01081 }; \
01082 STMPL \
01083 struct project_helper<I> { \
01084 TMPLVW static inline void \
01085 set_op (Modular(I)& y, const I& x) { \
01086 y= Modular(I) (x); } \
01087 static inline modular<modulus<I> > \
01088 op (const I& x) { \
01089 return modular<modulus<I> > (x); } \
01090 };
01091
01092 DECLARE_HELPER(unsigned char)
01093 DECLARE_HELPER(signed char)
01094 DECLARE_HELPER(unsigned short int)
01095 DECLARE_HELPER(signed short int)
01096 DECLARE_HELPER(unsigned int)
01097 DECLARE_HELPER(int)
01098 DECLARE_HELPER(unsigned long int)
01099 DECLARE_HELPER(long int)
01100 DECLARE_HELPER(unsigned long long int)
01101 DECLARE_HELPER(long long int)
01102 #undef DECLARE_HELPER
01103 #undef Modular
01104 #undef TMPLVW
01105
01106
01107
01108
01109
01110 #define DECLARE_HELPER(I) \
01111 template<typename V, typename W> \
01112 struct as_helper<modular<modulus<I,V>,W>, rational> { \
01113 static inline modular<modulus<I,V>,W> \
01114 cv (const rational& a) { \
01115 return as<modular<modulus<I,V>,W> > (numerator (a)) \
01116 / as<modular<modulus<I,V>,W> > (denominator (a)); } \
01117 };
01118
01119 DECLARE_HELPER(unsigned char)
01120 DECLARE_HELPER(signed char)
01121 DECLARE_HELPER(unsigned short int)
01122 DECLARE_HELPER(signed short int)
01123 DECLARE_HELPER(unsigned int)
01124 DECLARE_HELPER(int)
01125 DECLARE_HELPER(unsigned long int)
01126 DECLARE_HELPER(long int)
01127 DECLARE_HELPER(unsigned long long int)
01128 DECLARE_HELPER(long long int)
01129 #undef DECLARE_HELPER
01130
01131
01132
01133
01134
01135 #define DECLARE_HELPER(I) \
01136 template<typename V, typename W> \
01137 rational reconstruct (const modular<modulus<I,V>,W>& x) { \
01138 typedef signed_of_helper<I>::type J; \
01139 J n, d; \
01140 bool b= reconstruct (n, d, *x, *get_modulus (x)); \
01141 if (b) return rational (n) / rational (d); \
01142 ERROR ("rational reconstruction failed"); } \
01143 template<typename V, typename W> \
01144 bool is_reconstructible (const modular<modulus<I,V>,W>& x, \
01145 rational& r) { \
01146 typedef signed_of_helper<I>::type J; \
01147 J n, d; \
01148 bool b= reconstruct (n, d, *x, *get_modulus (x)); \
01149 if (!b) return false; \
01150 r= rational (n) / rational (d); \
01151 return true; }
01152
01153 DECLARE_HELPER(signed char)
01154 DECLARE_HELPER(signed short int)
01155 DECLARE_HELPER(int)
01156 DECLARE_HELPER(long int)
01157 DECLARE_HELPER(long long int)
01158 DECLARE_HELPER(unsigned char)
01159 DECLARE_HELPER(unsigned short int)
01160 DECLARE_HELPER(unsigned int)
01161 DECLARE_HELPER(unsigned long int)
01162 DECLARE_HELPER(unsigned long long int)
01163 #undef DECLARE_HELPER
01164
01165 }
01166 #endif //__MMX__MODULAR_INT__HPP