computing the number of partitions of an integer

[William Stein]

Hi,

Since we were recently discussing speed of computing
Bell numbers, maybe it would be interesting to discuss
computing the number of partitions of an integer as a sum
of other positive integers, i.e.,

sage: number_of_partitions(5)
7
sage: v = list(partitions(5)); v
[(1, 1, 1, 1, 1), (1, 1, 1, 2), (1, 2, 2), (1, 1, 3), (2, 3), (1, 4), (5,)]
sage: len(v)
7

Currently, I'm going through the "Mathematica tour" and
making a SAGE version of it. The Mathematica tour is
a free pdf download available here:

http://documents.wolfram.com/mathematica/Mathematica_V5_Book.zip

The first thing Mathematica can do that SAGE can't
is compute number_of_partitions(10^9) "in a few seconds -- a frontier
number theory calculation" (page 5). In SAGE it takes around
80-100 seconds to compute

sage: number_of_partitions(10^8, algorithm="pari")

Actually, on sage.math, Mathematica takes 67 seconds to
compute the number of partitions for 10^9 and about 10
seconds to do 10^8.

QUESTIONS: Why is Mathematica about 10 times faster than PARI
at this? What are the best ways to compute the number
of partitions of n? Is it a calculation involving fast arithmetic
with dense polynomials of large degree, which would be best done
using the upcoming FLINT library or NTL?

It's a little embarasing that the first thing Mathematica does
that SAGE doesn't do reasonably quickly in their tour is a
number theory calculation.

-- William

--
William Stein
Associate Professor of Mathematics
University of Washington
http://www.williamstein.org

[William Stein]

On 7/25/07, Timothy Clemans <timothy...@gmail.com> wrote:
> How do you in general find out how say Magma or Mathematica does something?

Short answer: it is often virtually impossible -- that's one of the main
reasons for the existence of SAGE. Read this page for the official
Mathematica stance on this sort of question:

http://reference.wolfram.com/mathematica/tutorial/WhyYouDoNotUsuallyNeedToKnowAboutInternals.html

-- William

[Mike Hansen]

The fastest way I know of to compute the number of integer partitions
is to use the Hardy-Ramanujan-Rademacher formula ( see Rademacher
series on http://en.wikipedia.org/wiki/Integer_partition ), and I'm
pretty sure this is what PARI uses. From modules/part.c:

* This program is basically the implementation of the script
*
* Psi(n, q) = local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b);
* (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c))
* L(n, q) = if(q==1,1,sum(h=1,q-1,if(gcd(h,q)>1,0,cos((g(h,q)-2*h*n)*Pi/q))))
* g(h, q) = if(q<3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2)))
* part(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q)))

I'm not sure where the bottleneck in PARI is since I can't imagine
Mathematica uses a different method to compute the number of
partitions.

--Mike

P.S. I'm going to push hard this weekend to get the combinatorics
stuff I've been working on to you.

[Alec Mihailovs]

From: "Mike Hansen" <mha...@gmail.com>

> I'm not sure where the bottleneck in PARI is since I can't imagine
> Mathematica uses a different method to compute the number of
> partitions.

I don't know what is used in the latest Mathematica version, but originally
NumberOfPartitions function in the Combinatorica package used the recursion
with pentagonal numbers, see

http://www.cs.uiowa.edu/~sriram/Combinatorica/NewCombinatorica.m

Alec

[Mike Hansen]

I just found this on
http://reference.wolfram.com/mathematica/note/SomeNotesOnInternalImplementation.html:

"PartitionsP[n] uses Euler's pentagonal formula for small n, and the
non-recursive Hardy-Ramanujan-Rademacher method for larger n."

--Mike

[Nick Alexander]

"William Stein" <wst...@gmail.com> writes:

> QUESTIONS: Why is Mathematica about 10 times faster than PARI
> at this? What are the best ways to compute the number
> of partitions of n? Is it a calculation involving fast arithmetic
> with dense polynomials of large degree, which would be best done
> using the upcoming FLINT library or NTL?

Please correct me if I'm crazy, but isn't there an asymptotic formula
due to Hardy and Rademacher that can evaluate $P(n)$ to a very high
accuracy very quickly? Surely both of these packages implement such a
convergent series approach, and it is possible that SAGE has faster
real arithmetic and could do this faster.

Again, I may be completely incorrect -- please let me know.

Nick

[Bill Hart]

Pari implements it, but incorrectly:

sage: number_of_partitions(5*10535+4,algorithm="pari")
132775694853416614410709359082850304628573507097
148711672236023178345995078715426574030466347126
367130008280402558755564770977624160764152691390
102001213796297899514590335375857238756973648770
246446295491295636766189177742381389268656730071
68971398574

But:

sage: number_of_partitions(5*10535+4)
132775694853416614410709359082850304628573507097
148711672236023178345995078715426574030466347126
367130008280402558755564770977624160764152691390
102001213796297899514590335375857238756973648770
246446295491295636766189177742381389268656730071
68971398575

At least GAP returns the right answer!! But after some time:

sage: number_of_partitions(5*10535+4)
---------------------------------------------------------------------------
<type 'exceptions.RuntimeError'> Traceback (most recent call
last)

/home/wbhart/flint/trunk/<ipython console> in <module>()

/home/was/s/local/lib/python2.5/site-packages/sage/combinat/
combinat.py in number_of_partitions(n, k, algorithm)
1245 if algorithm == 'gap':
1246 if k==None:
-> 1247 ans=gap.eval("NrPartitions(%s)"%(ZZ(n)))
1248 else:
1249 ans=gap.eval("NrPartitions(%s,%s)"%(ZZ(n),ZZ(k)))

/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/gap.py
in eval(self, x, newlines, strip)
298 if len(x) == 0 or x[len(x) - 1] != ';':
299 x += ';'
--> 300 s = Expect.eval(self, x)
301 if newlines:
302 return s

/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/
expect.py in eval(self, code, strip, **kwds)
519 code = code.strip()
520 try:
--> 521 return '\n'.join([self._eval_line(L, **kwds) for L
in code.split('\n') if L != ''])
522 except KeyboardInterrupt:
523 self._keyboard_interrupt()

/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/gap.py
in _eval_line(self, line, allow_use_file, wait_for_prompt)
482 return ''
483 else:
--> 484 raise RuntimeError, message
485
486 except KeyboardInterrupt:

<type 'exceptions.RuntimeError'>: Unexpected EOF from Gap executing
NrPartitions(52679);

Bill.

[Bill Hart]

There are two tricks I can see that are required to make this faster.

Firstly notice that L(n,q)*Psi(n,q) is relatively small once q gets
beyond a certain point. Thus L(n,q)*Psi(n,q) can be computed using
ordinary double precision for q greater than some very small limit
(depending on n). If one knew how fast this series converges (which
must have been worked out somewhere), one could even progressively
reduce the precision as the computation went, for even greater speed.
For example the following Pari script should compute all but the last
150 or so digits of P(10^9) quite quickly:

