### Solving polynomial equations with Groebner basis in $\mathbb{R}$

I am trying to solve the following system of polynomial equations:
$$x^2+y^2+z^2=4\
x^2+2y^2=5\
xz=1$$

I used the following command to get a Groebner basis

```
P.<x,y,z> = PolynomialRing(QQ,order='lex')
PList = [x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1]
I = ideal(PList)
B = I.groebner_basis(); B
[x + 2*z^3 - 3*z, y^2 - z^2 - 1, z^4 - 3/2*z^2 + 1/2]
```

Now I want to use B.subs() to plug the solutions for $z$ to solve for $x,y$. It worked well for $z=\pm 1$.

```
B.subs(z=1)
[x - 1, y^2 - 2, 0]
```

but not for the $1/\sqrt{2}$ since this is in polynomial ring $\mathbb{Q}[x,y,z]$. If I change the original QQ to RR, the numbers become all decimals.

```
B.subs(z=1/sqrt(2))
[x - 1.41421356237310, y^2 - 1.50000000000000, 0]
```

How can I still get exact solutions when using RR? I know I can compute this example by hands, but I'd like to know how to use Sage to solve this kind of problems.

Thank you for any help!