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...
Best regards,
Paul
tarte% svn diff src/lngamma.c src/gamma.c
Index: src/lngamma.c
===================================================================
--- src/lngamma.c (revision 8920)
+++ src/lngamma.c (working copy)
@@ -437,10 +437,13 @@
with respect to z0, and k = n, we get 1/(Pi*e)^(2n) ~ 2^(-w), i.e.,
k ~ w*log(2)/2/log(Pi*e) ~ 0.1616 * w.
However, since the series is more expensive to compute, the optimal
- value seems to be k ~ 4.5 * w experimentally. */
+ value seems to be k ~ 4.5 * w experimentally.
+ Note added February 12, 2014: for a target precision of 1000 bits,
+ gamma(pi^2) with k = 4.5*w gives m=55 and k=4639 which is about 60%
+ slower than k=1.5*w (m=69 and k=1540), thus we change for 1.5*w. */
mpfr_set_prec (s, 53);
mpfr_gamma_alpha (s, w);
- mpfr_set_ui_2exp (s, 9, -1, MPFR_RNDU);
+ mpfr_set_ui_2exp (s, 3, -1, MPFR_RNDU);
mpfr_mul_ui (s, s, w, MPFR_RNDU);
if (mpfr_cmp (z0, s) < 0)
{
Index: src/gamma.c
===================================================================
--- src/gamma.c (revision 8920)
+++ src/gamma.c (working copy)
@@ -258,6 +258,10 @@
mpfr_exp_t expxp;
MPFR_BLOCK_DECL (flags);
+ /* quick test for the default exponent range */
+ if (mpfr_get_emax () >= 1073741823UL && MPFR_GET_EXP(x) <= 25)
+ return mpfr_gamma_aux (gamma, x, rnd_mode);
+
/* 1/e rounded down to 53 bits */
#define EXPM1_STR "0.010111100010110101011000110110001011001110111100111"
mpfr_init2 (xp, 53);
--
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.