Psi(n, q) = local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b);

(sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c))
L(n, q) = if(q==1,1,sum(h=1,q-1,if(gcd(h,q)>1,0,cos((g(h,q)-2*h*n)*Pi/
q))))

g(h, q) = if(q<3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2)))

n=1000000004
\p36000
a1 = L(n,1)*Psi(n,1);
\p18000
a2 = L(n,2)*Psi(n,2);
\p12000
a3 = L(n,3)*Psi(n,3);
\p9000
a4 = L(n,4)*Psi(n,4);
\p8000
a5 = L(n,5)*Psi(n,5);
\p6000
a6 = L(n,6)*Psi(n,6);
\p6000
a7 = L(n,7)*Psi(n,7);
\p5000
a8 = L(n,8)*Psi(n,8);
\p4000
a9 = sum(y=9,14,L(n,y)*Psi(n,y));
\p3000
a10 = sum(y=15,17,L(n,y)*Psi(n,y));
\p2000
a11 = sum(y=18,34,L(n,y)*Psi(n,y));
\p1000
a12 = sum(y=35,300,L(n,y)*Psi(n,y));
round(a1+a2+a3+a4+a5+a6+a7+a8+a9+a10+a11+a12)

The second trick is that computing gcd(h,q) is costly. One should
factor q and then sieve for all h not coprime to q in short cache
friendly segments, especially for very large n.

But I don't believe Mathematica uses this Rademacher series thing
anyway. If you look at the inner most loop, that computation involving
frac must take at least 80 cycles in double precision. But for n =
10^9, that expression must get evaluated about 2*10^11 times. That's
about 1.6*10^13 cycles, which on sage.math has go to take about
10000s. Even if I'm out by a factor of 10 with the cycles, mathematica
clearly doesn't do this.

Perhaps if one had a fast way of evaluating the Dedekind sum, one
might have a chance.

Bill.

[Jonathan Bober]

On Thu, 2007-07-26 at 13:10 -0700, Bill Hart wrote:
>
> Perhaps if one had a fast way of evaluating the Dedekind sum, one
> might have a chance.
>
> Bill.
>

I think this gives a faster way to compute it:

Write the sum as

s(h,k) = sum_{j=1}^{k-1} j/k [hj/k - floor(hj/k) - 1/2]

(This isn't strictly correct in general, but in our case hj/k will never
be an integer, so we are ok.)

Then if we separate this into three different sums and use some simple
summation formulas, we get

s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).

(To compute the floor in the inner sum, you can just use integer
division.)

[Bill Hart]

On 26 Jul, 23:18, Jonathan Bober <jwbo...@gmail.com> wrote:

> s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
> (1/k) sum_{j=1}^{k-1} j*floor(hj/k).
>
> (To compute the floor in the inner sum, you can just use integer
> division.)

Yes, this will be faster. Of course integer division is not faster
than floating point division, but since we are doing a sum from j = 1
to k-1 we only need to know when floor(hj/k) increases by 1. This is
simply a matter of adding h each iteration and seeing if we have gone
over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
and so on.

This gets the inner computation down to about 8 cycles or so on
average. Not enough to beat Mathematica though, unless I made a
mistake in my back of the envelope computation somewhere.

Bill.

[Bill Hart]

OK, apparently the Dedekind sums should be done using an algorithm due
to Apostol. I haven't tracked it down yet, but this is clearly what we
need.

Bill.

[Bill Hart]

Here we go. Apparently g(h,q) = q*s(h,q) where s(h,q) is a dedekind
sum.

Then for h < q if r_0, r_1, ..., r_{n+1} are the remainders occurring
in the Euclidean algorithm when computing gcd(h,q) (of which there
should be about log q) then:

s(h,q) = 1/12 * sum_{i=1}^{n+1}((-1)^(i+1) (r_i^2 + r_{i-1}^2 + 1) /
(r_i * r_{i-1}) - ((-1)^n +1)/8)

This, with the other ideas I gave above, will definitely let us beat
Mathematica convincingly, as a simple back of the envelope calculation
shows. It will have the added benefit of allowing us to beat Magma at
something.

Bill.

[Bill Hart]

Alternatively you could just use the implementation here:

http://www.ark.in-berlin.de/part.c

which seems to only rely on mpfr and GMP. However, since the Pari
version was based on it, I suspect it may also need correction.

He does seem to adjust the precision as he goes, but I've no idea how
far above or below the required precision this is. However there is no
need to use a precision of 55 from a certain point on. That probably
forces it to use two doubles instead of one in the floating point
computations, which I think would be unnecesary. You can look on
mathworld for how well the Rademacher series converges.

It would probably be quicker to do a fresh implementation from scratch
given the application in mind. The code there may not be GPL'd also.

To check it works, one should at least look at the congruences modulo
5, 7, 11 and 13. The author of the mpfr version does seem to check
them mod 11 but for very small n and only for a very few values of n.

If the gcds prove to be a bottleneck, I have some code that is roughly
twice as fast as GMP for computing these. Certainly the code in the
mpfr version is suboptimal and needs some serious optimisation when
the terms in the dedekind sum fit into a single limb, which on a 64
bit machine is probably always, when n is of a reasonable size.

Another suggestion is that for sufficiently small values of h and or
k, it may be quicker to compute the dedekind sum directly rather than
use the gcd formula.

A lookup table would also be useful for the tail end of the gcd/
dedekind sum computation.

I'd be shocked if the computation for n = 10^9 couldn't be done in
under 10s on sage.math.

Bill.

[Pablo De Napoli]

I've tested it but seems to be buggy:: it works up to 156

./part 156
p(156)=
73232243759

