[ODE] FPU precission problems - suggestion

[Michal Bacik]

I'd like share my experience with handling FPU limited precission problems
with authors of ODE.

There're cases, where ODE potentially may divide by zero, or generate
overflow errors (I haven't experienced this with ODE yet, but looking at
code, the danger is undoubtedly there).

For example, function dCollideSpheres:
//...
dReal d = dDISTANCE (p1,p2);
if (d <= 0) {
   //...
}else {
   dReal d1 = dRecip (d);
   c->normal[0] = (p1[0]-p2[0])*d1;
   //...
}
//...
where dRecip is defined as simple divide:
#define dRecip(x) ((float)(1.0f/(x)))  /* reciprocal */

What's the problem: the "if" statement checks for zero or less, but what if
two spheres are positioned such that the distance is +1.0e-42 for example?
The 'if' gets after 'else' statement, because the value is still positive
number, although very small. But 1.0 / 1.0e-42 makes result 1.0e+42, which
is not representable by float precission, and you get float overflow -
infinity number, by which the function multiplies just one line below. And
we're done - your normal gets nicely infinite.

Then there's another problem - denormal FPU numbers, which I think are
non-zero numbers, but that small, that they generate divide-by-zero
exception.

Solution: definition of a "near-zero" value, e.g.:
#define MRG_ZERO 1e-8f        //marginal zero - thresh value, under which
number may be considered zero
and compare all values to be above this before performing divide. The code
will be more robust and stable.

It's also a good thing to use FPU exceptions for overflow and divide-by-zero
(at least during development), and watch all crashes resulting from similar
precission problems. And possibly it'd be wise to look at all divides at ODE
and make sure they're protected against similar problems. I've been writing
collision detection system for a while, so I know the problems being in this
area.

Hope my suggestions just help to make ODE better :)

- michal