/*============================================================================== | | File Name: | | Primpoly.c | | Description: | | Compute primitive polynomials of degree n modulo p. | | Functions: | | main Main driving routine. | | LEGAL | | Primpoly Version 13.0 - A Program for Computing Primitive Polynomials. | Copyright (C) 1999-2018 by Sean Erik O'Connor. All Rights Reserved. | | This program is free software: you can redistribute it and/or modify | it under the terms of the GNU General Public License as published by | the Free Software Foundation, either version 3 of the License, or | (at your option) any later version. | | This program is distributed in the hope that it will be useful, | but WITHOUT ANY WARRANTY; without even the implied warranty of | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | GNU General Public License for more details. | | You should have received a copy of the GNU General Public License | along with this program. If not, see <http://www.gnu.org/licenses/>. | | The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com | with the !DOT! replaced by . and the !AT! replaced by @ | ==============================================================================*/ /*------------------------------------------------------------------------------ | Include Files | ------------------------------------------------------------------------------*/ #include <stdio.h> #include <stdlib.h> #include <string.h> #include <math.h> #include "Primpoly.h" /*============================================================================== | main | ================================================================================ DESCRIPTION Program for finding a primitive polynomial of degree n modulo p for use in generating PN sequences and finite fields for error control coding. INPUT Run the program by typing $ Primpoly p n OUTPUT You will get an nth degree primitive polynomial modulo p. EXAMPLE CALLING SEQUENCE Let's find a primitive polynomial of degree 18, modulo the prime 11. After a few seconds, even on a slow PC, we get: $ primpoly -s 5 20 Primpoly Version 13.0 - A Program for Computing Primitive Polynomials. Copyright (C) 1999-2018 by Sean Erik O'Connor. All Rights Reserved. Primpoly comes with ABSOLUTELY NO WARRANTY; for details see the GNU General Public License. This is free software, and you are welcome to redistribute it under certain conditions; see the GNU General Public License for details. Factoring r = 23841857910156 into 2^2 3 11 13 41 71 521 9161 Total number of primitive polynomials = 1280348160000. Begin testing... Primitive polynomial modulo 5 of degree 20 x ^ 20 + x ^ 2 + 2 x + 3 +--------- Statistics ---------------------------- | | Total num. degree 20 polynomials mod 5 : 95367431640625 | Actually tested : 39 | Const. coeff. was primitive root : 16 | Free of linear factors : 5 | Irreducible or irred. to power : 3 | Had order r (x^r = integer) : 1 | Passed const. coeff. test : 1 | Had order m (x^m != integer) : 1 | +------------------------------------------------- METHOD Described in detail in my web pages at http://www.seanerikoconnor.freeservers.com http://seanerikoconnor.freeyellow.com (mirror site) under the topic Professional Interests->Generating Primitive Polynomials BUGS -------------------------------------------------------------------------------- | Function Call | ------------------------------------------------------------------------------*/ int main( int argc, char * argv[] ) { /*------------------------------------------------------------------------------ | Local Variables | ------------------------------------------------------------------------------*/ int n, /* Degree of the primitive polynomial f(x) to be found, n >= 2. */ p ; /* Modulo p arithmetic on the polynomial coefficients, p >= 2. */ bigint max_p_to_n = MAXPTON, /* Maximum value of p ^ n. */ max_num_poly, /* p ^ n, the number of polynomials to test for primitivity. */ num_poly = 0, /* Number of polynomials tested so far. */ r, /* The number (p ^ n - 1)/(p - 1). */ primes[ MAXNUMPRIMEFACTORS ], /* The distinct prime factors of r. */ prim_poly_count = 0, /* Counter for primitive polynomials found. */ num_prim_poly = 0 ; /* Total number of possible primitive polynomials. */ int count[ MAXNUMPRIMEFACTORS ], /* ... and their multiplicities. */ i, /* Prime index. */ prime_count, /* Primes are stored in array locations 0 through prime_count. */ f[ MAXDEGPOLY + 1 ], /* Coefficients of the polynomial f(x) which we test for primitivity. */ a = 0, /* Integer in the order r test. */ is_primitive_poly = NO, /* Equal to YES as soon as a primitive polynomial is found. */ num_free_of_linear_factors = 0,/* The number of polynomials which have no linear factors. */ num_const_coeff_prim_root = 0, /* Number of polynomials whose constant is a primitive root of p. */ num_passing_const_coeff_test = 0, /* Number of polynomials whose constant term passes a consistency check. */ num_irred_to_power = 0, /* Number of polynomials which are of the form irreducible poly to a power >= 1.*/ num_order_m = 0, num_order_r = 0, /* The number of polynomials which pass the order_m and order_r tests. */ stopTesting = NO, /* When to stop testing polynomials for primitivity. */ testPolynomialForPrimitivity = NO, /* Test a given input polnomial for primitivity? */ testPolynomial[ MAXDEGPOLY + 1 ], /* Coefficients of the test polynomial. */ listAllPrimitivePolynomials = NO, /* Print ALL primitive polynomials? */ printStatistics = NO, /* Print statistics? */ printHelp = NO, /* Print help information? */ selfCheck = NO, /* Do a self-check? Time consuming! */ /* x ^ n , ... , x ^ 2n-2 (mod f(x), p) */ power_table[ MAXDEGPOLY - 1 ] [ MAXDEGPOLY ] ; char outputFormat[ _MAX_PATH ] ; /* Formatting for printf's (used only when printing bigints) */ char * legalNotice = { "\n" "Primpoly Version 13.0 - A Program for Computing Primitive Polynomials.\n" "Copyright (C) 1999-2018 by Sean Erik O'Connor. All Rights Reserved.\n" "\n" "Primpoly comes with ABSOLUTELY NO WARRANTY; for details see the\n" "GNU General Public License. This is free software, and you are welcome\n" "to redistribute it under certain conditions; see the GNU General Public License\n" "for details.\n\n" } ; char * help = { "This program generates a primitive polynomial of degree n modulo p.\n\n" "Usage: primpoly p n\n\n" "Example: primpoly 2 4 \n" " generates the fourth degree polynomial\n\n" " x ^ 4 + x + 1, whose coefficients use modulo 2 arithmetic.\n\n" "Primitive polynomials find many uses in mathematics and communications \n" "engineering:\n" " * Generation of pseudonoise (PN) sequences for spread spectrum\n" " communications and chip fault testing.\n" " * Generation of CRC and Hamming codes.\n" " * Generation of Galois (finite) fields for use in decoding Reed-Solomon\n" " and BCH error correcting codes.\n\n" "Options:\n" " pp -c 2 4\n" " does an addtional time consuming double check on the primitivity.\n" " pp -s 2 4\n" " prints search statistics.\n" " pp -a 2 4\n" " lists ALL primitive polynomials of degree 4 modulo 2.\n" "\n\n" } ; /*------------------------------------------------------------------------------ | Function Body | ------------------------------------------------------------------------------*/ /* Show the legal notice first. */ printf( "%s", legalNotice ) ; /* Read the parameters p and n from the command line. Return with an error message if there are an incorrect number of inputs on the command line, or if p and n are out of bounds. */ parse_command_line( argc, argv, &testPolynomialForPrimitivity, &listAllPrimitivePolynomials, &printStatistics, &printHelp, &selfCheck, &p, &n, testPolynomial ) ; if (printHelp) { printf( "%s", help ) ; exit( 1 ) ; } if (p < 2) { printf( "ERROR: p must be 2 or more.\n\n" ) ; exit( 1 ) ; } if (n > MAXDEGPOLY || n < 2) { printf( "ERROR: n must be between 2 and %d\n\n", MAXDEGPOLY ) ; exit( 1 ) ; } /* Check to see if p is a prime. */ if (!is_almost_surely_prime( p )) { printf( "ERROR: p must be a prime number.\n\n" ) ; exit( 1 ) ; } /* n Return if p > MAXPTON because then we'll exceed the computer integer precision. */ if (n * log( (double) p ) > log( (double) (sbigint) max_p_to_n )) { printf( "ERROR: p to the nth power must be smaller than %llu\n\n", MAXPTON ) ; exit( 1 ) ; } /* n n p - 1 Compute p and r = ------. p - 1 */ max_num_poly = power( p, n ) ; r = (max_num_poly - 1) / (p - 1) ; /* Factor r into distinct primes. */ if (printStatistics) { sprintf( outputFormat, "%s%s%s", "\nFactoring r = ", bigintOutputFormat, " into\n " ) ; printf( outputFormat, r ) ; } prime_count = factor( r, primes, count ) ; if (printStatistics) { for (i = 0 ; i <= prime_count ; ++i) { if (count[ i ] == 1) { sprintf( outputFormat, "%s ", bigintOutputFormat ) ; printf( outputFormat, primes[ i ] ) ; } else { sprintf( outputFormat, "%s%s", bigintOutputFormat, "^%d " ) ; printf( outputFormat, primes[ i ], count[ i ] ) ; } } printf( "\n\n" ) ; } /* n Initialize f(x) to x + (-1). Then, when f(x) passes through function n next_trial_poly for the first time, it will have the correct value, x */ initial_trial_poly( f, n ) ; if (printStatistics || listAllPrimitivePolynomials) { sprintf( outputFormat, "%s%s%s", "Total number of primitive polynomials = ", bigintOutputFormat, ". Begin testing...\n\n" ) ; num_prim_poly = EulerPhi( power( p, n ) - 1 ) / n ; printf( outputFormat, num_prim_poly ) ; } /* Generate and test all possible n th degree, monic, modulo p polynomials f(x). A polynomial is primitive if passes all the tests successfully. */ do { next_trial_poly( f, n, p ) ; /* Try another polynomal. */ ++num_poly ; #ifdef DEBUG_PP_PRIMPOLY printf( "\nNext trial polynomial: " ) ; write_poly( f, n ) ; printf( "\n" ) ; #endif /* n 2n-2 Precompute the powers x , ..., x (mod f(x), p) for use in all later computations. */ construct_power_table( power_table, f, n, p ) ; /* Constant coefficient of f(x) * (-1)^n must be a primitive root of p. */ if (const_coeff_is_primitive_root( f, n, p )) { ++num_const_coeff_prim_root ; #ifdef DEBUG_PP_PRIMPOLY printf( "Coefficient of polynomial is primitive root.\n" ) ; #endif /* f(x) can't have any linear factors. */ if (!linear_factor( f, n, p )) { ++num_free_of_linear_factors ; #ifdef DEBUG_PP_PRIMPOLY printf( "Free of linear factors.\n" ) ; #endif /* f(x) can't have two or more distinct irreducible factors. */ if (!has_multi_irred_factors( power_table, n, p )) { ++num_irred_to_power ; #ifdef DEBUG_PP_PRIMPOLY printf( "Has one unique irreducible factor.\n" ) ; #endif /* x^r (mod f(x), p) = a must be an integer. */ if (order_r( power_table, n, p, r, &a )) { ++num_order_r ; #ifdef DEBUG_PP_PRIMPOLY printf( "Passes the order r test.\n" ) ; #endif /* Const coeff. of f(x)*(-1)^n must equal a mod p. */ if (const_coeff_test( f, n, p, a )) { ++num_passing_const_coeff_test ; #ifdef DEBUG_PP_PRIMPOLY printf( "Passes the constant coefficient test.\n" ) ; #endif /* x^m != integer for all m = r / q, q a prime divisor of r. */ if (order_m( power_table, n, p, r, primes, prime_count )) { ++num_order_m ; is_primitive_poly = YES ; #ifdef DEBUG_PP_PRIMPOLY printf( "Passes the order m tests.\n" ) ; #endif if (listAllPrimitivePolynomials) { printf( "\n\nPrimitive polynomial " ) ; sprintf( outputFormat, "%s of %s ", bigintOutputFormat, bigintOutputFormat ) ; printf( outputFormat, ++prim_poly_count, num_prim_poly ) ; printf( "modulo %d of degree %d\n\n", p, n ) ; write_poly( f, n ) ; printf( "\n\n" ) ; } } } /* end const coeff test */ } /* end order r */ } /* end can't determine if reducible */ } /* end no linear factors */ } /* end constant coefficient primitive. */ /* Stop when we've either checked all possible polynomials or we've not been asked to list all and found the first primtive one. */ stopTesting = (num_poly > max_num_poly) || (!listAllPrimitivePolynomials && is_primitive_poly) ; } while( !stopTesting ) ; printf( "\n\n" ) ; /* Report on success or failure. */ if (listAllPrimitivePolynomials) ; /* We're done */ else if (is_primitive_poly) { printf( "\n\nPrimitive polynomial modulo %d of degree %d\n\n", p, n ) ; write_poly( f, n ) ; printf( "\n\n" ) ; } else { printf( "Internal error: \n" "Tested all possible polynomials (%llu), but failed\n" "to find a primitive polynomial.\n" "Please let the author know by e-mail.\n", max_num_poly ) ; exit( 1 ) ; } /* Print the statistics of the primitivity tests. */ if (printStatistics) { printf( "+--------- Statistics -----------------------------------------------------------------\n" ) ; printf( "|\n" ) ; sprintf( outputFormat, "%s%s%s", "| Total num. degree %3d polynomials mod %3d : ", bigintOutputFormat, "\n" ) ; printf( outputFormat, n, p, max_num_poly ) ; sprintf( outputFormat, "%s%s%s", "| Actually tested : ", bigintOutputFormat, "\n" ) ; printf( outputFormat, num_poly ) ; printf( "| Const. coeff. was primitive root : %10d\n", num_const_coeff_prim_root ) ; printf( "| Free of linear factors : %10d\n", num_free_of_linear_factors ) ; printf( "| Irreducible or irred. to power : %10d\n", num_irred_to_power ) ; printf( "| Had order r (x^r = integer) : %10d\n", num_order_r ) ; printf( "| Passed const. coeff. test : %10d\n", num_passing_const_coeff_test ) ; printf( "| Had order m (x^m != integer) : %10d\n", num_order_m ) ; printf( "|\n" ) ; printf( "+--------------------------------------------------------------------------------------\n" ) ; } /* Confirm f(x) is primitive using a different, but extremely slow test for primitivity. Disabled when we list all primitive polynomials. */ if (selfCheck && !listAllPrimitivePolynomials) { printf( "\nConfirming polynomial is primitive with an independent check.\n" "Warning: You may wait an impossibly long time!\n\n" ) ; if (maximal_order( f, n, p )) printf( " -Polynomial is confirmed to be primitive.\n\n" ) ; else { printf( "Internal error: \n" "Primitive polynomial confirmation test failed.\n" "Please let the author know by e-mail.\n\n" ) ; return 1 ; } } return 0 ; } /* ========================== end of function main ======================== */
