Is there a double number better than 1.0 for representing the real number 1 + 1e-16?

zell 0
Viewed 33

It is well known floating-point computation is inexact. For example:

In [64]: 1+1e-16
Out[64]: 1.0

In this case, I would like to know if there is a double number that is better than 1.0 for representing the real number 1+1e-16 (being better means closer to the real number)? In other words, I try to figure out whether this inaccuracy is due to the problem itself or due to how this is solved.

rounding-error ieee-754 floating-point
2 Answers

The 53-bit significand precision gives from 15 to 17 significant decimal digits precision (2^−53 ≈ 1.11 × 10^−16). [...] 2^0 × (1 + 2^−52) ~= 1.0000000000000002, the smallest number > 1 wikipedia.

Explored the following 1+e-17 == 1+1e-16 == 1 != 1+e-15 so I think you are right.

#include <math.h>
#include <stdio.h>

int main(void) {
    double d[] = {
        1e-15,
        pow(2, -52),
        1e-16,
        1e-17
    };
    for(size_t i = 0; i < sizeof(d) / sizeof(*d); i++) {
        printf("%d: %.20e, %d\n", i, 1 + d[i], (1 + d[i]) == 1.0);
    }
}

and the output is:

0: 1.00000000000000111022e+00, 0
1: 1.00000000000000022204e+00, 0
2: 1.00000000000000000000e+00, 1
3: 1.00000000000000000000e+00, 1
Allan Wind0

Given that your original question seems to be in Python, here's a few batteries the language gives you to investigate:

from decimal import Decimal as D
from math import nextafter  # Python 3.9+
from sys import float_info

# just do the naive thing, as you note this will be 1.0
naive = 1 + 1e-16

# calculate next floating point number in a couple of different ways
next_via_eps = 1 + float_info.epsilon
next_via_math = nextafter(1, 2)

# make sure they are the same
assert next_via_eps == next_via_math

# print out hex encoded values, to show they only differ by their LSB
print(naive.hex(), next_via_eps.hex())

# use the decimal module to get an exact answer
exact = 1 + D(1e-16)
# note that decimal isn't always exact, but the default precision (28 decimal digits) is enough here

# calculate differences
approx_orig_err = abs(D(naive) - exact)
approx_next_err = abs(D(next_via_eps) - exact)

# make sure that the FPU is doing the right thing
assert approx_orig_err < approx_next_err

note that general use of float_info.epsilon requires a multiplication to shuffle the precision out, e.g.:

def demo(x):
    return x + x * float_info.epsilon

assert demo(256) == nextafter(256, 1000)

it doesn't need to appear when calculating next_via_eps because multiplying by 1 is an identity

Sam Mason0
Related