import MPFR_PREC_MAX from mpfr.h instead of hard coding it to the 32 bit limit

[williamstein]

The discussion below is besides the point. The main issue is that we define MPFR_PREC_MAX in Sage's sources instead of pulling it in from mpfr.h. We do hard code the 32 bit value, so on 64 bit boxen we limit the user to much lower precision than is actually technically feasible as pointed out below.

Cheers,

Michael

> >  CODE:
>
> >  s = pi.str(3000000*log(10,2))
> >  o = open('/Users/ericahls/Desktop/file.txt','w')
> >  o.write(str(s))
> >  o.close()
>
> >  --- Trying to get out to the farthest decimal point of PI I can.
>
> >  Error message:
>
> >  Traceback (most recent call last):    o.write(str(s))
> >   File "/Applications/sage/local/lib/python2.5/site-packages/sage/
> >  functions/functions.py", line 140, in str
> >     raise ValueError, "Number of bits must be at most 2^23."
> >  ValueError: Number of bits must be at most 2^23.
>
> >  ----If i Put 2000000 instead 0f 3000000 the equation works.  much over
> >  2 million the equation breaks dwon.
>
> I think that 2^23 is a bound in mpfr, and Sage uses mpfr to
> compute digits of pi.  I don't know if one can compute more than
> about 2^23 digits using mpfr.
>
> William


Yes we can. The issue was that MPFR used the stack instead of the heap
for certain operations [even when told not to use alloca] and would
smash it therefore with large number of digits. That has been fixed in
MPFR 2.3.1 (which we include) and all we need to do is to raise or
remove the limit in our code and do some testing. Care to open a
ticket?

 -- Mabshoff

CC: @zimmermann6 @kiwifb

Component: basic arithmetic

Author: Mike Hansen, Frédéric Chapoton

Branch/Commit: e49e83e

Reviewer: Paul Zimmermann, François Bissey

Issue created by migration from https://trac.sagemath.org/ticket/2567

[aghitza]

comment:1

As far as I can tell, there is still a limit, namely the maximum precision MPFR_PREC_MAX, which is 2^24 on my 32-bit machine.

I guess this means that we can relax the current 2^23 limit a little bit. However, the MPFR manual says: "Warning! MPFR needs to increase the precision internally, in order to provide accurate results (and in particular, correct rounding). Do not attempt to set the precision to any value near MPFR_PREC_MAX, otherwise MPFR will abort due to an assertion failure."

[fredrik-johansson]

comment:2

You could use the pi function in mpmath; as far as I know, it is limited only by available memory. I just verified that computing 100 million digits works on a 32-bit system. The last time someone compared, it was also about three times faster than MPFR (but probably less memory efficient).

Or you could perhaps use this code which is even faster: http://gmplib.org/pi-with-gmp.html

[sagetrac-mabshoff]

comment:3

The problem is not so much the computation of Pi itself, but that we limit the maximum amount of precision on can ask for when using MPFR. It seems the limit is good for 32 bits, but for 64 bits MPFR is also basically limited by memory only.

The solution to this ticket might be to import MPFR_PREC_MAX from mpfr.h or wherever it is defined and then the problem will magically go away.

Cheers,

Michael

[zimmermann6]

comment:4

The solution to this ticket might be to import MPFR_PREC_MAX from mpfr.h

right. It is defined there as 2³¹-1 on a 32-bit machine and 2⁶³-1 on a 64-bit machine.

[zimmermann6]

comment:5

You could use the pi function in mpmath; as far as I know, it is limited only by available memory. I just verified that computing 100 million digits works on a 32-bit system. The last time someone compared, it was also about three times faster than MPFR (but probably less memory efficient).

I am curious. Can you give real timings? Here is what I get on sage.math:

zimmerma@sage:~/mpfr-2.4.0/tests$ pwd
/home/zimmerma/mpfr-2.4.0/tests
zimmerma@sage:~/mpfr-2.4.0/tests$ time ./tconst_pi 332192809 0 0

real    42m21.420s
user    41m55.500s
sys     0m25.910s

This is without using the new FFT code we designed with Gaudry and Kruppa, which should give a
twofold speedup. Anyway if mpmath can compute 100 million digits in less than 15 minutes, I am
really impressed!

[fredrik-johansson]

comment:6

In mpmath on an Athlon 3700+ 2.21 GHz, 1 GB RAM,

10**6 digits took 5.96 seconds (4.77 calc, 1.19 convert)

10**7 digits took 109.45 seconds (82.16 calc, 27.28 convert)

10**8 digits took 2184.68 seconds (1634.65 calc, 550.02 convert)

I can't compare with MPFR on the same computer at the moment, due to network problems. (With an old version of sage, 3.0.2, %time str(pi.n(10**6*log(10.,2))) takes 43.06 s, but I don't trust that number).

This is the result reported by Ondrej a few months ago:

E.g. in sympy+gmpy:

In [3]: time a = pi.evalf(10**6)
CPU times: user 5.13 s, sys: 0.04 s, total: 5.17 s
Wall time: 5.19 s

