LEGAL
PrimpolySageMath Version 16.2 - A Program for Computing Primitive Polynomials using Sage Math.
Copyright (C) 1999-2023 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 seanerikoconnor!AT!gmail!DOT!com
with the !DOT! replaced by . and the !AT! replaced by @
def numPolynomials(p,n):
""""The total number of polynomials of degree n modulo p."""
return p ^ n
def numPrimitivePolynomials(p,n):
""""The number of primitive polynomials of degree n modulo p."""
return euler_phi( p ^ n - 1 ) / n
def probabilityRandomPolynomialIsPrimitive(p,n):
""""Probability that a polynomial of degree n modulo p chosen at random is primitive."""
return numPrimitivePolynomials( p, n ) / numPolynomials( p, n )
p = 5 ; n = 4
numPolynomials(p,n)
625
numPrimitivePolynomials(p,n)
48
N(probabilityRandomPolynomialIsPrimitive(p,n),digits=3)
0.0768
def isPrimRoot(a,p):
"""If the smallest n for which a^n = 1 is p-1, then a is a primitive root of p."""
if Mod(a,p) == 0:
return False
elif Mod(a,p).multiplicative_order() == p-1:
return True
else:
return False
isPrimRoot(3,p)
True
[Mod(3^i,p) for i in range(1,p)]
[3, 4, 2, 1]
[Mod(4^i,p) for i in range(1,p)]
[4, 1, 4, 1]
Mod(4,p).multiplicative_order()
2
R.<x>=PolynomialRing(GF(p))
R
Univariate Polynomial Ring in x over Finite Field of size 5
def numPolynomialsWithLinearFactors(p,n):
""""The number of polynomials of degree n modulo p which have one or more linear factors."""
if n >= p:
return p^n - ((p-1)^p) * p^(n-p)
else:
var('i') # Define i to be a symbolic variable.
return sum( -((-1)^i) * binomial( p, i ) * p^(n-i), i, 1, n)
def probRandomPolyHasLinearFactors(p,n):
"""" The probability that a polynomial of degree n modulo p has one or more linear factors."""
if n >= p:
return 1 - (1 - 1 / p)^p
else:
return numPolynomialsWithLinearFactors(p,n) / (p^n)
def polynomialHasLinearFactor(f):
"""" A test whether a polynomial of degree n modulo p has one or more linear factors."""
hasRoot=false
for i in range(0,p):
if f(i) == 0:
return True
return False
f = x^4 + x^2 + 2*x + 3 ; pretty_print(f)
numPolynomialsWithLinearFactors(p,n)
420
N(probRandomPolyHasLinearFactors(p,n),digits=3)
0.672
polynomialHasLinearFactor(f)
False
[f(k) for k in range(0,p)]
[3, 2, 2, 4, 3]
polynomialHasLinearFactor((x-2)*(x-3))
True
def isConstantPolynomial(f):
"""Is the polynomial constant?"""
# Get the polynomial coefficients and their number. For example, if
# f = x^4 + x^2 + 2*x + 3
# then f.list() = [3, 2, 1, 0, 1] and f.list()[0] = 3
coeffList = f.list()
numCoeff = len( coeffList )
# Polynomial has only a constant term.
if numCoeff == 1:
return True
else:
# Otherwise, do the powers x, x^2, ... have even one nonzero coefficient?
hasNonZeroPowers = False
for k in range(1, numCoeff):
if coeffList[ k ] != 0:
hasNonZeroPowers = True
break
# Nope, coefficients of the powers x, x^2, ... are all zero; polynomial is a constant.
return not hasNonZeroPowers
help(isConstantPolynomial)
Help on function isConstantPolynomial in module __main__: isConstantPolynomial(f) Is the polynomial constant?
isConstantPolynomial( f )
False
isConstantPolynomial(3)
True
Q.<z>=R.quotient( f )
Q
Univariate Quotient Polynomial Ring in z over Finite Field of size 5 with modulus x^4 + x^2 + 2*x + 3
def powersOfXModfp(f,n,p):
"""Compute the powers {z^0, z^p, z^2p, ..., z^(np)} modulo f(x) and p."""
# Locally define the quotient ring Q of polynomials modulo f(x) and p.
R.<x>=PolynomialRing(GF(p))
Q.<z>=R.quotient( f )
# Return the list of all powers.
return [z ^ (p * k) for k in range(0,n)]
listOfPowers = powersOfXModfp(f,n,p) ; pretty_print( listOfPowers )
def generateQIMatrix(f,n,p):
"""Generate the Q - I matrix. Rows of Q are coefficients of the powers z^0, ..., z^(np) modulo f(x), p"""
# Locally define the quotient ring Q of polynomials modulo f(x) and p.
R.<x>=PolynomialRing(GF(p))
Q.<z>=R.quotient( f )
# Compute the powers {z^0, z^p, z^2p, ..., z^(np)} modulo f(x) and p.
# e.g. for f = x^4 + x^2 + 2*x + 3 modulo p=5,
# listOfPowers = [1, 4*z^3 + 3*z^2 + 2*z, z^3 + 4*z^2 + 3, 3*z^3 + 4*z^2 + 4*z]
listOfPowers = powersOfXModfp(f,n,p)
# e.g. row 1 is listOfPowers[1].list() = [0, 2, 3, 4]
rows = [listOfPowers[k].list() for k in range(0,n)]
# Define a matrix space of n x n matrices modulo p.
M = MatrixSpace(GF(p),n,n)
# Load the rows.
m = M.matrix( rows )
# Q-I
return m - matrix.identity(n)
m = generateQIMatrix(f,n,p) ; pretty_print( m )
k = kernel(m) ; k
Vector space of degree 4 and dimension 1 over Finite Field of size 5 Basis matrix: [1 0 0 0]
m.left_nullity()
1
def hasTwoOrMoreDistinctIrreducibleFactors(f,n,p):
m=generateQIMatrix(f,n,p)
return m.left_nullity() >= 2
hasTwoOrMoreDistinctIrreducibleFactors(f,n,p)
False
r=(p^n-1)//(p-1) ; r
156
factorization = factor(r) ; pretty_print( factorization )
type( 156 / 4 )
<class 'sage.rings.rational.Rational'>
type( 156 // 4 )
<class 'sage.rings.integer.Integer'>
def distinctPrimeFactors(r):
# Factorization object. e.g. for r = 156, factor(r) = 2^2 * 3 * 13
factorization = factor(r)
# List of prime factors to powers. e.g. list(factor(r)) = [(2, 2), (3, 1), (13, 1)]
listOfPrimesToPowers = list(factorization)
# Pull out the prime factors only. e.g. [2, 3, 13]
distinctPrimeFactors = [listOfPrimesToPowers[k][0] \
for k in range(0,len(listOfPrimesToPowers))]
return distinctPrimeFactors
distinctPrimeFactors(r)
[2, 3, 13]
def isPrimitivePolynomial( f, n, p ):
"""Determine if the nth degree polynomial f(x) modulo p is primitive."""
isPrimitive = False
# Polynomial must be at least degree 2 modulo p >= 2
if p < 2 or n < 2:
return False
print( f"Passed step 1: p={p:d} >=2 and n = {n} >= 2")
a0=f[0]
if isPrimRoot( a0 * (-1)^n, p ):
print( f"Passed step 2: a0 = {a0} is a primitive root")
if not polynomialHasLinearFactor(f):
print( f"Passed step 3: f(x) = {f} has no linear factors" )
if not hasTwoOrMoreDistinctIrreducibleFactors(f,n,p):
print( f"Passed step 4: f(x) = {f} is either irreducible or irreducible ^ power" )
R.<x>=PolynomialRing(GF(p))
Q.<z>=R.quotient( f )
r=(p^n-1)//(p-1)
a=z^r
if isConstantPolynomial(a):
print( f"Passed step 5: a = {a} is constant" )
if a == Mod( a0 * (-1)^n, p):
print( f"Passed step 6: a = a0 (-1)^n mod p where a0 = {a0}" )
# Tentatively say it's primitive.
isPrimitive = True
primes = distinctPrimeFactors(r)
print( f"Begin step 7: Testing on prime factors p{0} ... p{len(primes)-1}" )
for k in range( 0 , len( primes )):
# Skip the test for kth prime pk if (p-1) | pk
if Mod( (p-1), primes[k] ) == 0:
print( f" Skip step 7 test since prime p{k} = {primes[k]} divides p-1 = {p-1}")
else:
# Oops! Failed the test for prime pk.
if isConstantPolynomial( z^(r//primes[k]) ):
print( f" Failed step 7 since x^m is an integer for prime p{k} = {primes[k]}")
isPrimitive = False
else:
print( f" Passed step 7 since x^m is not an integer for prime p{k} = \
{primes[k]}")
return isPrimitive
p = 5 ; n = 4
f = x^4 + x^2 + 2*x + 2 ; pretty_print( f )
isPrimitivePolynomial( f, n, p )
Passed step 1: p=5 >=2 and n = 4 >= 2 Passed step 2: a0 = 2 is a primitive root Passed step 3: f(x) = x^4 + x^2 + 2*x + 2 has no linear factors Passed step 4: f(x) = x^4 + x^2 + 2*x + 2 is either irreducible or irreducible ^ power Passed step 5: a = 2 is constant Passed step 6: a = a0 (-1)^n mod p where a0 = 2 Begin step 7: Testing on prime factors p0 ... p2 Skip step 7 test since prime p0 = 2 divides p-1 = 4 Passed step 7 since x^m is not an integer for prime p1 = 3 Passed step 7 since x^m is not an integer for prime p2 = 13
True
def generateAllNonzeroFieldElements( f, n, p ):
"""Generate all the nonzero finite field elements for GF(p^n) given polynomial f(x),
which is not necessarily primitive."""
# Define the quotient ring Q of polynomials modulo f(x) and p.
R.<x>=PolynomialRing(GF(p))
Q.<z>=R.quotient( f )
# A primitive polynomial f generates the nonzero elements of the field:
# z^m = 1 for m = p^n - 1 but z^m != 1 for 1 <= m <= p^n - 2
fieldElements = [z^m for m in range(1,p^n)]
return fieldElements
fieldElements = generateAllNonzeroFieldElements( f, n, p ) ; pretty_print( fieldElements )
len(fieldElements)
624
p^n-1
624
f = x^4 + x^2 + 2 ; pretty_print( f )
isPrimitivePolynomial( f, n, p )
Passed step 1: p=5 >=2 and n = 4 >= 2 Passed step 2: a0 = 2 is a primitive root Passed step 3: f(x) = x^4 + x^2 + 2 has no linear factors Passed step 4: f(x) = x^4 + x^2 + 2 is either irreducible or irreducible ^ power Passed step 5: a = 2 is constant Passed step 6: a = a0 (-1)^n mod p where a0 = 2 Begin step 7: Testing on prime factors p0 ... p2 Skip step 7 test since prime p0 = 2 divides p-1 = 4 Passed step 7 since x^m is not an integer for prime p1 = 3 Failed step 7 since x^m is an integer for prime p2 = 13
False
fieldElements = generateAllNonzeroFieldElements( f, n, p ); pretty_print( fieldElements )