On Wed, Feb 12, 2014 at 7:20 PM, Zimmermann Paul <paul.zimmerm...@inria.fr> wrote: > William, > > thank you for putting me in cc. > >> From: William Stein <wst...@gmail.com> >> Date: Wed, 12 Feb 2014 06:01:29 -0800 >> >> On Wed, Feb 12, 2014 at 4:55 AM, <ref...@uncg.edu> wrote: >> > Ah, I see what you mean. If that's the case then I understand. How does >> > one find out if this is true? >> >> 1. It is *NOT* true. Sage just directly calls the PARI C library, >> which is always loaded on startup of Sage. >> >> 2. Thanks for clarifying your original question. It's surprising that >> MPFR is a full *order of magnitude* slower than PARI at computing >> gamma on real input. It's pretty likely that when the line of code >> in question was written (prob in 2005 or 2006) this speed difference >> wasn't the case. I've cc'd Paul Zimmerman (author of MPFR) in case he >> has any comments about this: >> >> a = RealField(1000)(pi^2) >> b = ComplexField(1000)(pi^2) >> >> %timeit a.gamma() # SLOW -- uses mpfr -- "mpfr_gamma(x.value, >> self.value, (<RealField_class>self._parent).rnd)" >> %timeit b.gamma() # FAST -- uses pari -- >> "sage.libs.pari.all.pari.complex(self.real()._pari_(), >> self.imag()._pari_())" >> >> Maybe we should switch to PARI for this. At the least, you could add >> >> a.gamma(algorithm='pari') >> >> and at some point make that the default... >> >> -- William > > indeed mpfr_gamma is slow for input pi^2 at 1000 bits. I've investigated a bit > and with the following patch (with respect to the svn version of MPFR) I can > gain about 60%, i.e., the time decreases from about 2.9s on my computer to > about 1.8s for 1000 computations of gamma(pi^2) at 1000 bits. This is still > slower than Pari (which takes 506 us per loop). > > We could still save some time by caching the Bernoulli numbers (which take > 0.8s in the above 1.8s), but that would still be slower than Pari. Maybe the > algorithm used by Pari is better...
Tuning of the Stirling series is discussed briefly in my paper http://arxiv.org/abs/1310.3741 If you cache Bernoulli numbers, it becomes much faster to use a smaller argument reduction z -> z+k. My code uses k ~= 0.27*p where p is the precision in bits. Fredrik -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscr...@googlegroups.com. To post to this group, send email to sage-devel@googlegroups.com. Visit this group at http://groups.google.com/group/sage-devel. For more options, visit https://groups.google.com/groups/opt_out.
