On Tue, Apr 14, 2009 at 7:41 AM, William Stein <wst...@gmail.com> wrote:
> On Mon, Apr 13, 2009 at 10:13 PM, Nick Alexander <ncalexan...@gmail.com>
> wrote:
>>
>>
>> On 13-Apr-09, at 8:08 PM, Ondrej Certik wrote:
>>
>>>
>>>> Actually, it's not using fast_callable yet, but I do plan on changing
>>>> that (which should make things much faster). As for the zeta
>>>> function--is mpmath faster than pari? (If so, should we be using it
>>>> for zeta in general?)
>>>
>>> Probably not (yet). I just didn't know that Pari can do such nice
>>> graphs.
>>
>> I need to compute gamma functions and incomplete gamma functions.
>> mpmath is so slow that my doctests take hundreds of seconds (and I
>> always killed them) instead of 10-15 seconds using pari to compute
>> such things. The mpmath main loops are pure python, there's no way
>> it's going to be faster than even the most naive pari implementation.
Mpmath is slower than it needs to be for certain tasks (such as plotting)
because it doesn't use double precision numbers at all. The original
motivation for this approach was that the equivalent functionality for
double precision was already covered adequately by SciPy, but I've
gradually come to realize that this is a mistake because SciPy is
lacking a lot of functions (and many others are broken). There will
be better support for low precision calculations in future versions
of mpmath and where this can be used, it will likely provide an
order-of-magnitude speedup.
Also, the gamma and zeta functions in mpmath, for mpmath numbers,
are a little slower than they need to be because of naive algorithms
(this will improve to some extent in the future).
If you need repeated gamma function evaluations at high precision,
for small |Re(x)| (~ 0-100) and very small |Im(x)| (~= 0-3)
you may be interested in knowing that I'm working on a faster
algorithm in mpmath for this particular purpose. From a warm
cache, it currently computes gamma(3.7) slightly faster than PARI
at 50 digit precision and 3x faster than PARI at 500 digits (on my
system).
> Here is a quantitative comparison if the time of mpmath versus pari
> for computation of the Riemann zeta function. This uses the 100%
> pure Python version of mpmath, since Sage doesn't include gmpy by
> default, however I did try installing gmpy and that made no difference
> in timings.
It mostly makes a difference at huge precision. At low precision
there can be perhaps a 2x difference in some situations.
> SUMMARY: If you use Sage correctly, and are mainly interested in
> *double precision*, then Sage is between 30 and over 2400 (!) times
> faster than mpmath at evaluating the Riemann zeta function.
This is very much true. If you are interested in double precision, then
definitely, use double precision :-)
> For double precision reals, Sage is dramatically faster than
> PARI/mpmath "because" it uses GSL, which evidently has a very fast
> zeta implementation. For complex numbers it is about 30 times faster,
> because it uses PARI.
30 times sounds plausible. Part of the reason is that mpmath doesn't
use the Euler-Maclaurin formula (I presume PARI does), thus requiring
more exponentials, and part is the speed of the exponential function.
The former can obviously be fixed, and the latter could be fixed with
Mario Pernici's Taylor series code for gmpy. I don't think mpmath can
be quite as fast as PARI even with these changes, but I think
the factor could be 3 rather than 30.
FYI, I wrote a fairly straightforward pure-Python implementation of the
Hurwitz/Riemann zeta function for float/double and it appears to be
faster than PARI:
sage: a = RDF(1.23)
sage: timeit('zeta(a)')
625 loops, best of 3: 1.05 µs per loop
sage: zeta(a)
4.94152918822
sage: a = CDF(2.1,3.1)
sage: timeit('zeta(a)')
625 loops, best of 3: 261 µs per loop
sage: zeta(a)
0.809508074908 - 0.0982533523609*I
>>> ( paste source for hzeta() )
>>> from mpmath import timing
>>> timing(hzeta, 1.23) * 1e6 # * 1e6 for µs
17.935240373390116
>>> hzeta(1.23)
4.941529188217384
>>> timing(hzeta, 2.1+3.1j) * 1e6
47.687625102810216
>>> hzeta(2.1+3.1j)
(0.80950807490835608-0.09825335236086763j)
That is still 18 times slower than sage (6 times after I turned on psyco) in the
real case, but 261 / 48 = 5.5 times faster than sage in the complex case.
> I could not figure out how to make arbitrary precision numbers in mpmath:
Set the working precision with mp.dps = 100 (decimals) or mp.prec = 333 (bits).
Fredrik
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to
sage-devel-unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
