0

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()
Community
  • 1
  • 1
ocopa
  • 75
  • 8
  • First, I don't know how you got `random.randint` to work with these limits. It doesn't on my machine... But anyway I suspect this is related to float arithmetic accuracy. – Aguy Jul 27 '16 at 10:40
  • I think you are right, in which case my question is, how do i get round that, either in python or gmpy2 or some other library? – ocopa Jul 27 '16 at 13:00
  • I can get the anti-log by method described at https://en.wikipedia.org/wiki/Exponentiation_by_squaring, so for log 2 it is a matter of not going via natural log, but directly computing log base 2. – ocopa Jul 27 '16 at 13:23

3 Answers3

1

As long as you set the decimal module accuracy, the usual suggestion is to use Decimal datatype

from decimal import Decimal, getcontext

getcontext().prec = 1000

# Just a different method to get the random number:
x = Decimal(round(10**20 * (1 + 9 * random.random()))) 

x10 = Decimal.log10(x)
e10 = 10**x10

e10 - x
#outputs: Decimal('5.2E-978')

For different bases you may want to use the logarithmic formula:

x2 = Decimal.log10(x) / Decimal.log10(Decimal('2'))
e2 = 2**x2

e2 - x
#outputs: Decimal('3.9E-978')
Aguy
  • 7,851
  • 5
  • 31
  • 58
  • OK, thanks, this works for Decimal, log10, but what works in general for "various bases", ie different bases, any integer base other than natural log, eg 2, to recover x? – ocopa Jul 27 '16 at 22:54
  • @ocopa - see my edit regarding different bases. Just use log10(x)/log10(base). – Aguy Jul 28 '16 at 03:22
  • Yes, i knew that but i was trying to work out why i was getting the wrong answer from gmpy2, and what was the underlying code. It seems gmpy2 is using log base e behind the scene or maybe there is a float accuracy issue as we thought. i will mark this as answered after i post a couple of useful references, shortly. – ocopa Jul 28 '16 at 05:49
  • Useful Logarithm references https://en.wikipedia.org/wiki/Binary_logarithm answer by Dan04 at http://stackoverflow.com/questions/2882706/how-can-i-write-a-power-function-myself answer by log0 at http://stackoverflow.com/questions/3719631/log-to-the-base-2-in-python?rq=1 – ocopa Jul 28 '16 at 06:09
1

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
casevh
  • 11,093
  • 1
  • 24
  • 35
0

Aguy solved my problem, as acknowledged. I had not taken account of needing more than 15 digits of precision. This answer to another question covers that ground. gmpy2 log2 not accurate after 16 digits

Community
  • 1
  • 1
ocopa
  • 75
  • 8