1/*==============================================================================
  2| 
  3|  File Name:     
  4|
  5|     Primpoly.c
  6|
  7|  Description:   
  8|
  9|     Compute primitive polynomials of degree n modulo p.
 10| 
 11|  Functions:
 12|
 13|     main      Main driving routine.
 14|
 15|  LEGAL
 16|
 17|     Primpoly Version 16.4 - A Program for Computing Primitive Polynomials.
 18|     Copyright (C) 1999-2025 by Sean Erik O'Connor.  All Rights Reserved.
 19|
 20|     This program is free software: you can redistribute it and/or modify
 21|     it under the terms of the GNU General Public License as published by
 22|     the Free Software Foundation, either version 3 of the License, or
 23|     (at your option) any later version.
 24|
 25|     This program is distributed in the hope that it will be useful,
 26|     but WITHOUT ANY WARRANTY; without even the implied warranty of
 27|     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 28|     GNU General Public License for more details.
 29|
 30|     You should have received a copy of the GNU General Public License
 31|     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 32|     
 33|     The author's address is seanerikoconnor!AT!gmail!DOT!com
 34|     with the !DOT! replaced by . and the !AT! replaced by @
 35|
 36==============================================================================*/
 37
 38/*------------------------------------------------------------------------------
 39|                                Include Files                                 |
 40------------------------------------------------------------------------------*/
 41
 42#include <stdio.h>
 43#include <stdlib.h>
 44#include <string.h>
 45#include <math.h>
 46
 47#include "Primpoly.h"
 48
 49
 50/*==============================================================================
 51|                                    main                                      |
 52================================================================================
 53
 54DESCRIPTION
 55
 56     Program for finding a primitive polynomial of degree n modulo p for use in
 57     generating PN sequences and finite fields for error control coding.
 58
 59INPUT 
 60
 61     Run the program by typing
 62
 63         $ Primpoly p n
 64
 65OUTPUT
 66
 67     You will get an nth degree primitive polynomial modulo p.
 68
 69EXAMPLE CALLING SEQUENCE
 70
 71    Let's find a primitive polynomial of degree 18, modulo the prime 11.
 72    After a few seconds, even on a slow PC, we get:
 73
 74        $ primpoly -s 5 20
 75
 76        Primpoly Version 16.4 - A Program for Computing Primitive Polynomials.
 77        Copyright (C) 1999-2025 by Sean Erik O'Connor.  All Rights Reserved.
 78
 79        Primpoly comes with ABSOLUTELY NO WARRANTY; for details see the
 80        GNU General Public License.  This is free software, and you are welcome
 81        to redistribute it under certain conditions; see the GNU General Public License
 82        for details.
 83
 84
 85        Factoring r = 23841857910156 into
 86            2^2 3 11 13 41 71 521 9161
 87
 88        Total number of primitive polynomials = 1280348160000.  Begin testing...
 89
 90
 91        Primitive polynomial modulo 5 of degree 20
 92
 93         x ^ 20   +    x ^  2   +    2 x   +    3
 94
 95
 96        +--------- Statistics ----------------------------
 97        |
 98        | Total num. degree 20 polynomials mod  5 :  95367431640625
 99        | Actually tested :                          39
100        | Const. coeff. was primitive root :         16
101        | Free of linear factors :                    5
102        | Irreducible or irred. to power :            3
103        | Had order r (x^r = integer) :               1
104        | Passed const. coeff. test :                 1
105        | Had order m (x^m != integer) :              1
106        |
107        +-------------------------------------------------
108
109
110METHOD
111
112     Described in detail in my web pages at 
113           http://www.seanerikoconnor.freeservers.com 
114
115BUGS
116
117--------------------------------------------------------------------------------
118|                                Function Call                                 |
119------------------------------------------------------------------------------*/
120
121int main( int argc, char * argv[] )
122{
123
124/*------------------------------------------------------------------------------
125|                               Local Variables                                |
126------------------------------------------------------------------------------*/
127
128int 
129    n,      /*  Degree of the primitive polynomial f(x) to be found, n >= 2.  */
130
131    p ;     /*  Modulo p arithmetic on the polynomial coefficients, p >= 2.   */
132
133bigint
134    max_p_to_n = MAXPTON,          /* Maximum value of p ^ n.               */
135
136    max_num_poly,                  /* p ^ n, the number of polynomials to
137                                      test for primitivity.                 */
138
139    num_poly = 0,                  /* Number of polynomials tested so far.  */
140
141    r,                             /* The number (p ^ n - 1)/(p - 1).       */
142
143    primes[ MAXNUMPRIMEFACTORS ],  /* The distinct prime factors of r.      */
144    prim_poly_count = 0,           /* Counter for primitive polynomials found.        */
145    num_prim_poly = 0 ;            /* Total number of possible primitive polynomials. */
146
147int
148    count[ MAXNUMPRIMEFACTORS ],   /* ... and their multiplicities.         */
149    i,                             /* Prime index. */
150    prime_count,                   /* Primes are stored in array locations 0
151                                      through prime_count.                  */
152
153    f[ MAXDEGPOLY + 1 ],           /* Coefficients of the polynomial f(x)  
154                                      which we test for primitivity.        */
155
156    a = 0,                         /* Integer in the order r test.          */
157
158    is_primitive_poly = NO,        /* Equal to YES as soon as a primitive 
159                                      polynomial is found.                  */
160
161    num_free_of_linear_factors = 0,/* The number of polynomials which have  
162                                      no linear factors.                    */
163
164    num_const_coeff_prim_root = 0, /* Number of polynomials whose constant  
165                                      is a primitive root of p.             */
166
167    num_passing_const_coeff_test = 0, /* Number of polynomials whose constant
168                                         term passes a consistency check.   */
169    num_irred_to_power = 0,        /* Number of polynomials which are of the
170                                      form irreducible poly to a power >= 1.*/
171    num_order_m = 0,
172    num_order_r = 0,               /* The number of polynomials which pass
173                                      the order_m and order_r tests.                    */
174    stopTesting = NO,              /* When to stop testing polynomials for primitivity. */
175    testPolynomialForPrimitivity = NO, /* Test a given input polnomial for primitivity? */
176    testPolynomial[ MAXDEGPOLY + 1 ],  /* Coefficients of the test polynomial.          */
177    listAllPrimitivePolynomials  = NO, /* Print ALL primitive polynomials?              */
178    printStatistics              = NO, /* Print statistics?                             */
179    printHelp                    = NO, /* Print help information?                       */
180    selfCheck                    = NO, /* Do a self-check?  Time consuming!             */
181
182    /*  x ^ n , ... , x ^ 2n-2 (mod f(x), p) */
183    power_table[ MAXDEGPOLY - 1 ] [ MAXDEGPOLY ] ;
184
185
186char outputFormat[ _MAX_PATH ] ; /* Formatting for printf's (used only when printing bigints) */
187
188char * legalNotice = 
189{
190    "\n"
191    "Primpoly Version 16.4 - A Program for Computing Primitive Polynomials.\n"
192    "Copyright (C) 1999-2025 by Sean Erik O'Connor.  All Rights Reserved.\n"
193   "\n"
194    "Primpoly comes with ABSOLUTELY NO WARRANTY; for details see the\n"
195    "GNU General Public License.  This is free software, and you are welcome\n"
196    "to redistribute it under certain conditions; see the GNU General Public License\n"
197    "for details.\n\n"
198}  ;
199
200char * help =
201{
202     "This program generates a primitive polynomial of degree n modulo p.\n\n"
203         "Usage:    primpoly p n\n\n"
204         "Example:  primpoly 2 4 \n"
205	 "          generates the fourth degree polynomial\n\n" 
206         "          x ^ 4 + x + 1, whose coefficients use modulo 2 arithmetic.\n\n"
207
208     "Primitive polynomials find many uses in mathematics and communications \n"
209     "engineering:\n"
210     "   * Generation of pseudonoise (PN) sequences for spread spectrum\n"
211     "     communications and chip fault testing.\n"
212     "   * Generation of CRC and Hamming codes.\n"
213     "   * Generation of Galois (finite) fields for use in decoding Reed-Solomon\n"
214     "     and BCH error correcting codes.\n\n"
215
216     "Options:\n"
217     "   pp -c 2 4\n"
218     "       does an addtional time consuming double check on the primitivity.\n"
219     "   pp -s 2 4\n"
220     "       prints search statistics.\n"
221     "   pp -a 2 4\n"
222     "       lists ALL primitive polynomials of degree 4 modulo 2.\n"
223     "\n\n"
224} ;
225
226/*------------------------------------------------------------------------------
227|                                Function Body                                 |
228------------------------------------------------------------------------------*/
229
230
231/*  Show the legal notice first.  */
232
233printf(  "%s", legalNotice )  ;  
234
235
236/*  
237     Read the parameters p and n from the command line.  Return with an error
238     message if there are an incorrect number of inputs on the command line,
239     or if p and n are out of bounds.
240*/
241parse_command_line( argc, 
242                    argv,
243                    &testPolynomialForPrimitivity,
244                    &listAllPrimitivePolynomials,
245                    &printStatistics,
246                    &printHelp,
247                    &selfCheck,
248                    &p,
249                    &n,
250                    testPolynomial ) ;
251
252if (printHelp)
253{
254    printf(  "%s", help )  ;  
255
256    exit( 1 ) ;
257}
258
259if (p < 2)
260{
261    printf( "ERROR:  p must be 2 or more.\n\n" ) ;
262    exit( 1 ) ;
263}
264
265if (n > MAXDEGPOLY || n < 2)
266{
267    printf( "ERROR: n must be between 2 and %d\n\n", MAXDEGPOLY ) ;
268    exit( 1 ) ;
269}
270
271
272/*  Check to see if p is a prime. */
273if (!is_almost_surely_prime( p ))
274{
275    printf( "ERROR:  p must be a prime number.\n\n" ) ;
276    exit( 1 ) ;
277}
278
279
280
281/*
282               n
283    Return if p  > MAXPTON because then we'll exceed the computer integer precision.
284*/
285if (n * log( (double) p ) > log( (double) (sbigint) max_p_to_n ))
286{
287    printf( "ERROR:  p to the nth power must be smaller than %llu\n\n", MAXPTON ) ;
288    exit( 1 ) ;
289}
290
291
292
293/*
294                         n
295              n         p  - 1
296     Compute p  and r = ------.
297                        p - 1
298*/
299max_num_poly = power( p, n ) ;
300
301r = (max_num_poly - 1) / (p - 1) ;
302
303
304
305
306/*  Factor r into distinct primes. */
307if (printStatistics)
308{
309    sprintf( outputFormat, "%s%s%s", "\nFactoring r = ", bigintOutputFormat, " into\n    " ) ;
310    printf( outputFormat, r ) ;
311}
312
313prime_count = factor( r, primes, count ) ;
314
315if (printStatistics)
316{
317    for (i = 0 ;  i <= prime_count ;  ++i)
318    {
319        if (count[ i ] == 1)
320        {
321            sprintf( outputFormat, "%s ", bigintOutputFormat ) ;
322            printf(  outputFormat, primes[ i ] ) ;
323        }
324        else
325        {
326            sprintf( outputFormat, "%s%s", bigintOutputFormat, "^%d " ) ;
327            printf(  outputFormat, primes[ i ], count[ i ] ) ;
328        }
329    }
330    printf( "\n\n" ) ;
331}
332
333
334/*                       n
335     Initialize f(x) to x  + (-1).  Then, when f(x) passes through function 
336                                                                          n
337     next_trial_poly for the first time, it will have the correct value, x
338*/
339initial_trial_poly( f, n ) ;
340
341if (printStatistics || listAllPrimitivePolynomials)
342{
343    sprintf( outputFormat, "%s%s%s", "Total number of primitive polynomials = ", bigintOutputFormat, ".  Begin testing...\n\n" ) ;
344    num_prim_poly = EulerPhi( power( p, n ) - 1 ) / n ;
345    printf( outputFormat, num_prim_poly ) ;
346}
347
348
349
350/*
351     Generate and test all possible n th degree, monic, modulo p polynomials
352     f(x).  A polynomial is primitive if passes all the tests successfully.
353*/
354do {
355    next_trial_poly( f, n, p ) ;      /* Try another polynomal. */
356    ++num_poly ;
357
358    #ifdef DEBUG_PP_PRIMPOLY
359    printf( "\nNext trial polynomial:  " ) ;
360    write_poly( f, n ) ;
361    printf( "\n" ) ;
362    #endif
363
364    /*                         n         2n-2
365        Precompute the powers x ,  ..., x     (mod f(x), p)
366        for use in all later computations.
367    */
368    construct_power_table( power_table, f, n, p ) ;
369
370
371    /* Constant coefficient of f(x) * (-1)^n must be a primitive root of p. */
372    if (const_coeff_is_primitive_root( f, n, p ))
373    {
374        ++num_const_coeff_prim_root ;
375
376        #ifdef DEBUG_PP_PRIMPOLY
377        printf( "Coefficient of polynomial is primitive root.\n" ) ;
378        #endif
379
380        /* f(x) can't have any linear factors. */
381        if (!linear_factor( f, n, p ))
382        {
383            ++num_free_of_linear_factors ;
384
385            #ifdef DEBUG_PP_PRIMPOLY
386            printf( "Free of linear factors.\n" ) ;
387            #endif
388
389            /* f(x) can't have two or more distinct irreducible factors. */
390            if (!has_multi_irred_factors( power_table, n, p ))
391            {
392                ++num_irred_to_power ;
393
394                #ifdef DEBUG_PP_PRIMPOLY
395                printf( "Has one unique irreducible factor.\n" ) ;
396                #endif
397
398                /* x^r (mod f(x), p) = a must be an integer. */
399                if (order_r( power_table, n, p, r, &a ))
400                {
401                    ++num_order_r ;
402
403                     #ifdef DEBUG_PP_PRIMPOLY
404                     printf( "Passes the order r test.\n" ) ;
405                     #endif
406
407                     /*  Const coeff. of f(x)*(-1)^n must equal a mod p. */
408                     if (const_coeff_test( f, n, p, a ))
409                     {
410                         ++num_passing_const_coeff_test ;
411
412                         #ifdef DEBUG_PP_PRIMPOLY
413                         printf( "Passes the constant coefficient test.\n" ) ;
414                         #endif
415
416                         /*  x^m != integer for all m = r / q, q a prime divisor of r. */
417                         if (order_m( power_table, n, p, r, primes, prime_count ))
418                         {
419                             ++num_order_m ;
420                             is_primitive_poly = YES ;
421
422                             #ifdef DEBUG_PP_PRIMPOLY
423                             printf( "Passes the order m tests.\n" ) ;
424                             #endif
425
426                             if (listAllPrimitivePolynomials)
427                             {
428                                 printf( "\n\nPrimitive polynomial " ) ;
429                                 sprintf( outputFormat, "%s of %s ", bigintOutputFormat, bigintOutputFormat ) ;
430                                 printf(  outputFormat, ++prim_poly_count, num_prim_poly ) ;
431                                 printf( "modulo %d of degree %d\n\n", p, n ) ;
432                                 write_poly( f, n ) ;
433                                 printf( "\n\n" ) ;
434                             }
435                         }
436                     } /* end const coeff test */
437                } /* end order r */
438            } /* end can't determine if reducible */
439        } /* end no linear factors */
440    } /* end constant coefficient primitive. */
441
442    /* Stop when we've either checked all possible polynomials or 
443       we've not been asked to list all and found the first primtive one.  
444    */
445    stopTesting = (num_poly > max_num_poly) || 
446                  (!listAllPrimitivePolynomials && is_primitive_poly) ;
447
448} while( !stopTesting ) ;
449
450printf( "\n\n" ) ;
451
452
453/*
454     Report on success or failure.
455*/
456
457if (listAllPrimitivePolynomials)
458    ; /* We're done */
459else if (is_primitive_poly)
460{
461    printf( "\n\nPrimitive polynomial modulo %d of degree %d\n\n", 
462            p, n ) ;
463    write_poly( f, n ) ;
464    printf( "\n\n" ) ;
465}
466else {
467
468    printf( "Internal error:  \n"
469            "Tested all possible polynomials (%llu), but failed\n"
470            "to find a primitive polynomial.\n"
471            "Please let the author know by e-mail.\n",
472            max_num_poly ) ;
473    exit( 1 ) ;
474}
475
476
477/*  Print the statistics of the primitivity tests. */
478
479if (printStatistics)
480{
481    printf( "+--------- Statistics -----------------------------------------------------------------\n" ) ;
482    printf( "|\n" ) ;
483    sprintf( outputFormat, "%s%s%s", "| Total num. degree %3d polynomials mod %3d :    ", bigintOutputFormat, "\n" ) ;
484    printf( outputFormat, n, p, max_num_poly ) ;
485    sprintf( outputFormat, "%s%s%s", "| Actually tested :                              ", bigintOutputFormat, "\n" ) ;
486    printf( outputFormat,  num_poly ) ;
487    printf( "| Const. coeff. was primitive root :      %10d\n",  num_const_coeff_prim_root ) ;
488    printf( "| Free of linear factors :                %10d\n",  num_free_of_linear_factors ) ;
489    printf( "| Irreducible or irred. to power :        %10d\n",  num_irred_to_power ) ;
490    printf( "| Had order r (x^r = integer) :           %10d\n",  num_order_r ) ;
491    printf( "| Passed const. coeff. test :             %10d\n",  num_passing_const_coeff_test ) ;
492    printf( "| Had order m (x^m != integer) :          %10d\n",  num_order_m ) ;
493    printf( "|\n" ) ;
494    printf( "+--------------------------------------------------------------------------------------\n" ) ;
495}
496
497
498
499/*  Confirm f(x) is primitive using a different, but extremely slow test for 
500    primitivity.  Disabled when we list all primitive polynomials.
501*/
502if (selfCheck && !listAllPrimitivePolynomials)
503{
504
505    printf( "\nConfirming polynomial is primitive with an independent check.\n"
506            "Warning:  You may wait an impossibly long time!\n\n" ) ;
507
508    if (maximal_order( f, n, p ))
509        printf( "    -Polynomial is confirmed to be primitive.\n\n" ) ;
510    else
511    {
512        printf( "Internal error:  \n"
513                "Primitive polynomial confirmation test failed.\n"
514                "Please let the author know by e-mail.\n\n" ) ;
515        return 1 ;
516    }
517}
518
519return 0 ;
520
521} /* ========================== end of function main ======================== */