- 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。
(追記) sagemathのプログラムを追加しました。
問題 77. 整数を素数に分割
原文 Problem 77: Prime summations
問題の要約：自然数$n$を素数の和で表す組み合わせ初めて5000を超える$n$求めよ
今回は素数の和と言うことでとりあえず再帰関数で全探索しました。値があまり大きくないので何とか1分ルールはクリア。
from sympy import primerange
from copy import copy
def cutNumPrimes(L, n):
count = 0
for pr in Primes:
if len(L) > 0 and pr < L[-1]: continue # List should ascending order
L1 = copy(L)
L1.append(pr)
rem = n - pr # remaining numbers
if rem > 0:
count += cutNumPrimes(L1,rem) # search more
elif rem == 0:
count += 1
return count
Nmax = 100
Primes = list(primerange(1,Nmax+1))
for n in range(2,Nmax+1):
countSum = cutNumPrimes([],n)
if countSum > 5000:
break
print(f"Answer : {n}, {countSum}")
その後いろいろ探して見たところこの数列はオンライン整数列大辞典のA000607に以下のようなプログラムを発見しました。かなり速いです。
from sympy import primefactors
l = [1, 0]
for n in range(2, 101):
l.append(sum(sum(primefactors(k)) * l[n - k] for k in range(1, n + 1)) // n)
print(l)
これはどういうアルゴリズムなのか調べた結果なんとか見つけたのが、「Prime Partition(Wolfram MathWorld)とそこから参照されている「Euler Transform(Wolfram MathWorld)」の以下の式と説明。これを実装したのが上のプログラムかなと思いますが、下の式の意味は私の知識レベルでは理解できませんでした。誰か分かる方がいたらコメント頂ければ幸いです。
b_n=\frac{1}{n}[c_n+\sum_{k=1}^{n-1}c_k b_{n-k}]
if $a_n$=1 for n prime and $a_n$=0 for n composite, then $b_n$ is the number of partitions of n into prime parts (Sloane and Plouffe 1995, p. 21).
（開発環境：Google Colab）
sagemathのPartitionsを使う
sagemathのPartitionsにはparts_inというパラメータがあり使うパーツを指定できます。ここに素数のリストをprime_rangeで指定すると簡単に求まります。
n = 2
while Partitions(n, parts_in=prime_range(n + 1)).cardinality() < 5000:
n += 1
print(n)
（開発環境：CoCalc）