| 1 | <html> |
|---|
| 2 | <head> |
|---|
| 3 | <title> |
|---|
| 4 | A Tour of NTL: Examples: Modular Arithmetic </title> |
|---|
| 5 | </head> |
|---|
| 6 | |
|---|
| 7 | <body bgcolor="#fff9e6"> |
|---|
| 8 | <center> |
|---|
| 9 | <a href="tour-ex3.html"><img src="arrow1.gif" alt="[Previous]" align=bottom></a> |
|---|
| 10 | <a href="tour-examples.html"><img src="arrow2.gif" alt="[Up]" align=bottom></a> |
|---|
| 11 | <a href="tour-ex5.html"> <img src="arrow3.gif" alt="[Next]" align=bottom></a> |
|---|
| 12 | </center> |
|---|
| 13 | |
|---|
| 14 | <h1> |
|---|
| 15 | <p align=center> |
|---|
| 16 | A Tour of NTL: Examples: Modular Arithmetic |
|---|
| 17 | </p> |
|---|
| 18 | </h1> |
|---|
| 19 | |
|---|
| 20 | <p> <hr> <p> |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | NTL also supports modular integer arithmetic. |
|---|
| 24 | The class <tt>ZZ_p</tt> |
|---|
| 25 | represents the integers mod <tt>p</tt>. |
|---|
| 26 | Despite the notation, <tt>p</tt> need not in general be prime, |
|---|
| 27 | except in situations where this is mathematically required. |
|---|
| 28 | The classes <tt>vec_ZZ_p</tt>, <tt>mat_ZZ_p</tt>, |
|---|
| 29 | and <tt>ZZ_pX</tt> represent vectors, matrices, and polynomials |
|---|
| 30 | mod <tt>p</tt>, and work much the same way as the corresponding |
|---|
| 31 | classes for <tt>ZZ</tt>. |
|---|
| 32 | |
|---|
| 33 | <p> |
|---|
| 34 | Here is a program that reads a prime number <tt>p</tt>, |
|---|
| 35 | and a polynomial <tt>f</tt> modulo <tt>p</tt>, and factors it. |
|---|
| 36 | |
|---|
| 37 | <p> |
|---|
| 38 | <pre> |
|---|
| 39 | #include <NTL/ZZ_pXFactoring.h> |
|---|
| 40 | |
|---|
| 41 | int main() |
|---|
| 42 | { |
|---|
| 43 | ZZ p; |
|---|
| 44 | cin >> p; |
|---|
| 45 | ZZ_p::init(p); |
|---|
| 46 | |
|---|
| 47 | ZZ_pX f; |
|---|
| 48 | cin >> f; |
|---|
| 49 | |
|---|
| 50 | vec_pair_ZZ_pX_long factors; |
|---|
| 51 | |
|---|
| 52 | CanZass(factors, f); |
|---|
| 53 | |
|---|
| 54 | cout << factors << "\n"; |
|---|
| 55 | // calls "Cantor/Zassenhaus" algorithm |
|---|
| 56 | } |
|---|
| 57 | </pre> |
|---|
| 58 | |
|---|
| 59 | <p> |
|---|
| 60 | As a program is running, NTL keeps track of a "current modulus" |
|---|
| 61 | for the class <tt>ZZ_p</tt>, which can be initialized or changed |
|---|
| 62 | using <tt>ZZ_p::init</tt>. |
|---|
| 63 | This must be done before any variables are declared or |
|---|
| 64 | computations are done that depend on this modulus. |
|---|
| 65 | |
|---|
| 66 | <p> |
|---|
| 67 | Please note that for efficiency reasons, |
|---|
| 68 | NTL does not make any attempt to ensure that |
|---|
| 69 | variables declared under one modulus are not used |
|---|
| 70 | under a different one. |
|---|
| 71 | If that happens, the behavior of a program in this |
|---|
| 72 | case is completely unpredictable. |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | <p> <hr> <p> |
|---|
| 76 | |
|---|
| 77 | Here are two more examples that illustrate the <tt>ZZ_p</tt>-related |
|---|
| 78 | classes. |
|---|
| 79 | The first is a vector addition routine (already supplied by NTL): |
|---|
| 80 | |
|---|
| 81 | <p> |
|---|
| 82 | <pre> |
|---|
| 83 | #include <NTL/vec_ZZ_p.h> |
|---|
| 84 | |
|---|
| 85 | void add(vec_ZZ_p& x, const vec_ZZ_p& a, const vec_ZZ_p& b) |
|---|
| 86 | { |
|---|
| 87 | long n = a.length(); |
|---|
| 88 | if (b.length() != n) Error("vector add: dimension mismatch"); |
|---|
| 89 | |
|---|
| 90 | x.SetLength(n); |
|---|
| 91 | long i; |
|---|
| 92 | for (i = 0; i < n; i++) |
|---|
| 93 | add(x[i], a[i], b[i]); |
|---|
| 94 | } |
|---|
| 95 | </pre> |
|---|
| 96 | |
|---|
| 97 | <p> |
|---|
| 98 | |
|---|
| 99 | The second example is an inner product routine (also supplied by NTL): |
|---|
| 100 | |
|---|
| 101 | <p> |
|---|
| 102 | <pre> |
|---|
| 103 | #include <NTL/vec_ZZ_p.h> |
|---|
| 104 | |
|---|
| 105 | void InnerProduct(ZZ_p& x, const vec_ZZ_p& a, const vec_ZZ_p& b) |
|---|
| 106 | { |
|---|
| 107 | long n = min(a.length(), b.length()); |
|---|
| 108 | long i; |
|---|
| 109 | ZZ accum, t; |
|---|
| 110 | |
|---|
| 111 | accum = 0; |
|---|
| 112 | for (i = 0; i < n; i++) { |
|---|
| 113 | mul(t, rep(a[i]), rep(b[i])); |
|---|
| 114 | add(accum, accum, t); |
|---|
| 115 | } |
|---|
| 116 | |
|---|
| 117 | conv(x, accum); |
|---|
| 118 | } |
|---|
| 119 | </pre> |
|---|
| 120 | |
|---|
| 121 | <p> |
|---|
| 122 | This second example illustrates two things. |
|---|
| 123 | First, it illustrates the use of function <tt>rep</tt> which |
|---|
| 124 | returns a read-only reference to the representation of a <tt>ZZ_p</tt> |
|---|
| 125 | as a <tt>ZZ</tt> between <tt>0</tt> and <tt>p-1</tt>. |
|---|
| 126 | Second, it illustrates a useful algorithmic technique, |
|---|
| 127 | whereby one computes over <tt>ZZ</tt>, reducing mod <tt>p</tt> |
|---|
| 128 | only when necessary. |
|---|
| 129 | This reduces the number of divisions that need to be performed significantly, |
|---|
| 130 | leading to much faster execution. |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | <p> |
|---|
| 134 | The class <tt>ZZ_p</tt> supports all the basic arithmetic |
|---|
| 135 | operations in both operator and procedural form. |
|---|
| 136 | All of the basic operations support a "promotion logic", |
|---|
| 137 | promoting <tt>long</tt> to <tt>ZZ_p</tt>. |
|---|
| 138 | |
|---|
| 139 | <p> |
|---|
| 140 | Note that the class <tt>ZZ_p</tt> is mainly useful only |
|---|
| 141 | when you want to work with vectors, matrices, or polynomials |
|---|
| 142 | mod <tt>p</tt>. |
|---|
| 143 | If you just want to do some simple modular arithemtic, |
|---|
| 144 | it is probably easier to just work with <tt>ZZ</tt>s directly. |
|---|
| 145 | This is especially true if you want to work with many different |
|---|
| 146 | moduli: modulus switching is supported, but it is a bit awkward. |
|---|
| 147 | |
|---|
| 148 | <p> |
|---|
| 149 | The class <tt>ZZ_pX</tt> supports all the basic arithmetic |
|---|
| 150 | operations in both operator and procedural form. |
|---|
| 151 | All of the basic operations support a "promotion logic", |
|---|
| 152 | promoting both <tt>long</tt> and <tt>ZZ_p</tt> to <tt>ZZ_pX</tt>. |
|---|
| 153 | |
|---|
| 154 | <p> |
|---|
| 155 | See <a href="ZZ_p.txt"><tt>ZZ_p.txt</tt></a> for details on <tt>ZZ_p</tt>; |
|---|
| 156 | see <a href="ZZ_pX.txt"><tt>ZZ_pX.txt</tt></a> for details on <tt>ZZ_pX</tt>; |
|---|
| 157 | see <a href="ZZ_pXFactoring.txt"><tt>ZZ_pXFactoring.txt</tt></a> for details on |
|---|
| 158 | the routines for factoring polynomials over <tt>ZZ_p</tt>; |
|---|
| 159 | see <a href="vec_ZZ_p.txt"><tt>vec_ZZ_p.txt</tt></a> for details on <tt>vec_ZZ_p</tt>; |
|---|
| 160 | see <a href="mat_ZZ_p.txt"><tt>mat_ZZ_p.txt</tt></a> for details on <tt>mat_ZZ_p</tt>. |
|---|
| 161 | |
|---|
| 162 | <p> <hr> <p> |
|---|
| 163 | |
|---|
| 164 | There is a mechanism for saving and restoring a modulus, |
|---|
| 165 | which the following example illustrates. |
|---|
| 166 | This routine takes as input an integer polynomial |
|---|
| 167 | and a prime, and tests if the polynomial is irreducible modulo |
|---|
| 168 | the prime. |
|---|
| 169 | |
|---|
| 170 | <p> |
|---|
| 171 | <pre> |
|---|
| 172 | #include <NTL/ZZX.h> |
|---|
| 173 | #include <NTL/ZZ_pXFactoring.h> |
|---|
| 174 | |
|---|
| 175 | long IrredTestMod(const ZZX& f, const ZZ& p) |
|---|
| 176 | { |
|---|
| 177 | ZZ_pBak bak; // save current modulus in bak |
|---|
| 178 | bak.save(); |
|---|
| 179 | |
|---|
| 180 | ZZ_p::init(p); // set the current modulus to p |
|---|
| 181 | |
|---|
| 182 | return DetIrredTest(to_ZZ_pX(f)); |
|---|
| 183 | |
|---|
| 184 | // old modulus is restored automatically when bak is destroyed |
|---|
| 185 | // upon return |
|---|
| 186 | } |
|---|
| 187 | </pre> |
|---|
| 188 | |
|---|
| 189 | <p> |
|---|
| 190 | The modulus switching mechanism is actually quite a bit |
|---|
| 191 | more general and flexible than this example illustrates. |
|---|
| 192 | |
|---|
| 193 | <p> |
|---|
| 194 | The function <tt>to_ZZ_pX</tt> is yet another of NTL's many |
|---|
| 195 | conversion functions. |
|---|
| 196 | We could also have used the equivalent procedural form: |
|---|
| 197 | <pre> |
|---|
| 198 | ZZ_pX f1; |
|---|
| 199 | conv(f1, f); |
|---|
| 200 | return DetIrredTest(f1); |
|---|
| 201 | </pre> |
|---|
| 202 | |
|---|
| 203 | |
|---|
| 204 | |
|---|
| 205 | |
|---|
| 206 | |
|---|
| 207 | |
|---|
| 208 | <p> <hr> <p> |
|---|
| 209 | |
|---|
| 210 | Suppose in the above example that <tt>p</tt> is known in advance |
|---|
| 211 | to be a small, single-precision prime. |
|---|
| 212 | In this case, NTL provides a class <tt>zz_p</tt>, that |
|---|
| 213 | acts just like <tt>ZZ_p</tt>, |
|---|
| 214 | along with corresponding classes <tt>vec_zz_p</tt>, |
|---|
| 215 | <tt>mat_zz_p</tt>, and <tt>zz_pX</tt>. |
|---|
| 216 | The interfaces to all of the routines are generally identical |
|---|
| 217 | to those for <tt>ZZ_p</tt>. |
|---|
| 218 | However, the routines are much more efficient, in both time and space. |
|---|
| 219 | |
|---|
| 220 | <p> |
|---|
| 221 | For small primes, the routine in the previous example could be coded |
|---|
| 222 | as follows. |
|---|
| 223 | |
|---|
| 224 | |
|---|
| 225 | <p> |
|---|
| 226 | <pre> |
|---|
| 227 | #include <NTL/ZZX.h> |
|---|
| 228 | #include <NTL/lzz_pXFactoring.h> |
|---|
| 229 | |
|---|
| 230 | long IrredTestMod(const ZZX& f, long p) |
|---|
| 231 | { |
|---|
| 232 | zz_pBak bak; |
|---|
| 233 | bak.save(); |
|---|
| 234 | |
|---|
| 235 | zz_p::init(p); |
|---|
| 236 | |
|---|
| 237 | return DetIrredTest(to_zz_pX(f)); |
|---|
| 238 | } |
|---|
| 239 | </pre> |
|---|
| 240 | |
|---|
| 241 | <p> <hr> <p> |
|---|
| 242 | |
|---|
| 243 | The following is a routine (essentially the same as implemented in NTL) |
|---|
| 244 | for computing the GCD of polynomials with integer coefficients. |
|---|
| 245 | It uses a "modular" approach: the GCDs are computed modulo small |
|---|
| 246 | primes, and the results are combined using the Chinese Remainder Theorem (CRT). |
|---|
| 247 | The small primes are specially chosen "FFT primes", which are of |
|---|
| 248 | a special form that allows for particular fast polynomial arithmetic. |
|---|
| 249 | |
|---|
| 250 | <p> |
|---|
| 251 | <pre> |
|---|
| 252 | #include <NTL/ZZX.h> |
|---|
| 253 | |
|---|
| 254 | void GCD(ZZX& d, const ZZX& a, const ZZX& b) |
|---|
| 255 | { |
|---|
| 256 | if (a == 0) { |
|---|
| 257 | d = b; |
|---|
| 258 | if (LeadCoeff(d) < 0) negate(d, d); |
|---|
| 259 | return; |
|---|
| 260 | } |
|---|
| 261 | |
|---|
| 262 | if (b == 0) { |
|---|
| 263 | d = a; |
|---|
| 264 | if (LeadCoeff(d) < 0) negate(d, d); |
|---|
| 265 | return; |
|---|
| 266 | } |
|---|
| 267 | |
|---|
| 268 | ZZ c1, c2, c; |
|---|
| 269 | ZZX f1, f2; |
|---|
| 270 | |
|---|
| 271 | content(c1, a); |
|---|
| 272 | divide(f1, a, c1); |
|---|
| 273 | |
|---|
| 274 | content(c2, b); |
|---|
| 275 | divide(f2, b, c2); |
|---|
| 276 | |
|---|
| 277 | GCD(c, c1, c2); |
|---|
| 278 | |
|---|
| 279 | ZZ ld; |
|---|
| 280 | GCD(ld, LeadCoeff(f1), LeadCoeff(f2)); |
|---|
| 281 | |
|---|
| 282 | ZZX g, res; |
|---|
| 283 | |
|---|
| 284 | ZZ prod; |
|---|
| 285 | |
|---|
| 286 | zz_pBak bak; |
|---|
| 287 | bak.save(); |
|---|
| 288 | |
|---|
| 289 | long FirstTime = 1; |
|---|
| 290 | |
|---|
| 291 | long i; |
|---|
| 292 | for (i = 0; ;i++) { |
|---|
| 293 | zz_p::FFTInit(i); |
|---|
| 294 | long p = zz_p::modulus(); |
|---|
| 295 | |
|---|
| 296 | if (divide(LeadCoeff(f1), p) || divide(LeadCoeff(f2), p)) continue; |
|---|
| 297 | |
|---|
| 298 | zz_pX G, F1, F2; |
|---|
| 299 | zz_p LD; |
|---|
| 300 | |
|---|
| 301 | conv(F1, f1); |
|---|
| 302 | conv(F2, f2); |
|---|
| 303 | conv(LD, ld); |
|---|
| 304 | |
|---|
| 305 | GCD(G, F1, F2); |
|---|
| 306 | mul(G, G, LD); |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | if (deg(G) == 0) { |
|---|
| 310 | res = 1; |
|---|
| 311 | break; |
|---|
| 312 | } |
|---|
| 313 | |
|---|
| 314 | if (FirstTime || deg(G) < deg(g)) { |
|---|
| 315 | prod = 1; |
|---|
| 316 | g = 0; |
|---|
| 317 | FirstTime = 0; |
|---|
| 318 | } |
|---|
| 319 | else if (deg(G) > deg(g)) { |
|---|
| 320 | continue; |
|---|
| 321 | } |
|---|
| 322 | |
|---|
| 323 | if (!CRT(g, prod, G)) { |
|---|
| 324 | PrimitivePart(res, g); |
|---|
| 325 | if (divide(f1, res) && divide(f2, res)) |
|---|
| 326 | break; |
|---|
| 327 | } |
|---|
| 328 | |
|---|
| 329 | } |
|---|
| 330 | |
|---|
| 331 | mul(d, res, c); |
|---|
| 332 | if (LeadCoeff(d) < 0) negate(d, d); |
|---|
| 333 | } |
|---|
| 334 | </pre> |
|---|
| 335 | |
|---|
| 336 | |
|---|
| 337 | <p> |
|---|
| 338 | See <a href="lzz_p.txt"><tt>lzz_p.txt</tt></a> for details on <tt>zz_p</tt>; |
|---|
| 339 | see <a href="lzz_pX.txt"><tt>lzz_pX.txt</tt></a> for details on <tt>zz_pX</tt>; |
|---|
| 340 | see <a href="lzz_pXFactoring.txt"><tt>lzz_pXFactoring.txt</tt></a> for details on |
|---|
| 341 | the routines for factoring polynomials over <tt>zz_p</tt>; |
|---|
| 342 | see <a href="vec_lzz_p.txt"><tt>vec_lzz_p.txt</tt></a> for details on <tt>vec_zz_p</tt>; |
|---|
| 343 | see <a href="mat_lzz_p.txt"><tt>mat_lzz_p.txt</tt></a> for details on <tt>mat_zz_p</tt>. |
|---|
| 344 | |
|---|
| 345 | |
|---|
| 346 | |
|---|
| 347 | <p> <hr> <p> |
|---|
| 348 | |
|---|
| 349 | Arithmetic mod 2 is such an important special case that NTL |
|---|
| 350 | provides a class <tt>GF2</tt>, that |
|---|
| 351 | acts just like <tt>ZZ_p</tt> when <tt>p == 2</tt>, |
|---|
| 352 | along with corresponding classes <tt>vec_GF2</tt>, |
|---|
| 353 | <tt>mat_GF2</tt>, and <tt>GF2X</tt>. |
|---|
| 354 | The interfaces to all of the routines are generally identical |
|---|
| 355 | to those for <tt>ZZ_p</tt>. |
|---|
| 356 | However, the routines are much more efficient, in both time and space. |
|---|
| 357 | |
|---|
| 358 | <p> |
|---|
| 359 | |
|---|
| 360 | This example illustrates the <tt>GF2X</tt> and <tt>mat_GF2</tt> |
|---|
| 361 | classes with a simple routine to test if a polynomial over GF(2) |
|---|
| 362 | is irreducible using linear algebra. |
|---|
| 363 | NTL's built-in irreducibility test is to be preferred, however. |
|---|
| 364 | |
|---|
| 365 | <pre> |
|---|
| 366 | |
|---|
| 367 | #include <NTL/GF2X.h> |
|---|
| 368 | #include <NTL/mat_GF2.h> |
|---|
| 369 | |
|---|
| 370 | long MatIrredTest(const GF2X& f) |
|---|
| 371 | { |
|---|
| 372 | long n = deg(f); |
|---|
| 373 | |
|---|
| 374 | if (n <= 0) return 0; |
|---|
| 375 | if (n == 1) return 1; |
|---|
| 376 | |
|---|
| 377 | if (GCD(f, diff(f)) != 1) return 0; |
|---|
| 378 | |
|---|
| 379 | mat_GF2 M; |
|---|
| 380 | |
|---|
| 381 | M.SetDims(n, n); |
|---|
| 382 | |
|---|
| 383 | GF2X x_squared = GF2X(2, 1); |
|---|
| 384 | |
|---|
| 385 | GF2X g; |
|---|
| 386 | g = 1; |
|---|
| 387 | |
|---|
| 388 | for (long i = 0; i < n; i++) { |
|---|
| 389 | VectorCopy(M[i], g, n); |
|---|
| 390 | M[i][i] += 1; |
|---|
| 391 | g = (g * x_squared) % f; |
|---|
| 392 | } |
|---|
| 393 | |
|---|
| 394 | long rank = gauss(M); |
|---|
| 395 | |
|---|
| 396 | if (rank == n-1) |
|---|
| 397 | return 1; |
|---|
| 398 | else |
|---|
| 399 | return 0; |
|---|
| 400 | } |
|---|
| 401 | </pre> |
|---|
| 402 | |
|---|
| 403 | <p> |
|---|
| 404 | Note that the statement |
|---|
| 405 | <pre> |
|---|
| 406 | g = (g * x_squared) % f; |
|---|
| 407 | </pre> |
|---|
| 408 | could be replace d by the more efficient code sequence |
|---|
| 409 | <pre> |
|---|
| 410 | MulByXMod(g, g, f); |
|---|
| 411 | MulByXMod(g, g, f); |
|---|
| 412 | </pre> |
|---|
| 413 | but this would not significantly impact the overall |
|---|
| 414 | running time, since it is the Gaussian elimination that |
|---|
| 415 | dominates the running time. |
|---|
| 416 | |
|---|
| 417 | <p> |
|---|
| 418 | See <a href="GF2.txt"><tt>GF2.txt</tt></a> for details on <tt>GF2</tt>; |
|---|
| 419 | see <a href="GF2X.txt"><tt>GF2X.txt</tt></a> for details on <tt>GF2X</tt>; |
|---|
| 420 | see <a href="GF2XFactoring.txt"><tt>GF2XFactoring.txt</tt></a> for details on |
|---|
| 421 | the routines for factoring polynomials over <tt>GF2</tt>; |
|---|
| 422 | see <a href="vec_GF2.txt"><tt>vec_GF2.txt</tt></a> for details on <tt>vec_GF2</tt>; |
|---|
| 423 | see <a href="mat_GF2.txt"><tt>mat_GF2.txt</tt></a> for details on <tt>mat_GF2</tt>. |
|---|
| 424 | |
|---|
| 425 | <p> |
|---|
| 426 | |
|---|
| 427 | <center> |
|---|
| 428 | <a href="tour-ex3.html"><img src="arrow1.gif" alt="[Previous]" align=bottom></a> |
|---|
| 429 | <a href="tour-examples.html"><img src="arrow2.gif" alt="[Up]" align=bottom></a> |
|---|
| 430 | <a href="tour-ex5.html"> <img src="arrow3.gif" alt="[Next]" align=bottom></a> |
|---|
| 431 | </center> |
|---|
| 432 | |
|---|
| 433 | </body> |
|---|
| 434 | </html> |
|---|