but for 157 gives a floating point exception error
(and a gdb tracing says it is in the line 176 of
the source code

r2 = r0 % r1;

in function g

I've compiled it using

gcc part.c -g -lmpfr -lm -DTEST_CODE -o part

and my mpfr library is 2.2.1_p5 and gmp is 4.2.1

The code theres seems indeed to be GPL-ed, it has a copyright notice that says:

"(C) 2002 Ralf Stephan (ra...@ark.in-berlin.de).
* See part.pdf on the same path for a summary of the math.
* Distributed under GPL (see gnu.org). Version 2002-12."

The pdf is:

http://www.ark.in-berlin.de/part.pdf

Pablo

[William Stein]

On 7/27/07, Pablo De Napoli <pde...@gmail.com> wrote:
>
> I've tested it but seems to be buggy:: it works up to 156
>
> ./part 156
> p(156)=
> 73232243759
>
> but for 157 gives a floating point exception error
> (and a gdb tracing says it is in the line 176 of
> the source code
>
> r2 = r0 % r1;
>
> in function g

I've attached a slightly modified version of part.c that seems to work
(it doesn't
crash and gives correct answers in all the cases I tested). I just changed the
code in g slighlty so if r1 is 0, then r2 is also set to 0. I compile
it using

gcc part.c -O3 -g -lmpfr -lm -o part -I/home/was/sage/local/include
-L/home/was/sage/local/lib/ -lgmp

where /home/was/sage is the path to my SAGE install.

Interestingly, when I time part.c it is not much better than CVS pari
(see below).

Bill Hart wrote:
>This, with the other ideas I gave above, will definitely let us beat
>Mathematica convincingly, as a simple back of the envelope calculation
>shows. It will have the added benefit of allowing us to beat Magma at
>something.

Magma is already terrible at computing NumberOfPartitions, though
not as bad as GAP (which is even worse):

mwas@ubuntu:~$ magma
Magma V2.13-5 Fri Jul 27 2007 10:54:07 on ubuntu [Seed = 1720324423]
Type ? for help. Type <Ctrl>-D to quit.
> time k := NumberOfPartitions(10^5);
Time: 3.070

With the latest CVS pari:
? gettime; n=numbpart(10^5); gettime/1000.0
%1 = 0.01200000000000000000000000000
? gettime; n=numbpart(10^6); gettime/1000.0
%2 = 0.2080000000000000000000000000
? gettime; n=numbpart(10^7); gettime/1000.0
%3 = 6.973000000000000000000000000
? gettime; n=numbpart(10^8); gettime/1000.0

With part.c:

I use time ./part 100000 > out and record system + user time:
10^6 0.04 seconds
10^7 4.584 seconds
10^8 47.483 seconds

So, to clarify or summarize, Bill is there one thing that we could
change in part.c that would
speed it up? I realized you sort of answered this before, but the
sequence of previous
emails yesterday was rather long and may have got confused.

Thanks,
william

[William Stein]

---------- Forwarded message ----------
From: William Hart <har...@yahoo.com>
Date: Jul 27, 2007 11:58 AM
Subject: Re: [sage-devel] Re: computing the number of partitions of an integer
To: William Stein <wst...@gmail.com>

Some things to try (in order of importance):

* how long does it take to compute the first few (say
100 odd) terms. It should be nearly instant. If not
the precision is held too high for too long.

* much of the time is surely taken up by GMP in the
s(h,k) computation. One shouldn't use GMP for that. It
really needs to be optimally fast. I'd code a C
version, though to make any difference you probably
need to use GMP's longlong.h to get 128 bit
arithmetic. Obviously you still need to use GMP
outside that range. Frankly I can't see where the time
is going in there. Maybe I've just underestimated how
long it takes to do one of these, or underestimated
how many of them need to be done. Time how long each
s(h, k) is taking when the series has tailed off a bit
(say under 150 digits for n= 10^9). My bet is all the
time is being taken up by an inefficient
implementation of these.

* time how much time is spent in gcds. I have code
that will double the speed of those if they are a
bottleneck. But I sincerely doubt it.

Once you know where the bottleneck is, you can look
more carefully at that section of code and I might
have some more ideas.

BTW google groups will not accept postings from this
email address.

Bill.

[Pablo De Napoli]

Your version did work for me!
Pablo

[William Stein]

On 7/27/07, Pablo De Napoli <pde...@gmail.com> wrote:
>
> Your version did work for me!
> Pablo

Great. Is there any chance you could work on any of the optimization or
timing ideas Bill Hart suggested in the previous post? I won't be able to
work on this myself for a while.

-- William

[Ralf Stephan]

On Jul 27, 2:16 am, Bill Hart <goodwillh...@googlemail.com> wrote:
> Alternatively you could just use the implementation here:
>
> http://www.ark.in-berlin.de/part.c

In fact, that was my attempt of a C implementation of this Pari
script:

* Psi(n, q) = local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b);
* (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c))
* L(n, q) =
if(q==1,1,sum(h=1,q-1,if(gcd(h,q)>1,0,cos((g(h,q)-2*h*n)*Pi/q))))
* g(h, q) = if(q<3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2)))
* part(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q)))

which is also of my design. The C version then served as base for the
implementation of numbpart() within Pari proper. So, they all use the
same
algorithm (modulo some changes by Karim Belabas etc).

Just to get the history straight,
Ralf Stephan

[Jonathan Bober]

I've been working on a from-scratch implementation (attached). Right now
it runs a faster than Ralf Stephan's part.c, but not as fast as we would
like. (And it seems to work, although I can't guarantee that right
now.)

