Re: [sage-devel] Re: Issue With Implementation of Gamma Function

[Zimmermann Paul]

   Hi Jori,

> Date: Mon, 17 Feb 2014 14:56:56 +0200 (EET)
> From: Jori Mantysalo 
>=20
> On Mon, 17 Feb 2014, Zimmermann Paul wrote:
>=20
> > On my computer, computing gamma(Pi^2) to 1 bits takes about 1.4s (f=
or the
> > first computation, when Bernoulli numbers are not cached) instead of 6.=
1s with
> > Pari/GP.
>=20
> Out of curiosity: How this compares to, for example, Mathematica?
>=20
> --=20
> Jori M=C3=A4ntysalo

I don't have access to Mathematica, but according to [1], Mathematica 8.0 t=
akes
about 1.7 seconds for the first computation at 3000 decimal digits, which i=
s
comparable to 1 bits (however this is for complex input, I don't know i=
f
Mathematica is faster for real numbers).

Subsequent computations on my machine take 0.083s each with the development
version of MPFR (at 1 bits).

Paul

[1] http://fredrikj.net/blog/2013/02/timings-for-the-complex-gamma-function=
-in-arb/

--=20
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 e=
mail 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.

Re: [sage-devel] Re: Issue With Implementation of Gamma Function

[Zimmermann Paul]

a last update: MPFR now uses the Von Staudt=E2=80=93Clausen theorem to comp=
ute
Bernoulli numbers (thanks to Fredrik Johansson for pointing out this formul=
a).
On my computer, computing gamma(Pi^2) to 1 bits takes about 1.4s (for t=
he
first computation, when Bernoulli numbers are not cached) instead of 6.1s w=
ith
Pari/GP.

Paul Zimmermann

--=20
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 e=
mail 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.

Re: [sage-devel] Re: Issue With Implementation of Gamma Function

[Zimmermann Paul]

an update on this: we now cache Bernoulli numbers in MPFR, and the time to
compute 1000 times gamma(pi^2) at precision 1000 bits has decreased to 0.4s
on my computer (was 2.9s before). In comparison, Pari/GP (which also does
cache Bernoulli numbers) takes 0.5s.

Best regards,
Paul Zimmermann

-- 
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.

Re: [sage-devel] Re: Issue With Implementation of Gamma Function

[Zimmermann Paul]

   Dear Jori,

> And reason is of course clear, as Fredrik Johansson wrote "If you cache 
> Bernoulli numbers, - -".

in fact there is another reason: the MPFR code computes the Bernoulli numbers
exactly, as integers B(2n)*(2n+1)!, whereas Pari/GP computes a floating-point
approximation. For 1000-bit precision with input pi^2, and the parameters of
Pari/GP, this requires computing Bernoulli numbers of 3800 bits.

We should compute floating-point approximations of the Bernoulli numbers in
MPFR too, but this will require redoing the error analysis, which is non
trivial.

Paul

-- 
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.

Re: [sage-devel] Re: Issue With Implementation of Gamma Function

[Zimmermann Paul]

   William,

thank you for putting me in cc.

