[ODE] What is the ODE approach to denormals?

[Jon Watte]

I tracked down an annoying problem this week-end (there's a code 
diff at the bottom of this message). The cause was that the 
internals of ODE ended up with denormal numbers, and those 
denormals were then fed to the bottom part of the division 
operator, resulting in infinity, and later, 0 * result was added 
to a variable. As everyone knows, 0 * infinity == infinity. This 
led to my capsule-based player (sitting inside a GeomTransform 
to put him straight in my Y-up world) walking straight through a 
box (just hanging in a slider joint in the world).

There are two problems with denormals:

a) they cause FPUs to run up to 1000x slower than they do with 
   regular floating point values.

b) they cause incorrectness and/or surprising results, such as 
   players walking through boxes.

There are a few ways of dealing with denormals:

1) whenever data is entered, or persistent intermediate values 
   are generated, detect denormals and set them to 0.

2) every so often (every frame, or every 10 frames, or whatever), 
   examine the entire state of the simulation and nuke any denormal 
   values to 0.

3) add a jitter noise to all simulation inputs, making denormals 
   extremely unlikely to occur.

4) wherever data is examined, if denormals matter, do the right 
   thing, which may include branching, or treating as 0, or using 
   epsilons in the code.

In digital audio, where I first fought denormals, a) is the big 
problem, and it's usually dealt with using 3) (or, sometimes, 2) 
for performance reasons.

In physics simulation, we're concerned both about a) and b). 
However, I don't think either of solutions 1) through 3) are 
really that suited for simulation; in fact, I tried solution 1) 
first, and it didn't actually fix the problem, because denormals 
were generated as part of the simulation, not as an input.

Thus, I believe that all internal physics code should be written 
to properly deal with denormals. This is, sadly, quite tricky, and 
a bit of an unknown art, which has shrouded the whole issue and 
clouded it with "just add epsilons" webs of sorcery.

You can change code to do things like:

  if (fabsf(x) < EPSILON) {
    x = 0;
  }

However, this may cause extra branching if your compiler isn't 
smart enough to use fcomi and fcmov. Evenso, it adds clutter to 
the code. Luckily, in the specific case I ran into, there was a 
comparision already, that could be suitably changed, leading to 
a zero run-time cost of the fix.

Okay, enough soap-boxing. The specific code diff is available at 
http://www.mindcontrol.org/~hplus/ode/ (the diff date is 040426). 
This fixes capsules-in-transforms-walking-through-suspended-boxes. 
It does it by changing a comparision:

  if (v[i] > 0) {

into:

  if (v[i] > XXX) {

where "XXX" depends on dSINGLE/dDOUBLE, according to the following 
comment:

-    if (v[i] > 0) {
+    //  Denormals are a problem, because we divide by v[i], and then 
+    //  multiply that by 0. Alas, infinity times 0 is infinity (!)
+    //  We also use v2[i], which is v[i] squared. Here's how the epsilons 
+    //  are chosen:
+    //  float epsilon = 1.175494e-038 (smallest non-denormal number)
+    //  double epsilon = 2.225074e-308 (smallest non-denormal number)
+    //  For single precision, choose an epsilon such that v[i] squared is 
+    //  not a denormal; this is for performance.
+    //  For double precision, choose an epsilon such that v[i] is not a 
+    //  denormal; this is for correctness.
+#if defined( dSINGLE )
+    if (v[i] > dReal(1e-19))
+#else
+    if (v[i] > dReal(1e-307))
+#endif
+    {



Is this the right approach? Should we document it somewhere?

Cheers,

			/ h+