On my Core Duo 1.8 ghz , it computes p(10^8) in about 17 seconds,
compared to about 70 seconds for the part.c previously posted. However,
it took about 270 seconds for p(10^9). (I don't have Mathematica, so I
can't give that comparison.) On the other hand, I don't know how much
faster sage.math is than my laptop, but since it is 64 bit, it might run
the code much faster.

I think that there is still a good amount of optimization that can be
done to make this faster. Some things that might a lot help include
better error estimates for the tail end of the series (if such estimates
exist) and, in general, going over everything carefully to try to make
sure that no unneeded precision is ever used. (Once the code decides
that no more than 53 bits of precision are needed, it switches over to
computing with C doubles, and the rest of the computation finishes
"instantly".)

Note that this is C++ code, but it could be switched to pure C quite
easily, since is doesn't actually use any real C++.

[William Stein]

On 7/28/07, Jonathan Bober <jwb...@gmail.com> wrote:
> I've been working on a from-scratch implementation (attached). Right now
> it runs a faster than Ralf Stephan's part.c, but not as fast as we would
> like. (And it seems to work, although I can't guarantee that right
> now.)

Great work! Many many thanks. Could you resend (to me at least)
your part.cc but with a copyright statement (e.g., a GPL header)? Thanks.
I'd like to include it in SAGE-2.7.2 (asap) as an optional non-default
way to compute number_of_partitions.

This will put your code under revision control and making it easier for
other people to contribute to. Once we're confident that
your program is producing right answers, we can switch it over
to be the default number_of_partitions functions.

> On my Core Duo 1.8 ghz , it computes p(10^8) in about 17 seconds,
> compared to about 70 seconds for the part.c previously posted. However,
> it took about 270 seconds for p(10^9). (I don't have Mathematica, so I
> can't give that comparison.) On the other hand, I don't know how much
> faster sage.math is than my laptop, but since it is 64 bit, it might run
> the code much faster.

The wall times for Mathematica 6.0 doing this calculation on my 2.33Ghz laptop
under 32-bit linux (via Parallels) are:
10^8 12 seconds
10^9 95 seconds

On sage.math again (via mathematica 5.2):
10^8 10.88 seconds (wall time)
10^9 69.80 seconds (wall time)

The times for your code on the same system are:
10^7 wall: 1.2 seconds; cpu: 1.096
10^8 wall:16 seconds; cpu: 10.705
10^9 wall: 213.383; cpu: 145.673 seconds

On sage.math with your program:
10^7 cpu time: 1.124 seconds (cpu = user + sys; mostly user
in each case)
10^8 cpu time: 13.465 seconds
10^9 cpu time: 160.932 seconds

So you're basically within a factor of 2 of Mathematica, and I think you
easily now have the fastest open source implementation of counting
the number of partitions of an integer.

> I think that there is still a good amount of optimization that can be
> done to make this faster. Some things that might a lot help include
> better error estimates for the tail end of the series (if such estimates
> exist) and, in general, going over everything carefully to try to make
> sure that no unneeded precision is ever used. (Once the code decides
> that no more than 53 bits of precision are needed, it switches over to
> computing with C doubles, and the rest of the computation finishes
> "instantly".)
>
> Note that this is C++ code, but it could be switched to pure C quite
> easily, since is doesn't actually use any real C++.

I'm completely fine with C++ code.

-- William

[William Stein]

On 7/28/07, Jonathan Bober <jwb...@gmail.com> wrote:

> I've been working on a from-scratch implementation (attached). Right now
> it runs a faster than Ralf Stephan's part.c, but not as fast as we would
> like. (And it seems to work, although I can't guarantee that right
> now.)

Jonathan,

I've now included your code in SAGE (for 2.7.2) and made it available
(not by default) so it can be used as follows:
sage: time v=number_of_partitions(10^7, algorithm='bobber')
CPU times: user 0.86 s, sys: 0.00 s, total: 0.87 s
Wall time: 1.21

Your code agrees with Mathematica for 10^7 and 10^8, by the way:
sage: s=mathematica.eval('PartitionsP[10^7]')
sage: s = ''.join(s.replace('>','').split())
sage: s == str(v)
True
sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 13.89
sage: time v=number_of_partitions(10^8, algorithm='bobber')
CPU times: user 11.94 s, sys: 0.10 s, total: 12.04 s
Wall time: 16.44
sage: s = ''.join(s.replace('>','').split())
sage: s == str(v)
True

I want to release sage-2.7.2 soon, so please let me know if including part.cc
with it is OK.

-- William

[Jonathan Bober]

This is fine, except that my last name is Bober, not Bobber.

(Random question, since it just caused problems for me again - it there
a way to sign up for sage-devel with a non gmail email address?)

[William Stein]

On 7/28/07, Jonathan Bober <jwb...@gmail.com> wrote:
>
> This is fine, except that my last name is Bober, not Bobber.

Oops, I'll fix that. Thanks.

> (Random question, since it just caused problems for me again - it there
> a way to sign up for sage-devel with a non gmail email address?)

I don't know. I don't think so. Let me know if this isn't true.

William

[Michael Abshoff]

William Stein wrote:

Hello,

>
>> (Random question, since it just caused problems for me again - it there
>> a way to sign up for sage-devel with a non gmail email address?)
>
> I don't know. I don't think so. Let me know if this isn't true.

In my experience the participation in google groups doesn't requires a
gmail email address. I have always uses a non-gmail address. What exactly
is your problem?

Cheers,

Michael

>
>
> William
>
> >
>

[Jonathan Bober]

Attached is a somewhat improved version.

(I don't want to flood the mailing list with attachments, so I won't
send this file as an attachment to the list again after this.)

The main improvement is that the code is better commented and more
readable now. The secondary improvement is an increase in speed - I
think it's about 15% faster, or something like that, at computing
p(10^9) than the previous code. (Although I'm not exactly sure what the
phrase "15% faster" means, and I was having some weird issues with my
laptop when I timed the old code, so I'm not sure exactly how fast it
ran.) There are no improvements in automated testing, so I'm no more
confident of this correctness of this code than I was of the old code,
but it certainly seems to work.

I'm confident that this can be sped up more, but I'm not sure that we
can reach as low as 10 seconds to compute p(10^9). It certainly might be
possible, though.

[Pablo De Napoli]

great work, Jonathan!

I've tested, and I've found the following problems:

1) part 1 hangs

2) compiling with -Wall gives this warning

part.cc:865: warning: unused variable 'temp2'

3) part without arguments returns 42

Pablo

[William Stein]

On 7/29/07, Pablo De Napoli <pde...@gmail.com> wrote:
> great work, Jonathan!
>
> I've tested, and I've found the following problems:
>
> 1) part 1 hangs

Yep. In the SAGE rapper for part, I just do that
as a special case (i.e., dont' call part for n <= 1)

> 2) compiling with -Wall gives this warning
>
> part.cc:865: warning: unused variable 'temp2'
>
> 3) part without arguments returns 42

Hmm. This reminds me of the Hitchhiker's guide to the universe....

Anyway, if you look at the source code, you'll see that
the default input is 10 for some reason. This

William

[Bill Hart]

Wow!! Excellent work indeed.

In fact on 64 bit X86 systems you could actually use the 128 bit long
doubles to give you a little bit more precision (I believe it only
gives you 80 bits including exponent and sign, so probably 64 bit
mantissa).

It would be interesting to see the time for Mathematica on a 32 bit
X86 machine, since this would tell us if that is what they do.

Certainly I think you are right that any remaining optimization would
be in making sure it uses no unnecessary precision. Here is a page
giving information about the remainder of the series after N terms:

http://mathworld.wolfram.com/PartitionFunctionP.html (see eqn 26).

Also in Pari, I noted that the computation of the Psi function could
dramatically slow the whole computation. I was surprised to find it
figured in the runtime. It may be worthwhile checking if this is
slowing things down at all. It should be computed relatively quickly
if implemented correctly. The main issue was again using the minimum
possible precision.

Bill.

> part.cc
> 22KDownload

[Bill Hart]

On 31 Jul, 01:24, Bill Hart <goodwillh...@googlemail.com> wrote:
> It would be interesting to see the time for Mathematica on a 32 bit
> X86 machine, since this would tell us if that is what they do.

Doh! I should have read William's timings more carefully. He gives the
times for a 32 bit machine. So I guess Mathematica doesn't use 80 bit
long doubles on a 64 bit X86 then. Still it is an option for us.

Once there is a stable version of the new code which seems to give
correct results, I'll take a closer look and see if I can spot any
obvious speed improvements. I can't promise anything. I suspect my
fundamental mistake was not realising that you still needed quite a
bit of multi-precision code for quite a few terms. In fact now that I
think about it, I don't see why I thought you could compute all the
s(h,k)'s using single limb arithmetic.

