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