I have the following system of equations,
1+x+y+z==0, 1+xy+yz+xz==0
which I want to solve in the extension field of GF(2) for example. There is a parametric solution of these equations in terms of the parameter s as x=1+s, y=1+ω s, z=1+ω2s where s is the parameter and ω2+ω+1=0. How can I modify the following for Sage to be able to output parametric solutions like this one?
R.<x,y,z> = PolynomialRing(GF(4))
I = R.ideal([1 + x + y + z, 1 + x*y + y*z + x*z])
I.variety()