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 : 39100 | Const. coeff. was primitive root : 16101 | Free of linear factors : 5102 | Irreducible or irred. to power : 3103 | Had order r (x^r = integer) : 1104 | Passed const. coeff. test : 1105 | Had order m (x^m != integer) : 1106 |107 +-------------------------------------------------108 109 110 METHOD111 112 Described in detail in my web pages at 113 http://www.seanerikoconnor.freeservers.com114 115 BUGS116 117 --------------------------------------------------------------------------------118 | Function Call |119 ------------------------------------------------------------------------------*/120 121 int main( int argc, char * argv[] )
122 {
123 124 /*------------------------------------------------------------------------------125 | Local Variables |126 ------------------------------------------------------------------------------*/127 128 int129 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 to137 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 int148 count[ MAXNUMPRIMEFACTORS ], /* ... and their multiplicities. */149 i, /* Prime index. */150 prime_count, /* Primes are stored in array locations 0151 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 constant168 term passes a consistency check. */169 num_irred_to_power = 0, /* Number of polynomials which are of the170 form irreducible poly to a power >= 1.*/171 num_order_m = 0,
172 num_order_r = 0, /* The number of polynomials which pass173 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 error238 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 n283 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 n295 n p - 1296 Compute p and r = ------.297 p - 1298 */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 else325 {
326 sprintf( outputFormat, "%s%s", bigintOutputFormat, "^%d " ) ;
327 printf( outputFormat, primes[ i ], count[ i ] ) ;
328 }
329 }
330 printf( "\n\n" ) ;
331 }
332 333 334 /* n335 Initialize f(x) to x + (-1). Then, when f(x) passes through function 336 n337 next_trial_poly for the first time, it will have the correct value, x338 */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 polynomials352 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_PRIMPOLY359 printf( "\nNext trial polynomial: " ) ;
360 write_poly( f, n ) ;
361 printf( "\n" ) ;
362 #endif363 364 /* n 2n-2365 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_PRIMPOLY377 printf( "Coefficient of polynomial is primitive root.\n" ) ;
378 #endif379 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_PRIMPOLY386 printf( "Free of linear factors.\n" ) ;
387 #endif388 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_PRIMPOLY395 printf( "Has one unique irreducible factor.\n" ) ;
396 #endif397 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_PRIMPOLY404 printf( "Passes the order r test.\n" ) ;
405 #endif406 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_PRIMPOLY413 printf( "Passes the constant coefficient test.\n" ) ;
414 #endif415 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_PRIMPOLY423 printf( "Passes the order m tests.\n" ) ;
424 #endif425 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 elseif (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 else511 {
512 printf( "Internal error: \n"513 "Primitive polynomial confirmation test failed.\n"514 "Please let the author know by e-mail.\n\n" ) ;
515 return1 ;
516 }
517 }
518 519 return0 ;
520 521 } /* ========================== end of function main ======================== */