> From: William Stein 
> Date: Wed, 12 Feb 2014 06:01:29 -0800
> 
> On Wed, Feb 12, 2014 at 4:55 AM,   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, (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.0100001011010101100011011000101100111000111"
   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.

[sage-devel] Re: real literals and IEEE-754

[Zimmermann Paul]

   Nils,

> Date: Mon, 6 Jan 2014 13:18:07 -0800 (PST)
> From: Nils Bruin 
> 
> On Monday, January 6, 2014 12:46:19 PM UTC-8, Zimmermann Paul wrote:
> 
> > What about getting rid of real literals? 
> >
> 
> I think they exist mainly to let
> 
> RealField(200)(1e-20)
> 
> work in the first place: The parser needs to generate code to instantiate 
> the constant 1e-20 somewhere and that constant should remember its original 
> string representation to allow convenient creation of an element in 
> RealField(200).

but RealField(200)("1e-20") already does the job:

sage: RealField(200)("1e-20")
1.00e-20

> Consider for instance
> 
> sage: RealField(200)(a+(1e-800))
> 9.9994515327145420957165172950370278739244710772e-21
> 
> (adding this constant doesn't change the 53-bit representation of a but it 
> does ensure the resulting number is not considered a RealLiteral any more).

this demonstrates another inconsistency:

sage: a=1e-20
sage: RealField(200)(a)
1.00e-20
sage: RealField(200)(a+0)
9.9994515327145420957165172950370278739244710772e-21
sage: RealField(200)(1*a)
9.9994515327145420957165172950370278739244710772e-21

> So back to your question: I expect that real literals can't be removed 
> because it would make creating high precision constants too much of a pain.

I don't think so. We should educate the user:

* if she/he writes 1e-20, she/he should be aware that (like in other languages)
  1e-20 is not exactly representable in binary, thus will be approximated to
  the current precision

* if she/he wants to define a high precision constant, use "1e-20"

Paul

-- 
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.

[sage-devel] real literals and IEEE-754

[Zimmermann Paul]

   Hi,

[since I am not subscribed to sage-devel, please keep me in cc]

the concept of real literals, which was intended (as far as I understand)
to keep exactly track of inputs like "1e-20", leads to the following:

sage: a=RealField(53)(1e-20)
sage: Reals(200)(a)
1.00e-20
sage: Reals(200)(a.exact_rational())
9.9994515327145420957165172950370278739244710772e-21

I consider this to be a bug.

In the C language, when one writes "double a = 1e-20; printf ("%.20e\n", a);",
everybody who is fluent with IEEE-754 2008 knows that 10^-20 is not exactly
representable in the binary64 format, thus nobody is surprised to see the
output 9.99945153e-21.

Moreover real literals seem to be created only with 53-bit precision:

sage: Reals(200)(RealField(52)(1e-20))
1.956165483594623726717277804672348392078969e-20
sage: Reals(200)(RealField(53)(1e-20))
1.00e-20
sage: Reals(200)(RealField(54)(1e-20))
1.20384909906835972161728642085058275023e-20

What about getting rid of real literals?

Paul



-- 
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.

Re: [sage-devel] imag(CC(infinity)) is 0?

[Zimmermann Paul]

   Peter,

> I think a case could be made for having two versions of the current RR: one 
> like the current one (more like a model of the extended real line) and one 
> where overflow or division by zero raises an exception instead of returning 
> +/- infinity (more like a model of the usual real numbers).  I can think of 
> at least two reasons for having the latter:

note that MPFR internally raises exceptions, and it would be easy to provide
Sage functions to query/set/clear the exception flags.

Paul

-- 
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.

Re: [sage-devel] imag(CC(infinity)) is 0?

[Zimmermann Paul]

   William,

[please forward to sage-devel since I'm not sure I'm allowed to post there]

> The implementation of RR and CC in Sage are a very direct wrapping of
> MPFR, which is the most well-thought out efficient implementation of
> floating point real numbers I've ever seen.  It is worth visiting
> http://www.mpfr.org/mpfr-current/mpfr.html and searching for
> "infinity".

thank you for the pointer. However for questions whether Inf should be in RR
or not, the MPFR documentation does not say much. About this topic, I
recommend to look at the discussions around the P1788 standard for interval
arithmetic : .

Best regards,
Paul

-- 
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.

[sage-devel] Re: Build Error in gf2x package on AMD Athlon XP 2600+

[Zimmermann Paul]

   Jean-Pierre,

[please forward to sage-devel, since I'm not subscribed, thus my mail will
 be rejected]

thank you for forwarding us that message:

> Date: Fri, 20 Sep 2013 11:36:26 -0700 (PDT)
> From: Jean-Pierre Flori 
> Cc: Zimmermann Paul 
> 
> [1:multipart/alternative Hide]
> [1/1:text/plain Hide]
> Hey,
> 
> I think you should report this to the gf2x dev as well... you never know.
> I've CC'ed Paul Zimmermann in case he cares.

yes we care!

> Best,
> JP
> 
> On Friday, September 20, 2013 3:01:13 PM UTC+2, Andrew Fiori wrote:
> >
> > This is likely a build error bug on all systems which do not support SSE2 
> > but where building is done using the current version of gcc. It may also 
> > manifest as a runtime error on binary builds installed on systems which do 
> > not support SSE2 (but I have not tested that).
> >
> > I am submitting this here rather than to bug tracker for several reasons:
> >- It is unlikely anyone cares enough to fix this properly upstream 
> > (which is where the bug really is).
> >- The build error only effects building from source on systems that are 
> > out of date and will becoming even less common going forward.
> >
> > I am mentioning it at all, because it will help anyone who comes across 
> > this error in the future. The same type of error may cause somewhat 
> > confusing runtime errors in other packages (bug tracker searches suggest 
> > this type of issue has come up before in other packages).
> >
> > The bug:
> >   The bug is somewhere in the gf2x build scripts "src/config/acinclude.m4" 
> > or "src/configure" somehow it decides whether or not the system supports 
> > sse2 (look around line 75 of acinclude.m4. It appears to do this by 
> > checking if gcc can compile something that uses sse2. It turns out that it 
> > can... even if the system you are on cannot. Consequently it tries to build 
> > everything with sse2. This causes a problem when it tries to tune itself 
> > "src/src/tune-lowlevel.pl". Around line 70-78 it tries to run the 
> > compiled code. Most of these don't require sse2 really, and so work, some 
> > do (mul3t, mul3k, mul3k2) these fail, it then can't pick which one to use, 
> > and terminates building.
> >If there are other places in sage where similar sse2 detection is used, 
> > but which are not executed at build time (most things are not), then sage 
> > will encounter unexpected runtime bugs when pieces of code that can benefit 
> > from sse2 are used.
> >
> >
> > The solution:
> >gf2x doesn't require sse2 to work. The easiest way to build it on a 
> > system where the above bug is encountered is to modify:
> >spkg-install:38 to become:
> >./configure --enable-sse2=no --prefix="$SAGE_LOCAL"
> > You don't really want to do this in the main branch though, as this slows 
> > things down for the 99+% of people who don't have an archaic system they 
> > are trying to run sage on.
> > A better solution is to add an if statement at this place and environment 
> > variable globally to disable sse2 when it doesn't exist.
> > Alternatively one could invest the time and figure out then fix gf2x build 
> > system.
> >
> > - Andrew

I'm puzzled because in the version of GF2X shipped with Sage 5.11,
src/config/acinclude.m4 defines SSE2_EXAMPLE on lines 75-79, then
it is used in AC_RUN_IFELSE (not AC_COMPILE_IFELSE) on lines 90 and 94.
Thus I would say it is a bug in gcc. Please can you send us the corresponding
config.log file from gf2x?

Best regards,
Paul

-- 
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.

[sage-devel] Annoying inconsistency of type in doctests

[Zimmermann Paul]

this is now ticket http://trac.sagemath.org/sage_trac/ticket/14466

Paul

-- 
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?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.