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?


1 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
        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
  • 2
    $\begingroup$ I just came across this question that I asked three years ago and saw the above answer with an upvote. This answer is misleading unfortunately. I am not using this site anymore so I can't write comments or change the original question, but I would like to point out that it is not about single or double precision. The problem is that calculating a Sturm sequence with remainder sequences is inherently numerically unstable and this problem persists with any finite floating point precision. $\endgroup$
    – user17241
    Sep 21, 2021 at 19:35
  • 1
    $\begingroup$ There are still many visible artifacts when rendering implicit surfaces with this algorithm in double precision and there will still be artifacts when rendering with even greater floating point precision. Hence, the above question is asking for different methods / algorithms to compute a Sturm sequence that don't have these numerical problems. $\endgroup$
    – user17241
    Sep 21, 2021 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.