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