Precision limits of numpy floats?

Instead of adding 1, double the denominator. You can safely assume that it's some binary number. Here is one simple method:

one = np.float64(1.0) two = np.float64(2.0) n = one bits = 0 while one + n != one: bits += 1 n /= two 

You start with bits = 0 because otherwise you would get the count of bits that took you past resolution.

In the end, you get bits = 53, which is the number of bits in an IEEE-754 encoded 64-bit floating point number.

That means that for any number, which is encoded in what is effectively binary scientific notation, the ULP (unit of least precision) is approximately n * 2**-53. Specifically, where n is the number rounded to its highest bit. You won't be able to resolve smaller relative changes in a float.

Bonus: Running the above code for the other floating point types gives:

float16 (half): 11 bits float32 (single): 24 bits float64 (double): 53 bits float96 (sometimes longdouble): 80 bits float128 (when available): 113 bits 

You can modify the code above to work for any target number:

target = np.float16(0.0004883) one = np.float16(1.0) two = np.float16(2.0) n = two**(np.floor(np.log2(target)) - one) bits = 0 while target + n != target: bits += 1 n /= two 

The result (ULP) is given by n * 2 since the loop stops after you lose resolution. This is the same reason we start with bits = 0. In this case:

>>> n * two 5e-07 

You can entirely short circuit the computation if you know the number of bits in the mantissa up front. So for float16, where bits = 11, you can do

>>> two**(np.floor(np.log2(target)) - np.float16(bits)) 5e-07 

Read more here:

• https://en.wikipedia.org/wiki/IEEE_754
• https://en.wikipedia.org/wiki/Unit_in_the_last_place

My other answer provides the theory behind what you are actually asking, but needs some non-trivial interpretation. Here is the missing step:

Given some integer i, you can write

1 / i - 1 / (i + 1) = (i + 1 - i) / (i * (i + 1)) = 1 / (i * (i + 1)) = 1 / (i**2 + i) 

To find an i such that 1 / (i**2 + i) falls below ULP of 1 / i in some binary representation, you can use my other answer pretty directly.

The ULP of 1 / i is given by

ulp = 2**(floor(log2(1 / i)) - (bits + 1)) 

You can try to find an i such that

1 / (i**2 + i) < 2**(floor(log2(1 / i)) - (bits + 1)) 1 / (i**2 + i) < 2**floor(log2(1 / i)) / 2**(bits + 1) 2**(bits + 1) < (i**2 + i) * 2**floor(log2(1 / i)) 

It's not trivial as written because of the floor operation, and Wolfram Alpha runs out of time. Since I am cheap and don't want to buy Mathematica, let's just approximate:

2**(bits + 1) < (i**2 + i) * 2**floor(log2(1 / i)) 2**(bits + 1) < (i**2 + i) / i 2**(bits + 1) < i + 1 

You may be off by one or so, but you should see that around i = 2**(bits + 1) - 1, the difference stops being resolvable. And indeed for the 11 bit mantissa of float16, we see:

>>> np.float16(1 / (2**12 - 1)) - np.float16(1 / (2**12)) 0.0 

The actual number is a tiny bit less here (remember that approximation where we took away the floor):

>>> np.float16(1 / (2**12 - 5)) - np.float16(1 / (2**12 - 4)) 0.0 >>> np.float16(1 / (2**12 - 6)) - np.float16(1 / (2**12 - 5)) 2.4e-07 

As you noted in your comments, i is

>>> 2**12 - 6 4090 

You can compute the exact values for all the other floating point types in a similar manner. But that is indeed left as an exercise for the reader.