Cython mpmath performance

I'm presently working on a Cython-based backend for mpmath, with the immediate goal to make mpmath competitive as a component of Sage. The Cython code currently depends on utility functions in the Sage library but should eventually be possible to compile for standalone use or together with SymPy's future Cython backend (linking directly against GMP/MPIR).

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.

Report from Sage Days 15

Sage Days 15 is now over. I'm still in Seattle as I write this; I depart early tomorrow. It's been a great week in almost every way. I didn't take photos, but see William Stein's photos and video recordings.

Edit: David Joyner also has photos:
http://sage.math.washington.edu/home/wdj/sagedays/sd15/

People

My prejudices were confirmed insofar as that the Sage developers are an exceptionally friendly bunch. Everybody was helpful and had a good sense of humor. And of course, these people are all very good at what they work on. It's both humbling and inspiring to be surrounded by people who are far more knowledgeable than yourself about almost everything.

There were a couple of talks daily (except for the last two days), most of which were very good; not to mention all the interesting informal discussions. Personally, I thought it was particularly interesting to get Bill Hart's (MPIR and FLINT guy) perspective of his work; he also provided useful technical advice. William's presentation of the general direction of the Sage project and the Cython talk by Robert Bradshaw and Craig Citro were also quite informative.

I gave a short talk about mpmath on Monday (slides). I'm not sure if the talk itself was good, but several people seemed to be interested in my work. I've received questions and comments about mpmath all week, and the general attitute towards using it in Sage seems positive.

Cython

On the purely technical side, by far the coolest thing about SD15 was the chance to finally learn Cython. I've seen Cython in action before, but never actually used it myself. My impressions are exclusively positive. I have plenty of experience with Python and a decent amount of experience with C, and Cython truly combines the best (alternatively, removes the worst) of both worlds. The biggest surprise (and of course this is Cython's selling point) is how simple the interfacing between high level and low level code becomes, and the fact that it is all very robust.

It's exiciting to see that there are several active projects around that attempt to speed up Python. The nice thing about Cython is that it doesn't give you "half the speed of C" or "maybe nearly the speed of C, 3 years from now" -- it gives the real deal, -O3 C, and it works right now. I hope Cython will be part of the standard Python distribution within a year or two from now, and thereby make it even easier to build Cython extensions, especially on Windows.

On a tangential note, it's good to see that progress is being made on the native Windows port of Sage. (I dual boot between Windows and Ubuntu, and vmware Sage on Windows is plainly inconvenient.) The pragmatic attitude of various Sage developers towards issues such as support for proprietary platforms gives the impression of a project in good hands.

Coding

I offered right before SD15 to look at speeding up the prime counting function in Sage by implementing the asymptotically fast Meissel-Lehmer-Lagarias-Miller-Odlyzko method. Doing this from scratch turned out to be unnecessary when Victor Miller emailed me the C implementation of the algorithm which he wrote two decades ago! The code compiles, and is amazingly fast on my laptop; unfortunately, it has some issues (such as not working for small n). It would be possible to wrap it in Sage as-is, and just use a different algorithm for small n. But what's probably an even better solution is that Victor himself (he showed up at SD) is now working on reimplementing his algorithm in Sage.

Most of my productive coding time (I also spent quite a lot of time just poking around at the Sage codebase -- which I've never really looked at before) went into wrapping mpmath for Sage. The purpose of wrapping mpmath is mainly to improve the support for special functions in Sage, firstly by just making the mpmath function library available and secondly by further improving the algorithms in mpmath. (I am getting funded via William's research grant to work on the latter this summer.)

As a first step, I made mpmath able to use sage.Integer instead of Python long. This was mostly trivial, the only problem being a coercion-related bug in sage.Integer that had to be fixed (William helped me write the patch, my first for Sage!).

When all mpmath tests finally passed, they unfortunately turned out to run 40% times slower than pure Python and 2.6x (!) slower than gmpy mpmath. The slowdown was in fact expected since sage.Integer was known to be substantially slower than gmpy.mpz. I therefore helped Robert Bradshaw identify the performance bottlenecks (they were mostly due either to unnecessary coercions or unnecessary memory reallocation), and he promptly patched them up. A nice side effect is that many parts of Sage should be faster now, since much code depends on the Integer type.

After I also implemented a few exra helper functions in Cython, mpmath on top of Sage now runs just a few percent slower than gmpy mpmath. Much of the remaining difference is probably due to inefficiencies in my own Cython code.

Finally, I also wrote wrapper code in Cython that permits Sage to call mpmath functions with Sage objects as input and getting Sage RealNumber or ComplexNumber back, all with reasonably low conversion overhead. This should now make it trivial to fix the state of numerical special functions in Sage (wrapping an mpmath function is literally a one-liner). I'm just waiting for Sage 4.0 with the new symbolics code to come out before I submit a patch and start adding wrapped functions; I want to ensure that the symbolics and numerics integrate well.

The state of mpmath is that it is still a bit slow compared to some numerical functions in Sage, but it's not at all bad. The long term plan is to speed up mpmath further by providing a pure-Cython backend, which along with certain algorithmic improvements should make it roughly as fast as MPFR and Pari for cases where it's currently slower.

Travel

The flight over here was reasonably convenient considering the circa 20-hour duration with transfers accounted for (the repeated security checks were really the biggest hassle). Hopefully the flight back tomorrow will proceed as smoothly.

Seattle, at least the University District, is a fairly nice city (not least because of the trees everywhere). Except for the trip to Microsoft Research on Tuesday, all activities took place
at the University of Washington campus. The campus alone is huge, so there has been plenty to see. I've lived very conveniently at the College Inn, right next to both the campus and The Ave.

I did make one excursion to downtown Seattle. I arrived early last Friday and with nothing on the schedule that day, I had time to meet up with Huy Pham (fellow Doomer and online acquintance) and we went to see Star Trek in IMAX (at the Pacific Science Center), which was well worth the admission.

As a last remark, I'm very grateful to William and the other Sage Days 15 organizers for their good job and for inviting me (and paying for my trip!). I've had a great time, had the opportunity meet interesting people, learned many things, and now I'm sad it's over.