It is the multi-precision stuff that is slowing it down, no doubt.
mpfr has a 15x overhead over ordinary double precision, even at 53
bits, or so I have read. I guess there is a lot of branching to ensure
the accuracy of arithmetic. Whilst that is needed for many
applications, it probably isn't here. Sadly there don't seem to be any
decent open source alternatives for when that accuracy is not
required. I have a similar problem in some code I am currently
writing. I need precisely quad precision, so mpfr is out of the
question.

Bill.

[didier deshommes]

2007/7/30, Bill Hart <goodwi...@googlemail.com>:

> I have a similar problem in some code I am currently
> writing. I need precisely quad precision, so mpfr is out of the
> question.

Hi Bill,
You might want to consider Yozo Hida's quaddouble C/C++ package here:
http://www.cs.berkeley.edu/~yozo/

There is also a wrapper for it in SAGE.

didier

[Bill Hart]

Hi Didier,

Thanks. I also just found:

http://www.nongnu.org/hpalib/

which fascinates me. Has anyone used it?

Bill.

On 31 Jul, 01:46, "didier deshommes" <dfdes...@gmail.com> wrote:
> 2007/7/30, Bill Hart <goodwillh...@googlemail.com>:

[boo...@u.washington.edu]

Last update, 2005?

[Bill Hart]

Yes, I've ended up using Hida and Bailey's quad-double package. Very
cool.

But the license just says not to use the LB name to promote any
derived product. Am I right in assuming this is GPL compatible, i.e.
because they are a public institution everything is automatically
public domain?

Bill.

[Pablo De Napoli]

Certainly no, not everything from a public institution is in the
public domain. This
should be analyzed case by case.
In case of doubt it would be better to ask the author.

Pablo

[Pablo De Napoli]

In this case, the license does not says that it is in the public domain
(but that it is a copyrighted work!), but you can use it as
"AS IS".
I think that the only condition that is imposed to us is
to include the declaimer.

Pablo

[Bill Hart]

Oh well, I don't understand all this licensing stuff. So do you
understand this license? What about derived works. Does that mean it
is not possible to modify this library and redistribute the modified
version?

In particular, as a C++ library it is no use to me unmodified. There
are also some functions I'd like to add, and the calling convention is
the opposite to what I use. So basically I'd have to modify it to make
it useful.

Bill.

On 31 Jul, 20:35, "Pablo De Napoli" <pden...@gmail.com> wrote:
> In this case, the license does not says that it is in the public domain
> (but that it is a copyrighted work!), but you can use it as
> "AS IS".
> I think that the only condition that is imposed to us is
> to include the declaimer.
>
> Pablo
>

> On 7/31/07, Pablo De Napoli <pden...@gmail.com> wrote:
>
> > Certainly no, not everything from a public institution is in the
> > public domain. This
> > should be analyzed case by case.
> > In case of doubt it would be better to ask the author.
>
> > Pablo
>

[Jonathan Bober]

On Mon, 2007-07-30 at 17:24 -0700, Bill Hart wrote:
> Wow!! Excellent work indeed.
>
> In fact on 64 bit X86 systems you could actually use the 128 bit long
> doubles to give you a little bit more precision (I believe it only
> gives you 80 bits including exponent and sign, so probably 64 bit
> mantissa).
>

Actually, even on my 32 bit core duo, the long double type seems to give
64 bits of precision, so using it might help a little. Do you have any
idea how to check at run/compile time what the precision of a double or
a long double is?

> It would be interesting to see the time for Mathematica on a 32 bit
> X86 machine, since this would tell us if that is what they do.
>
> Certainly I think you are right that any remaining optimization would
> be in making sure it uses no unnecessary precision. Here is a page
> giving information about the remainder of the series after N terms:
>
> http://mathworld.wolfram.com/PartitionFunctionP.html (see eqn 26).
>

That's what I'm using, but I took it straight from Rademacher's paper. I
think that what is needed is a careful analysis of how many terms (N)
will need to be computed to get the remainder to be less than, say, r
< .5. Then we know that the error in each term will have to be less than
(.5 - r)/N to guarantee that when we round, we will get the right
answer.

> Also in Pari, I noted that the computation of the Psi function could
> dramatically slow the whole computation. I was surprised to find it
> figured in the runtime. It may be worthwhile checking if this is
> slowing things down at all. It should be computed relatively quickly
> if implemented correctly. The main issue was again using the minimum
> possible precision.
>

I don't think that this is slowing things down much, but I could be
wrong. I think that most time is spent in the actual computation of
A(n,k), (this is the notation I use for the sum of exponentials -- I
think that is the notation that Mathworld uses also.) If I set s(h,k) to
return 1 every time, without computing anything, I only save about 20
seconds on a 2m 30s computation of p(10^9), I think.

[Mike Hansen]

From COPYING,

"Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met"

--Mike

On 7/31/07, Bill Hart <goodwi...@googlemail.com> wrote:
>

[Bill Hart]

Ah, that's better. Excellent. I feel much happier with this library
now.

Thanks.

Bill.

On 31 Jul, 21:38, "Mike Hansen" <mhan...@gmail.com> wrote:
> From COPYING,
>
> "Redistribution and use in source and binary forms, with or without
> modification, are permitted provided that the following conditions are
> met"
>
> --Mike
>

[Bill Hart]

I believe that the IEEE standard guarantees you 80 bits (though it's
only 64 bits of mantissa or something like that). The trouble is, you
aren't guaranteed the IEEE standard.

I've spent much time researching this, but either I didn't look at the
right websites, or this stuff isn't documented well.

At any rate, does autoconf allow you to check this? Unfortunately I
don't know about such things. I'm an admitted ignoramous when it comes
to anything to do with runtime/compile time checking.

Bill.

[Bill Hart]

I do highly recommend this quad double library by the way. And they've
implemented all manor of transcendental functions too!! The quad-
doubles would give you 206 bits, even on your machine.

Bill.

On 31 Jul, 21:33, Jonathan Bober <jwbo...@gmail.com> wrote:

[Bill Hart]

Or probably 212 actually. :-)

[Alec Mihailovs]

> Actually, even on my 32 bit core duo, the long double type seems to give
> 64 bits of precision, so using it might help a little. Do you have any
> idea how to check at run/compile time what the precision of a double or
> a long double is?

Being an assembler programmer, I can say definitely that all FPU registers
on x86 are 80-bit and all compilers that I know compile long double as
80-bit numbers.

Alec

[Alec Mihailovs]

From: "Alec Mihailovs" <al...@mihailovs.com>

> Being an assembler programmer, I can say definitely that all FPU registers
> on x86 are 80-bit and all compilers that I know compile long double as
> 80-bit numbers.

From other point of view, 80-bit real gives 64-bit precision in usual sense
(mantissa), while a 64-bit double has 53-bit precision in usual sense.

Alec

[William Stein]

Hi,

Just to extend this thread some more, a few remarks.

