Power Operator Fails with Decimal Exponent

Greetings,

I’m trying to calculate 2**64.5 in Python. The answer I get is 26087635650665566208.00000000 when it should be 26087635650665564424.6991436125… Is this a limitation with decimal exponents?

Thanks,
Bassam Abdul-Baki

No, it is a limitation of the floating point data type, which has only 64 bits.

We can see the same issue with floating point calculations even when the
exponent is a whole number:

# print all 34 digits of the float calculation
print('%.34g' % (3.0**70))
# -> prints 2503155504993241614088572898377728

print(3**70)  # uses exact integer calculation
# -> prints 2503155504993241601315571986085849

When you calculate 2**64.5 the calculation is done using binary floats with 64 bits, which is only enough storage for about 17 decimal places. If you compare the two values, you will see that the float value is acurate to 16 digits:

26087635650665566208.0000000000...
26087635650665564424.6991436125...

You can use the Decimal type instead, which allows you to choose how many decimal digits to store. By default it will give you the answer you want to eight decimal places (or 28 significant figures):

from decimal import Decimal
2**Decimal('64.5')
# returns Decimal('26087635650665564424.69914361')

But Decimal is slower than float, and it has its own traps and gotchas too.

At the cost of more memory and more time, we can calculate the answer to 100 significant figures:

import decimal
with decimal.localcontext() as ctx:
    ctx.prec = 100
    2**Decimal('64.5')

which will return:

Decimal('26087635650665564424.69914361250501673776655257918571715746415135307033663760525399712144022679277348')
4 Likes

Further to the discussion, we can see what resolution is available around the value we want.

We want the closest possible value to the mathematically exact √(2**129), or 2**64.5. To four decimal places, that is

26087635650665564424.6991

but if we try entering that as a 64-bit float, the closest number available is just

26087635650665566208

Around that value, each float is 4096 units apart from the next closest float!

x = 2**64.5
print(x - math.nextafter(x, 0))
print(math.nextafter(x, 2*x) - x)

will print 4096.0 twice. This tells us that the difference between two floats close to our target value 2**64.5 is 4096. Numbers less than that apart are impossible to represent as floats.

The exact value we want is 26087635650665564424.6991… but the three closest values that are available with 64 bits in that area are

exists-  26087635650665562112
wanted-  26087635650665564424.6991...
exists-  26087635650665566208  # what we actually get
exists-  26087635650665570304    

As you can see, the value that we actually get is the closest possible float.

How wrong is the answer that we get? It is off by about 1784, which seems like a lot, but it’s only 0.000000000000007% wrong! So it’s actually pretty close.

4 Likes

Thanks for the explanation! All I really care about is that the integer part is correct and one decimal place is shown. Speed is not an issue, so the Decimal format will work.

Thanks again!
Bassam