00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__MODULAR_INTEGER__HPP
00014 #define __MMX__MODULAR_INTEGER__HPP
00015 #include <numerix/rational.hpp>
00016 #include <numerix/modular.hpp>
00017 #include <numerix/modular_int.hpp>
00018
00019 namespace mmx {
00020
00021 #define TMPL template<typename C, typename M>
00022
00023
00024
00025
00026
00027
00028
00029 struct modulus_normalization_integer_naive {
00030
00031 template<typename C> static inline bool
00032 normalize (C& p) {
00033 if (p < 0) p = -p;
00034 return true;
00035 }
00036 };
00037
00038 template<typename V>
00039 struct modulus_add_integer_naive : public V {
00040
00041 TMPL static inline void
00042 neg_mod (C& dest, const M& m) {
00043 if (dest != 0) dest = m.p - dest; }
00044
00045 TMPL static inline void
00046 neg_mod (C& dest, const M& m, C& carry) {
00047 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00048 if (dest != 0 || carry != 0) { dest= m.p - dest - carry; carry= 1; } }
00049
00050 TMPL static inline void
00051 neg_mod (C& dest, const C& s, const M& m) {
00052 if (s != 0) dest = m.p - s; else dest = s;
00053 }
00054
00055 TMPL static inline void
00056 neg_mod (C& dest, const C& s, const M& m, C& carry) {
00057 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00058 if (s != 0 || carry != 0) { dest= m.p - s - carry; carry= 1; }
00059 else dest= 0; }
00060
00061 TMPL static inline void
00062 add_mod (C& dest, const C& s, const M& m) {
00063 dest += s;
00064 if (dest >= m.p) dest -= m.p; }
00065
00066 TMPL static inline void
00067 add_mod (C& dest, const C& s, const M& m, C& carry) {
00068 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00069 dest += s + carry;
00070 if (dest >= m.p) { dest -= m.p; carry= 1; } else carry= 0; }
00071
00072 TMPL static inline void
00073 add_mod (C& dest, const C& s1, const C& s2, const M& m) {
00074 dest = s1;
00075 add_mod (dest, s2, m); }
00076
00077 TMPL static inline void
00078 add_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
00079 dest = s1;
00080 add_mod (dest, s2, m, carry); }
00081
00082 template<typename C> static inline void
00083 sub_mod_core (C& dest, const C& s, const C& p) {
00084 if (dest < s) dest += p;
00085 dest -= s; }
00086
00087 template<typename C> static inline void
00088 sub_mod_core (C& dest, const C& s, const C& p, C& carry) {
00089 VERIFY (carry == 0 || carry == 1, "unexpected large carry");
00090 C t = s + carry;
00091 if (t == p || dest < t) { dest += p; carry= 1; } else carry= 0;
00092 dest -= t; }
00093
00094 TMPL static inline void
00095 sub_mod (C& dest, const C& s, const M& m) {
00096 sub_mod_core (dest, s, m.p); }
00097
00098 TMPL static inline void
00099 sub_mod (C& dest, const C& s, const M& m, C& carry) {
00100 sub_mod_core (dest, s, m.p, carry); }
00101
00102 TMPL static inline void
00103 sub_mod (C& dest, const C& s1, const C& s2, const M& m) {
00104 dest = s1;
00105 sub_mod (dest, s2, m); }
00106
00107 TMPL static inline void
00108 sub_mod (C& dest, const C& s1, const C& s2, const M& m, C& carry) {
00109 dest = s1;
00110 sub_mod (dest, s2, m, carry); }
00111 };
00112
00113 template<typename V>
00114 struct modulus_inv_integer_naive : public V {
00115
00116 TMPL static inline bool
00117 is_invertible_mod (const C& s, const M& m) {
00118 return abs (gcd (s, m.p)) == 1; }
00119
00120 TMPL static inline void
00121 inv_mod (C& a, const M& m) {
00122 if (m.p == 0) {
00123 if (a == 1 || a == -1)
00124 return;
00125 else
00126 ERROR ("inv_mod: argument is not invertible");
00127 }
00128 a = invert_modulo (a, m.p);
00129 if (a == 0 && m.p != 1)
00130 ERROR ("inv_mod: argument is not invertible"); }
00131
00132 TMPL static inline void
00133 inv_mod (C& dest, const C& s, const M& m) {
00134 dest = s;
00135 inv_mod (dest, m); }
00136 };
00137
00138 template<typename V>
00139 struct modulus_encoding_integer_naive : public V {
00140 TMPL static inline void
00141 encode_mod (C& dest, const C& s, const M& m) {
00142 if (s < 0 && m.p != 0) {
00143 dest = - s;
00144 V::reduce_mod (dest, m);
00145 dest = m.p - dest;
00146 }
00147 else
00148 dest = s;
00149 V::reduce_mod (dest, m);
00150 }
00151 TMPL static inline void
00152 decode_mod (C& dest, const C& s, const M& m) {
00153 (void) m;
00154 dest = s; }
00155 };
00156
00157 struct modulus_integer_naive:
00158 public modulus_encoding_integer_naive<
00159 modulus_div_naive<
00160 modulus_inv_integer_naive<
00161 modulus_mul_naive<
00162 modulus_add_integer_naive<
00163 modulus_reduction_naive<
00164 modulus_normalization_integer_naive> > > > > > {};
00165
00166
00167
00168
00169
00170 STMPL
00171 struct modulus_variant_helper<integer> {
00172 typedef modulus_integer_naive MV; };
00173
00174
00175
00176
00177
00178 #define TMPLVW template<typename V, typename W>
00179 #define Modular(I) modular<modulus<I,V>,W>
00180
00181 TMPLVW
00182 struct lift_helper<Modular(integer)> {
00183 template<typename R> static inline void
00184 set_op (R&y, const Modular(integer)& x) {
00185 set_as (y, *x); }
00186 static inline integer
00187 op (const Modular(integer)& x) {
00188 return as<integer> (*x); }
00189 };
00190
00191 #undef Modular
00192 #define Modular(I) modular<modulus<I,V>,W>
00193
00194 STMPL
00195 struct project_helper<integer> {
00196 template<typename I, typename V, typename W> static inline void
00197 set_op (Modular(I)& y, const integer& x) {
00198 integer p= *Modular(I)::get_modulus ();
00199 y= (x >= 0) ? Modular(I) (as<I> (rem (x, p)))
00200 : - Modular(I) (as<I> (rem (-x, p))); }
00201 static inline modular<modulus<integer> >
00202 op (const integer& x) {
00203 return modular<modulus<integer> > (x); }
00204 };
00205
00206 #undef Modular
00207 #undef TMPLVW
00208
00209
00210
00211
00212
00213 template<typename C, typename V, typename W>
00214 struct as_helper<modular<modulus<C,V>,W>, rational> {
00215 static inline modular<modulus<C,V>,W>
00216 cv (const rational& a) {
00217 return as<modular<modulus<C,V>,W> > (numerator (a))
00218 / as<modular<modulus<C,V>,W> > (denominator (a)); }
00219 };
00220
00221
00222
00223
00224
00225 #undef TMPL
00226 #define TMPL template<typename V, typename W>
00227 #define Modular modular<modulus<integer,V>,W>
00228 UNARY_RETURN_TYPE(TMPL,reconstruct_op,Modular,rational);
00229
00230 template<typename V, typename W>
00231 rational
00232 reconstruct (const modular<modulus<integer,V>,W>& x) {
00233 integer u= *x;
00234 integer m= *get_modulus (x);
00235 integer N= floor_sqrt (m >> 1);
00236 integer D= N, n, d;
00237 if (reconstruct (n, d, u, m, N, D))
00238 return rational (n) / rational (d);
00239 ERROR ("rational reconstruction failed"); }
00240
00241 template<typename V, typename W>
00242 bool
00243 is_reconstructible (const modular<modulus<integer,V>,W>& x, rational& r) {
00244 integer u= *x;
00245 integer m= *get_modulus (x);
00246 integer N= floor_sqrt (m >> 1);
00247 integer D= N, n, d;
00248 if (!reconstruct (n, d, u, m, N, D)) return false;
00249 r= rational (n) / rational (d);
00250 return true; }
00251
00252 #undef Modular
00253 #undef TMPL
00254 }
00255 #endif //__MMX__MODULAR_INTEGER__HPP