Finding Primitive Polynomials

LEGAL

PrimpolySageMath Version 16.2 - A Program for Computing Primitive Polynomials using Sage Math.
Copyright (C) 1999-2025 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 @

This notebook is called PrimpolySageMath.ipynb

Run it using

jupyter lab PrimpolySageMath.ipynb

or first launch the SageMath app by clicking on it, or from the command line with

sage --notebook

But then you need to make sure to select SageMath as the kernel:

Kernel ➤ Change Kernel ➤ Select Kernel = SageMath 9.5

When you get done, please create an HTML version of this notebook for viewing on the web by running

jupyter nbconvert PrimpolySageMath.ipynb --to html

I urge you to do Kernel->Restart Kernel and Clear All Outputs , followed by Kernel->Restart Kernel and Run All Cells

Counting the number of primitive polynomials of degree $n$ modulo $p$. The probability that a polynomial chosen at random is primitive.

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 )

Let's do a low degree example with a small modulus.

p = 5 ; n = 4

numPolynomials(p,n)

numPrimitivePolynomials(p,n)

N(probabilityRandomPolynomialIsPrimitive(p,n),digits=3)

Testing whether $a$ is a primitive root of $p$

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)

Since 3 is a primitive root of 5, it generates the nonzero elements of GF(5).

[Mod(3^i,p) for i in range(1,p)]

But 4 is not a primitive root of 5 since its order is only 2 instead of 4.

[Mod(4^i,p) for i in range(1,p)]

Mod(4,p).multiplicative_order()

Define the ring of polynomials with coefficients over $GF(p)$ with generator $x$.

R.<x>=PolynomialRing(GF(p))

R

Testing if polynomials have linear factors.

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

An example polynomial of low degree.

f = x^4 + x^2 + 2*x + 3 ; pretty_print(f)

numPolynomialsWithLinearFactors(p,n)

N(probRandomPolyHasLinearFactors(p,n),digits=3)

polynomialHasLinearFactor(f)

Confirm that f has no linear factors.

[f(k) for k in range(0,p)]

Confirm that this polynomial with two roots, i.e. two linear factors, does indeed pass the one or more linear factors test.

polynomialHasLinearFactor((x-2)*(x-3))

Testing if a polynomial is constant

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)

isConstantPolynomial( f )

isConstantPolynomial(3)

Define the quotient ring $Q$ of polynomials modulo $f(x)$ with generator $z$.

We previously defined the variable x to denote the generator for the polynomial ring R with coefficients in $GF(p)$. Sage is smart enough to remember that.

Q.<z>=R.quotient( f )

Q

Compute the powers $z^n (mod f(x), p)$

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 )

Irreducibility Testing

Compute $\mathbf{Q-I}$ where the $n$ rows of matrix $\mathbf{Q}$ are the coefficients of the powers $z^0, z^p, ..., z^{np} (mod f(x), p)$

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 )

This is the left kernel (nullspace), i.e. all vectors w satisfying $\mathbf{w A = 0}$

k = kernel(m) ; k

m.left_nullity()

def hasTwoOrMoreDistinctIrreducibleFactors(f,n,p):
    m=generateQIMatrix(f,n,p)
    return m.left_nullity() >= 2

hasTwoOrMoreDistinctIrreducibleFactors(f,n,p)

Factor the expression r into distinct primes.

If you interrupt the kernel while factoring, you can see that SageMath uses [https://pari.math.u-bordeaux.fr/dochtml/html/Arithmetic_functions.html#se:factorint](PARI's integer factorization) which uses the following methods:

• Factors the integer n into a product of pseudoprimes (see ispseudoprime), using a combination of the Shanks SQUFOF and Pollard Rho method (with modifications due to Brent), Lenstra's ECM (with modifications by Montgomery), and MPQS (the latter adapted from the LiDIA code with the kind permission of the LiDIA maintainers), as well as a search for pure powers.

r=(p^n-1)//(p-1) ; r

factorization = factor(r) ; pretty_print( factorization )

# p = 2 ; n = 929 ; r=(p^n-1)//(p-1) ; pretty_print( r )
# factorization = factor(r) ; pretty_print( factorization )

NOTE: Use // for integer divide instead of / or else we'll get a rational number type for the answer even if $(p-1) | (p^n-1)$.

type( 156 / 4 )

type( 156 // 4 )

Distinct prime factors of r

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)

Testing whether a polynomial $f(x)$ modulo $p$ is primitive.

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

Examples: Primitive and non-primitive polynomials as generators of field elements.

A primitive polynomial generates all of the $p^n - 1$ nonzero field elements.

p = 5 ; n = 4

f = x^4 + x^2 + 2*x + 2 ; pretty_print( f )

isPrimitivePolynomial( f, n, p )

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)

p^n-1

A non-primitive polynomial won't generate all the nonzero field elements. Note the repetitions in the polynomal list below.

f = x^4 + x^2 + 2 ; pretty_print( f )

isPrimitivePolynomial( f, n, p )

fieldElements = generateAllNonzeroFieldElements( f, n, p ); pretty_print( fieldElements )

Finding Conway polynomials.

Conway polynomials and the related pseudo-Conway polynomials are a type of primitive polynomial with additional nice properties.

https://doc.sagemath.org/html/en/reference/finite_rings/sage/rings/finite_rings/conway_polynomials.html

from sage.rings.finite_rings.conway_polynomials import PseudoConwayLattice

def findConwayPolynomial( p, n ):
    """Conway polynomials and the related pseudo-Conway polynomials are a type of primitive polynomial with additional nice properties."""
    try:
        poly = conway_polynomial(p,n)
        print( f"Conway polynomial of degree n = {n} modulo p = {p}\n{poly}" )
        return poly
    except Exception as e:
        print( f"Couldn't find a Conway polynomial of degree n = {n} modulo p = {p}:  {e}")
        return None

def findPseudoConwayPolynomial( p, n ):        
    """Conway polynomials and the related pseudo-Conway polynomials are a type of primitive polynomial with additional nice properties."""
    PCL = PseudoConwayLattice(p, use_database=False)
    try:
        poly=PCL.polynomial(n)
        print( f"pseudo-Conway polynomial of degree n = {n} modulo p = {p}\n{poly}" )
        return poly
    except Exception as e:
        print( f"Couldn't find a pseudo-Conway polynomial of degree n = {n} modulo p = {p}:  {e}")
        return None

p = 2 ; n = 600

poly = findConwayPolynomial( p, n )

%time poly = findPseudoConwayPolynomial( p, n )

%time isPrimitivePolynomial( poly, n, p )

Try a smaller degree polynomial to see if we can find a Conway polynomial.

n = 100

%time poly = findConwayPolynomial( p, n )

isPrimitivePolynomial( poly, n, p )

This size of polynomal was running for several hours, so I aborted the computation.

p = 2 ; n = 929 %time poly = findPseudoConwayPolynomial( p, n )