Sage 3.1.1:

sage: time a = pi.n(digits=10**6)
CPU times: user 14.06 s, sys: 0.06 s, total: 14.12 s
Wall time: 14.34 s 

This is without using the new FFT code we designed with Gaudry and Kruppa, which should give a twofold speedup. Anyway if mpmath can compute 100 million digits in less than 15 minutes, I am really impressed!

Mpmath relies directly on multiplication of GMP mpz's. If it is faster than MPFR, that is entirely due to using a better formula. Before using the Chudnovsky series, mpmath used AGM which has better theoretical complexity but was 3x slower up to at least 1M digits.

[sagetrac-mabshoff]

comment:7

This discussion is very interesting, but it should not happen in trac, but on some mailing list since there people actually tend to see it. Trac's audience is rather limited and the comment section isn't meant for discussions :)

I have changed the ticket to reflect Paul's pointer about pulling in MPFR_PREC_MAX from mpfr.h

Cheers,

Michael

[zimmermann6]

comment:8

Mpmath relies directly on multiplication of GMP mpz's. If it is faster than MPFR, that is entirely due to using a better formula. Before using the Chudnovsky series, mpmath used AGM which has better theoretical complexity but was 3x slower up to at least 1M digits.

yes in a previous version MPFR did use the Chudnovsky series, but it only gives a fixed number of terms per iteration, whereas the current AGM-based code doubles the accuracy at each iteration, thus is asymptotically better. Also when the division in GMP is really O(M(n)) the current MPFR code should be much faster. However we should use the Chudnovsky series for small precision and the AGM for large precision.

[zimmermann6]

comment:9

I guess this ticket is now obsolete.

[mwhansen]

Author: Mike Hansen

[mwhansen]

comment:10

I've attached a patch which removes our hard coded value.

[zimmermann6]

Reviewer: Paul Zimmermann

[zimmermann6]

comment:11

I've attached a patch which removes our hard coded value.

Mike, I guess this has to be rebased due to #8157. Or #8157 should be undone since your patch is
better (in particular on 64-bit machines we should be able to compute with more that 2³¹ bits).

Paul

[mwhansen]

comment:12

Oops, I knew that I had seen that ticket in the title a few days ago, but wasn't able to find it again :-) I'll rebase this one.

The warning in the comment seems not to apply now.

[zimmermann6]

comment:13

Mike,

The warning in the comment seems not to apply now.

which warning do you mean?

[mwhansen]

comment:14

Sorry, I knew after posting that I should have been more specific. There's a comment in the source code that said that things totally break if we use the value from mpfr.h.

[mwhansen]

comment:15

I've rebased this and posted a new patch. This needs to be tested on a 32-bit machine.

[zimmermann6]

comment:16

while trying to review that ticket, I get with the input in the description:

sage: s = pi.str(3000000*log(10,2))
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)

/localdisk/tmp/sage-4.4.4/<ipython console> in <module>()

/localdisk/tmp/sage-4.4.4/local/lib/python2.6/site-packages/sage/structure/element.so in sage.structure.element.Element.__getattr__ (sage/structure/element.c:2632)()

/localdisk/tmp/sage-4.4.4/local/lib/python2.6/site-packages/sage/structure/parent.so in sage.structure.parent.getattr_from_other_class (sage/structure/parent.c:2835)()

/localdisk/tmp/sage-4.4.4/local/lib/python2.6/site-packages/sage/structure/parent.so in sage.structure.parent.raise_attribute_error (sage/structure/parent.c:2602)()

AttributeError: 'sage.symbolic.expression.Expression' object has no attribute 'str'

Is that normal?

Paul

[mwhansen]

comment:33

Paul, does MPFR_PREC_MAX now reflect the 2^31 - 257 limit in 3.1.0? For 64-bit, it seems to be set at 2^63.

[zimmermann6]

comment:34

in 3.1.x MPFR_PREC_MAX is 2⁶³-1 on 64-bit computers. In the development version of MPFR it is set to 2³¹-257 in 32-bit and 2⁶³-257 in 64-bit mode.

Paul

[zimmermann6]

comment:35

any progress on this ticket?

Paul

[fchapoton]

Branch: u/chapoton/2567

[fchapoton]

Commit: e49e83e

[fchapoton]

New commits:

e49e83e

trac 2567 MPFR_PREC_MAX

[fchapoton]

comment:42

Well, I made a branch and this seems to work. Anybody interested ?

[kiwifb]

Changed reviewer from Paul Zimmermann to Paul Zimmermann, François Bissey

[kiwifb]

comment:43

Well, it looks good to me. I cannot see clear objections in the earlier ticket comments so let's see what happens when we include it for real.

The author field may need updating though.

[fchapoton]

comment:44

voilà, voilà.

[fchapoton]

Changed author from Mike Hansen to Mike Hansen, Frédéric Chapoton

[kiwifb]

comment:45

Comme une lettre à la poste :)

[vbraun]

Changed branch from u/chapoton/2567 to e49e83e