# Ray tracing implicit surfaces using Sturm's theorem

I am trying to render implicit surfaces given by a polynomial equation $$F(x,y,z) = 0$$ using Sturm's theorem.

Plugging the parametric ray equation $$r(t) = o + td$$, where $$o\in\mathbb{R}^3$$ is the ray origin and $$d\in\mathbb{R}^3\backslash\{0\}$$ is the ray direction, into $$F$$ yields a polynomial equation $$p(t) = 0$$ for the intersections of the ray with the surface. I noticed that the straightforward implementation of calculating a Sturm sequence as

$$p_0 = p, \\ p_1 = p', \\ p_{i} = -\text{rem}(p_{i-2},p_{i-1}) \quad\forall i\ge2$$

is numerically unstable. When rendering a torus for example, the equation for one particular ray in my program was

$$p(t) := t^4 + 0.00329589844t^2 - 660.766052t + 9627.17676 = 0.$$

In single precision, the associated Sturm sequence is computed as

$$p_0(t) = t^4 + 0.00329589844t^2 - 660.766052t + 9627.17676,\\ p_1(t) = 4t^3+0.00659179688t-660.766052,\\ p_2(t) = -0.001647949220t^2+495.5745390t-9627.17676,\\ p_3(t) = -3.617114498\cdot 10^{11}t+7.027166964\cdot 10^{12},\\ p_4(t) = 0,$$

which indicates that (up to multiplicity) $$p$$ has exactly one real root. But in fact $$p$$ has no real roots, which can easily be verified with a computer algebra system. If the same computation is done with higher precision, it turns out that $$p_4(t)\approx 0.000039>0$$. In this case the correct number of roots is calculated.

Is there anything that can be done to avoid this problem when calculating a Sturm sequence?