Bug fix: gcfactor(x*%i) => lisp error
[maxima.git] / tests / rtest_numth.mac
blob0f6536399bbfbf90b88d363dbbbcd90b6a65979f
1 (kill(all), 0);
2 0;
4 /* (Z/pZ)* p prime */
5 (p : 2^61-1, fs : ifactors(p - 1), g : zn_primroot(p, fs));
6 37;
8 zn_primroot_p(power_mod(g, 7, p), p, fs);
9 false;
11 zn_primroot_p(power_mod(g, 17, p), p, fs);
12 true;
14 is(zn_order(g, p, fs) = totient(p));
15 true;
17 is(zn_order(power_mod(g, 7, p), p, fs) = zn_order(g, p, fs));
18 false;
20 is(zn_order(power_mod(g, 17, p), p, fs) = zn_order(g, p, fs));
21 true;
23 (a : power_mod(g, 1234567890, p), zn_log(a, g, p, fs));
24 1234567890;
26 (g : 3, ord_3 : zn_order(3,p,fs), fs_ord_3 : ifactors(ord_3), (p-1)/ord_3);
29 (a : power_mod(g, 1234567890, p), zn_log(a, g, p, fs));
30 1234567890;
32 /* (Z/nZ)* n composite */
33 n : 22;
34 22;
36 g : zn_primroot(n);
39 zn_primroot_p(power_mod(g, 2, n), n);
40 false;
42 zn_primroot_p(power_mod(g, 3, n), n);
43 true;
45 zn_order(power_mod(g, 2, n), n);
48 zn_order(g, n);
49 10;
51 a : power_mod(g, 8, n);
54 zn_log(a, g, n);
57 (g : 5, a : power_mod(g, 4, n), zn_log(a, g, n));
60 /* CRT */
61 mods : [1009, 1013, 1019];
62 [1009, 1013, 1019];
64 x : 374599943;
65 374599943;
67 rems : map(lambda([z], mod(x, z)), mods);
68 [621, 647, 258];
70 chinese(rems, mods);
71 374599943;
73 (remvalue(p,fs,g,a,n,mods,x,rems), 0);
76 (kill(all), modulus_copy : modulus, modulus : false, 
77 gf_coeff_limit_copy : gf_coeff_limit, gf_coeff_limit : 256, 0); 
80 /* Tests adapted from contrib/gf/gf_test.mac */
82 ( gf_set_data(123127, x^5+2*x+1), gf_infolist() );
83 [123127,x^5+2*x+1,x+4,28298700309333316062584407,28298700309333316062584406];
85 ( gf_set_data(7, x^10+5*x^2+x+5), gf_infolist() );
86 [7,x^10+5*x^2+x+5,x,282475249,282475248];
88 gf_exp(gf_primitive(), gf_index(x^9+3*x^6+x^5+2*x^2+6));
89 x^9+3*x^6+x^5+2*x^2+6;
91 gf_minimal_poly(x^9+3*x^6+x^5+2*x^2+6);
92 z^10+5*z^9+6*z^8+5*z^6+3*z^5+4*z^4+z^3+2*z^2+4*z+3;
94 ( gf_set_data(19), gf_infolist() );
95 [19,x,2,19,18];
97 gf_exp(gf_primitive(), gf_index(17));
98 17;
100 ( gf_set_data(10000000019), gf_infolist() );
101 [10000000019,x,2,10000000019,10000000018];
103 gf_exp(gf_primitive() ,gf_index(3));
106 ( gf_set_data(32717, x^11+x^5+x^2+x+1), gf_infolist() );
107 [32717,x^11+x^5+x^2+x+1,x+2,
108         45973568360012658888852552517205008541393124962933,
109         45973568360012658888852552517205008541393124962932];
111 ( gf_set_data(211, x^17+2*x^2+1), gf_infolist() );
112 [211,x^17+2*x^2+1,x+6,3256879871129217157345956711624826081171,
113         3256879871129217157345956711624826081170];
115 ( gf_set_data(2, x^20+x^3+x^2+x+1), gf_infolist() );
116 [2,x^20+x^3+x^2+x+1,x^2+x,1048576,1048575];
118 gf_exp(gf_primitive(), gf_index(x^10+1));
119 x^10+1;
121 gf_minimal_poly(x^10+1);
122 z^20+z^13+z^12+z^5+z^4+z^3+1;
124 ( gf_set_data(3, x^91+x^35+x+1), gf_infolist() );
125 [3,x^91+x^35+x+1,x,26183890704263137277674192438430182020124347,
126         26183890704263137277674192438430182020124346];
128 /* Tests adapted from contrib/gf/gf_hard_test.mac */
130 ( gf_set_data(7, x^61+x^4+1), gf_infolist() );
131 [7,x^61+x^4+1,x+3,3556153025177363557255317383565515512407041673852007,
132         3556153025177363557255317383565515512407041673852006];
134 ( gf_set_data(197, x^24-x^8+2), gf_infolist() );
135 [197,x^24+196*x^8+2,x+19,
136         11673186598630578538556565100133681446610566511878526881,
137         11673186598630578538556565100133681446610566511878526880];
139 ( gf_set_data(5, x^84+x^41+x^2+1), gf_infolist() );
140 [5,x^84+x^41+x^2+1,x^2+1,
141         51698788284564229679463043254372678347863256931304931640625,
142         51698788284564229679463043254372678347863256931304931640624];
144 ( gf_set_data(2, x^102+x^29+1), gf_infolist() );
145 [2,x^102+x^29+1,x+1,5070602400912917605986812821504,
146         5070602400912917605986812821503];
148 ( gf_set_data(5, x^61+x^15+1), gf_infolist() );
149 [5,x^61+x^15+1,x+4,4336808689942017736029811203479766845703125,
150         4336808689942017736029811203479766845703124];
152 ( gf_set_data(8796519617, x^8+3*x^6+x+1), gf_infolist() );
153 [8796519617,x^8+3*x^6+x+1,x+9,
154         35849822058178726610670969179311817327626124038937602048832281182665519944803841,
155         35849822058178726610670969179311817327626124038937602048832281182665519944803840];
157 ( gf_set_data(3, x^120+x^4-1), gf_infolist() );
158 [3,x^120+x^4+2,x^3+x^2+2,
159         1797010299914431210413179829509605039731475627537851106401,
160         1797010299914431210413179829509605039731475627537851106400];
162 ( gf_set_data(18659817111137), gf_infolist() );
163 [18659817111137,x,3,18659817111137,18659817111136];
165 gf_log(7);
166 15962290024269;
168 /* Examples adapted from contrib/gf/gf_manual.pdf */
170 ( gf_set_data(2, x^4+x+1), gf_infolist() );
171 [2,x^4+x+1,x,16,15];
173 (a : x^3+x^2+1, b : x^3+x+1, 0);
176 gf_add(a, b);
177 x^2+x;
179 gf_mult(a, b);
180 x^2+x;
182 gf_inv(b);
183 x^2+1;
185 gf_div(a, b);
186 x^3+x^2;
188 gf_mult(a, gf_inv(b));
189 x^3+x^2;
191 gf_exp(a, 10);
192 x^2+x+1;
194 gf_exp(a, 15);
197 gf_primitive();
200 gf_index(a);
203 ev(a = gf_exp(x, 13)), pred;
204 true;
206 (gf_make_logs(), 0);
209 gf_logs[10];
212 gf_n2p(10);
213 x^3+x;
215 gf_index(x^3+x);
218 (a : x^2+x+1, b : x^3+x^2+1, 0);
221 gf_log(a, b);
224 gf_primitive_p(x^3+x+1);
225 true;
227 gf_primitive_p(x^2+x);
228 false;
230 gf_order(x^2+x);
233 gf_order(x^3+x+1);
236 (a : x^3+x+1, 0);
239 p : gf_minimal_poly(a);
240 z^4+z^3+1;
242 p : subst(a, z, p);
243 (x^3+x+1)^4+(x^3+x+1)^3+1;
245 gf_eval(p);
248 ( gf_set_data(7, x^4+3*x^3+5*x^2+6*x+2), gf_infolist() );
249 [7,x^4+3*x^3+5*x^2+6*x+2,x+4,2401,2400];
251 (char : 7, exp : 4, g : 3*x^3+3*x^2+6, 0);
254 gf_primitive_p(g);
255 true;
257 a : makelist(gf_log(x+i, g),i, 0, 6);
258 [1028,1935,2054,1008,379,1780,223];
260 d : 1702;
261 1702;
263 ord : char^exp - 1;
264 2400;
266 c : makelist(mod(a[i] + d, ord), i, 1, 7);
267 [330,1237,1356,310,2081,1082,1925];
269 m : [1,0,1,1,0,0,1];
270 [1,0,1,1,0,0,1];
272 c : mod(sum(m[i] * c[i], i, 1, 7), ord);
273 1521;
275 r : mod(c - exp*d, ord);
276 1913;
278 u : gf_exp(g, r);
279 x^3+3*x^2+2*x+5;
281 s : u + gf_reduction();
282 x^4+4*x^3+8*x^2+8*x+7;
284 gf_factor(s);
285 x*(x+2)*(x+3)*(x+6);
287 ( gf_set_data(2, x^4+x+1), gf_infolist() );
288 [2,x^4+x+1,x,16,15];
290 (m : matrix([x^3+x^2+x, x^3, x^3+x^2], [x^2, x^3+x^2+1, x^3+x+1], [x^2+x, x^3+x^2+x+1, x^2]), 0);
293 mi : gf_invert_by_lu(m);
294 matrix([x^2, x^3+x^2+x+1, x^3], [x^3+x+1, x^2+x, x^3+1], [x, 0, x^3+x+1]);
296 gf_matmult(m, mi);
297 matrix([1,0,0], [0,1,0], [0,0,1]);
299 ( gf_set_data(2, x^10+x^3+1), gf_infolist() );
300 [2,x^10+x^3+1,x,1024,1023];
302 elem : gf_normal();
303 x^7;
305 m : gf_normal_basis(elem);
306 matrix([0,0,0,1,1,0,1,1,0,0],[0,0,1,1,0,1,1,0,0,0],
307        [1,1,1,1,1,1,1,1,1,1],[0,0,0,1,1,0,0,0,0,0],
308        [0,0,0,0,1,1,0,0,0,0],[0,1,1,1,0,1,1,1,0,0],
309        [0,0,0,0,0,1,1,0,0,0],[0,0,0,0,1,0,1,0,1,1],
310        [0,0,0,0,1,1,0,1,1,0],[0,0,0,0,0,1,0,0,0,0]);
312 gf_normal_basis_rep(elem, mi : gf_invert_by_lu(m));
313 [1,0,0,0,0,0,0,0,0,0];
315 (a : x^9+x^7+x^6+x^5+x^3+x^2+x, 0);
318 gf_normal_basis_rep(a, mi);
319 [1,1,1,0,1,0,1,1,1,0];
321 gf_normal_basis_rep(gf_exp(a, 2), mi);
322 [0,1,1,1,0,1,0,1,1,1];
324 ( gf_set_data(2, x^20+x^3+1), gf_infolist() );
325 [2,x^20+x^3+1,x,1048576,1048575];
327 (a : x^15+x^5+1, 0);
330 gf_index(a);
331 720548;
333 gf_exp(a, 3^12);
334 x^17+x^16+x^13+x^12+x^11+x^3+x^2+x;
336 /* some new tests */
338 ( gf_set_data(2, x^12+x^3+1), gf_infolist() );
339 [2,x^12+x^3+1,x+1,4096,4095];
341 gf_log(gf_exp(x+1, 1234));
342 1234;
344 ( gf_set_data(8796519617, x^2+3), gf_infolist() );
345 [8796519617,x^2+3,x+9,77378757372265826689,77378757372265826688];
347 gf_log(gf_exp(x+9, 1234567890));
348 1234567890;
350 ( gf_set_data(2, 4), gf_infolist() );
351 [2,x^4+x+1,x,16,15];
353 ( prod : (z - 0), 
354   for i:1 thru 15 do prod : prod*(z - gf_exp(x,i)), 
355   block([modulus:2], polymod(remainder(prod, x^4+x+1))) );
356 z^16+z;
358 fs : gf_factor(x^16-x, 2);
359 x*(x+1)*(x^2+x+1)*(x^4+x+1)*(x^4+x^3+1)*(x^4+x^3+x^2+x+1);
361 (gf_minimal_set(2, x^17), 0);
364 apply(gf_mult, args(fs));
365 x^16+x;
367 ( gf_set_data(13, x^7+3*x+2), gf_infolist() );
368 [13,x^7+3*x+2,x,62748517,62748516];
370 p : 9*x^6+12*x^5+5*x^4+x^3+8*x^2+2*x;
371 9*x^6+12*x^5+5*x^4+x^3+8*x^2+2*x;
373 gf_normal_p(p);
374 true;
376 m : gf_normal_basis(p);
377 matrix([9,7,10,4,4,2,3],[12,1,8,5,9,2,2],[5,12,8,9,10,3,5],
378    [1,5,6,8,11,8,0],[8,4,11,6,12,6,5],[2,2,12,11,1,5,6],[0,6,10,2,2,8,5]);
380 gf_normal_basis_rep(p, mi : gf_invert_by_lu(m));
381 [1,0,0,0,0,0,0];
383 coeffs : gf_normal_basis_rep(x^2+4*x+7, mi);
384 [8,8,7,2,5,1,2];
386 ( basis : map(gf_l2p, args(transpose(m))), 
387   apply(gf_add, map(gf_mult, coeffs, basis)) );
388 x^2+4*x+7;
390 ( gf_set_data(7, 4), gf_infolist() );
391 [7,x^4+x+1,x+5,2401,2400];
393 ( gf_set_data(7, x^4+x^2+3), gf_infolist() );
394 [7,x^4+x^2+3,x+1,2401,2400];
396 ( gf_set_data(2, x^8+1), gf_infolist() );
397 [2,x^8+1,unknown,256,128];
399 gf_order();
400 128;
402 gf_inv(x);
403 x^7;
405 gf_inv(x+1);
406 false;
408 gf_gcdex(x+1, gf_reduction());
409 [1, 0, x+1];
411 ( gf_set_data(3, x^8+1), gf_infolist() );
412 [3,x^8+1,unknown,6561,6400];
414 gf_order();
415 6400;
417 ( gf_set_data(3, x^8+x+1), gf_infolist() );
418 [3,x^8+x+1,unknown,6561,4368];
420 gf_order();
421 4368;
423 ( gf_set_data(13, x^8+2), gf_infolist() );
424 [13,x^8+2,x+2,815730721,815730720];
426 (a : x^7+x+1, b : x^3+3*x^2+9*x+7, 0);
429 gf_gcd(a, b);
430 x^2+2*x+7;
432 gf_gcdex(a, b);
433 [7, 6*x^4+8*x^3+3*x, x^2+2*x+7];
435 gf_primitive_poly(7, 8);
436 x^8+x+3;
438 gf_primitive_poly_p(x^8+x+3, 7);
439 true;
441 gf_primitive_poly_p(gf_irreducible(7, 8), 7);
442 true;
444 gf_primitive_poly(2, 8);
445 x^8+x^4+x^3+x^2+1;
447 gf_primitive_poly_p(x^8+x^4+x^3+x^2+1, 2);
448 true;
450 gf_primitive_poly_p(gf_irreducible(2, 8), 2);
451 false;
453 block([count:0, end:2*3^4], 
454   for n:3^4 thru end do 
455     if gf_primitive_poly_p(p:gf_n2p(n, 3), 3) then count:count+1, 
456   count);
459 /* lu decomposition: */
461 ( gf_set_data(2, x^8+x^4+x^3+x+1), gf_infolist() );
462 [2,x^8+x^4+x^3+x+1,x+1,256,255];
464 m : matrix([2,3,1,1], [1,2,3,1], [1,1,2,3], [3,1,1,2]);
465 matrix([2,3,1,1], [1,2,3,1], [1,1,2,3], [3,1,1,2]);
467 mp : matrixmap(gf_n2p, m);
468 matrix([x,x+1,1,1], [1,x,x+1,1], [1,1,x,x+1], [x+1,1,1,x]);
470 mpi : gf_invert_by_lu(mp);
471 matrix([x^3+x^2+x, x^3+x+1,   x^3+x^2+1, x^3+1    ],
472        [x^3+1,     x^3+x^2+x, x^3+x+1,   x^3+x^2+1],
473        [x^3+x^2+1, x^3+1,     x^3+x^2+x, x^3+x+1  ],
474        [x^3+x+1,   x^3+x^2+1, x^3+1,     x^3+x^2+x]);
476 mi : matrixmap(gf_p2n, mpi);
477 matrix([14,11,13,9], [9,14,11,13], [13,9,14,11], [11,13,9,14]);
479 ( ef_set_data(x), ef_infolist() );
480 [x,3,256,255];
482 is( ef_invert_by_lu(m) = mi );
483 true;
485 ef_unset();
486 true;
488 ( p : 17, n : 4, gf_set_data(p,n), 0);
491 (elem : 6*x^3+13*x^2+4*x+2, gf_normal_p(elem));
492 true;
494 m : gf_normal_basis(elem);
495 matrix([6,7,11,10],[13,4,13,4],[4,1,13,16],[2,2,2,2]);
497 mi : gf_invert_by_lu(m);
498 matrix([5,1,16,15],[14,16,13,15],[12,1,1,15],[3,16,4,15]);
500 is(gf_matmult(m, mi) = diagmatrix(n, 1));
501 true;
503 is(zn_invert_by_lu(m, p) = mi);
504 true;
506 mod(determinant(m), p);
509 mod(determinant_by_lu(m, 'generalring), p);
512 gf_determinant(m);
515 zn_determinant(m, p);
518 /* extension fields: AES mix columns operation */
520 ( gf_set_data(2, 8), gf_infolist() );
521 [2,x^8+x^4+x^3+x+1,x+1,256,255];
523 (ef_minimal_set(x^4+1), 0);
526 ef_irreducible_p(x^4+1);
527 false;
529 (ibase : obase : 16, 0);
532 m : matrix([0d4,0e0,0b8, 1e], [0bf,0b4, 41, 27], [ 5d, 52, 11, 98], [ 30,0ae,0f1,0e5]);
533 matrix([0D4,0E0,0B8,1E],[0BF,0B4,41,27],[5D,52,11,98],[30,0AE,0F1,0E5]);
535 c : ef_l2p(reverse(flatten(args(col(m, 1)))));
536 30*x^3+5D*x^2+0BF*x+0D4;
538 p3 : 3*x^3+x^2+x+2;
539 3*x^3+x^2+x+2;
541 ef_add(p3, c);
542 33*x^3+5C*x^2+0BE*x+0D6;
544 ef_mult(p3, c);
545 0E5*x^3+81*x^2+66*x+4;
547 i3 : ef_inv(p3);
548 0B*x^3+0D*x^2+9*x+0E;
550 ef_div(1, p3);
551 0B*x^3+0D*x^2+9*x+0E;
553 ef_exp(p3, -1);
554 0B*x^3+0D*x^2+9*x+0E;
556 ef_mult(p3, i3);
559 (ibase : obase : 10., 0);
562 /* this time use lookup tables : */
564 (gf_make_logs(), 0);
567 ord : gf_order();
568 255;
570 ( gf_coeff_mult(a, b) := 
571     if a = 0 or b = 0 then 0
572     else gf_powers[ ?mod(gf_logs[a] + gf_logs[b], ord) ],
574   gf_coeff_inv(a) := 
575     if a = 0 then 0
576     else gf_powers[ord - gf_logs[a]],
578   gf_coeff_add : ?logxor,
583 (ibase : obase : 16, 0);
586 ef_add(p3, c);
587 33*x^3+5C*x^2+0BE*x+0D6;
589 ef_mult(p3, c);
590 0E5*x^3+81*x^2+66*x+4;
592 i3 : ef_inv(p3);
593 0B*x^3+0D*x^2+9*x+0E;
595 ef_div(1, p3);
596 0B*x^3+0D*x^2+9*x+0E;
598 ef_exp(p3, -1);
599 0B*x^3+0D*x^2+9*x+0E;
601 ef_mult(p3, i3);
604 (ibase : obase : 10., 0);
607 /* 
608 extension fields:
610 examples taken from
611 Efficient Software Implementations of Large Finite Fields
612 by LUO, BOWERS, OPREA, XU
615 ( gf_set_data(2, x^8+x^4+x^3+x^2+1), gf_infolist() );
616 [2,x^8+x^4+x^3+x^2+1,x,256,255];
618 ef_irreducible_p(x^4+x^2+6*x+1);
619 true;
621 ef_irreducible_p(x^6+x^2+x+32);
622 true;
624 ef_irreducible_p(x^8+x^3+x+9);
625 true;
627 ef_irreducible_p(x^10+x^3+x+32);
628 true;
630 ef_irreducible_p(x^12+x^3+x+2);
631 true;
633 ef_irreducible_p(x^14+x^3+x+33);
634 true;
636 ef_irreducible_p(x^16+x^3+x+6);
637 true;
639 ( gf_set_data(2, x^16+x^12+x^3+x+1), gf_infolist() );
640 [2,x^16+x^12+x^3+x+1,x,65536,65535];
642 ef_irreducible_p(x^2+x+8192);
643 true;
645 ef_irreducible_p(x^3+x+1);
646 true;
648 ef_irreducible_p(x^4+x^2+2*x+1);
649 true;
651 ef_irreducible_p(x^5+x^2+1);
652 true;
654 ef_irreducible_p(x^6+x^3+8192);
655 true;
657 ef_irreducible_p(x^7+x+1);
658 true;
660 ef_irreducible_p(x^8+x^3+x+8);
661 true;
663 (remvalue(a,b,c,d,p,n,g,m,mi,ord,elem,r,s,u,prod,mp,mpi,fs,fs_ord,basis,coeffs,p3,i3), 
664   modulus : modulus_copy, gf_coeff_limit : gf_coeff_limit_copy, 0);
667 gf_unset();
668 true;
670 /* Bug #3427: gcfactor(3+4*%i) => error */
672 gcfactor(3+4*%i);
673 (2+%i)^2;
675 /* Bug #3265: gcfactor(0) -> division by zero */
677 gcfactor(0);
680 /* Bug #2839: gcfactor(9) => 9 */
682 block ([g : gcfactor (9)], [op (g), args (g)]);
683 ["^", [3, 2]];
685 block ([g : gcfactor (343)], [op (g), args (g)]);
686 ["^", [7, 3]];
688 /* A bogus comparison caused various things like gcfactor(x*%i)
689  * to signal lisp errors
690  */
692 gcfactor (x * %i);
693 x * %i;
695 gcfactor (%i / 2);
696 %i / 2;