🐺

プログラムを書き始める前に、数式でどのように表されるのかを考える〜ゴールドバッハ予想を例に〜

に公開
2025/12/16

アドベントカレンダー(自称)vol.13

※2025年12月16日修正

皆さん、こんにちは。yumyum116 です。SWE転職を目指す者です。
アドベントカレンダーに便乗して、1日1記事投稿に挑戦するとともに、執筆活動の習慣化にも挑戦してみます。

この記事は、線形探索で総当りに解答を調べるアプローチしか思いつけない自分が綴る、効率的に解を求める方法を思いつくアプローチを身につけるシリーズです。
自分と同じようなところで躓いている方、プログラミング初学者の方の参考になれば嬉しいです。

今回のテーマはゴールドバッハ予想です。

ゴールドバッハの予想を用いて、問題を解いてみる

ゴールドバッハ予想とは、次の命題です。

2より大きな偶数は、2つの素数の和で表せる。

"命題" とはいったものの、実はこの命題は未解決問題の一つです。
1742年に提唱され、今なお証明されていないものの、計算機による検証では、4 * 10^{18} まで成立することが、2012年に Oliveira e Silva によって確認されています(参考)。

この計算機による検証結果に基づき、ゴールドバッハ予想を用いて問題を解いてみる、というのが本記事の主旨です。

さて、次の問題を解くことを考えてみます。

4以上の偶数 N について、2つの素数の和で表し、標準出力に出力してください。
ただし、答えが複数個ある場合には、素数の組合せのうち、その積が最も大きくなるような組合せを出力すること。

まずは、私の思考アプローチをご紹介します。

  1. 素数リストを作成する
  2. 1 で作成した素数リストを用いて、2つの素数の和が N となる素数の組を線形探索により探索する
  3. 2 で見つけた素数の組合せを、配列に格納する。
  4. 3 で作成した配列の中から、積が最大となる組合せのインデックスを調べる
  5. 4 で調べたインデックスの位置にある配列の要素を出力する

上記のアプローチの場合、時間計算量は O(n^2) になります。

プログラムで表すと、次のように表すことができます。

goldbach_original.py
def goldbach_conjecture(n):
    prime = [2]
    for i in range(3, n + 1, 2):
        for j in range(3, i):
            if i % j == 0:
                break
        else:
            prime.append(i)

    prime_combi = []
    for i in range(len(prime)):
        if n - prime[i] in prime:
            prime_combi.append([prime[i], n - prime[i]])
    if not prime_combi:
            return None

    max_res = -1
    max_idx = 0
    for j in range(len(prime_combi)):
        prod = prime_combi[j][0] * prime_combi[j][1]
        if prod > max_res:      // 判定式で値の更新有無を確認しておらず、積が最大となるか判定できないため、修正しました。
            max_res = prod
            max_idx = j

    return prime_combi[max_idx]

n = int(input())
result = goldbach_conjecture(n)
if result is None:
	print("No combination found")
else:
	for prime in result:
		print(prime)

上記のプログラムでも動くものの、N の値が大きくなるについて、TLEとなる場合があります(実際に、テストケースによっては TLE になりました)。
本プログラムのアプローチを再度見直すと、とくに次のステップの計算量を削減できそうです。

  1. 2 で見つけた素数の組合せを、配列に格納する。
  2. 3 で作成した配列の中から、積が最大となる組合せのインデックスを調べる

さて、そもそも求めたい値がどのような値であったかを振り返ります。
そう、2つの素数について、その和が N となるような組合せの中で、その積が最大となる素数の組合せを求める でした。

言わずもがな、上記の課題要件に当てはまらない素数の組合せについては、計算する必要も、その組合せを配列に保存しておく必要もありません。

それでは、どのようにして「積が最大となる素数の組合せ」を求めればよいのでしょうか。

ここで、数式を用いて考えてみます。

2つの素数 p, q を用いて、

N = p + q

と表せるものとする。

この時、2つの素数の積は次のように表せる。

p * q = p * (N - p) = -p^2 + Nq

これは、下に凸の二次関数であるため、最大値は

p = \frac{N}{2}

である。

したがって、\frac{N}{2} に最も近い素数同士が積を最大化する。

この考えに基づくと、\frac{N}{2} から外側に向かって探索するアプローチが一番計算量を抑えられそうである、という算段がつきます。

さて、見直すべきステップの見通しが立ったところで、アプローチを次のように修正してみます。

  1. 素数リストを作成する
  2. 1 で作成した素数リストを用いて、2つの素数の和が N となる素数の組を線形探索により探索する。なお、探索は \frac{N}{2} から始めるものとする。
  3. 2 の探索により見つかった素数の組を戻り値として返す

ステップもスッキリしましたね。
なお、時間計算量は、O(N * log(log(N))) + O(N) = O(N * log(log(N))) となります。

プログラムは、次のように記述することができます。

goldbach.py
def goldbach_conjecture(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False

    for i in range(2, int(n ** 0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

    for p in range(n // 2, 1, -1):
        q = n - p
        if is_prime[p] and is_prime[q]:
            return p, q

n = int(input())
p, q = goldbach_conjecture(n)
print(p)
print(q)

プログラムもかなりスッキリしました。


プログラムの高速化は、計算量が最小となるようなアルゴリズムの知識やI/Oの最適化だけではなく、求めるべき解答を、数式によって表す数学的なアプローチもあることを学びました。

とくに、今回の問題は特別な名前がついた性質や定理等を使うわけではないため、数学を専門に学んでいない方であっても、十分に対応できると思います。

今後、問題を解くアプローチを検討する上で、数式を用いた解法の最適化も手段の一つとして頭に入れておくと、発想の枠が一気に拡がることと思います。

記事内に誤謬等ある場合には、修正いたします。その際は、ご連絡いただけますと幸いです。

それでは、また。

GitHubで編集を提案

Discussion