00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_MODULAR_POLYNOMIAL_HPP
00014 #define __MMX_MODULAR_POLYNOMIAL_HPP
00015 #include <numerix/modular.hpp>
00016 #include <algebramix/polynomial.hpp>
00017 #include <algebramix/quotient_polynomial.hpp>
00018
00019 namespace mmx {
00020
00021 #define TMPL template<typename C, typename M>
00022 #define Polynomial polynomial<C,V>
00023
00024
00025
00026
00027
00028 template<typename V>
00029 struct modulus_polynomial_inv_naive : public V {
00030
00031 TMPL static inline bool
00032 is_invertible_mod (const C& s, const M& m) {
00033 return abs (gcd (s, m.p)) == 1; }
00034
00035 TMPL static inline void
00036 inv_mod (C& a, const M& m) {
00037 if (m.p == 0) {
00038 if (deg (a) == 0)
00039 a= 1 / a [0];
00040 else
00041 ERROR ("inv_mod: argument is not invertible");
00042 }
00043 a = invert_modulo (a, m.p);
00044 if (a == 0)
00045 ERROR ("inv_mod: argument is not invertible"); }
00046
00047 TMPL static inline void
00048 inv_mod (C& dest, const C& s, const M& m) {
00049 dest = s;
00050 inv_mod (dest, m); }
00051 };
00052
00053 struct modulus_polynomial_naive:
00054 public modulus_encoding_naive<
00055 modulus_div_naive<
00056 modulus_polynomial_inv_naive<
00057 modulus_mul_naive<
00058 modulus_add_naive<
00059 modulus_reduction_naive<
00060 modulus_normalization_naive> > > > > > {};
00061
00062
00063
00064
00065
00066 template<typename C, typename V>
00067 struct modulus_variant_helper<polynomial<C,V> > {
00068 typedef modulus_polynomial_naive MV; };
00069
00070
00071
00072
00073
00074 template<typename C, typename PV, typename MV, typename MW> inline syntactic
00075 flatten (const modular<modulus<polynomial<C, PV>, MW>, MV>& c) {
00076 generic x= polynomial<C, PV>::get_variable_name ();
00077 generic a= gen (GEN_PRIME, x);
00078 if (x == "x") a= "a";
00079 if (x == "y") a= "b";
00080 if (x == "z") a= "c";
00081 return flatten (*c, as_syntactic (a));
00082 }
00083
00084 template<typename C, typename PV, typename MV> inline syntactic
00085 flatten (const modulus<polynomial<C, PV>, MV>& c) {
00086 generic x= polynomial<C, PV>::get_variable_name ();
00087 generic a= gen (GEN_PRIME, x);
00088 if (x == "x") a= "a";
00089 if (x == "y") a= "b";
00090 if (x == "z") a= "c";
00091 return flatten (*c, as_syntactic (a));
00092 }
00093
00094
00095
00096
00097
00098 template<typename X>
00099 struct modulus_polynomial_reduction_preinverse : public X {
00100 template<typename C,typename V, typename M> static inline void
00101 reduce_mod (Polynomial& dest, const M& m) {
00102 if (is_exact_zero (m.p)) return;
00103 int d= degree (m.p);
00104 if (deg (dest) < (d << 1)) {
00105 Polynomial carry;
00106 carry = rshiftz (rshiftz (dest, d) * m.q, d);
00107 dest = range (dest - carry * m.p, 0, d);
00108 }
00109 else dest= rem (dest, m.p); }
00110
00111 template<typename C,typename V, typename M> static inline void
00112 reduce_mod (Polynomial& dest, const M& m, Polynomial& carry) {
00113 int d= degree (m.p);
00114 if (is_exact_zero (m.p) != 0)
00115 carry= 0;
00116 else if (deg (dest) < (d << 1)) {
00117 carry = rshiftz (rshiftz (dest, d) * m.q, d);
00118 dest = range (dest - carry * m.p, 0, d);
00119 }
00120 else {
00121 carry= quo (dest, m.p);
00122 dest = rem (dest, m.p); } }
00123
00124 template<typename C,typename V, typename M> static inline void
00125 reduce_mod (Polynomial& dest, const Polynomial& s, const M& m) {
00126 dest= s; reduce_mod (dest, m); }
00127
00128 template<typename C,typename V, typename M> static inline void
00129 reduce_mod (Polynomial& dest, const Polynomial& s,
00130 const M& m, Polynomial& carry) {
00131 dest= s; reduce_mod (dest, m, carry); }
00132 };
00133
00134 template<typename X,typename W=modulus_mul_naive<X> >
00135 struct modulus_polynomial_mul_preinverse : public X {
00136
00137 template<typename C, typename V, typename M> static inline void
00138 mul_mod (Polynomial& dest, const Polynomial& s, const M& m) {
00139 if (is_exact_zero (*m)) dest *= s;
00140 else {
00141 int d= degree (*m);
00142 Polynomial a= dest * s;
00143 Polynomial b= rshiftz (rshiftz (a, d) * m.q, d);
00144 dest= a - b * (* m);
00145 }
00146 }
00147
00148 template<typename C, typename V, typename M> static inline void
00149 mul_mod (Polynomial& dest, const Polynomial& s, const M& m,
00150 Polynomial& carry) {
00151 if (is_exact_zero (*m)) {
00152 dest = dest * s + carry;
00153 carry = 0;
00154 }
00155 else {
00156 int d= degree (*m);
00157 Polynomial a= dest * s + carry;
00158 carry= rshiftz (rshiftz (a, d) * m.q, d);
00159 dest= a - carry * (* m);
00160 }
00161 }
00162
00163 template<typename C, typename V, typename M> static inline void
00164 mul_mod (Polynomial& dest, const Polynomial& s1, const Polynomial& s2,
00165 const M& m) {
00166 dest = s1;
00167 mul_mod (dest, s2, m); }
00168
00169 template<typename C, typename V, typename M> static inline void
00170 mul_mod (Polynomial& dest, const Polynomial& s1, const Polynomial& s2,
00171 const M& m, Polynomial& carry) {
00172 dest = s1;
00173 mul_mod (dest, s2, m, carry); }
00174
00175 };
00176
00177 template<typename V=modulus_polynomial_naive>
00178 struct modulus_polynomial_preinverse:
00179 public modulus_encoding_naive<
00180 modulus_div_naive<
00181 modulus_polynomial_inv_naive<
00182 modulus_polynomial_mul_preinverse<
00183 modulus_add_naive<
00184 modulus_polynomial_reduction_preinverse<
00185 modulus_normalization_naive> > > > > > {};
00186
00187 template<typename C, typename V, typename W>
00188 class modulus <Polynomial, modulus_polynomial_preinverse<W> > {
00189 MMX_ALLOCATORS
00190 public:
00191 Polynomial p;
00192 Polynomial q;
00193
00194 typedef Polynomial base;
00195 typedef modulus_polynomial_preinverse<W> variant;
00196 typedef modulus <Polynomial, modulus_polynomial_preinverse<W> > self_type;
00197
00198 inline Polynomial operator * () const { return p; }
00199
00200 inline modulus () : p (0), q (0) {}
00201
00202 inline modulus (const Polynomial& p2): p (p2) {
00203 q= is_exact_zero (p) ? p : invert_hi (p); }
00204
00205 inline modulus (const self_type& x) : p (x.p), q (x.q) {}
00206
00207 inline self_type&
00208 operator = (const self_type& x) {
00209 p = x.p;
00210 q = x.q;
00211 return *this; }
00212
00213 inline bool operator == (const self_type& a) const {
00214 return p == a.p; }
00215
00216 inline bool operator != (const self_type& a) const {
00217 return p != a.p; }
00218 };
00219
00220
00221
00222
00223
00224 template<typename X,typename W=modulus_mul_naive<X> >
00225 struct modulus_polynomial_mul_power_of_the_variable : public X {
00226
00227 template<typename C, typename V, typename M> static inline void
00228 mul_mod (Polynomial& dest, const Polynomial& s, const M& m) {
00229 if (is_exact_zero (*m)) dest *= s;
00230 else {
00231 int d= degree (*m);
00232 Polynomial a= dest * s;
00233 dest = range (a, 0, d);
00234 }
00235 }
00236
00237 template<typename C, typename V, typename M> static inline void
00238 mul_mod (Polynomial& dest, const Polynomial& s, const M& m,
00239 Polynomial& carry) {
00240 if (is_exact_zero (*m)) {
00241 dest = dest * s + carry;
00242 carry = 0;
00243 }
00244 else {
00245 int d= degree (*m);
00246 Polynomial a= dest * s + carry;
00247 carry= range (a, d, 2*d-1);
00248 dest= range (a, 0, d);
00249 }
00250 }
00251
00252 template<typename C, typename V, typename M> static inline void
00253 mul_mod (Polynomial& dest, const Polynomial& s1, const Polynomial& s2,
00254 const M& m) {
00255 dest = s1;
00256 mul_mod (dest, s2, m); }
00257
00258 template<typename C, typename V, typename M> static inline void
00259 mul_mod (Polynomial& dest, const Polynomial& s1, const Polynomial& s2,
00260 const M& m, Polynomial& carry) {
00261 dest = s1;
00262 mul_mod (dest, s2, m, carry); }
00263
00264 };
00265
00266 template<typename V=modulus_polynomial_naive>
00267 struct modulus_polynomial_power_of_the_variable:
00268 public modulus_encoding_naive<
00269 modulus_div_naive<
00270 modulus_polynomial_inv_naive<
00271 modulus_polynomial_mul_power_of_the_variable<
00272 modulus_add_naive<
00273 modulus_reduction_naive<
00274 modulus_normalization_naive> > > > > > {};
00275
00276 template<typename C, typename V, typename W>
00277 class modulus <Polynomial, modulus_polynomial_power_of_the_variable<W> > {
00278 MMX_ALLOCATORS
00279 public:
00280 Polynomial p;
00281
00282 typedef Polynomial base;
00283 typedef modulus_polynomial_power_of_the_variable<W> variant;
00284 typedef modulus <Polynomial, modulus_polynomial_power_of_the_variable<W> >
00285 self_type;
00286
00287 inline Polynomial operator * () const { return p; }
00288
00289 inline modulus () : p (0) {}
00290
00291 inline modulus (const Polynomial& p2): p (p2) {
00292 if (is_exact_zero (p)) return;
00293 ASSERT (is_exact_zero (range (p, 0, deg (p))),
00294 "Modulus must be a power of the variable");
00295 }
00296
00297 inline modulus (const self_type& x) : p (x.p) {}
00298
00299 inline self_type&
00300 operator = (const self_type& x) {
00301 p = x.p;
00302 return *this; }
00303
00304 inline bool operator == (const self_type& a) const {
00305 return p == a.p; }
00306
00307 inline bool operator != (const self_type& a) const {
00308 return p != a.p; }
00309 };
00310
00311
00312
00313
00314
00315 #undef TMPL
00316 #define TMPL template<typename C,typename V,typename MoV,typename MaV>
00317 #define Quotient quotient<polynomial<C,V>,polynomial<C,V> >
00318 #define Modular modular<modulus<polynomial<C,V>, MoV>,MaV>
00319
00320 TMPL
00321 struct as_helper<Modular,Quotient> {
00322 static Modular cv (const Quotient& x) {
00323 return as<Modular> (numerator (x)) / as<Modular> (denominator (x)); }
00324 };
00325
00326
00327
00328
00329
00330 UNARY_RETURN_TYPE(TMPL,reconstruct_op,Modular,Quotient);
00331
00332 TMPL Quotient
00333 reconstruct (const modular<modulus<polynomial<C,V>, MoV>,MaV>& x) {
00334 polynomial<C,V> num, den, q= *get_modulus (x);
00335 nat k= N(q) >> 1;
00336 bool b= reconstruct (*x, q, k, num, den);
00337 ASSERT (b, "rational reconstruction failed");
00338 return Quotient (num, den); }
00339
00340 TMPL bool
00341 is_reconstructible (const modular<modulus<polynomial<C,V>, MoV>,MaV>& x,
00342 Quotient& y) {
00343 polynomial<C,V> num, den, q= *get_modulus (x);
00344 nat k= N(q) >> 1;
00345 if (!reconstruct (*x, q, k, num, den)) return false;
00346 y= Quotient (num, den);
00347 return true; }
00348
00349 #undef Quotient
00350 #undef Modular
00351
00352 #undef TMPL
00353 #undef Polynomial
00354 }
00355 #endif // __MMX_MODULAR_POLYNOMIAL_HPP
