```  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.
19 |
20 |     This program is free software: you can redistribute it and/or modify
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.
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"
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 ======================== */
```