Compute the number of ways how you can express a number as a sum [closed]

For a given natural number n, compute the number of ways how one can express n as a sum of positive natural numbers. For example, for n = 1 there is just one way. For n = 4 we have 5 possible ways:

1 + 1 + 1 + 1
1 + 1 + 2
1 + 3
2 + 2
4


The ordering of terms in a sum doesn't matter, so 1 + 3 and 3 + 1 are considered the same. Your program should be effective enough to return the result for n = 1000 in order of minutes.

The winner will the solution with most up-votes in (at least) 2 weeks.

• For people who remember having already done it on Project Euler: Problem 76 May 13, 2013 at 7:49
• Why is this a popcon rather than code-golf? Jun 27, 2016 at 13:13
• @ppperry Honestly I don't remember exactly. In general I favor nice code over golfed code, so this is probably the reason. I guess [fastest-code] or [fastest-algorithm] might also be good for the question.
– Petr
Jun 27, 2016 at 14:02

An efficient Haskell version using a generating function with Memoization. Solve instantaneously the problem for n = 1000. No precision issues even with large numbers:

time ./partitions 10000
36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144
./partitions 10000  17.02s user 0.00s system 100% cpu 17.021 total


The code

generalizedPentagonal = map pentagonal $[1..] >>= \i->[i,-i] where pentagonal n = (3*n^2-n)quot2 partitions = (map partitions' [0..] !!) where partitions' 0 = 1 partitions' n = sum . alternate . map partitions . takeWhile (>=0) . map (n-)$ generalizedPentagonal
alternate = zipWith (*) (concat . repeat $[1,1,-1,-1])  Code-golf version p=(map a[0..]!!)where a 0=1;a n=sum.zipWith(*)(concat.repeat$[1,1,-1,-1]).map p.takeWhile(>=0).map(n-).map(\n->(3*n^2-n)quot2)$[1..]>>= \i->[i,-i]  Efficiency improvements Replace the memoization with an Array instead of a list: import Data.Array [...] bnd = (0,10000) partitions = (array bnd [(i, partitions' i) | i <- range bnd] !) [...]  Time for partitions 10000 drop from 17s to 960ms Results Interpreted: partitions 1000 24061467864032622473692149727991 (0.11 secs, 21833064 bytes)  Compiled: time ./partitions 1000 24061467864032622473692149727991 ./partitions 1000 0.08s user 0.00s system 99% cpu 0.084 total  Useful links • Nice. You can shorten a few things, for instance you don't need where. concat[[i,-i]|i<-[1..]] is more concisely [1..]>>= \i->[i,-i] (where you might as well directly build in the pentagonal call). May 13, 2013 at 23:10 Mathematica 25 11 PartitionsP  Usage PartitionsP[4] (* 5 *)  Calculating the result for n = 1000 is instantaneous Timing[PartitionsP@1000] (* {0., 24061467864032622473692149727991} *)  • I was going to propose this solution but thought it was a bit too ready-made. May 12, 2013 at 17:59 • @DavidCarraher All APL solutions seems like that to me :D May 12, 2013 at 18:00 • I know just what you mean. May 12, 2013 at 18:01 JavaScript var partition = function(n) { var sigma=[], p=[1], i,j; for (i=1; i<=n; i++) for (j=i; j<=n; j+=i) sigma[j] = (sigma[j]||0) + i; for (i=1; i<=n; i++) { p[i] = 0; for (j=1; j<=i; j++) p[i] += p[i-j] * sigma[j]; p[i] /= i; } return p[n]; }  This computes the partition function A000041 as the Euler transform of the sequence A000012 (1,1,1,...). The intermediate sequence sigma is A000203, the sum-of-divisors function. For n=1000 it takes approx 0.05 seconds in Node on my 2GHz machine, although since it's working with JS numbers it runs into precision issues. A slightly more complex version (online demo, it's a bit long to paste) adds a naïve big integer implementation and takes approx 0.5s in Node. JavaScript, 105 103 for(n=prompt(a=[b=1]);b<=n;a[b++]=s)for(s=0,k=1,j=b-1;j>-1;j-=k&1?k:k/2)s+=++k&2?a[j]:-a[j];alert(a[n])  Calculation takes less than 250ms for n=1000. The result isn't precise, because JS doesn't have arbitrary precision numbers, but the order of the result is correct. GolfScript, table recurrence [[1]]\({..[[{1$,1$,-=}%0@0=+]zip{{+}*}:^%]\+}*0=^  Assumes input as an integer on the stack: leaves the result on the stack. Online demo counting partitions of 10. This isn't as efficient as the other answers (it computes the partitions of 1000 in approximately 30 seconds), but the mathematical background is much simpler. Instead of considering just the number of partitions of n, consider the number of partitions of n into k parts: p(n, k). Each such partition either contains a 1 or doesn't. The number of partitions which contain a 1 is easily seen to be the same as the number of partitions of n-1 into k-1 parts. For the partitions which don't contain a 1, we can subtract 1 from each part without getting any zeroes, so their number is equal to the number of partitions of n-k into k parts. So we have the recurrence p(n, k) = p(n-1, k-1) + p(n-k, k). This program builds the table of p(n, k) starting from p(1, 1) = 1, and then at the end adds up the partitions of n into 1..n parts to get the result. Dissection: # Initialise the table p(x,k) with the row for x=1 # Note that we will prepend new rows to the table, so it is built backwards [[1]] # Repeat n-1 times \({ # We need a couple of copies for the two parts of the sum in the recurrence .. # Collect in an array... [ # Collect in an array... [ # First an array of p(x-k, k) for k = 1 to x-1 # x is the length of the current table plus one # We iterate over the rows of the table { # Stack: table table row # Row x-k is of length x-k, so k = x - length(row) # But because the row is indexed from 0, we want to find k - 1 # = length(table) - length(row) 1$,1$,- = }% # Stack: table table [p(x-1,1) ... p(1,x-1)] # (Although actually it will be shorter, only until about p(x/2, x/2)) 0@0=+ ] # Stack: table [[p(x-1,1) ... p(1,x-1)] [0 p(x-1,1) ... p(x-1,x-1)]] # Vector sum. Note that we store the array-sum operation in ^ for later use. zip{{+}*}:^% ] # Prepend to the table \+ }* # Take the first row of the table and sum it 0=^  (NB If this were then I would save 1 char by repeating the loop n times and taking the second row of the table at the end). Python2, 110 (Bruteforce O(n^n) implementation) A bruteforce implementation in python which takes the input of stdin: from itertools import* n=input() print len(set(frozenset(p)for p in product(range(n+1),repeat=n)if sum(p)==n))  Computing the result for n=8 takes 10 seconds on my machine. Anything over that is way over a few minutes. GolfScript 2,\.,{[2+.2&(0@,{.2%!)/+}/]}%{0\{~.4$,<*~)3$=*+}/+}+*-1=  Just because we need a GolfScript entry. Thanks to Peter for several ideas. Example: 1000 2,\.,{[2+.2&(0@,{.2%!)/+}/]}%{0\{~.4$,<*~)3\$=*+}/+}+*-1=
# => 24061467864032622473692149727991

• Why ..+ near the start? Doesn't .) give you enough pentagonal numbers? May 13, 2013 at 14:52
• @PeterTaylor It does ;-) May 13, 2013 at 14:56
• Calculating the pentagonal numbers as partial sums of A026741 gives a considerable simplification: replace ),...%1> with ,{[2+.2&(0@,{.2%!)/+}/]}% May 13, 2013 at 15:07
• @PeterTaylor Pushing it inside the loop shortens even further ... but unfortunately makes it slower also. May 13, 2013 at 15:31
• I was thinking along similar lines to your change to the second half, but using negative indexing into arrays. Your reversal is nicer, though, because it lets you extract the first rather than the last element for the final result. May 13, 2013 at 15:37

This depends on a BigNumber implementation which supports the following functions: plus, minus, times, div, sqrt, isZero, lt, gt, round). It also uses config to control the precision used in divisions and to set the rounding mode.

var partition = function(n) {
// Hardy-Ramanujan estimate to set the precision with appropriate margin
BigNumber.config((5 + 1.115 * Math.sqrt(n) + Math.log(n) / Math.log(100))|0, 6);

var zero = new BigNumber(0),
one = new BigNumber(1),
// \sum_{i=0}^\infty 2^{i+1} i!^2 / {2i+1}!
PI = function(){ var t = new BigNumber(2), s=t, i=1; while (!t.isZero()) s = s.plus(t = t.times(i).div((i+ ++i))); return s }(),
A = new BigNumber(n).minus(one.div(24)),
B = A.times(2).div(3).sqrt().times(PI),
ABsqrt12 = A.times(B).times(new BigNumber(12).sqrt()),
gcd = function(x, y) { return y ? gcd(y, x % y) : x },
genexp = function(x, a, b, c, d) {
// \sum_{i=0}^\infty ax^{4i}/(4i)! + bx^{4i+1}/(4i+1)! + cx^{4i+2}/(4i+2)! + dx^{4i+3}/(4i+3)!
var res = zero, z = one, i;
for (i=0; !z.isZero(); i+=4) {
res = res.plus(z.times(a)); z = z.times(x).div(i+1);
res = res.plus(z.times(b)); z = z.times(x).div(i+2);
res = res.plus(z.times(c)); z = z.times(x).div(i+3);
res = res.plus(z.times(d)); z = z.times(x).div(i+4);
}
return res;
},
p = zero, q = 0, C, L, Psi_, h, k, s, max_L = zero;
for (;;) {
// Adaptive precision calculation for performance
if (++q > 1) BigNumber.config((5 + 1.115 * Math.sqrt(n) / q + Math.log(n) / Math.log(100))|0);

L = zero;
for (h = 0; h < q; h++) {
if (gcd(h,q) > 1) continue;
for (k=s=0; k < q; k++) s+= (2*(h*k %q) - q) * k - 4*h*n;
// NB The %(4*q*q) is critical for performance
L = L.plus(genexp(new BigNumber(s % (4*q*q)).div(2*q*q).times(PI),1,0,-1,0));
}
if (L.gt(max_L)) max_L = L;

C = B.div(q);
Psi_ = genexp(C,C,-1,C,-1).times(new BigNumber(q).sqrt());
p = p.plus(L.times(Psi_));
if (Psi_.times(max_L).abs().lt(ABsqrt12)) break;
}
return p.div(ABsqrt12).round();
};


It computes partition(10000) in about 0.67 seconds using Node on a 3.5GHz PC.