Benchmark using floating point math

Question

I wrote some code in Python, to benchmark mathematical operations. It does not do anything in particular besides just calculating a lot of floating point numbers.

To clarify: I need 128-bits precision floating point numbers.

from __future__ import division
from numpy import arctan, sin, absolute, log10
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 = mpf(x + 1)
 y = mpf(y + x * x)
 z = mpf(z + sin(y))
 t = mpf(t + absolute(z))
 u = mpf(u + log10(t))
print("--- %s seconds ---" % (time.time() - start_time))
print x
print y
print z
print t
print u

This took approximately 87 secs to print results. How can I improve this code? My friend in fortran wrote similar program and it took only 3.14 secs for him to print results in 128-bit precision.

Also in general I want to improve speed of my code. Instead 87 secs I want to reach times much closer to fortran.

Answer 1

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