Project Euler Problem 78 (C++和Python代码实现和解析)

Problem 78 : Coin partitions

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

OOOOO
OOOO O
OOO OO
OOO O O
OO OO O
OO O O O
O O O O O
Find the least value of n for which p(n) is divisible by one million.

1. 欧拉项目第78道题 硬币划分(堆)

让p(n)表示n个硬币可以被分成成堆的不同方式的数目。例如，五个硬币可以用七种不同的方式分开成堆，所以p(5)=7。

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

求n的最小值，其中p(n)可被100万整除。

2. 求解分析

这个p(n)果然来头不小，用普通的生成函数也很费时！找到了这个，Computing the Partitions of n，简单朴实的网页里给出了这样的一个公式。

其中k从1开始，直到迭代到无法迭代，即n<0。经过观察，k(3k+1)/2是一个等差数列和的形式，其通项为3k+2，同理后面的那一个通项为3k+1。

我们就不必用普通的生成函数来求解这道题了，直接使用公式来求解。

3. C++ 代码实现

C++采用了和求解分析一样的方法。a_k是一个等差数列的通项 a_k，b_k也是一个等差数列的通项 b_k。为了避免频繁地申请内存，还是使用了数组 P[max_n + 1]。

C++代码

#include <iostream>
#include <cassert>

using namespace std;

class PE0078
{
private:
     static const int max_n = 100000;

     int coinPartitions(int n, int P[]);

public:
     int get_n_ForPnDivbyOneMillion();
};

int PE0078::coinPartitions(int n, int P[]) 
{
    int result = 0, k = 1, coeff = 1; // (-1)^0 = 1
    int a_k = k * (3 * k + 1) / 2, b_k = k * (3 * k - 1) / 2;

    while (n >= a_k)
    {
        result += coeff * (P[n - a_k] + P[n - b_k]);
        a_k    += 3 * k + 2;
        b_k    += 3 * k + 1;

        coeff *= -1;
        k      += 1;
    }

    // # (n - a_k) < 0, check whether (n - b_k) > 0
    result += (n >= b_k) ? coeff * P[n - b_k] : 0;

    return result % 1000000;  // % one million
}
    
int PE0078::get_n_ForPnDivbyOneMillion()
{
     int P[max_n + 1];

     P[0] = 1;
     P[1] = 1;

     int n = 2;

     while (true)
     {
          P[n] = coinPartitions(n, P);

          if (P[n] == 0)
          {
               break;
          }
          else
          {
               n++;
          }
     }

     assert(P[5] == 7);

     return n;
}

int main()
{
     PE0078 pe0078;

     int n = pe0078.get_n_ForPnDivbyOneMillion();

    cout << "Found the least value of " << n << " for which p(";
    cout << n << ") is divisible by one million." << endl;

    return 0;
}

4. Python 代码实现

Python采用了和求解分析一样的方法，列表Pn_list用来保存Pn序列的每一项，但列表ab_list的每一项是元组(a_k, b_k)，而且我们提前把列表ab_list的前200项的值计算出来了。

Python 代码

import time
    
class CoinPartitions(object):
    def __init__(self):
        """
        Pn_list is used to save values of p(n)
        ab_list is used to save (k*(3*k+1)//2, k*(3*k-1)//2)  
        """
        self.Pn_list = [1, 1]    # Pn_list[0] = Pn_list[1] = 1
        self.ab_list = [ (k*(3*k+1)//2, k*(3*k-1)//2) \
            for k in range(1, 200) ]

    def performCoinPartitions(self, n):
        """
        implement p(n), but return (Pn % 1000000)
        """
        Pn, k, coeff = 0, 1, 1  # (-1)^0 = 1
        (a_k, b_k)   = self.ab_list[0]

        while n >= a_k:
            Pn += coeff * (self.Pn_list[n - a_k] + self.Pn_list[n - b_k])   
            coeff, k   = -coeff, k+1   # coeff *= -1, k += 1
            (a_k, b_k) = self.ab_list[k-1]

        # (n - a_k) < 0, check whether (n - b_k) > 0
        Pn += coeff * self.Pn_list[n - b_k] if n >= b_k else 0

        return Pn % 1000000  # Pn is divisible by one million

    def get_n_ForPnDivbyOneMillion(self):
        """
        find the least value n for which p(n) is divisible by one million
        """
        n = 2
        while True:
            Pn = self.performCoinPartitions(n)
            if 0 == Pn:
                break
            else:
                self.Pn_list += [ Pn ]
                n += 1
        return n
        
def main():
    start = time.process_time()

    cp = CoinPartitions()

    n = cp.get_n_ForPnDivbyOneMillion()

    print("Found the least value of %d for which p(%d" %(n, n), end='')
    print(") is divisible by one million.")

    end = time.process_time()
    print('CPU processing time : %.5f' %(end-start))

if  __name__ == '__main__':
    main()

P.S 本文参考了 https://blog.youkuaiyun.com/u013975830/article/details/55105696 的方法。