Using Sage numbers in mpmath

January 27, 2010

I’ve written a very basic mpmath context for computing “natively” with Sage’s real and complex numbers. It can be tried out by upgrading mpmath in Sage to the svn trunk, applying this patch this patch to fix the helper code in Sage, and then importing this file.

Some examples:

------------------------------------------------------------------
| Sage Version 4.3.1, Release Date: 2010-01-20                   |
| Type notebook() for the GUI, and license() for information.    |
------------------------------------------------------------------
sage: from mpsage import SageContext
sage: from mpmath import mp, timing
sage: mp.prec = 150
sage: sc = SageContext(150)
sage: 
sage: print sc.quad(lambda t: sc.exp(-t**2), [0,sc.inf])
0.88622692545275801364908374167057259139877471
sage: print mp.quad(lambda t: mp.exp(-t**2), [0,mp.inf])
0.88622692545275801364908374167057259139877473
sage: timing(sc.quad, lambda t: sc.exp(-t**2), [0,sc.inf])
0.40903496742248535
sage: timing(mp.quad, lambda t: mp.exp(-t**2), [0,mp.inf])
0.43610286712646484
sage: 
sage: print sc.zeta(4+5*I)
0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897*I
sage: print mp.zeta(4+5*I, method='euler-maclaurin')
(0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897j)
sage: timing(sc.zeta, 4+5*I)
0.015111517906188966
sage: timing(mp.zeta, 4+5*I, method='euler-maclaurin')
0.028736209869384764
sage: 
sage: print sc.zeta(1+100000000*I)
1.8349308953278466065175598876062914152382527 + 1.0691847904236128969174476736978194200591565*I
sage: print mp.zeta(1+100000000*I)
(1.8349308953278466065175598876062914152563156 + 1.069184790423612896917447673697819420203942j)
sage: timing(sc.zeta, 1+100000000*I)
1.0952680110931396
sage: timing(mp.zeta, 1+100000000*I)
2.0537140369415283

Unfortunately, there are some issues (besides the fact that some methods are missing, so not everything works yet).

This context doesn’t provide variable precision, so the user has to manually allocate sufficient extra precision to compensate for rounding errors.

I first tried to do it, but variable precision is very inconvenient to implement using Sage’s way of managing precision. There is no direct way to perform operations with a given target precision (independent of the inputs), and the second best option (which is to perform coercions to the target precision everywhere) is very slow, besides losing accuracy. The only way to implement variable precision in a reasonable way is to bypass Sage’s RealNumber and ComplexNumber (at least their public interfaces) and wrap MPFR directly using a precision model similar to what MPFR and mpmath uses, where the precision of the result is always specified and independent of the inputs.

Secondly, there is not that much of a speedup (in the examples above, the speedup is at most about 2x). This is mainly due to the fact that the context uses wrapper classes for RealNumber and ComplexNumber, with all interface and conversion code written in Python. So the overhead is about the same as for the corresponding code in vanilla mpmath (where it accounts for about 50% of the total overhead). The reason RealNumber and ComplexNumber can’t be used directly, even in a fixed-precision setting, is that mpmath in many places multiplies numbers by floats (usually exact float literals like 0.5), and Sage always coerces to the number with less precision. This could presumably be fixed by replacing all float and complex constants in mpmath, but I’m not in the mood to do that right now.

There should be a significant performance benefit if direct-from-Cython classes and conversion methods were used. It’s also very important to optimize certain core functions; for example, quadrature would benefit greatly from a fast fdot implementation.

All things considered, I’m probably not going to continue much more down this particular path. It’s better to write a fully compatible Cython context with classes designed directly for mpmath. Trying to wrap Sage’s numbers did however help identify a few problems with the existing interfaces, so I might extend the work on this context a little more just to find more such issues.

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor