Which is the next double?

[merse]

numeric_limits<double>::epsilon() gives the error of (double)1.0
in the sense that 1+epsilon is the next double value which is represented.
But how can we calculate the error of arbitrary double number?

[arkoenig]

Given a number x, you want to find epsilon with the following two properties:

1) x + epsilon != x

2) If 0 <= y < epsilon, x + y == x.

It should not be terribly hard to find the value of epsilon by binary search.

[merse]

Yes, but its very slow for me.
But an estimate is good for me.
We can use the known representation of double numbers.

[arkoenig]

If you use the known representation of double numbers, your program will work only on machines that use that representation. If you use the binary-search technique I suggested, it will run anywhere.

Admittedly the binary search is slower. What is it about this particular operation that you need to be able to do quickly?

[merse]

I want to use error handling, which means every number has a mean value and an absolute error, which can be a real error or error due to the representation.

[merse]

I think that basically the double representation is the same on every platform. Or the only variable is the number of used bytes. But this is not a serious problem. So, can anyone tell me a function which can be used as an estimate?

[merse]

I have tried the following code:

#include <iostream>
#include <limits>

using namespace std;

int main()
{
  for(double x = 1e-10; x < 1e11; x *= 10)
  {
    cout << x;

    double s = 0;

    for(double y = x; y == x; y = x+(s+=x*1e-20)) {} 
    
    cout << " " << s/x << endl; 
  }

  cout << "epsilon(): " << numeric_limits<double>::epsilon() << endl;
 
  return 0;
}

The output is:

1e-010 6.46e-017
1e-009 1.034e-016
1e-008 8.28e-017
1e-007 6.63e-017
1e-006 1.059e-016
1e-005 8.48e-017
0.0001 6.78e-017
0.001 1.085e-016
0.01 8.67e-017
0.1 6.95e-017
1 1.111e-016
10 8.89e-017
100 7.11e-017
1000 5.69e-017
10000 9.1e-017
100000 7.28e-017
1e+006 5.83e-017
1e+007 9.32e-017
1e+008 7.46e-017
1e+009 5.97e-017
1e+010 9.55e-017
epsilon(): 2.22045e-016

Why is the calculated value for 1 is only the half of epsilon()?

I think that a good estimation is error x = x*1e-16;