Implement the partition function

Question

I am trying to write my own version of the PartitionsP function. Here is my code:

P[1, 1] = 1;
P[n_, k_] := With[{x = n - k},
   If[x == 0, 1,
    If[x <= k, Sum[P[x, i], {i, 1, x}], Sum[P[x, i], {i, 1, k}]]]
   ];
P[n_] := Sum[P[n, i], {i, 1, n}];

I am using the fact that the number of partitions of n into k parts is equal to the number of partitions whose largest part is k. This code gives me the right answers. However, it is very inefficient and takes a long time for anything larger than $n = 50$. How can I improve the efficiency of the code? Thanks in advance.

Answer 1

For example, you can use memoization

P[n_, k_] := 
  P[n, k] = 
   With[{x = n - k}, 
    If[x == 0, 1, If[x <= k, Sum[P[x, i], {i, 1, x}], Sum[P[x, i], {i, 1, k}]]]
   ];

And as a second step, you can get rid of the Ifs:

ClearAll[Q];
Q[n_, n_] = 1;
Q[n_, k_] := Q[n, k] = With[{x = n - k}, Sum[Q[x, i], {i, 1, Min[x, k]}]];
Q[n_] := Sum[Q[n, i], {i, 1, n}];

Comparing the freshly defined Q to the original P, I get for example

a = P[50]; // AbsoluteTiming // First
b = Q[50]; // AbsoluteTiming // First
a == b

22.9355

0.00772

True

Answer 2

Rademacher like formula:

 ded[h_, k_] := If[k < 3, 0, Sum[m*(FractionalPart[h*m/k] - 1/2), {m, 1, k - 1}]]; 
 omega[k_, n_] := If[k == 1, 1, Sum[If[GCD[h, k] == 1, Cos[(ded[h, k] - 2*h*n)*Pi/k], 0], {h, 1, k - 1}]];
 part[n_] := Round[1/(Pi*2*Sqrt[2])*Sum[omega[k, n]*Sqrt[k] * D[Exp[Pi/k*Sqrt[2/3*(nn - 1/24)]]/Sqrt[nn - 1/24], nn] /. nn -> n, {k, 1, Max[5, Sqrt[n]/2]}]];

 part[1000] // AbsoluteTiming
 (* {0.0156423, 24061467864032622473692149727991} *)

 PartitionsP[1000] // AbsoluteTiming
 (* {0.0298074, 24061467864032622473692149727991} *)

Observation:

c*Sqrt[n] terms suffice to compute part(n) exactly by forcing the error < 1/2 and rounding to the nearest integer. Interesting is a constant c. With c = 0.5 is the difference between rounded and real value within the required bounds (but for example with c = 0.3 not!).

 c = 0.5;
 partfloat[n_] := 1/(Pi*2*Sqrt[2]) * Sum[omega[k, n]*Sqrt[k] * 
 D[Exp[Pi/k*Sqrt[2/3*(nn - 1/24)]]/Sqrt[nn - 1/24], nn] /. nn -> n, {k, 1, 
 Max[5, c*Sqrt[n]]}];
 dif = Table[N[partfloat[n] - PartitionsP[n], 100], {n, 1, 1000}];
 Min[dif]
 Max[dif]
 (* -0.1924581822178408 *)
 (* 0.19346720647243795 *)

 ListPlot[dif]

 c = 0.3;
 dif = Table[N[partfloat[n] - PartitionsP[n], 100], {n, 1, 1000}];
 ListPlot[dif]

Answer 3

As Henrik said, use memorization (Functions That Remember Values They Have Found)

Clear[P]

P[1, 1] = 1;
P[n_, k_] := 
  P[n, k] = With[{x = n - k}, 
    If[x == 0, 1, 
     If[x <= k, Sum[P[x, i], {i, 1, x}], Sum[P[x, i], {i, 1, k}]]]];
P[n_] := Sum[P[n, i], {i, 1, n}];

However, then generate a sequence

seq = P /@ Range[10];

and find the corresponding sequence function

FindSequenceFunction[seq, n]

(* PartitionsP[n] *)

EDIT: Checking equivalence over a wider range

(P /@ Range[200]) == PartitionsP@Range[200]

(* True *)

where I have made use of the fact that PartitionsP is Listable

Attributes[PartitionsP]

(* {Listable, Protected} *)