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 ======================== */
