【例題 １】 N以下の正の整数で2つの素数の和で表される数の個数を求めよ
まず例として$N=20$のとき素数の数を数えると8個あります。
import sympy as sp
N = 20
primes = list(sp.primerange(N+1))
print(len(primes), primes)
# 8 [2, 3, 5, 7, 11, 13, 17, 19]
この中から重複を許して2個足し合わせてできる数を求めます。
from itertools import combinations_with_replacement as combr
pp = set([sum(ps) for ps in combr(primes,2) if sum(ps) <= N ])
print(f"p2 = {pp}")
#p2 = {4, 5, 6, 7, 8, 9, 10, 12, 13, 14, 15, 16, 18, 19, 20}
これだけだと分かりづらいので奇数と偶数に分けます。
even, odd = set(list(range(2,N+1,2))), set(list(range(1,N+1,2)))
print(f" even : {sorted(list(pp&even))}")
print(f" odd : {sorted(list(pp&odd))}")
# even : [4, 6, 8, 10, 12, 14, 16, 18, 20]
# odd : [5, 7, 9, 13, 15, 19]
これをよく見ると以下のことが分かります。
- $4$以上の偶数はすべて含まれている
- 奇数はもとの素数に$2$を足したものが含まれている
1.はゴールドバッハの予想 (Wikipedia)が正しければいつも真となります。
すべての 2 よりも大きな偶数は2つの素数の和として表すことができる。このとき、2つの素数は同じであってもよい。
2.は奇数になるのは元の素数のうち$2$以外の奇素数と$2$を足した場合だけなので、その数は元の8個から$2$と$2$を足して$20$を超える$19$を除いて$8-2=6$となるわけです。
元の$N$以下の素数の数はすべてリストしなくてもsympyのprimepiを使うと求められるので。例題 １の答えは以下のコードで求められます。
import sympy as sp
np1 = sp.primepi(N) # number of primes [1,n]
np2 = N // 2 - 1 + (np1 - 1) # Even (>=4) + Odd (>=5)
if N > 2 and (sp.isprime(N) or sp.isprime(N-1)): np2 -= 1 # +2 makes out of range
print(f"np2 = {np2}")
# np2 = 15
【例題 ２】 N以下の正の整数で3つの素数の和で表される数の個数を求めよ
上のプログラムのcombrの第二引数を３にすれば良いので。
pp = set([sum(ps) for ps in combr(primes,3) if sum(ps) <= N ])
allN = set(list(range(1,N+1)))
print(f"p3 = {pp}")
print(f" not in p3 : {sorted(list(allN-pp))}")
# p3 = {6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}
# not in p3 : [1, 2, 3, 4, 5]
こうなると含まれない数の方が少ないので見てみると$[1,2,3,4,5]$となります。これは２つの素数の和で「４以上の偶数」はすべて含まれていたので、これに$2$を足すと「６以上の偶数」、$3$を足すと「７以上の奇数」はすべて含まれることからも導かれます。したがって例題 ２の答えは以下の簡単なコードで求められます。
p3 = max(n - 5, 0)
この考え方はProject Euler Problem 543 Prime-Sum Numbersを解くのに役に立ちます。
（開発環境：Google Colab）