Problem 78 of Project Euler has been in my scope for a long time, since it is the first exercise on the list which is solved by less than 5000 people at the time of writing this post. The problem reads

Let p(n) represent the number of different ways in whichncoins can be separated into piles. For example, five coins can separated into piles in exactly seven different ways, so p(5)=7.

OOOOOOOOO OOOO OOOOO O OOO OO OOO O O OO O O O O

Find the least value ofnfor which p(n) is divisible by one million.

To me this was a rather easy problem to solve, and now you might ask why since it doesn’t seem straight forward. The answer to that lies in the fact that I was reading the problem description a good while ago since as mentioned earlier it is the first problem which has been solved by less than 5000 people. A day or so later I listened to the second episode of Strongly Connect Components where Bruce Reznick mentioned that his prime research is on integer partitions and he mentions the exact example given in the problem description.

So before starting on the problem I had a very good idea where I should look for a possible solution. A quick search on Google turned up the Wikipedia page on integer partitions where the section on generating functions seems to provide the methodology for solving this problem.

## Generating functions

The generating function for p(n) is given as

Which according to wikipedia can be rewritten as

p(n) = p(n – 1) + p(k – 2) – p(k – 5) – p(k – 7) + p(k – 12) + p(k – 15) – p(k – 22)…

where p(0) = 1 and p(n) = 0 for n < 0.

The sequence to use are the generalized pentagonal numbers which are given on the form f(k) = k(3k-1)/2 for both negative and positive k. This can be written as a positive sequence of integers m=1,2,3,…. such that

where the division is an integer division.

The signs of the function follows the pattern +,+,-,-,+,+,-,-…. so that is pretty easy to replicate.

So if we know all previous integer divisions we can also calculate the next one.

I don’t say I completely understand the underlying theory of the integer partition formula, so I wont go into details about it. But you can start at wikipedia for info.

## Avoiding BigInteger

My first implementation of the code for the generating function included the use of the BigInteger class. However that becomes cumbersome. So we want to eliminate that. However, since we are not interested in finding the actual number of partitions but just the first one divisible by 1.000.000 and since we are summing up things we only need to save the last 7 digits. So we just take modulo of 1.000.000 before we store the result.

## C# implementation of the code

Since we now have all we need in order to implement the code, I don’t think I need to blabber on. Here is the code

List<int> p = new List<int>(); p.Add(1); int n = 1; while(true){ int i = 0; int penta = 1; p.Add(0); while (penta <= n){ int sign = (i % 4 > 1) ? -1 : 1; p[n] += sign * p[n - penta]; p[n] %= 1000000; i++; int j = (i % 2 == 0) ? i / 2 + 1 : -(i / 2 + 1); penta = j * (3 * j - 1) / 2; } if (p[n] == 0) break; n++; }

## Wrapping up

Now that we have the code lets execute it and get the solution

n = 55374 is the least value for which p(n) is divisible by 1.000.000 Solution took 299,6262 ms

I can execute the code in just under 300ms, which I guess is pretty okay for this problem. However if you have more efficient methods I would be most interested in hearing about it.

As usual you can find the source code here.

The image used for the post header is shared under the creative commons license by Tax Brackets.

Thank you very much, this problem i wasn’t able to solve for a long time. Your description helped me a lot. And than a I could code it myself (with some problems).

How efficient is it to implement the following recurrence relation?

P(n,k)=P(n-1,k-1)+P(n-k,k)

P(n,k) means number of partitions of size k

We can sum P(n,k) from k=1 to k=n

P(n,k)=0 for k>n, P(n,n)=1, and P(n,0)=0

Source: http://mathworld.wolfram.com/PartitionFunctionP.html

My C++ implementation ran in just 30 ms. The difference is that my inner loop does much less work than yours. It’s faster to generate the generalized pentagon numbers ahead of time — there are only ~1600 less than a million.

I generated the generalized pentagon numbers as follows:

int32_t iN, iNext1 = 0, iNext2 = 0, iOff1 = 1, iOff2 = 2;

