Floating-point and ieee-754

Question

I wrote such a code to find the grid step(difference between two adjacent machine-representable numbers, machine epsilon the difference between the minimum machine-representable number greater than one and one). How correct is this program?

#include <stdio.h>
int main(void)
{
    long double x,eps=1.0,a; 
    scanf("%Le",&x);
    if(x>=1.8446744e+19) eps=x;
    do
    {
        a=eps;
        eps/=2.0;
    }
    while( x+eps>x);
    printf("Grid step: %Le",(long double)a);
    return 0;
}

Answer 1

There are a number of things I'd suggest looking at.

One is that there's no difference between double and long double in Visual C++, but that's a minor thing. Still, it's worth knowing that you're getting 64-bit IEEE binary floating point format.

Read Eric Postpischil's comment. The step you're looking for is always a power of 2 but it might be greater than, less than or equal to 1. You can start out at 1 and use two while loops. The first will halve the step until adding the step to the original number leaves it unchanged.

    long double eps = 1.0;
    while (original + eps != original)
        eps *= 0.5;

Multiplying by 0.5 is exactly the same as dividing by 2, only faster on most machines.

The second loop will double the step until adding it to the original does make a difference. Same pattern as above, changing the != to == and 0.5 to 2.0. The second loop will always run at least once, if only to back out the last change to eps and get to the power of 2 that does make a difference.

An edge case to watch out for: If the original is 0 then the step is going to underflow to 0 during the first loop. The simplest way to handle this is to include <float.h> and test for the original being zero. If so, return LDBL_MIN as the result before the first loop.