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.

## No comments:

Post a Comment