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+$\omega$ s, z=1+$\omega^2$s where $s$ is the parameter and $\omega^2+\omega+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()