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.4 - A Program for Computing Primitive Polynomials.
25| Copyright (C) 1999-2025 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 seanerikoconnor!AT!gmail!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
60DESCRIPTION
61
62 Return the initial monic polynomial in the sequence for f(x).
63
64INPUT
65
66 f (int *) Monic polynomial f(x).
67 n (int, n >= 1) Degree of f(x).
68
69RETURNS
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
77METHOD
78
79BUGS
80
81 None.
82
83--------------------------------------------------------------------------------
84| Function Call |
85------------------------------------------------------------------------------*/
86
87void initial_trial_poly( int * f, int n )
88{
89
90/*------------------------------------------------------------------------------
91| Local Variables |
92------------------------------------------------------------------------------*/
93
94int i ;
95
96/*------------------------------------------------------------------------------
97| Function Body |
98------------------------------------------------------------------------------*/
99
100f[ 0 ] = -1 ;
101
102for (i = 1 ; i <= n-1 ; ++i)
103
104 f[ i ] = 0 ;
105
106f[ n ] = 1 ;
107
108} /* ================= end of function initial_trial_poly ==================== */
109
110
111
112/*==============================================================================
113| next_trial_poly |
114================================================================================
115
116DESCRIPTION
117
118 Return the next monic polynomial in the sequence after f(x).
119
120INPUT
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
126RETURNS
127
128 f (int *) Overwrites f(x) with the next polynomial
129 after it in the sequence (explained below).
130EXAMPLE
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
139METHOD
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
149BUGS
150
151 None.
152
153--------------------------------------------------------------------------------
154| Function Call |
155------------------------------------------------------------------------------*/
156
157void
158 next_trial_poly( int * f, int n, int p )
159{
160
161/*------------------------------------------------------------------------------
162| Local Variables |
163------------------------------------------------------------------------------*/
164
165int
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
179for (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
195DESCRIPTION
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
202INPUT
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
209RETURNS
210
211 YES if we pass the test.
212 NO otherwise
213
214EXAMPLE
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
222METHOD
223 n
224 We test if (a - (-1) a ) mod p = 0.
225 0
226BUGS
227
228 None.
229
230--------------------------------------------------------------------------------
231| Function Call |
232------------------------------------------------------------------------------*/
233
234int const_coeff_test( int * f, int n, int p, int a )
235{
236
237/*------------------------------------------------------------------------------
238| Local Variables |
239------------------------------------------------------------------------------*/
240
241int constantCoeff = f[ 0 ] ;
242
243/*------------------------------------------------------------------------------
244| Function Body |
245------------------------------------------------------------------------------*/
246
247if (n % 2 != 0)
248 constantCoeff = -constantCoeff ; /* (-1)^n < 0 for odd n. */
249
250return( (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
260DESCRIPTION
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
268INPUT
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
274EXAMPLE
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
286METHOD
287
288 n
289 We test if (-1) a mod p is a primitive root mod p.
290 0
291BUGS
292
293 None.
294
295--------------------------------------------------------------------------------
296| Function Call |
297------------------------------------------------------------------------------*/
298
299int const_coeff_is_primitive_root( int * f, int n, int p )
300{
301
302/*------------------------------------------------------------------------------
303| Local Variables |
304------------------------------------------------------------------------------*/
305
306int constantCoeff = f[ 0 ] ;
307
308/*------------------------------------------------------------------------------
309| Function Body |
310------------------------------------------------------------------------------*/
311
312/* (-1)^n < 0 for odd n. */
313if (n % 2 != 0)
314 constantCoeff = -constantCoeff ;
315
316return 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
326DESCRIPTION
327
328 Make the test p | (p - 1).
329 i
330
331INPUT
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
339RETURNS
340
341 YES if the test succeeds for some i.
342 NO otherwise
343
344EXAMPLE
345
346 Suppose i = 0, primes[ 0 ] = 2 and p = 5. Return YES since 2 | 5-1.
347
348METHOD
349
350 Test if (p - 1) mod p = 0 for all i.
351 i
352BUGS
353
354 None.
355
356--------------------------------------------------------------------------------
357| Function Call |
358------------------------------------------------------------------------------*/
359
360int
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*/
372if ( (bigint)( p-1 ) < primes[ i ] )
373 return NO ;
374else
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
386DESCRIPTION
387
388 Find out if the monic polynomial f( x ) has two or more distinct
389 irreducible factors.
390
391INPUT
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
397RETURNS
398
399 1 if the monic polynomial f( x ) has two or more distinct rreducible factors
400 0 otherwise.
401
402EXAMPLE
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
422METHOD
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
431BUGS
432
433 None.
434
435--------------------------------------------------------------------------------
436| Function Call |
437------------------------------------------------------------------------------*/
438
439int 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
489DESCRIPTION
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
500INPUT
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
512RETURNS
513
514
515EXAMPLE
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
536METHOD
537
538 Modified from ART OF COMPUTER PROGRAMMING, Vol. 2, 2nd ed., Donald E. Knuth,
539 Addison-Wesley.
540
541BUGS
542
543 None.
544
545--------------------------------------------------------------------------------
546| Function Call |
547------------------------------------------------------------------------------*/
548
549void
550 generate_Q_matrix( int ** Q, int power_table[][ MAXDEGPOLY ], int n, int p )
551{
552
553/*------------------------------------------------------------------------------
554| Local Variables |
555------------------------------------------------------------------------------*/
556
557int xp[ MAXDEGPOLY ] ; /* x^p (mod f(x),p) */
558int q [ MAXDEGPOLY ] ; /* Current row of Q matrix. */
559
560int row = 0 ;
561
562
563/*------------------------------------------------------------------------------
564| Function Body |
565------------------------------------------------------------------------------*/
566
567/* Check for invalid inputs. */
568if (n < 2 || p < 2 || Q == (int **)0)
569{
570 return ;
571}
572
573
574
575/* Row 0 of Q = 1. */
576Q[ 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*/
583x_to_power( (bigint) p, xp, power_table, n, p ) ;
584
585memcpy( q, xp, n * sizeof( int ) ) ;
586memcpy( 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*/
594for (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
603for (row = 0 ; row < n ; ++row)
604{
605 Q[ row ][ row ] = mod( Q[ row ][ row ] - 1, p ) ;
606}
607
608return ;
609
610} /* ================== end of function generate_Q_matrix =================== */
611
612
613
614
615/*==============================================================================
616| find_nullity |
617================================================================================
618
619DESCRIPTION
620
621
622INPUT
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
628RETURNS
629
630 Nullity of Q. However if the nullity is 2 or more, we just return 2.
631
632EXAMPLE
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
681METHOD
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
717BUGS
718
719
720--------------------------------------------------------------------------------
721| Function Call |
722------------------------------------------------------------------------------*/
723
724int find_nullity( int ** Q, int n, int p )
725{
726
727int colFlag[ MAXDEGPOLY ] ; /* Is -1 if the column has no pivotal element. */
728int nullity = 0 ;
729int row, col ;
730int r ;
731int found = NO ;
732int pivotCol = -1 ;
733int q = 0 ;
734int t = 0 ;
735
736/* Initialize column flags to -1. */
737memset( colFlag, -1, n * sizeof( int ) ) ;
738
739/* Sweep through each row. */
740for (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
802return nullity ;
803
804} /* ===================== end of function find_nullity ===================== */
805
806
807#if 0
808int test
809 /* n 2n-2
810 Precompute the powers x , ..., x (mod f(x), p)
811 for use in all later computations.
812 */
813n = 4 ; p = 5 ;
814f[0] = 2 ; f[1] = 3 ; f[2] = 3 ; f[3] = 3 ; f[4] = 1 ;
815
816construct_power_table( power_table, f, n, p ) ;
817if (has_multi_irred_factors( power_table, n, p ) == 1)
818 printf( "Pass\n" ) ;
819else
820 printf( "Fail\n" ) ;
821
822#endif