Disclaimer: I maintain the gmpy2 library which also provides arbitrary-precision arithmetic.
gmpy2
includes fast integer arithmetic based on the GMP library. If gmpy2
is not installed, mpmath
uses Python's native integer type(s) for computation. If gmpy2
is installed, mpmath
uses the faster gmpy2.mpz
integer type. You didn't mention if you have gmpy2
(or gmpy
, an older version) installed, so I ran one test with your original code, mpmath
, and without gmpy2
. The running time was ~91 seconds. With gmpy2
installed, the running time was ~63 seconds.
The rest of the examples assume gmpy2
is installed.
Your code include superfluous calls to mpf()
. Since x
is already an mpf
, the result of x+1
will also be an mpf
so the mpf()
call is not needed. If I remove those calls, the running time drops to ~56 seconds.
You import from both numpy
and mpmath
. mpmath.sin
and mpmath.log10
replace the functions imported from numpy
. But you are still using numpy.absolute
. Numpy can be very fast but only when uses types it natively supports. If I remove numpy
and change absolute
to abs
, the running time drops to ~48 seconds.
Here is the code with all the changes listed from above:
from __future__ import division
from mpmath import *
import time
imax = 1000001
x = mpf(0)
y = mpf(0)
z = mpf(0)
t = mpf(0)
u = mpf(0)
i = 1
mp.prec = 128
start_time = time.time()
while i < 1000001:
i += 1
x = x + 1
y = y + x * x
z = z + sin(y)
t = t + abs(z)
u = u + log10(t)
print("--- %s seconds ---" % (time.time() - start_time))
print x
print y
print z
print t
print u
I don't see any other significant improvements for an mpmath
based approach.
gmpy2
also provides arbitrary-precision floating point based on the MPFR library. If I use gmpy2
instead, the running time is reduced to ~18 seconds.
Here is the code for a gmpy2
based solution:
from __future__ import division
import gmpy2
from gmpy2 import mpfr, sin, log10
import time
imax = 1000001
x = mpfr(0)
y = mpfr(0)
z = mpfr(0)
t = mpfr(0)
u = mpfr(0)
i = 1
gmpy2.get_context().precision = 128
start_time = time.time()
while i < 1000001:
i += 1
x = x + 1
y = y + x * x
z = z + sin(y)
t = t + abs(z)
u = u + log10(t)
print("--- %s seconds ---" % (time.time() - start_time))
print x
print y
print z
print t
print u