I want to compute all solutions over $\mathbb{Z}_9$ in parameters $a,b,c$ for the equation
$$(ax+b\sqrt2+c)^2=1$$
where $x$ is such that $(x+\sqrt2)^2=2(x+\sqrt2)$.
I'm first defining the polynomial ring over $\mathbb{Z}_9$ in variables $a,b,c,x,y$, then factoring by the ideal generated by
$$y^2-2, (x+y)^2-2(x+y),$$
but then I don't know which command to use in order to get the solutions. I have tried "solve" and "variety", but they do not seem to work. The code up to this point is just
R.<x,y,a,b,c> = PolynomialRing(IntegerModRing(9),order='lex')
J= R.ideal(x^2-2,(x+y)^2-2*(x+y))
S=R.quotient(J)
Which function should I use?