When using log2() in gmpy2 it does not seem to be accurate after 16 digits. It seems to work fine at 15 digits but after that the answer is not correct using mpz(mpfr(2) ** mpfr(x)). Do I need change the precision? I thought python by itself would be accurate up to 53 digits.
Additionally, is there a way in gmpy2 to use a logarithm operation in bases besides 10 and 2? For example, base 8 or 16.
The standard Python float type is accurate to 53 bits which is roughly 16 decimal digits. gmpy2 uses a default precision of 53 bits. If you want more accurate results, you will need to increase the precision.
>>> import gmpy2
>>> from gmpy2 import mpz,mpfr,log2
>>> a=12345678901234567890
>>> gmpy2.get_context().precision=70
>>> mpz(2**log2(a))
mpz(12345678901234567890L)
To calculate a logarithm in a different, just use
>>> gmpy2.log(x)/gmpy2.log(base)
Update
Recovering an exact integer result from a sequence of floating point calculations is generally not possible. Depending on the actual calculations, you can increase the precision until you get "close enough".
Let's look at the impact of precision. Note that a is 57 bits long so it cannot be exactly represented with 53 bits of floating point precision.
>>> a=123543221556677776
>>> a.bit_length()
57
>>> gmpy2.get_context().precision=53
>>> mpfr(a);2**log2(a)
mpfr('1.2354322155667778e+17')
mpfr('1.2354322155667752e+17')
Since conversion of a binary floating point number to decimal can introduce a conversion error, lets look at the results in binary.
>>> mpfr(a).digits(2);(2**log2(a)).digits(2)
('11011011011101001111001111100101101001011000011001001', 57, 53)
('11011011011101001111001111100101101001011000010111001', 57, 53)
Let's trying increasing the precision to 57 bits.
>>> gmpy2.get_context().precision=57
>>> mpfr(a).digits(2);(2**log2(a)).digits(2)
('110110110111010011110011111001011010010110000110010010000', 57, 57)
('110110110111010011110011111001011010010110000110010011000', 57, 57)
Notice more bits are correct but there is still an error. Let's try 64 bits.
>>> gmpy2.get_context().precision=64
>>> mpfr(a);2**log2(a)
mpfr('123543221556677776.0',64)
mpfr('123543221556677775.953',64)
>>> mpfr(a).digits(2);(2**log2(a)).digits(2)
('1101101101110100111100111110010110100101100001100100100000000000', 57, 64)
('1101101101110100111100111110010110100101100001100100011111111010', 57, 64)
The large number of trailing 1's is roughly equivalent to trailing 9's in decimal.
Once you get "close enough", you can convert to an integer which will round the result to the expected value.
Why isn't 57 bits sufficient? The MPFR library that is used by gmpy2 does perform correct rounding. There is still a small error. Let's also look at the results using the floating point values immediately above and below the correctly rounded value.
>>> gmpy2.get_context().precision=57
>>> b=log2(a)
>>> 2**gmpy2.next_below(b);2**log2(a);2**gmpy2.next_above(b)
mpfr('123543221556677746.0',57)
mpfr('123543221556677784.0',57)
mpfr('123543221556677822.0',57)
Notice that even a small change in b causes a much larger change in 2**b.
Update 2
Floating point arithmetic is only an approximation to the mathematical properties of real numbers. Some numbers are rational (they can be written as a fraction) but most numbers are irrational (they can never be written exactly as a fraction). Floating point arithmetic actually uses a rational approximation to a number.
I've skipped some of the details in the following - I assume all numbers are between 0 and 1.
With binary floating point (what most computers use), the denominator of the rational approximation must be a power of 2. Numbers like 1/2 or 1/4 can be represented exactly. Decimal floating point uses rational approximations that have a denominator that is a power of 10. Numbers like 1/2, '1/4', '1/5', and 1/20 can all be represented exactly. Neither can represent 1/3 exactly. A base-6 implementation of floating point arithmetic can represent 1/2 and 1/3 exactly but not 1/10. The precision of a particular format just specifies the maximum size of the numerator. There will always be some rational numbers that cannot be represented exactly by a given base.
Since irrational numbers can't be written as a rational number, they can not be represented exactly by a given base. Since logarithm and exponential functions almost always result in irrational values, the calculations are almost never exact. By increasing the precision, you can usually get "close enough" but you can never get exact.
There are programs that work symbolically - they remember that a is log2(n) and when you do 2**a, the exact value of a is returned. See SymPy.
