So far I've implemented a real type (mpf replacement) and a complex type (mpc replacement); addition, multiplication, and the real exponential function; just enough to do some basic benchmarking. Below is a comparison between standard Python mpmath (running on top of sage.Integer), the new Cython-based types, and Sage's RealNumber / ComplexNumber types. I used the inputs x = sqrt(3)-1, y = sqrt(5)-1, and for the complex cases z = x+yi, w = y+xi. The precision ranges between 53 bits (IEEE double compatible) and 3333 bits (≈1000 decimal digits).

Addition, x+y

mpmath cython sage

53 8.72 µs 979 ns 707 ns

100 8.81 µs 1.01 µs 733 ns

333 10.1 µs 1.14 µs 737 ns

3333 10.6 µs 1.59 µs 1.15 µs

Multiplication, x*y

53 9.04 µs 968 ns 728 ns

100 9.56 µs 1.14 µs 901 ns

333 11.3 µs 1.59 µs 1.2 µs

3333 35.7 µs 25.4 µs 17.8 µs

Real exponential function, exp(x)

53 56.2 µs 9.41 µs 16 µs

100 52.8 µs 11.6 µs 11.3 µs

333 122 µs 25.4 µs 26.9 µs

3333 1.38 ms 831 µs 1.14 ms

Complex addition, z+w

53 14.7 µs 1.45 µs 978 ns

100 14.5 µs 1.7 µs 1 µs

333 16.3 µs 1.85 µs 992 ns

3333 17.7 µs 3.2 µs 1.73 µs

Complex multiplication, z*w

53 35.8 µs 2.72 µs 1.68 µs

100 34.5 µs 2.76 µs 2.08 µs

333 43.2 µs 4.38 µs 3.61 µs

3333 139 µs 98.6 µs 70 µs

A few observations can be made. First off, basic arithmetic with number instances is an order of magnitude faster when fully Cythonized. This isn't surprising because something as simple as an addition of two mpmath.mpf instances has to go through some 50 lines of Python code to check types, check for special cases, and round the result (creating lots of temporary objects, etc).

Reimplemented in Cython, arithmetic comes well within half the speed of Sage's numbers. I believe MPFR (which Sage uses) does arithmetic faster because it works directly with the limb data and uses some tricks to speed up the rounding, whereas my implementation uses the mpz interface straightforwardly. Some difference also comes from the fact that MPFR uses machine-precision integers for exponents, whereas I'm using mpz_t to allow arbitrary-size exponents. In any case, the results are very good.

At very low precision, memory allocations account for probably half the time (for both my Cython-based types and Sage's RealNumber). Thus there is some room for improvement later on by implementing a freelist.

The best news (for special functions) is that my exponential function is as fast as MPFR's, so the same should be true for other functions when I get there. The new exp (which uses a slightly different algorithm from the one currently in mpmath) is actually slightly faster than MPFR, although it should be said that performance for transcendental functions depends heavily on tuning, the inputs, and possibly other factors (I have no idea why Sage's RealNumber.exp in my benchmark is much slower at 53 bits than at 100 bits, for example), so this is not a conclusive result.

The code itself is currently too messy and ad-hoc to share publicly, sorry.

## No comments:

Post a Comment