How are Galois fields represented in SymPy? I couldn't find any documentation for this online, but SymPy contains a module called "galoistools", so I thought I should give it a try. I tried the following experiment:
from sympy import *
x = symbols("x")
A = [LC(Poly(i*x, modulus=8) * Poly(j*x, modulus=8)) for i in range(1, 8) for j in range(1, i+1)]
B = [LC(Poly(i*x, domain=GF(8)) * Poly(j*x, domain=GF(8))) for i in range(1, 8) for j in range(1, i+1)]
However, the resulting lists A and B are identical, so I'm obviously misunderstanding how this is supposed to be used. I'm trying to work in GF(8), i.e. GF(2^3), which is not the same as computing modulo 8.
At present SymPy does not have support for finite fields other than Z/pZ. The existing class GF(n) is misleadingly named; it actually implements Z/nZ as you observed.
However, using the low-level routines already present in SymPy's galoistools module, one can create a class for general finite fields GF(p^n) and for polynomials over such a field: see this answer where these classes are implemented (for the purpose of computing an interpolating polynomial, but they can be used for other things too). This is just a minimal class; it does not interface with advanced polynomial manipulation methods that are implemented in SymPy.