1 /*============================================================================== 2 | 3 | File Name: 4 | 5 | ppHelperFunc.c 6 | 7 | Description: 8 | 9 | Higher level helper functions which don't belong elsewhere. 10 | 11 | Functions: 12 | 13 | initial_trial_poly 14 | next_trial_poly 15 | const_coeff_test 16 | const_coeff_is_primitive_root 17 | skip_test 18 | has_multi_irred_factors 19 | generate_Q_matrix 20 | find_nullity 21 | 22 | LEGAL 23 | 24 | Primpoly Version 16.1 - A Program for Computing Primitive Polynomials. 25 | Copyright (C) 1999-2021 by Sean Erik O'Connor. All Rights Reserved. 26 | 27 | This program is free software: you can redistribute it and/or modify 28 | it under the terms of the GNU General Public License as published by 29 | the Free Software Foundation, either version 3 of the License, or 30 | (at your option) any later version. 31 | 32 | This program is distributed in the hope that it will be useful, 33 | but WITHOUT ANY WARRANTY; without even the implied warranty of 34 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 35 | GNU General Public License for more details. 36 | 37 | You should have received a copy of the GNU General Public License 38 | along with this program. If not, see <http://www.gnu.org/licenses/>. 39 | 40 | The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com 41 | with !DOT! replaced by . and the !AT! replaced by @ 42 | 43 ==============================================================================*/ 44 45 /*------------------------------------------------------------------------------ 46 | Include Files | 47 ------------------------------------------------------------------------------*/ 48 49 #include <stdio.h> 50 #include <stdlib.h> 51 #include <memory.h> 52 53 #include "Primpoly.h" 54 55 56 /*============================================================================== 57 | initial_trial_poly | 58 ================================================================================ 59 60 DESCRIPTION 61 62 Return the initial monic polynomial in the sequence for f(x). 63 64 INPUT 65 66 f (int *) Monic polynomial f(x). 67 n (int, n >= 1) Degree of f(x). 68 69 RETURNS 70 n 71 f (int *) Sets f( x ) = x - 1 72 73 EXAMPLE 74 4 75 Let n = 4. Set f(x) = x - 1. 76 77 METHOD 78 79 BUGS 80 81 None. 82 83 -------------------------------------------------------------------------------- 84 | Function Call | 85 ------------------------------------------------------------------------------*/ 86 87 void initial_trial_poly( int * f, int n ) 88 { 89 90 /*------------------------------------------------------------------------------ 91 | Local Variables | 92 ------------------------------------------------------------------------------*/ 93 94 int i ; 95 96 /*------------------------------------------------------------------------------ 97 | Function Body | 98 ------------------------------------------------------------------------------*/ 99 100 f[ 0 ] = -1 ; 101 102 for (i = 1 ; i <= n-1 ; ++i) 103 104 f[ i ] = 0 ; 105 106 f[ n ] = 1 ; 107 108 } /* ================= end of function initial_trial_poly ==================== */ 109 110 111 112 /*============================================================================== 113 | next_trial_poly | 114 ================================================================================ 115 116 DESCRIPTION 117 118 Return the next monic polynomial in the sequence after f(x). 119 120 INPUT 121 122 f (int *) Monic polynomial f(x). 123 n (int, n >= 1) Degree of monic polynomial f(x). 124 p (int, p >= 2) Modulo p coefficient arithmetic. 125 126 RETURNS 127 128 f (int *) Overwrites f(x) with the next polynomial 129 after it in the sequence (explained below). 130 EXAMPLE 131 3 132 Let n = 3 and p = 5. Suppose f(x) = x + 4 x + 4. As a mod p number, 133 this is 1 0 4 4. Adding 1 gives 1 0 4 5. We reduce modulo 134 5 and propagate the carry to get the number 1 0 5 0. Propagating 135 the carry again and reducing gives 1 1 0 0. The next polynomial 136 3 2 137 after f(x) is x + x . 138 139 METHOD 140 141 Think of the polynomial coefficients as the digits of a number written 142 in base p. The "next" polynomial is the one you would get by adding 1 143 to this number in multiple precision arithmetic. Our intention is to 144 run through all possible monic polynomials modulo p. 145 146 Propagate carries in digits 1 through n-2 when any digit exceeds p. No 147 carries take place in the n-1 st digit because our polynomial is monic. 148 149 BUGS 150 151 None. 152 153 -------------------------------------------------------------------------------- 154 | Function Call | 155 ------------------------------------------------------------------------------*/ 156 157 void 158 next_trial_poly( int * f, int n, int p ) 159 { 160 161 /*------------------------------------------------------------------------------ 162 | Local Variables | 163 ------------------------------------------------------------------------------*/ 164 165 int 166 digit_num ; /* Loop counter and digit number. */ 167 168 /*------------------------------------------------------------------------------ 169 | Function Body | 170 ------------------------------------------------------------------------------*/ 171 172 ++f[ 0 ] ; /* Add 1, i.e. increment the coefficient of the x term. */ 173 174 /* 175 Sweep through the number from right to left, propagating carries. Skip 176 the constant and the nth degree terms. 177 */ 178 179 for (digit_num = 0 ; digit_num <= n - 2 ; ++digit_num) 180 { 181 if (f[ digit_num ] == p) /* Propagate carry to next digit. */ 182 { 183 f[ digit_num ] = 0 ; 184 ++f[ digit_num + 1 ] ; 185 } 186 } 187 188 } /* ================= end of function next_trial_poly ====================== */ 189 190 191 /*============================================================================== 192 | const_coeff_test | 193 ================================================================================ 194 195 DESCRIPTION 196 197 n 198 Test if a = (-1) a (mod p) where a is the constant coefficient 199 0 0 200 of polynomial f(x) and a is the number from the order_r() test. 201 202 INPUT 203 204 f (int *) nth degree monic mod p polynomial f(x). 205 n (int) Its degree. 206 p (int) Modulus for coefficient arithmetic. 207 a (int) 208 209 RETURNS 210 211 YES if we pass the test. 212 NO otherwise 213 214 EXAMPLE 215 11 3 216 Let p = 5, n = 11, f( x ) = x + 2 x + 1, and a = 4. 217 Then return YES since 218 219 11 220 (-1) * 1 (mod 5) = -1 (mod 5) = 4 (mod 5). 221 222 METHOD 223 n 224 We test if (a - (-1) a ) mod p = 0. 225 0 226 BUGS 227 228 None. 229 230 -------------------------------------------------------------------------------- 231 | Function Call | 232 ------------------------------------------------------------------------------*/ 233 234 int const_coeff_test( int * f, int n, int p, int a ) 235 { 236 237 /*------------------------------------------------------------------------------ 238 | Local Variables | 239 ------------------------------------------------------------------------------*/ 240 241 int constantCoeff = f[ 0 ] ; 242 243 /*------------------------------------------------------------------------------ 244 | Function Body | 245 ------------------------------------------------------------------------------*/ 246 247 if (n % 2 != 0) 248 constantCoeff = -constantCoeff ; /* (-1)^n < 0 for odd n. */ 249 250 return( (mod( a - constantCoeff, p ) == 0) ? YES : NO ) ; 251 252 } /* =================== end of function const_coeff_test ================ */ 253 254 255 256 /*============================================================================== 257 | const_coeff_is_primitive_root | 258 ================================================================================ 259 260 DESCRIPTION 261 262 n 263 Test if (-1) a (mod p) is a primitive root of the prime p where 264 0 265 a is the constant term of the polynomial f(x). 266 0 267 268 INPUT 269 270 f (int *) nth degree monic mod p polynomial f(x). 271 n (int) Its degree. 272 p (int) Modulus for coefficient arithmetic. 273 274 EXAMPLE 275 11 3 276 Let p = 7, n = 11, f( x ) = x + 2 x + 4. 277 Then return YES since 278 279 11 280 (-1) * 4 (mod 7) = -4 (mod 7) = 3 (mod 7), and 3 is a 281 1 2 3 282 primitive root of the prime 7 since 3 = 3, 3 = 2, 3 = 6, 283 4 5 7-1 284 3 = 4, 3 = 5 (mod 7); all not equal to 1 until 3 = 1 (mod 7). 285 286 METHOD 287 288 n 289 We test if (-1) a mod p is a primitive root mod p. 290 0 291 BUGS 292 293 None. 294 295 -------------------------------------------------------------------------------- 296 | Function Call | 297 ------------------------------------------------------------------------------*/ 298 299 int const_coeff_is_primitive_root( int * f, int n, int p ) 300 { 301 302 /*------------------------------------------------------------------------------ 303 | Local Variables | 304 ------------------------------------------------------------------------------*/ 305 306 int constantCoeff = f[ 0 ] ; 307 308 /*------------------------------------------------------------------------------ 309 | Function Body | 310 ------------------------------------------------------------------------------*/ 311 312 /* (-1)^n < 0 for odd n. */ 313 if (n % 2 != 0) 314 constantCoeff = -constantCoeff ; 315 316 return is_primitive_root( mod( constantCoeff, p ), p ) ; 317 318 } /* ============= of function const_coeff_is_primitive_root ================ */ 319 320 321 322 /*============================================================================== 323 | skip_test | 324 ================================================================================ 325 326 DESCRIPTION 327 328 Make the test p | (p - 1). 329 i 330 331 INPUT 332 333 i (int) ith prime divisor of r. 334 335 primes (bigint *) Array of distict prime factors of r. primes[ i ] = p 336 i 337 p (int) p >= 2. 338 339 RETURNS 340 341 YES if the test succeeds for some i. 342 NO otherwise 343 344 EXAMPLE 345 346 Suppose i = 0, primes[ 0 ] = 2 and p = 5. Return YES since 2 | 5-1. 347 348 METHOD 349 350 Test if (p - 1) mod p = 0 for all i. 351 i 352 BUGS 353 354 None. 355 356 -------------------------------------------------------------------------------- 357 | Function Call | 358 ------------------------------------------------------------------------------*/ 359 360 int 361 skip_test( int i, bigint * primes, int p ) 362 { 363 364 /*------------------------------------------------------------------------------ 365 | Function Body | 366 ------------------------------------------------------------------------------*/ 367 368 /* 369 p cannot divide the smaller number (p-1) 370 i 371 */ 372 if ( (bigint)( p-1 ) < primes[ i ] ) 373 return NO ; 374 else 375 return( ((bigint)(p-1) % primes[ i ] == 0) ? YES : NO ) ; 376 377 } /* ====================== end of function skip_test ======================= */ 378 379 380 381 382 /*============================================================================== 383 | has_multi_irred_factors | 384 ================================================================================ 385 386 DESCRIPTION 387 388 Find out if the monic polynomial f( x ) has two or more distinct 389 irreducible factors. 390 391 INPUT 392 393 power_table (int **) x ^ k (mod f(x), p) for n <= k <= 2n-2, f monic. 394 n (int, n >= 1) Degree of monic polynomial f(x). 395 p (int, p >= 2) Modulo p coefficient arithmetic. 396 397 RETURNS 398 399 1 if the monic polynomial f( x ) has two or more distinct rreducible factors 400 0 otherwise. 401 402 EXAMPLE 403 404 Let n = 4, p = 5 405 406 4 407 f( x ) = x + 2 is irreducible, so has one distinct factor. 408 409 4 3 2 4 410 f( x ) = x + 4x + x + 4x + 1 = (x + 1) has one distinct factor. 411 412 3 2 413 f( x ) = x + 3 = (x + 3x + 4)(x + 2) has two distinct irreducible factors. 414 415 4 3 2 2 416 f( x ) = x + 3x + 3x + 3x + 2 = (x + 1) (x + 2) (x + 3) has 3 distinct 417 irreducible factors. 418 419 4 420 f(x) = x + 4 = (x+1)(x+2)(x+3)(x+4) has 4 distinct irreducible factors. 421 422 METHOD 423 424 Berlekamp's method for factoring polynomials over GF( p ), modified to test 425 for irreducibility only. 426 427 See my notes; I skip the polynomial GCD step which ensures polynomials 428 are square-free due to time constraints, but this requires a proof that 429 the method is still valid. 430 431 BUGS 432 433 None. 434 435 -------------------------------------------------------------------------------- 436 | Function Call | 437 ------------------------------------------------------------------------------*/ 438 439 int has_multi_irred_factors( int power_table[][ MAXDEGPOLY ], int n, int p ) 440 { 441 int ** Q ; 442 int row ; 443 int nullity = 0 ; 444 445 446 /* Allocate space for the Q matrix. */ 447 Q = (int **) calloc( n, sizeof( int * ) ) ; 448 449 for (row = 0 ; row < n ; ++row) 450 { 451 Q[ row ] = (int *)calloc( n, sizeof( int ) ) ; 452 } 453 454 455 /* Generate the Q-I matrix. */ 456 generate_Q_matrix( Q, power_table, n, p ) ; 457 458 459 /* Find nullity of Q-I */ 460 nullity = find_nullity( Q, n, p ) ; 461 462 /* Free space for the Q matrix. */ 463 for (row = 0 ; row < n ; ++row) 464 { 465 free( Q[ row ] ) ; 466 } 467 468 free( Q ) ; 469 470 471 /* If nullity >= 2, f( x ) is a reducible polynomial modulo p since it has */ 472 /* two or more distinct irreducible factors. */ 473 /* e */ 474 /* Nullity of 1 implies f( x ) = p( x ) for some power e >= 1 so we cannot */ 475 /* determine reducibility. */ 476 if (nullity >= 2) 477 return 1 ; 478 479 return 0 ; 480 481 } /* =============== end of function has_multi_irred_factors ================ */ 482 483 484 485 /*============================================================================== 486 | generate_Q_matrix | 487 ================================================================================ 488 489 DESCRIPTION 490 491 492 Generate n x n matrix Q - I, where rows of Q are the powers, 493 494 p 2p (n-1) p 495 1, x (mod f(x),p), x (mod f(x), p), ... , x (mod f(x), p) 496 497 for monic polynomial f(x). 498 499 500 INPUT 501 502 Q (int **) Memory is allocated for this matrix already and 503 all entries are 0. 504 505 k 506 power_table (int **) x (mod f(x), p) for n <= k <= 2n-2, f monic. 507 508 n (int, n >= 1) Degree of monic polynomial f(x). 509 510 p (int, p >= 2) Modulo p coefficient arithmetic. 511 512 RETURNS 513 514 515 EXAMPLE 516 517 4 518 f(x) = x + 2, n = 4, p = 5 519 520 ( 1 ) ( 1 ) ( 1 0 0 0 ) 521 ( 5 ) ( ) ( ) 522 Q = ( x ) ( 3x ) ( 0 3 0 0 ) 523 ( 10 ) = ( 2 ) = ( ) 524 ( x ) ( 4x ) ( 0 0 4 0 ) 525 ( 15 ) ( 3 ) ( ) 526 ( x ) ( 2x ) ( 0 0 0 2 ) 527 528 2 3 529 1 x x x 530 531 ( 0 0 0 0 ) 532 ( 0 2 0 0 ) 533 Q - I = ( 0 0 3 0 ) 534 ( 0 0 0 1 ) 535 536 METHOD 537 538 Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, 539 Addison-Wesley. 540 541 BUGS 542 543 None. 544 545 -------------------------------------------------------------------------------- 546 | Function Call | 547 ------------------------------------------------------------------------------*/ 548 549 void 550 generate_Q_matrix( int ** Q, int power_table[][ MAXDEGPOLY ], int n, int p ) 551 { 552 553 /*------------------------------------------------------------------------------ 554 | Local Variables | 555 ------------------------------------------------------------------------------*/ 556 557 int xp[ MAXDEGPOLY ] ; /* x^p (mod f(x),p) */ 558 int q [ MAXDEGPOLY ] ; /* Current row of Q matrix. */ 559 560 int row = 0 ; 561 562 563 /*------------------------------------------------------------------------------ 564 | Function Body | 565 ------------------------------------------------------------------------------*/ 566 567 /* Check for invalid inputs. */ 568 if (n < 2 || p < 2 || Q == (int **)0) 569 { 570 return ; 571 } 572 573 574 575 /* Row 0 of Q = 1. */ 576 Q[ 0 ][ 0 ] = 1 ; 577 578 579 /* p 580 Find x (mod f(x),p) save the value into q(x) and 581 set row 1 of Q to this value. 582 */ 583 x_to_power( (bigint) p, xp, power_table, n, p ) ; 584 585 memcpy( q, xp, n * sizeof( int ) ) ; 586 memcpy( Q[ 1 ], q, n * sizeof( int ) ) ; 587 588 589 /* pk 590 Row k of Q = x (mod f(x), p) 2 <= k <= n-1, computed by 591 p 592 multiplying each previous row by x (mod f(x),p). 593 */ 594 for (row = 2 ; row <= n-1 ; ++row) 595 { 596 product( q, xp, power_table, n, p ) ; 597 memcpy( Q[ row ], q, n * sizeof( int ) ) ; 598 } 599 600 601 /* Subtract Q - I */ 602 603 for (row = 0 ; row < n ; ++row) 604 { 605 Q[ row ][ row ] = mod( Q[ row ][ row ] - 1, p ) ; 606 } 607 608 return ; 609 610 } /* ================== end of function generate_Q_matrix =================== */ 611 612 613 614 615 /*============================================================================== 616 | find_nullity | 617 ================================================================================ 618 619 DESCRIPTION 620 621 622 INPUT 623 624 Q (int **) Matrix of integers mod p. 625 n (int, n >= 1) Degree of monic polynomial f(x). 626 p (int, p >= 2) Modulo p coefficient arithmetic. 627 628 RETURNS 629 630 Nullity of Q. However if the nullity is 2 or more, we just return 2. 631 632 EXAMPLE 633 634 Let p = 5 and n = 3. We use the facts that -1 = 4 (mod 5), 1/2 = 3, -1/2 = 2, 635 1/3 = 2, -1/3 = 3, 1/4 = 4, -1/4 = 1. 636 637 Consider the matrix 638 639 ( 2 3 4 ) 640 Q = ( 0 2 1 ) 641 ( 3 3 3 ) 642 643 Begin with row 1. No pivotal columns have been selected yet. Scan row 1 and 644 pick column 1 as the pivotal column because it contains a nonzero element. 645 646 Normalizing column 1 by multiplying by -1/pivot = -1/2 = 2 gives 647 648 ( 4 3 4 ) 649 ( 0 2 1 ) 650 ( 1 3 3 ) 651 652 Now perform column reduction on column 2 by multiplying the pivotal column 1 653 by 3 (the column 2 element in the current row) and adding back to row 2. 654 655 ( 4 0 4 ) 656 ( 0 2 1 ) 657 ( 1 1 3 ) 658 659 Column reduce column 3 by multiplying the pivotal column by 4 and adding back to row 3, 660 661 ( 4 0 0 ) 662 ( 0 2 1 ) 663 ( 1 1 2 ) 664 665 For row 2, pick column 2 as the pivotal column, normalize it and reduce columns 1, then 3, 666 667 ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) 668 ( 0 2 1 ) => ( 0 4 1 ) => ( 0 4 1 ) => ( 0 4 0 ) 669 ( 1 1 2 ) ( 1 2 2 ) ( 1 2 2 ) ( 1 2 4 ) 670 norm. c.r. 1 c.r. 3 671 672 For row 3, we must pick column 3 as pivotal column because we've used up columns 1 and 2, 673 674 ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) ( 4 0 0 ) 675 ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) => ( 0 4 0 ) 676 ( 1 2 4 ) ( 1 2 4 ) ( 1 2 4 ) ( 0 0 4 ) 677 norm. c.r. 1 c.r. 2 678 679 The nullity is zero, since we were always able to find a pivot in each row. 680 681 METHOD 682 683 Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth, Addison-Wesley. 684 685 We combine operations of normalization of columns, 686 687 ( c ) ( C ) 688 ( c ) ( C ) 689 ( . ) ( C ) 690 ( . . . q . . . ) row ================> ( . . . -1 . . . ) row 691 ( c ) ( C ) 692 ( c ) column times ( C ) 693 ( c ) -1/a modulo p ( C ) 694 pivotCol pivotCol 695 696 and column reduction, 697 698 ( C b ) ( C B ) 699 ( C b ) ( C B ) 700 ( C b ) ( C B ) 701 ( . . . -1 . . .e . . . ) row ================> ( . . . -1 . . . 0 . . . ) 702 ( C b ) ( C B ) 703 ( C b ) pivotCol times ( C B ) 704 ( C b ) e added back to col ( C B ) 705 pivotCol col col 706 707 708 to reduce the matrix to a form in which columns having pivots are zero until 709 the pivotal row. 710 711 The column operations don't change the left nullspace of the 712 matrix. 713 714 The matrix rank is the number of pivotal rows since they are linearly 715 independent. The nullity is then the number of non-pivotal rows. 716 717 BUGS 718 719 720 -------------------------------------------------------------------------------- 721 | Function Call | 722 ------------------------------------------------------------------------------*/ 723 724 int find_nullity( int ** Q, int n, int p ) 725 { 726 727 int colFlag[ MAXDEGPOLY ] ; /* Is -1 if the column has no pivotal element. */ 728 int nullity = 0 ; 729 int row, col ; 730 int r ; 731 int found = NO ; 732 int pivotCol = -1 ; 733 int q = 0 ; 734 int t = 0 ; 735 736 /* Initialize column flags to -1. */ 737 memset( colFlag, -1, n * sizeof( int ) ) ; 738 739 /* Sweep through each row. */ 740 for (row = 0 ; row < n ; ++row) 741 { 742 /* Search for a pivot in this row: a non-zero element 743 in a column which had no previous pivot. 744 */ 745 found = NO ; 746 for (col = 0 ; col < n ; ++col) 747 { 748 if (Q[ row ][ col ] > 0 && colFlag[ col ] < 0) 749 { 750 found = YES ; 751 pivotCol = col ; 752 break ; 753 } 754 } 755 756 /* No pivot; increase nullity by 1. */ 757 if (found == NO) 758 { 759 /* Early out. */ 760 if (++nullity >= 2) 761 { 762 return nullity ; 763 } 764 } 765 /* Found a pivot, q. */ 766 else 767 { 768 q = Q[ row ][ pivotCol ] ; 769 770 /* Compute -1/q (mod p) */ 771 t = mod( -inverse_mod_p( q, p ), p ) ; 772 773 /* Normalize the pivotal column. */ 774 for (r = 0 ; r < n ; ++r) 775 { 776 Q[ r ][ pivotCol ] = mod( t * Q[ r ][ pivotCol ], p ) ; 777 } 778 779 /* Do column reduction: Add C times the pivotal column to the other 780 columns where C = element in the other column at current row. */ 781 for (col = 0 ; col < n ; ++col) 782 { 783 if (col != pivotCol) 784 { 785 q = Q[ row ][ col ] ; 786 787 for (r = 0 ; r < n ; ++r) 788 { 789 t = mod( q * Q[ r ][ pivotCol ], p ) ; 790 Q[ r ][ col ] = mod( t + Q[ r ][ col ], p ) ; 791 } 792 } 793 } 794 795 /* Record the presence of a pivot in this column. */ 796 colFlag[ pivotCol ] = row ; 797 798 } /* found a pivot */ 799 800 } /* end for row */ 801 802 return nullity ; 803 804 } /* ===================== end of function find_nullity ===================== */ 805 806 807 #if 0 808 int test 809 /* n 2n-2 810 Precompute the powers x , ..., x (mod f(x), p) 811 for use in all later computations. 812 */ 813 n = 4 ; p = 5 ; 814 f[0] = 2 ; f[1] = 3 ; f[2] = 3 ; f[3] = 3 ; f[4] = 1 ; 815 816 construct_power_table( power_table, f, n, p ) ; 817 if (has_multi_irred_factors( power_table, n, p ) == 1) 818 printf( "Pass\n" ) ; 819 else 820 printf( "Fail\n" ) ; 821 822 #endif
