floating point number imprecision while iterating

Question

I have a function that computes a point in 3d spaced based on a value in range [0, 1]. The problem I'm facing is, that a binary floating-point number cannot represent exactly 1.

The mathematical expression that is evaluated in the function is able to compute a value for t=1.0, but the value will never be accepted by the function because it checks if for the range before computing.

curves_error curves_bezier(curves_PointList* list, curves_Point* dest, curves_float t) {
    /* ... */
    if (t < 0 || t > 1)
        return curves_invalid_args;
    /* ... */
    return curves_no_error;
}

How can I, with this function, compute the 3d point at t=1.0? I heard something about an ELLIPSIS some time ago that I think had to do with such an issue, but I'm not sure.

Thanks

EDIT: Ok, I'm sorry. I assumed a float cannot represent exactly 1, because of the issue I'm facing. The problem may be because I was doing an iteration like this:

for (t=0; t <= 1.0; t += 0.1) {
    curves_error error = curves_bezier(points, point, t);
    if (error != curves_no_error)
        printf("Error with t = %f.\n", t);
    else
        printf("t = %f is ok.\n", t);
}

Answer 1

for (t=0; t <= 1.0; t += 0.1) {

your problem is that a binary floating point number cannot exactly represent 0.1.

The closest 32-bit single precision IEEE754 floating point number is 0.100000001490116119384765625 and the closest 64-bit double precision one 0.1000000000000000055511151231257827021181583404541015625. If the arithmetic is performed strictly at 32-bit precision, the result of adding 0.1f ten times to 0 is

1.00000011920928955078125

If intermediate computations are performed at greater precision than float has, it could result in exactly 1.0 or even slightly smaller numbers.

To fix your problem, in this case you could use

for(k = 0; k <= 10; ++k) {
    t = k*0.1;

because 10 * 0.1f is exactly 1.0.

Another option is to use a small tolerance in your curves_bezier function,

if (t > 1 && t < 1 + epsilon) {
    t = 1;
}

for a suitably small epsilon, maybe float epsilon = 1e-6;.

Answer 2

binary floating-point number cannot represent exactly 1

Proof that it can can be found here.

Most accurate representation = 1.0E0

There could be problems with

1. fractional numbers that would have infinite fractional digits in radix of 2
2. numbers that are too small to exactly represent without losing precision
3. numbers that are too large to represent without losing precision.

But 1.0 is none of them!

However 0.1 is a problem case, violating point number 1, look at this:

Most accurate representation = 1.00000001490116119384765625E-1

So if you add up 0.1 ten times, you will get 1.00000001490116119384765625E-0 which is greater than 1.0.

(examples are in IEEE754 single precision 32 bit floating point numbers)

Possible solution:

int i;
for (i=0; i <= 10; i++) {
    t=i/10.0;
    curves_error error = curves_bezier(points, point, t);
    if (error != curves_no_error) {
        printf("Error with t = %f.\n", t);
    }
    else {
        printf("t = %f is ok.\n", t);
    }
}

This way, the error of the binary format does not get summed up!

(Note: I used extra curly braces for the if and else statements. Do that, you'll thank yourself one day.)

Answer 3

When comparing floating point numbers you should check if they are close enough not exactly equal, for the reasons mentioned in other answers, something like:

#define EPSILON 0.000001f
#define FEQUAL(a,b) (fabs((a) - (b)) < EPSILON)