Ray tracing implicit surfaces using Sturm's theorem

Question

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?

Answer 1

This is a rather late reply but it might interest others who come up to this question. Your question assumes that your ray tracing platform does not support double precision. This is unlikely.

If you are using CUDA or OpenCL they both support double precision vectors.

If you are using Vulkan things are a little more complicated due to specification saying: "The precision of double-precision instructions is at least that of single precision."

Though not a ray tracing platform OpenGL also supports double precision vectors from OpenGL 4.0 and on.

If you are tracing on CPU, chances are you can define a double precision vector data type, if it is not already defined.

---

If, in the unlikely case that your ray tracing platform does not support double precision, or you are obliged to use single precision in any way, the only thing that I can think of is using pseudo-remainder sequences for rem- i.e. euclidean division. Particularly you might want to check out Subresultant pseudo-remainder sequence from the wikipedia section in the link. I am adding the pseudocode of the algorithm in here:

r_0 = a
r_1 = b
for (i:= 1;r_i not equal to 0; i := i+1)
    d_i := degree_of(r_{i-1}) - degree_of(r_i)
    g_i := leading_coefficient_of(r_i)
    if i = 1 then
        b_1 := (-1)^{d_1 + 1}
        p_1 := -1
    else
        b_i := -g_{i-1} (p_i)^d_i
        p_i := ((-g_{i-1})^d_{i-1}) / (p_{i-1}^{d_{i-1}-1}
    end if
    r_i := rem((g_i)^{d_i+1} r_{i-1}, r_i) / b_i
end for