Getting underestimated floating point factors

Question

Most physics engines support doing an object trajectory trace that will return a factor between 0.0 and 1.0 representing how far along an object's trajectory it will first hit an object.

The problem I'm concerned about is cases where moving the object by that factor of its trajectory results in its position being past the boundary it was supposed to hit and stop at (due to floating point rounding).

For instance, I created a C program that tried random cases until it ran into this issue, and came up with this example (although I've experienced ones with much less extreme movement, so it's not specific to large floating points):

float start     = 4884361.0f;
float wall      = 37961976.0f;
float end       = 1398674432.0f;
float time      = (wall - start) / (end - start);
float new_pos   = start + time * (end - start);
printf("We hit %f, but moving left us at %f.\n", wall, new_pos);

And this case prints out: We hit 37961976.000000, but moving left us at 37961980.000000.

So the position moved beyond the wall position and now the object is stuck inside the wall.

Is there a way to generate the factor or perform the factor multiplication such that the floating point error will always undershoot the actual value for all possible values?

Answer 1

The calculated value is the next (or nearly the next) floating point number. We are at the limits of float precision. To insure an answer is at or to one side of the expected answer, there are a number of approaches

1) Higher intermediate precision: (should come up with the right answer far more often)

float start     = 4884361.0f;
float wall      = 37961976.0f;
float end       = 1398674432.0f;
double time     = ((double)wall - start) / ((double)end - start);
float new_pos   = start + time * ((double)end - start);

2) Logical: (this absolutely will work)

if (new_pos > wall) new_pos = wall;

3) Use a slightly lower time value: (a gentle hack)

float new_pos = start + nextafterf(time,0.0f) * (end - start);

4) Change FP rounding mode to round toward zero: (this may have a large impact though)

fesetround(FE_TOWARDZERO); 

5) A simple factor:

static const float factor = 0.99999;
float new_pos = start + factor*time * (end - start);

Lots of pros & cons per approach.