2つの素数の和で表される数（ゴールドバッハの予想を仮定して）

【例題 １】 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]

これをよく見ると以下のことが分かります。

1. $4$以上の偶数はすべて含まれている
2. 奇数はもとの素数に$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）