【Project Euler】Problem 78: コインの分割

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題 78.コインの分割

原文 Problem 78: Coin partitions

問題の要約：コイン$n$個を分割する方法の数$p(n)$が$10^6$の倍数になる最小の$n$を求めよ

分割問題3部作の最後です。「Problem 76: 整数の分割」と同じなのですが数が全然大きいのでまた違ったアプローチが必要になります。とても100番以下の問題とは思えませんね。

とりあえずProblem 76で紹介したnpartitionで単純に探して見ます。5分以上かかりますが答えは見つかります。

from sympy.ntheory import npartitions
from itertools import count
for n in count():
  p = npartitions(n)
  if p % 10**6 == 0:
    break
print(f"Answer: {n}")

これを超える方法なんてあるのと思いますが探すとあるんですね。「分割数(Wikipedia)」の「分割数の母函数」セクションで以下のような漸化式が紹介されてます。

\begin{align}
&p(k) = p(k − 1) + p(k − 2) − p(k - 5) − p(k − 7) + p(k − 12) + p(k − 15) − p(k − 22) − ...  \\
&ここで1,2,5,7,12,15,22 は\frac{n(3n − 1)}{3} の形の一般五角数でnの符号が交互(+/-)\\
&また各項の符号は +, +, −, −, +, +, ... と続く
\end{align}

これをPythonで実装したのがこちらです。

from itertools import count
from functools import lru_cache
@lru_cache(maxsize=None)
def pentnum(m):
  k = [1, -1][m%2] * (m // 2 + 1)
  return (k * (3 * k - 1)) // 2

p = [1]
for n in count(1):
  p.append(0)
  for m in count(0):
    penta = pentnum(m)
    if penta > n: break
    p[n] += [1,1,-1,-1][m%4] * p[n - penta]
    p[n] %= 10**6
  if (p[n] == 0): break
print(f"Answer : {n}")

（開発環境：Google Colab）