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