(1) quaddouble has been included in SAGE for several months
now, thanks to the hard work of Didier Deshommes and Robert
Bradshaw.

sage: RQDF
Real Quad Double Field
sage: RQDF(2).sin()
0.909297426825681695396019865911744842702254971447890268378973011

(2) Great work on figuring out the problem with mpfr in SAGE!
I have fixed this for sage-2.7.3.

(3) Here are the latest timings of SAGE versus Mathematica
on (a) my Intel OS X laptop, and
(b) on an unloaded sage.math (64-bit opteron), both with
mathematica 5.2. This is after rebuilding mpfr with -O2.

On sage.math:
SAGE:
sage: time n = number_of_partitions(10^7)
CPU times: user 0.73 s, sys: 0.00 s, total: 0.73 s
Wall time: 0.73
sage: time n = number_of_partitions(10^8)
CPU times: user 8.42 s, sys: 0.00 s, total: 8.42 s
Wall time: 8.42
sage: time n = number_of_partitions(10^9)
CPU times: user 105.82 s, sys: 0.00 s, total: 105.82 s
Wall time: 105.81

Mathematica:
sage: time s=mathematica.eval('PartitionsP[10^7]')

CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s

Wall time: 1.90

sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s

Wall time: 9.91
sage: time s=mathematica.eval('PartitionsP[10^9]')
CPU times: user 0.03 s, sys: 0.01 s, total: 0.04 s
Wall time: 70.43

Interesting, SAGE seems better at Mathematica for
smaller input.)

On my OS X core 2 duo 2.33Ghz laptop:

SAGE:
sage: time n=number_of_partitions(10^7)
CPU times: user 0.62 s, sys: 0.00 s, total: 0.62 s
Wall time: 0.62
sage: time n=number_of_partitions(10^8)
CPU times: user 6.99 s, sys: 0.02 s, total: 7.00 s
Wall time: 7.03
sage: time n=number_of_partitions(10^9)
CPU times: user 94.71 s, sys: 0.28 s, total: 94.99 s
Wall time: 95.56

MATHEMATICA:
sage: time s=mathematica.eval('PartitionsP[10^7]')

CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s

Wall time: 8.65

sage: time s=mathematica.eval('PartitionsP[10^8]')

CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 48.08
sage: time s=mathematica.eval('PartitionsP[10^9]')
CPU times: user 0.04 s, sys: 0.00 s, total: 0.05 s
Wall time: 350.30

Yep -- Mathematica 5.2 interestingly totally sucks at
computing the number of partitions on an Intel OSX
machine... and SAGE rocks.

-- William

[Justin C. Walker]

On Jul 31, 2007, at 18:36 , William Stein wrote:

>
> Hi,
>
> Just to extend this thread some more, a few remarks.

While we're extending this, here's my $0.02 canadian (:-})

Checking partition count computation:

On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
Mathematica 5.2 (PartitionsP[10^9]: 95.5115 s
Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
Jon Bober's code (i386: 'jb 1000000000'): 160 s

I wonder why our Core 2 Duo timings are so different?

All of the executables are "i386" binaries (32-bit x86).

One oddity: Jon Bober's code produced a value ending in
'30457526857797923685688339', while Mathematica and SAGE produced one
ending in '30457526831036667979062760', so they differ.

Justin

--
Justin C. Walker, Curmudgeon at Large
Institute for the Absorption of Federal Funds
-----------
If it weren't for carbon-14, I wouldn't date at all.
-----------

[William Stein]

On 7/31/07, Justin C. Walker <jus...@mac.com> wrote:
> While we're extending this, here's my $0.02 canadian (:-})

> > Yep -- Mathematica 5.2 interestingly totally sucks at
> > computing the number of partitions on an Intel OSX
> > machine... and SAGE rocks.
>
> Checking partition count computation:
>
> On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
> Mathematica 5.2 (PartitionsP[10^9]: 95.5115 s
> Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
> Jon Bober's code (i386: 'jb 1000000000'): 160 s
>
> I wonder why our Core 2 Duo timings are so different?

SAGE 2.7.2.1 is different than what I just tested, since I
had rebuilt mpfr with optimization turned on, which definitely
improves the timings to "95seconds". I am surprised
that you get 95.5115s out of mathematica, given that
it takes 350 seconds with Mathematica to do the same
calculation on my laptop. Let's race Mathematica's
at the workshop tomorrow :-).

> All of the executables are "i386" binaries (32-bit x86).
>
> One oddity: Jon Bober's code produced a value ending in
> '30457526857797923685688339', while Mathematica and SAGE produced one
> ending in '30457526831036667979062760', so they differ.

You probably have an old version of his code? I think he
emailed newer versions to me, which are what's in SAGE,
and which he might not have posted to the list.

-- William

[Jonathan Bober]

On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:
>
> On Jul 31, 2007, at 18:36 , William Stein wrote:
>
> >
> > Hi,
> >
> > Just to extend this thread some more, a few remarks.
>
> While we're extending this, here's my $0.02 canadian (:-})

[...]

