All the solutions of a polynomial system over a finite ring

edit

In fact, my custom implementation of variety() method from other answer works fine here, provided two things:

• a custom implementation of roots() method that supports polynomials with the content non-coprime to the ring modulus;
• a correct implementation of interreduced_basis()method avoiding the bug https://trac.sagemath.org/ticket/33982

It finds all 1152 solutions within a few seconds.

import itertools
from sage.rings.finite_rings.integer_mod_ring import crt as mycrt

def myinterreduced_basis(J):
    B = list(J.gens())
    R = J.ring()
    while True:
        updated = False
        for i in range(len(B)):
            r = R.ideal(list(set(B)-{B[i]})).reduce(B[i])
            if r!=B[i]:
                B[i] = r
                updated = True
        if not updated:
            return B
        B = [b for b in B if b!=0]

def myroots(p):
    R = p.base_ring()
    if not isinstance(R,sage.rings.abc.IntegerModRing):
        yield from p.roots(multiplicities=False)
        return

    m = R.zero().modulus()
    pz = p.change_ring(ZZ)
    g = gcd(gcd(pz.coefficients()),m)
    mg = m//g

    roots = (pz/g).roots(ring=Integers(mg),multiplicities=False)  # roots modulo mg = m/g

    # embedding roots modulo mg to modulo m=mg*g
    for rr in itertools.product( *( [Mod(r*q^valuation(mg,q),q^k) for r in range(q^valuation(g,q))] for q,k in factor(m) ) ):
        r2 = mycrt(rr)
        yield from (r2 + ZZ(r) for r in roots)

In fact, computing the interreduced basis of G immediately eliminates quadratic terms, and so the problem can be reduced solving linear equations right away as explained in the comments. This way it takes just a fraction of a second to get all the solutions.