Implementation of the Randomized Cubic Primality Test in Python

This post provides a complete, runnable Python implementation of the randomized cubic primality test using explicit Frobenius computation in a cubic ring, as described in my article “A Randomized Cubic Primality Test Using Explicit Frobenius Computation”.

The test operates in the ring

$$
R_N = (\mathbb{Z}/N\mathbb{Z})[x]/(f(x)),\qquad f(x)=x^3-qx-q,
$$

where $N$ is the integer under test and $q$ is a carefully chosen auxiliary prime of the form

$$
q=m^2+m+7,\qquad m\ge 0.
$$

Let $\alpha$ be the image of $x$ in $R_N$. Every element of $R_N$ can be written uniquely as

$$
c_0 + c_1\alpha + c_2\alpha^2,
$$

and in the code it is represented by the coefficient triple $(c_0,c_1,c_2)$ with coefficients reduced modulo $N$.

The implementation uses the gmpy2 library for fast arbitrary precision arithmetic and primality testing. It includes

* validity checks for a round $(N,q)$
* ring arithmetic in $R_N$, multiplication and exponentiation
* the core cubic test comparing the computed Frobenius image $\alpha^N$ to the predicted conjugate
* a simple random search function to hunt for mismatches, useful when looking for potential pseudoprimes.

Python

import random
import gmpy2
from gmpy2 import mpz


def aa(m):
    m = mpz(m)
    return m * m + m + 7


def is_square(n):
    return bool(gmpy2.is_square(mpz(n)))


def is_cube(n):
    n = mpz(n)
    if n < 0:
        return False
    root, exact = gmpy2.iroot(n, 3)
    return bool(exact)


def ispseudoprime(n):
    return bool(gmpy2.is_prime(mpz(n)))


def inv_mod(a, n):
    a = mpz(a) % n
    inv = gmpy2.invert(a, n)
    if inv == 0:
        raise ZeroDivisionError("no modular inverse")
    return inv


def mul(u, v, q, N):
    c0 = (u[0] * v[0]) % N
    c1 = (u[0] * v[1] + u[1] * v[0]) % N
    c2 = (u[0] * v[2] + u[1] * v[1] + u[2] * v[0]) % N
    c3 = (u[1] * v[2] + u[2] * v[1]) % N
    c4 = (u[2] * v[2]) % N

    dq = (c4 * q) % N
    c2 = (c2 + dq) % N
    c1 = (c1 + dq) % N

    dq = (c3 * q) % N
    c1 = (c1 + dq) % N
    c0 = (c0 + dq) % N

    return (c0, c1, c2)


def pow_ring(base, exp, q, N):
    res = (mpz(1), mpz(0), mpz(0))
    exp = mpz(exp)
    while exp > 0:
        if exp & 1:
            res = mul(res, base, q, N)
        base = mul(base, base, q, N)
        exp >>= 1
    return res


def valid(N, q):
    N = mpz(N)
    q = mpz(q)

    if N <= 2:
        return False
    if (N & 1) == 0:
        return False
    if is_square(N):
        return False
    if is_cube(N):
        return False

    if not gmpy2.is_prime(q):
        return False

    D = 4 * q - 27
    if not gmpy2.is_square(D):
        return False

    if (q - 1) % 3 != 0:
        return False

    if gmpy2.gcd(N, q) != 1:
        return False

    a = gmpy2.isqrt(D)
    if gmpy2.gcd(a, N) != 1:
        return False

    c = gmpy2.powmod(N % q, (q - 1) // 3, q)
    if c <= 1:
        return False

    return True


def cubictest(N, q):
    N = mpz(N)
    q = mpz(q)

    a = gmpy2.isqrt(4 * q - 27)
    d = a * q

    inv6 = inv_mod(6, q)
    s1 = ((-a - 3) * inv6) % q
    s2 = ((a - 3) * inv6) % q

    c = gmpy2.powmod(N % q, (q - 1) // 3, q)

    inv_dN = inv_mod(d, N)
    C2 = (3 * q * inv_dN) % N
    C1 = ((-d - 9 * q) * inv_mod(2 * d, N)) % N
    C0 = (-2 * q * q * inv_dN) % N

    beta = (C0, C1, C2)
    gamma = ((-C0) % N, (-C1 - 1) % N, (-C2) % N)

    alpha = (mpz(0), mpz(1), mpz(0))
    frob = pow_ring(alpha, N, q, N)

    if c == s1:
        predicted = gamma
    elif c == s2:
        predicted = beta
    else:
        return False

    return frob == predicted


def search(M):
    for _ in range(M):
        N = mpz(random.randrange(10**9))
        rm = random.randrange(10**3)
        q = aa(rm)

        if not valid(N, q):
            continue

        ct = cubictest(N, q)
        psp = ispseudoprime(N)

        if ct != psp:
            print("N =", N, "q =", q, "cubictest =", int(ct), "ispseudoprime =", int(psp))
            return (N, q, rm, ct, psp)

    return None




if __name__ == "__main__":
    print("sanity valid 1997,7 =", valid(1997, 7))
    print("sanity cubictest 1997,7 =", cubictest(1997, 7))
    print("search result =", search(1000000))