do

{

iNext1 += iOff1;

iOff1 += 3;

Pents.push_back(iNext1);

iNext2 += iOff2;

iOff2 += 3;

Pents.push_back(iNext2);

}while (Pents.back() < iMod);

And this was my main loop. Note that I used bitwise and to distinguish the sign, and only did one modulo operation at the end (this cut my runtime by a good bit).

while (true)

{

iSum = 0;

for (iN = 0; Pents[iN] <= Parts.size(); ++iN)

{

if (iN & 2u)

iSum -= Parts[ Parts.size() – Pents[iN] ];

else

iSum += Parts[ Parts.size() – Pents[iN] ];

}

iSum %= iMod;

if (iSum == 0)

break;

Parts.push_back(iSum);

}

[…] air’s memory and marches slowly around p(37000). I looked up the solution online and found this blog that explains everything, apparently there is a well developed math theory called integer […]

There is a good Paper about p(n) with different solotions + source code.

A solution for this problem:

[latex]private static long p(long n)

{

// trivial special cases

if(n < 0) return 0;

if(n == 0) return 1;

// important things

// K in 1/2*3k^2+k

//

// we need at least one extra K to stop loops.

long maxK = (long) ((Math.sqrt(1+24*(double)n)-1)/6) + 1;

long pentNums[] = new long[2*(int)maxK];

long parts[] = new long[(int)n+1]; // From 0 to n means n+1 elements.

long partialSum = 0;

int i, j;

int sign = 1;

parts[0] = 1;

// generate all the pentagonal numbers we need.

for(i = 1; i <= maxK; i++)

{

pentNums[2*i -2] = (3*i*i-i)/2; // the -2 in the array index is to zero-base

pentNums[2*i+1 -2] = (3*i*i+i)/2;

}

// calculate p(i) for i <= n.

for(i = 1; i <= n; i++)

{

partialSum = 0;

// sum p(n-1) + p(n-2) – p(n-5) – …

for(j = 1; pentNums[j-1] <= i; j++)

{

// determine j´s residue class mod 4. Signs come in pairs: +, +, -, -, ..

sign = (int)(((float)j/4 – (float)(j/4))*4);

switch(sign)

{

case 1:

partialSum += parts[i-(int)pentNums[j-1]];

partialSum %= 1000000;

break;

case 2:

partialSum += parts[i-(int)pentNums[j-1]];

partialSum %= 1000000;

break;

case 3:

partialSum -= parts[i-(int)pentNums[j-1]];

partialSum %= 1000000;

break;

case 0:

partialSum -= parts[i-(int)pentNums[j-1]];

partialSum %= 1000000;

break;

}

}

parts[i] = partialSum;

if(partialSum == 0)

{

return i;

}

}

return partialSum;

}[/latex]

Time: 93 ms

The least value of n for which p(n) is divisible by one million.

36325300925435785930832331577396761646715836173633893227071086460709268608053489541731404543537668438991170680745272159154493740615385823202158167635276250554555342115855424598920159035413044811245082197335097953570911884252410730174907784762924663654000000

http://blog.dreamshire.com/project-euler-78-solution/

p(n) = p(n – 1)+ p(k – 2) – p(k – 5) – p(k – 7) + p(k – 12) + p(k – 15) – p(k – 22)…needs to be changed to

p(k) = p(k − 1)+ p(k − 2) − p(k − 5) − p(k − 7) + p(k − 12) + p(k − 15) − p(k − 22) − …I think in there is a mistyp in formula for generalized pentagonal numbers: not “if k mod 2 = 0”, but “if m mod 2 = 0” instead.

when i ran the code, for example, p(100) is a negative number, should it be a positive one ? can u explain this ? tks for your great work!!

Hey!

You might be interested in my solution.

It’s well explained: https://godovich.com/2017/08/18/coin-partitions-pe78/

for calculating the next pentagonal number instead of those two lines only one conditional

statement will do the job;

pent += ((i%2)?i/2:i)+1;