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.1 - A Program for Computing Primitive Polynomials.
 18 |     Copyright (C) 1999-2021 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 artificer!AT!seanerikoconnor!DOT!freeservers!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 
 54 DESCRIPTION
 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 
 59 INPUT 
 60 
 61      Run the program by typing
 62 
 63          $ Primpoly p n
 64 
 65 OUTPUT
 66 
 67      You will get an nth degree primitive polynomial modulo p.
 68 
 69 EXAMPLE 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.1 - A Program for Computing Primitive Polynomials.
 77         Copyright (C) 1999-2021 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 
110 METHOD
111 
112      Described in detail in my web pages at 
113            http://www.seanerikoconnor.freeservers.com 
114 
115 BUGS
116 
117 --------------------------------------------------------------------------------
118 |                                Function Call                                 |
119 ------------------------------------------------------------------------------*/
120 
121 int main( int argc, char * argv[] )
122 {
123 
124 /*------------------------------------------------------------------------------
125 |                               Local Variables                                |
126 ------------------------------------------------------------------------------*/
127 
128 int
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 
133 bigint
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 
147 int
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 
186 char outputFormat[ _MAX_PATH ] ; /* Formatting for printf's (used only when printing bigints) */
187 
188 char * legalNotice =
189 {
190     "\n"
191     "Primpoly Version 16.1 - A Program for Computing Primitive Polynomials.\n"
192     "Copyright (C) 1999-2021 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 
200 char * 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 
233 printf(  "%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 */
241 parse_command_line( argc,
242                     argv,
243                     &testPolynomialForPrimitivity,
244                     &listAllPrimitivePolynomials,
245                     &printStatistics,
246                     &printHelp,
247                     &selfCheck,
248                     &p,
249                     &n,
250                     testPolynomial ) ;
251 
252 if (printHelp)
253 {
254     printf(  "%s", help )  ;
255 
256     exit( 1 ) ;
257 }
258 
259 if (p < 2)
260 {
261     printf( "ERROR:  p must be 2 or more.\n\n" ) ;
262     exit( 1 ) ;
263 }
264 
265 if (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. */
273 if (!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 */
285 if (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 */
299 max_num_poly = power( p, n ) ;
300 
301 r = (max_num_poly - 1) / (p - 1) ;
302 
303 
304 
305 
306 /*  Factor r into distinct primes. */
307 if (printStatistics)
308 {
309     sprintf( outputFormat, "%s%s%s", "\nFactoring r = ", bigintOutputFormat, " into\n    " ) ;
310     printf( outputFormat, r ) ;
311 }
312 
313 prime_count = factor( r, primes, count ) ;
314 
315 if (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 */
339 initial_trial_poly( f, n ) ;
340 
341 if (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 */
354 do {
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 
450 printf( "\n\n" ) ;
451 
452 
453 /*
454      Report on success or failure.
455 */
456 
457 if (listAllPrimitivePolynomials)
458     ; /* We're done */
459 else 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 }
466 else {
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 
479 if (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 */
502 if (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 
519 return 0 ;
520 
521 } /* ========================== end of function main ======================== */