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?


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.