Member Avatar for 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?

  • 2 Contributors
  • 6 Replies
  • 196 Views
  • 1 Day Discussion Span
  • Latest Post Latest Post by merse

Dani AI

Generated

As 's experiment shows, the spacing between adjacent IEEE‑754 doubles depends on the exponent. The standard constant often quoted — numeric_limits<double>::epsilon() — is the gap at 1.0 (2^-52). The smallest increment that will actually change a rounded result is commonly half that (2^-53), because ordinary round‑to‑nearest uses a half‑ulp threshold (this is why measured roughly 0.5*epsilon). 's binary‑search idea is portable but slow; two practical alternatives are shown below: one exact (portable and simple) and one fast (uses the floating‑point layout).

#include <cmath>
#include <limits>

double ulp_nextafter(double x) {
  if (!std::isfinite(x)) return std::numeric_limits<double>::quiet_NaN();
  double nx = std::nextafter(x, std::numeric_limits<double>::infinity());
  return nx - x;                    // exact distance to next representable double
}

double ulp_frexp(double x) {
  x = std::fabs(x);
  if (x == 0.0) return std::numeric_limits<double>::denorm_min();
  int exp;
  std::frexp(x, &exp);              // x = m * 2^exp, 0.5 <= m < 1
  return std::ldexp(1.0, exp - std::numeric_limits<double>::digits);
}

Notes and practical tips:

  • ulp_nextafter gives the exact spacing including subnormals and handles edge cases. It is the most robust replacement for a binary search.
  • ulp_frexp is very fast and gives the spacing for normal numbers via the exponent; it assumes standard binary floating format (IEEE‑754) and needs special handling for zero/subnormals.
  • For rounding/error bounds use the unit roundoff u = 0.5epsilon (≈2^-53) as the typical relative rounding bound; for absolute tolerances use ulp(x) or `kulp(x)` (k = 1..3) rather than a fixed decimal fraction like 1e-16.
  • Be careful with infinities/NaNs and when comparing values near zero — subnormals make the spacing much larger relatively.

Recommended Answers

All 6 Replies

Member Avatar for 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.

Member Avatar for merse

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

Member Avatar for 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?

Member Avatar for 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.

Member Avatar for 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?

Member Avatar for 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;

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.