Search code examples
pythonlogarithmgmpy

Python - How to avoid discrepancy of base y**log base y of x, in gmpy2


The Following code examples my problem which does not occur between 10 power 10 and 10 power 11, but does for the example given in code and above it.

I can't see where in my code I am not properly handling the retrieval of the original value. May be I have just missed something simple.

I need to be sure that I can recover x from log x for various bases. Rather than rely on a library function such as gmpy2, is there any reverse anti-log algorithm which guarantees that for say 2**log2(x) it will give x.

I can see how to directly develop a log, but not how to get back, eg, Taylor series needs a lot of terms... How can I write a power function myself? and @dan04 reply. Code follows.

from gmpy2 import gcd, floor, next_prime, is_prime    
from gmpy2 import factorial, sqrt, exp, log,log2,log10,exp2,exp10    
from gmpy2 import mpz, mpq, mpfr, mpc, f_mod, c_mod,lgamma    
from time import clock    
import random    
from decimal import getcontext
x=getcontext().prec=1000 #also tried 56, 28
print(getcontext())

def rint():#check accuracy of exp(log(x))
    e=exp(1)
    l2=log(2)
    l10=log(10)
    #x=random.randint(10**20,10**21) --replaced with an actual value on next line
    x=481945878080003762113
    # logs to different bases
    x2=log2(x)
    x10=log10(x)
    xe=log(x)
    # logs back to base e
    x2e=xe/l2
    x10e=xe/l10
    #
    e2=round(2**x2)
    e10=round(10**x10)
    ex=round(e**xe)
    #
    ex2e=round(2**x2e)
    ex10e=round(10**x10e)
    error=5*x-(e2+e10+ex+ex2e+ex10e)
    print(x,"error sum",error)
    #print(x,x2,x10,xe)
    #print(x2e,x10e)
    print(e2,e10,ex)
    print(ex2e,ex10e)
 rint()

Solution

  • Note: I maintain the gmpy2 library.

    In your example, you are using getcontext() from the decimal module. You are not changing the precision used by gmpy2. Since the default precision of gmpy2 is 53 bits and your value of x requires 69 bits, it is expected that you have an error.

    Here is a corrected version of your example that illustrates how the accumulated error changes as you increase the precision.

    import gmpy2
    
    def rint(n):
        gmpy2.get_context().precision = n
        # check accuracy of exp(log(x))
        e = gmpy2.exp(1)
        l2 = gmpy2.log(2)
        l10 = gmpy2.log(10)
        x = 481945878080003762113
        # logs to different bases
        x2 = gmpy2.log2(x)
        x10 = gmpy2.log10(x)
        xe = gmpy2.log(x)
        # logs back to base e
        x2e = xe/l2
        x10e = xe/l10
        #
        e2 = round(2**x2)
        e10 = round(10**x10)
        ex = round(e**xe)
        #
        ex2e = round(2**x2e)
        ex10e = round(10**x10e)
        error = 5 * x - (e2 + e10 + ex + ex2e + ex10e)
        print("precision", n, "value", x, "error sum", error)
    
    for n in range(65, 81):
        rint(n)
    

    And here are the results.

    precision 65 value 481945878080003762113 error sum 1061
    precision 66 value 481945878080003762113 error sum 525
    precision 67 value 481945878080003762113 error sum -219
    precision 68 value 481945878080003762113 error sum 181
    precision 69 value 481945878080003762113 error sum -79
    precision 70 value 481945878080003762113 error sum 50
    precision 71 value 481945878080003762113 error sum -15
    precision 72 value 481945878080003762113 error sum -14
    precision 73 value 481945878080003762113 error sum 0
    precision 74 value 481945878080003762113 error sum -2
    precision 75 value 481945878080003762113 error sum 1
    precision 76 value 481945878080003762113 error sum 0
    precision 77 value 481945878080003762113 error sum 0
    precision 78 value 481945878080003762113 error sum 0
    precision 79 value 481945878080003762113 error sum 0
    precision 80 value 481945878080003762113 error sum 0