> Checking partition count computation:
>
> On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
> Mathematica 5.2 (PartitionsP[10^9]: 95.5115 s
> Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
> Jon Bober's code (i386: 'jb 1000000000'): 160 s
>
> I wonder why our Core 2 Duo timings are so different?
>

That is puzzling. Are you sure that you have the latest version of the
code? The source file should have 'Version .3' at the top. You can get
the latest version, or at least that latest version minus some memory
leak fixes, via hg_sage.pull(), or from

http://www.math.lsa.umich.edu/~bober/partitions_c.cc

The last version I posted to the list was significantly slower than what
I have now (especially since I think I accidently had significant
improvements commented out in the last code I posted to the list.)

> All of the executables are "i386" binaries (32-bit x86).
>
> One oddity: Jon Bober's code produced a value ending in
> '30457526857797923685688339', while Mathematica and SAGE produced one
> ending in '30457526831036667979062760', so they differ.
>

I am fairly confident that the answer ends in '339', since I have seen
it so many times while testing my code. In fact, this is what is in the
graphic on page 5 (page 27 of the pdf) of the Mathematica book.

Wikipedia doesn't 'know' any congruences that p(10^9) should satisfy, so
I can't use that to check the value. I'll try to figure out what the
right answer should be.

If your willing, try computing p(10^9 + 139) with all three. The answer
should be divisible by 5, 7, and 11. (We should have been using this as
a test case all along.)

Incidentally, I while I was writing this email, some tests I was running
just finished. The latest version of my code and the current version of
sage both give the answer ending in '339'.

If Mathematica and SAGE (ie PARI?) both give the same wrong answer, then
there might be something fishy going on.

[Jonathan Bober]

On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:

[...]

> Checking partition count computation:
>
> On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
> Mathematica 5.2 (PartitionsP[10^9]: 95.5115 s
> Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
> Jon Bober's code (i386: 'jb 1000000000'): 160 s
>
> I wonder why our Core 2 Duo timings are so different?
>
> All of the executables are "i386" binaries (32-bit x86).
>
> One oddity: Jon Bober's code produced a value ending in
> '30457526857797923685688339', while Mathematica and SAGE produced one
> ending in '30457526831036667979062760', so they differ.
>

Addendum to last email:

It looks like the 2.7.2.1 is using the new code by default, so the
'sage' time above should be the time for the newest version of the code,
so I'm going to assume that the other timing is for an old version of my
code. (I running with algorithm='pari' right now and it is taking a
really long time.)

While I'm assuming things, I'm also going to assume that you must have
mixed up the output values, and thus Mathematica must be producing the
wrong answer, in which case its no longer fair to compare running times
to Mathematica. (Of course, Mathematica does have the number ending in
339 in the Mathematica book.)

[Jonathan Bober]

Thanks. I should probably update the code to make use of long doubles,
then.

Incidentally, answering my own previous question float.h defines

FLT_RADIX (the base of the floating point representation)

and

FLT_MANT_DIG
DBL_MANT_DIG
LDBL_MANT_DIG

which give the precision.

[Justin C. Walker]

On Jul 31, 2007, at 22:54 , Jonathan Bober wrote:
> On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:
>> On Jul 31, 2007, at 18:36 , William Stein wrote:

[snip]

>> On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
>> Mathematica 5.2 (PartitionsP[10^9]: 95.5115 s
>> Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
>> Jon Bober's code (i386: 'jb 1000000000'): 160 s
>>
>> I wonder why our Core 2 Duo timings are so different?
> That is puzzling. Are you sure that you have the latest version of the
> code? The source file should have 'Version .3' at the top.

Yup. Double-checked the file; and re-ran the test. Essentially the
same time. I'll check this tomorrow, along with the Mathematica oddity.

> You can get
> the latest version, or at least that latest version minus some memory
> leak fixes, via hg_sage.pull(), or from
>
> http://www.math.lsa.umich.edu/~bober/partitions_c.cc
>
> The last version I posted to the list was significantly slower than
> what
> I have now (especially since I think I accidently had significant
> improvements commented out in the last code I posted to the list.)

That was the the "Jon Bober" timing above. I'll try the latest one
later.

Justin

--
Justin C. Walker, Curmudgeon-At-Large

Institute for the Absorption of Federal Funds
--------

Men are from Earth.
Women are from Earth.
Deal with it.
--------

[William Stein]

On 8/1/07, Justin C. Walker <jus...@mac.com> wrote:
> > The last version I posted to the list was significantly slower than
> > what
> > I have now (especially since I think I accidently had significant
> > improvements commented out in the last code I posted to the list.)
>
> That was the the "Jon Bober" timing above. I'll try the latest one
> later.

Two comments for this thread:

(1) The PARI in SAGE is the latest stable release, which means
number_of_partitions(n, algorithm='pari') is *wrong* even on
some fairly small values of n.

(2) If you do
sage: hg_sage.apply('http://sage.math.washington.edu/home/was/sage.hg')
you'll get an update to my latest code (not 100% doc tested!) which
has the latest version of jon's code, and for which
number_of_partitions(n)
uses Jon's code by default.

William

[William Stein]

Hi,

Regarding the discrepancy in timings between me
and Justin Walker using OS X *intel* Mathematica,
it turns out I was running Mathematica on my
laptop under OS X via Rosetta, so those times should
be ignored. The timings I've posted under Linux are
all fine. SO, currently SAGE and Mathematica
on OSX Intel MacBook Pro 2.33Ghz take almost
exactly the same amount of time (95 seconds) to
compute p(10^9).

On 64-bit Linux Mathematica is still faster than
SAGE, so probably there is a 64-bit only precision
trick that is relevant...

-- William

[Justin C. Walker]

On Jul 31, 2007, at 10:54 PM, Jonathan Bober wrote:
> On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:
>> On Jul 31, 2007, at 18:36 , William Stein wrote:

[snip]

> That is puzzling. Are you sure that you have the latest version of the
> code?

I downloaded the .3 source and ran it on my 2.33 GHz Core 2 Duo. I
compiled the program with several "-O" settings, and ran them with
the argument '1000000000':

Opt none: 2m11s
Opt 3: 2m8 s
Opt 2: 2m10s
Opt 1: 2m8 s
Opt s: 2m6 s
Opt z: 2m7 s
Opt fast: 2m6 s

Not a lot of difference, but all at least 126 seconds. Any
suggestions to figure out why this seems so far off from your (and
other's) experience?

Also, FWIW, if I run this through sage ("time number_of_partitions
(10^9)"), the time is ~134s.

Justin

--
Justin C. Walker, Curmudgeon-at-Large
() The ASCII Ribbon Campaign
/\ Help Cure HTML Email

[Justin C. Walker]

On Aug 1, 2007, at 12:29 PM, Justin C. Walker wrote:

>
>
> On Jul 31, 2007, at 10:54 PM, Jonathan Bober wrote:
>> On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:
>>> On Jul 31, 2007, at 18:36 , William Stein wrote:
> [snip]
>> That is puzzling. Are you sure that you have the latest version of
>> the
>> code?
>
> I downloaded the .3 source and ran it on my 2.33 GHz Core 2 Duo. I
> compiled the program with several "-O" settings, and ran them with
> the argument '1000000000':

Adding to the background noise:

On a 2.33GHz Core 2 Duo (Mac OS X, 10.4.10):
Mathematica 6.0: 83s
Mathematica 6.0.1: 77s

Justin

--
Justin C. Walker, Curmudgeon-At-Large
Institute for the Enhancement of the Director's Income
--------
Experience is what you get
when you don't get what you want.
--------

[Jonathan Bober]

On Wed, 2007-08-01 at 12:29 -0700, Justin C. Walker wrote:
>
> On Jul 31, 2007, at 10:54 PM, Jonathan Bober wrote:
> > On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote:
> >> On Jul 31, 2007, at 18:36 , William Stein wrote:
> [snip]
> > That is puzzling. Are you sure that you have the latest version of the
> > code?
>
> I downloaded the .3 source and ran it on my 2.33 GHz Core 2 Duo. I
> compiled the program with several "-O" settings, and ran them with
> the argument '1000000000':
>
> Opt none: 2m11s
> Opt 3: 2m8 s
> Opt 2: 2m10s
> Opt 1: 2m8 s
> Opt s: 2m6 s
> Opt z: 2m7 s
> Opt fast: 2m6 s
>
> Not a lot of difference, but all at least 126 seconds. Any
> suggestions to figure out why this seems so far off from your (and
> other's) experience?
>
> Also, FWIW, if I run this through sage ("time number_of_partitions
> (10^9)"), the time is ~134s.
>
> Justin
>

Since the time with the self-compiled code is so similar to the time
with the sage-compiled code, I would guess that you are linking to the
sage mpfr library, which wasn't compiled with -O2 (or at least linking
to another version of mpfr that wasn't compiled with optimizations.)

Also, did you double check the answers? If my code is giving the wrong
answer on your computer, there could be some processor specific bug
somewhere.

[Jonathan Bober]

On Tue, 2007-07-31 at 14:24 -0700, Bill Hart wrote:
> I do highly recommend this quad double library by the way. And they've
> implemented all manor of transcendental functions too!! The quad-
> doubles would give you 206 bits, even on your machine.
>
> Bill.

> URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
> -~----------~----~----~----~------~----~------~--~---

Have you found the quad double library to be faster than mpfr, or just
more convenient? I tried using it in the partition counting code, and it
actually slowed things down when I used it for all computations between
200 and 64 bits. Alternately, if I just use it between 200 and, say, 180
bits, it gives almost no change in speed.

[Jonathan Bober]

Incidentally, I should have mentioned here that I submitted a patch for
version .4, and also updated it at

http://www.math.lsa.umich.edu/~bober/partitions_c.cc

It uses long doubles now when then precision is small enough (and then,
later, just doubles like before), and the speedup is significant.

[Justin C. Walker]

On Aug 1, 2007, at 1:26 PM, Jonathan Bober wrote:
> On Wed, 2007-08-01 at 12:29 -0700, Justin C. Walker wrote:

[snip]

> Since the time with the self-compiled code is so similar to the time
> with the sage-compiled code, I would guess that you are linking to the
> sage mpfr library, which wasn't compiled with -O2 (or at least linking
> to another version of mpfr that wasn't compiled with optimizations.)

Oops. You're right. I forgot about this part of the discussion.
I'll get an 'mpfr' library compiled with more optimization.

Has anyone looked at the differences between the various optimization
choices? On Mac OS X, there is "-Ox", for x in {0,1,2,3,s,z}, "-O",
and "-fast".

> Also, did you double check the answers? If my code is giving the wrong
> answer on your computer, there could be some processor specific bug
> somewhere.

I'm verifying this now, but I think that my claim (that Mathematica
and sage agreed and computed a different value than your code) is a
typo. It should have said 'sage', which at the time was using Pari
(and according to Stein, was generating incorrect results).

Yup: all versions of Mathematica (5.2, 6.0, 6.0.1) agree with you :-}

Justin

--
Justin C. Walker, Curmudgeon-At-Large
Director

Institute for the Enhancement of the Director's Income
--------

"Weaseling out of things is what separates us from the animals.
Well, except the weasel."
- Homer J Simpson
--------

[Alec Mihailovs]

From: "Jonathan Bober" <jwb...@gmail.com>

>
> Incidentally, I should have mentioned here that I submitted a patch for
> version .4, and also updated it at
>
> http://www.math.lsa.umich.edu/~bober/partitions_c.cc
>
> It uses long doubles now when then precision is small enough (and then,
> later, just doubles like before), and the speedup is significant.

By the way, I just installed SAGE 2.7.2 (from source) on my Ubuntu 7.04
system, with AMD Athlon 64 X2 5200+ and tested the version included there
(using algorithm='bober'). It has the following feature - for input 2^32 it
gave the ValueError saying that the input should be less than 2^64-2
(evaluated).

Alec

[William Stein]

On 8/2/07, Alec Mihailovs <al...@mihailovs.com> wrote:
> > It uses long doubles now when then precision is small enough (and then,
> > later, just doubles like before), and the speedup is significant.
>
> By the way, I just installed SAGE 2.7.2 (from source) on my Ubuntu 7.04
> system, with AMD Athlon 64 X2 5200+ and tested the version included there
> (using algorithm='bober'). It has the following feature - for input 2^32 it
> gave the ValueError saying that the input should be less than 2^64-2
> (evaluated).

Thanks. I've fixed this for sage >= 2.7.3 in the attached patch.

William

[Bill Hart]

That's interesting. I haven't done any speed comparisons. But it seems
impossible that a carefully coded fixed-precision library should not
be faster than a multi-precision library.

Did you use double-doubles when you got below 106 bits?

I *have* made the assumption here that these guys know how to code
efficiently. At some point, you are right, I should do some timings.

At any rate, for my application, performance seems fine. It seems to
be within a factor of two of what I could get with hand crafted
purpose written code to do the same job. But it is definitely more
convenient.

Bill.

On 1 Aug, 21:30, Jonathan Bober <jwbo...@gmail.com> wrote:
> On Tue, 2007-07-31 at 14:24 -0700, Bill Hart wrote:
> > I do highly recommend this quad double library by the way. And they've
> > implemented all manor of transcendental functions too!! The quad-
> > doubles would give you 206 bits, even on your machine.
>
> > Bill.

> > URLs:http://sage.scipy.org/sage/andhttp://modular.math.washington.edu/sage/

[Alec Mihailovs]

From: "William Stein" <wst...@gmail.com>

> Thanks. I've fixed this for sage >= 2.7.3 in the attached patch.

Just upgraded SAGE to 2.7.3 and number_of_partitions works twice as fast as
in 2.7.2 with algorithm='bober' for 10^9 - it took 112 sec in 2.7.2 and only
56 sec in 2.7.3. This is in Ubuntu 7.04 on Dual Core 64-bit 2.6GHz Athlon
with 2 GB RAM. By the way, to compare with other systems, SAGE 2.7.2
installation took about 53 minutes.

I noticed the following interesting thing with time measuring: assigning
number_of_partitions to a variable (that shouldn't take time at all), took
about 10% or more of total time. For example,

sage: time number_of_partitions(10^7)

took about 0.8 sec while

sage: time a=number_of_partitions(10^7)

took slightly more than 1 sec. For 10^8 the corresponding times were 5 sec
and 5.45 sec.

That seems odd, isn't it?

Alec

[Jonathan Bober]

Since you've extended this thread some more, I'll mention that the
current version of the code that I am working on runs a little bit more
than twice as fast as that version. It's posted at

http://www.math.lsa.umich.edu/~bober/partitions_c.cc

if anyone is interested in testing it. (But don't waste too much time
testing right now, because it is likely to change.)

It's not quite ready for inclusion, but, eventually, maybe Bill Hart's
prediction about it only taking 10 seconds to compute p(10^9) on
sage.math will come true.

[Alec Mihailovs]

From: "Jonathan Bober" <jwb...@gmail.com>

> It's not quite ready for inclusion, but, eventually, maybe Bill Hart's
> prediction about it only taking 10 seconds to compute p(10^9) on
> sage.math will come true.

That's fantastic! Great work!

Alec