tests/rtest_numth.mac

   1 (kill(all), 0);
   2 0;
   3
   4 /* (Z/pZ)* p prime */
   5 (p : 2^61-1, fs : ifactors(p - 1), g : zn_primroot(p, fs));
   6 37;
   7
   8 zn_primroot_p(power_mod(g, 7, p), p, fs);
   9 false;
  10
  11 zn_primroot_p(power_mod(g, 17, p), p, fs);
  12 true;
  13
  14 is(zn_order(g, p, fs) = totient(p));
  15 true;
  16
  17 is(zn_order(power_mod(g, 7, p), p, fs) = zn_order(g, p, fs));
  18 false;
  19
  20 is(zn_order(power_mod(g, 17, p), p, fs) = zn_order(g, p, fs));
  21 true;
  22
  23 (a : power_mod(g, 1234567890, p), zn_log(a, g, p, fs));
  24 1234567890;
  25
  26 (g : 3, ord_3 : zn_order(3,p,fs), fs_ord_3 : ifactors(ord_3), (p-1)/ord_3);
  27 9;
  28
  29 (a : power_mod(g, 1234567890, p), zn_log(a, g, p, fs));
  30 1234567890;
  31
  32 /* (Z/nZ)* n composite */
  33 n : 22;
  34 22;
  35
  36 g : zn_primroot(n);
  37 7;
  38
  39 zn_primroot_p(power_mod(g, 2, n), n);
  40 false;
  41
  42 zn_primroot_p(power_mod(g, 3, n), n);
  43 true;
  44
  45 zn_order(power_mod(g, 2, n), n);
  46 5;
  47
  48 zn_order(g, n);
  49 10;
  50
  51 a : power_mod(g, 8, n);
  52 9;
  53
  54 zn_log(a, g, n);
  55 8;
  56
  57 (g : 5, a : power_mod(g, 4, n), zn_log(a, g, n));
  58 4;
  59
  60 /* CRT */
  61 mods : [1009, 1013, 1019];
  62 [1009, 1013, 1019];
  63
  64 x : 374599943;
  65 374599943;
  66
  67 rems : map(lambda([z], mod(x, z)), mods);
  68 [621, 647, 258];
  69
  70 chinese(rems, mods);
  71 374599943;
  72
  73 (remvalue(p,fs,g,a,n,mods,x,rems), 0);
  74 0;
  75
  76 (kill(all), modulus_copy : modulus, modulus : false,
  77 gf_coeff_limit_copy : gf_coeff_limit, gf_coeff_limit : 256, 0);
  78 0;
  79
  80 /* Tests adapted from contrib/gf/gf_test.mac */
  81
  82 ( gf_set_data(123127, x^5+2*x+1), gf_infolist() );
  83 [123127,x^5+2*x+1,x+4,28298700309333316062584407,28298700309333316062584406];
  84
  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];
  87
  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;
  90
  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;
  93
  94 ( gf_set_data(19), gf_infolist() );
  95 [19,x,2,19,18];
  96
  97 gf_exp(gf_primitive(), gf_index(17));
  98 17;
  99
 100 ( gf_set_data(10000000019), gf_infolist() );
 101 [10000000019,x,2,10000000019,10000000018];
 102
 103 gf_exp(gf_primitive() ,gf_index(3));
 104 3;
 105
 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];
 110
 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];
 114
 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];
 117
 118 gf_exp(gf_primitive(), gf_index(x^10+1));
 119 x^10+1;
 120
 121 gf_minimal_poly(x^10+1);
 122 z^20+z^13+z^12+z^5+z^4+z^3+1;
 123
 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];
 127
 128 /* Tests adapted from contrib/gf/gf_hard_test.mac */
 129
 130 ( gf_set_data(7, x^61+x^4+1), gf_infolist() );
 131 [7,x^61+x^4+1,x+3,3556153025177363557255317383565515512407041673852007,
 132         3556153025177363557255317383565515512407041673852006];
 133
 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];
 138
 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];
 143
 144 ( gf_set_data(2, x^102+x^29+1), gf_infolist() );
 145 [2,x^102+x^29+1,x+1,5070602400912917605986812821504,
 146         5070602400912917605986812821503];
 147
 148 ( gf_set_data(5, x^61+x^15+1), gf_infolist() );
 149 [5,x^61+x^15+1,x+4,4336808689942017736029811203479766845703125,
 150         4336808689942017736029811203479766845703124];
 151
 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];
 156
 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];
 161
 162 ( gf_set_data(18659817111137), gf_infolist() );
 163 [18659817111137,x,3,18659817111137,18659817111136];
 164
 165 gf_log(7);
 166 15962290024269;
 167
 168 /* Examples adapted from contrib/gf/gf_manual.pdf */
 169
 170 ( gf_set_data(2, x^4+x+1), gf_infolist() );
 171 [2,x^4+x+1,x,16,15];
 172
 173 (a : x^3+x^2+1, b : x^3+x+1, 0);
 174 0;
 175
 176 gf_add(a, b);
 177 x^2+x;
 178
 179 gf_mult(a, b);
 180 x^2+x;
 181
 182 gf_inv(b);
 183 x^2+1;
 184
 185 gf_div(a, b);
 186 x^3+x^2;
 187
 188 gf_mult(a, gf_inv(b));
 189 x^3+x^2;
 190
 191 gf_exp(a, 10);
 192 x^2+x+1;
 193
 194 gf_exp(a, 15);
 195 1;
 196
 197 gf_primitive();
 198 x;
 199
 200 gf_index(a);
 201 13;
 202
 203 ev(a = gf_exp(x, 13)), pred;
 204 true;
 205
 206 (gf_make_logs(), 0);
 207 0;
 208
 209 gf_logs[10];
 210 9;
 211
 212 gf_n2p(10);
 213 x^3+x;
 214
 215 gf_index(x^3+x);
 216 9;
 217
 218 (a : x^2+x+1, b : x^3+x^2+1, 0);
 219 0;
 220
 221 gf_log(a, b);
 222 10;
 223
 224 gf_primitive_p(x^3+x+1);
 225 true;
 226
 227 gf_primitive_p(x^2+x);
 228 false;
 229
 230 gf_order(x^2+x);
 231 3;
 232
 233 gf_order(x^3+x+1);
 234 15;
 235
 236 (a : x^3+x+1, 0);
 237 0;
 238
 239 p : gf_minimal_poly(a);
 240 z^4+z^3+1;
 241
 242 p : subst(a, z, p);
 243 (x^3+x+1)^4+(x^3+x+1)^3+1;
 244
 245 gf_eval(p);
 246 0;
 247
 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];
 250
 251 (char : 7, exp : 4, g : 3*x^3+3*x^2+6, 0);
 252 0;
 253
 254 gf_primitive_p(g);
 255 true;
 256
 257 a : makelist(gf_log(x+i, g),i, 0, 6);
 258 [1028,1935,2054,1008,379,1780,223];
 259
 260 d : 1702;
 261 1702;
 262
 263 ord : char^exp - 1;
 264 2400;
 265
 266 c : makelist(mod(a[i] + d, ord), i, 1, 7);
 267 [330,1237,1356,310,2081,1082,1925];
 268
 269 m : [1,0,1,1,0,0,1];
 270 [1,0,1,1,0,0,1];
 271
 272 c : mod(sum(m[i] * c[i], i, 1, 7), ord);
 273 1521;
 274
 275 r : mod(c - exp*d, ord);
 276 1913;
 277
 278 u : gf_exp(g, r);
 279 x^3+3*x^2+2*x+5;
 280
 281 s : u + gf_reduction();
 282 x^4+4*x^3+8*x^2+8*x+7;
 283
 284 gf_factor(s);
 285 x*(x+2)*(x+3)*(x+6);
 286
 287 ( gf_set_data(2, x^4+x+1), gf_infolist() );
 288 [2,x^4+x+1,x,16,15];
 289
 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);
 291 0;
 292
 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]);
 295
 296 gf_matmult(m, mi);
 297 matrix([1,0,0], [0,1,0], [0,0,1]);
 298
 299 ( gf_set_data(2, x^10+x^3+1), gf_infolist() );
 300 [2,x^10+x^3+1,x,1024,1023];
 301
 302 elem : gf_normal();
 303 x^7;
 304
 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]);
 311
 312 gf_normal_basis_rep(elem, mi : gf_invert_by_lu(m));
 313 [1,0,0,0,0,0,0,0,0,0];
 314
 315 (a : x^9+x^7+x^6+x^5+x^3+x^2+x, 0);
 316 0;
 317
 318 gf_normal_basis_rep(a, mi);
 319 [1,1,1,0,1,0,1,1,1,0];
 320
 321 gf_normal_basis_rep(gf_exp(a, 2), mi);
 322 [0,1,1,1,0,1,0,1,1,1];
 323
 324 ( gf_set_data(2, x^20+x^3+1), gf_infolist() );
 325 [2,x^20+x^3+1,x,1048576,1048575];
 326
 327 (a : x^15+x^5+1, 0);
 328 0;
 329
 330 gf_index(a);
 331 720548;
 332
 333 gf_exp(a, 3^12);
 334 x^17+x^16+x^13+x^12+x^11+x^3+x^2+x;
 335
 336 /* some new tests */
 337
 338 ( gf_set_data(2, x^12+x^3+1), gf_infolist() );
 339 [2,x^12+x^3+1,x+1,4096,4095];
 340
 341 gf_log(gf_exp(x+1, 1234));
 342 1234;
 343
 344 ( gf_set_data(8796519617, x^2+3), gf_infolist() );
 345 [8796519617,x^2+3,x+9,77378757372265826689,77378757372265826688];
 346
 347 gf_log(gf_exp(x+9, 1234567890));
 348 1234567890;
 349
 350 ( gf_set_data(2, 4), gf_infolist() );
 351 [2,x^4+x+1,x,16,15];
 352
 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;
 357
 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);
 360
 361 (gf_minimal_set(2, x^17), 0);
 362 0;
 363
 364 apply(gf_mult, args(fs));
 365 x^16+x;
 366
 367 ( gf_set_data(13, x^7+3*x+2), gf_infolist() );
 368 [13,x^7+3*x+2,x,62748517,62748516];
 369
 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;
 372
 373 gf_normal_p(p);
 374 true;
 375
 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]);
 379
 380 gf_normal_basis_rep(p, mi : gf_invert_by_lu(m));
 381 [1,0,0,0,0,0,0];
 382
 383 coeffs : gf_normal_basis_rep(x^2+4*x+7, mi);
 384 [8,8,7,2,5,1,2];
 385
 386 ( basis : map(gf_l2p, args(transpose(m))),
 387   apply(gf_add, map(gf_mult, coeffs, basis)) );
 388 x^2+4*x+7;
 389
 390 ( gf_set_data(7, 4), gf_infolist() );
 391 [7,x^4+x+1,x+5,2401,2400];
 392
 393 ( gf_set_data(7, x^4+x^2+3), gf_infolist() );
 394 [7,x^4+x^2+3,x+1,2401,2400];
 395
 396 ( gf_set_data(2, x^8+1), gf_infolist() );
 397 [2,x^8+1,unknown,256,128];
 398
 399 gf_order();
 400 128;
 401
 402 gf_inv(x);
 403 x^7;
 404
 405 gf_inv(x+1);
 406 false;
 407
 408 gf_gcdex(x+1, gf_reduction());
 409 [1, 0, x+1];
 410
 411 ( gf_set_data(3, x^8+1), gf_infolist() );
 412 [3,x^8+1,unknown,6561,6400];
 413
 414 gf_order();
 415 6400;
 416
 417 ( gf_set_data(3, x^8+x+1), gf_infolist() );
 418 [3,x^8+x+1,unknown,6561,4368];
 419
 420 gf_order();
 421 4368;
 422
 423 ( gf_set_data(13, x^8+2), gf_infolist() );
 424 [13,x^8+2,x+2,815730721,815730720];
 425
 426 (a : x^7+x+1, b : x^3+3*x^2+9*x+7, 0);
 427 0;
 428
 429 gf_gcd(a, b);
 430 x^2+2*x+7;
 431
 432 gf_gcdex(a, b);
 433 [7, 6*x^4+8*x^3+3*x, x^2+2*x+7];
 434
 435 gf_primitive_poly(7, 8);
 436 x^8+x+3;
 437
 438 gf_primitive_poly_p(x^8+x+3, 7);
 439 true;
 440
 441 gf_primitive_poly_p(gf_irreducible(7, 8), 7);
 442 true;
 443
 444 gf_primitive_poly(2, 8);
 445 x^8+x^4+x^3+x^2+1;
 446
 447 gf_primitive_poly_p(x^8+x^4+x^3+x^2+1, 2);
 448 true;
 449
 450 gf_primitive_poly_p(gf_irreducible(2, 8), 2);
 451 false;
 452
 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);
 457 8;
 458
 459 /* lu decomposition: */
 460
 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];
 463
 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]);
 466
 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]);
 469
 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]);
 475
 476 mi : matrixmap(gf_p2n, mpi);
 477 matrix([14,11,13,9], [9,14,11,13], [13,9,14,11], [11,13,9,14]);
 478
 479 ( ef_set_data(x), ef_infolist() );
 480 [x,3,256,255];
 481
 482 is( ef_invert_by_lu(m) = mi );
 483 true;
 484
 485 ef_unset();
 486 true;
 487
 488 ( p : 17, n : 4, gf_set_data(p,n), 0);
 489 0;
 490
 491 (elem : 6*x^3+13*x^2+4*x+2, gf_normal_p(elem));
 492 true;
 493
 494 m : gf_normal_basis(elem);
 495 matrix([6,7,11,10],[13,4,13,4],[4,1,13,16],[2,2,2,2]);
 496
 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]);
 499
 500 is(gf_matmult(m, mi) = diagmatrix(n, 1));
 501 true;
 502
 503 is(zn_invert_by_lu(m, p) = mi);
 504 true;
 505
 506 mod(determinant(m), p);
 507 3;
 508
 509 mod(determinant_by_lu(m, 'generalring), p);
 510 3;
 511
 512 gf_determinant(m);
 513 3;
 514
 515 zn_determinant(m, p);
 516 3;
 517
 518 /* extension fields: AES mix columns operation */
 519
 520 ( gf_set_data(2, 8), gf_infolist() );
 521 [2,x^8+x^4+x^3+x+1,x+1,256,255];
 522
 523 (ef_minimal_set(x^4+1), 0);
 524 0;
 525
 526 ef_irreducible_p(x^4+1);
 527 false;
 528
 529 (ibase : obase : 16, 0);
 530 0;
 531
 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]);
 534
 535 c : ef_l2p(reverse(flatten(args(col(m, 1)))));
 536 30*x^3+5D*x^2+0BF*x+0D4;
 537
 538 p3 : 3*x^3+x^2+x+2;
 539 3*x^3+x^2+x+2;
 540
 541 ef_add(p3, c);
 542 33*x^3+5C*x^2+0BE*x+0D6;
 543
 544 ef_mult(p3, c);
 545 0E5*x^3+81*x^2+66*x+4;
 546
 547 i3 : ef_inv(p3);
 548 0B*x^3+0D*x^2+9*x+0E;
 549
 550 ef_div(1, p3);
 551 0B*x^3+0D*x^2+9*x+0E;
 552
 553 ef_exp(p3, -1);
 554 0B*x^3+0D*x^2+9*x+0E;
 555
 556 ef_mult(p3, i3);
 557 1;
 558
 559 (ibase : obase : 10., 0);
 560 0;
 561
 562 /* this time use lookup tables : */
 563
 564 (gf_make_logs(), 0);
 565 0;
 566
 567 ord : gf_order();
 568 255;
 569
 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) ],
 573
 574   gf_coeff_inv(a) :=
 575     if a = 0 then 0
 576     else gf_powers[ord - gf_logs[a]],
 577
 578   gf_coeff_add : ?logxor,
 579
 580 0);
 581 0;
 582
 583 (ibase : obase : 16, 0);
 584 0;
 585
 586 ef_add(p3, c);
 587 33*x^3+5C*x^2+0BE*x+0D6;
 588
 589 ef_mult(p3, c);
 590 0E5*x^3+81*x^2+66*x+4;
 591
 592 i3 : ef_inv(p3);
 593 0B*x^3+0D*x^2+9*x+0E;
 594
 595 ef_div(1, p3);
 596 0B*x^3+0D*x^2+9*x+0E;
 597
 598 ef_exp(p3, -1);
 599 0B*x^3+0D*x^2+9*x+0E;
 600
 601 ef_mult(p3, i3);
 602 1;
 603
 604 (ibase : obase : 10., 0);
 605 0;
 606
 607 /*
 608 extension fields:
 609
 610 examples taken from
 611 Efficient Software Implementations of Large Finite Fields
 612 by LUO, BOWERS, OPREA, XU
 613 */
 614
 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];
 617
 618 ef_irreducible_p(x^4+x^2+6*x+1);
 619 true;
 620
 621 ef_irreducible_p(x^6+x^2+x+32);
 622 true;
 623
 624 ef_irreducible_p(x^8+x^3+x+9);
 625 true;
 626
 627 ef_irreducible_p(x^10+x^3+x+32);
 628 true;
 629
 630 ef_irreducible_p(x^12+x^3+x+2);
 631 true;
 632
 633 ef_irreducible_p(x^14+x^3+x+33);
 634 true;
 635
 636 ef_irreducible_p(x^16+x^3+x+6);
 637 true;
 638
 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];
 641
 642 ef_irreducible_p(x^2+x+8192);
 643 true;
 644
 645 ef_irreducible_p(x^3+x+1);
 646 true;
 647
 648 ef_irreducible_p(x^4+x^2+2*x+1);
 649 true;
 650
 651 ef_irreducible_p(x^5+x^2+1);
 652 true;
 653
 654 ef_irreducible_p(x^6+x^3+8192);
 655 true;
 656
 657 ef_irreducible_p(x^7+x+1);
 658 true;
 659
 660 ef_irreducible_p(x^8+x^3+x+8);
 661 true;
 662
 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);
 665 0;
 666
 667 gf_unset();
 668 true;
 669
 670 /* Bug #3427: gcfactor(3+4*%i) => error */
 671
 672 gcfactor(3+4*%i);
 673 (2+%i)^2;
 674
 675 /* Bug #3265: gcfactor(0) -> division by zero */
 676
 677 gcfactor(0);
 678 0;
 679
 680 /* Bug #2839: gcfactor(9) => 9 */
 681
 682 block ([g : gcfactor (9)], [op (g), args (g)]);
 683 ["^", [3, 2]];
 684
 685 block ([g : gcfactor (343)], [op (g), args (g)]);
 686 ["^", [7, 3]];
 687
 688 /* A bogus comparison caused various things like gcfactor(x*%i)
 689  * to signal lisp errors
 690  */
 691
 692 gcfactor (x * %i);
 693 x * %i;
 694
 695 gcfactor (%i / 2);
 696 %i / 2;
 697
 698 /* Exhaustive testing of zn_determinant and zn_invert_by_lu,
 699  * inspired by:
 700  * https://stackoverflow.com/questions/68541899/maxima-wrongly-says-a-matrix-is-non-invertible
 701  *
 702  * Enumerate all matrices of size 2x2 and 3x3 and calculate
 703  * zn_determinant and zn_invert_by_lu with modulus = 2.
 704  * Repeat with modulus = 3 but process only 2x2 matrices.
 705  * There are enough 3x3 matrices with modulus = 3 to cause
 706  * trouble with some Lisp implementations; too many arguments.
 707  */
 708
 709 (enumerate_matrices_modulo (p, n) :=
 710   (enumerate_lists_modulo (p, n^2),
 711    map (lambda ([L], reshape_matrix_from_list (L, n)), %%)),
 712
 713  enumerate_lists_modulo (p, nn) :=
 714    if nn = 1
 715      then makelist ([i], i, 0, p - 1)
 716      else tree_reduce ('append,
 717                        map (lambda ([L], makelist (cons (i, L), i, 0, p - 1)),
 718                             enumerate_lists_modulo (p, nn - 1))),
 719
 720  reshape_matrix_from_list (L, n) :=
 721    apply ('matrix, makelist (makelist (L[1 + (i - 1)*n + (j - 1)], j, 1, n), i, 1, n)),
 722
 723  test_zn_determinant (p, M) :=
 724    block ([det_M: determinant (M),
 725            zn_det_M: zn_determinant (M, p)],
 726            if zn_det_M = mod (det_M, p)
 727              then true
 728              else failed ('zn_determinant (M, p))),
 729
 730  test_zn_invert_by_lu (p, M) :=
 731    block ([zn_inv_M: zn_invert_by_lu (M, p)],
 732           if zn_inv_M = false
 733             then if mod (determinant (M), p) = 0
 734                    then true
 735                    else failed_a ('zn_invert_by_lu (M, p))
 736             else block ([M_Minv, Minv_M,
 737                          scalarmatrixp: false],
 738                         M_Minv: mod (M . zn_inv_M, p),
 739                         Minv_M: mod (zn_inv_M . M, p),
 740                         if identity_matrix_p (M_Minv) and identity_matrix_p (Minv_M)
 741                           then true
 742                           else failed_b ('zn_invert_by_lu (M, p)))),
 743
 744  identity_matrix_p (M) := is (M = ident (length (M))),
 745
 746  0);
 747 0;
 748
 749 (enumerate_matrices_modulo (2, 1),
 750  map (lambda ([M], test_zn_determinant (2, M)), %%),
 751  sublist (%%, lambda ([x], x # true)));
 752 [];
 753
 754 (enumerate_matrices_modulo (2, 2),
 755  map (lambda ([M], test_zn_determinant (2, M)), %%),
 756  sublist (%%, lambda ([x], x # true)));
 757 [];
 758
 759 (enumerate_matrices_modulo (2, 3),
 760  map (lambda ([M], test_zn_determinant (2, M)), %%),
 761  sublist (%%, lambda ([x], x # true)));
 762 [];
 763
 764 (enumerate_matrices_modulo (3, 1),
 765  map (lambda ([M], test_zn_determinant (3, M)), %%),
 766  sublist (%%, lambda ([x], x # true)));
 767 [];
 768
 769 (enumerate_matrices_modulo (3, 2),
 770  map (lambda ([M], test_zn_determinant (3, M)), %%),
 771  sublist (%%, lambda ([x], x # true)));
 772 [];
 773
 774 /* time consuming, skip it for now
 775 (enumerate_matrices_modulo (3, 3),
 776  map (lambda ([M], test_zn_determinant (3, M)), %%),
 777  sublist (%%, lambda ([x], x # true)));
 778 [];
 779  */
 780
 781 (enumerate_matrices_modulo (2, 1),
 782  map (lambda ([M], test_zn_invert_by_lu (2, M)), %%),
 783  sublist (%%, lambda ([x], x # true)));
 784 [];
 785
 786 (enumerate_matrices_modulo (2, 2),
 787  map (lambda ([M], test_zn_invert_by_lu (2, M)), %%),
 788  sublist (%%, lambda ([x], x # true)));
 789 [];
 790
 791 (enumerate_matrices_modulo (2, 3),
 792  map (lambda ([M], test_zn_invert_by_lu (2, M)), %%),
 793  sublist (%%, lambda ([x], x # true)));
 794 [];
 795
 796 (enumerate_matrices_modulo (3, 1),
 797  map (lambda ([M], test_zn_invert_by_lu (3, M)), %%),
 798  sublist (%%, lambda ([x], x # true)));
 799 [];
 800
 801 (enumerate_matrices_modulo (3, 2),
 802  map (lambda ([M], test_zn_invert_by_lu (3, M)), %%),
 803  sublist (%%, lambda ([x], x # true)));
 804 [];
 805
 806 /* time consuming, skip it for now
 807 (enumerate_matrices_modulo (3, 3),
 808  map (lambda ([M], test_zn_invert_by_lu (3, M)), %%),
 809  sublist (%%, lambda ([x], x # true)));
 810 [];
 811  */