· 開発

Project Euler 3: 割れるだけ割って最大の素因数を求める

小さい数から割れるだけ割っていく試し割り法で、12桁の整数の最大の素因数を求めます。最後に残る値がなぜ素数なのかの証明、√n の反復回数の上限、64ビット整数の範囲、C++17とPython 3の実装まで扱います。

Project Eulerアルゴリズム数学素因数分解

ネタバレ注意 · この記事には解き方と答えがそのまま書かれています。先に自分で解いてみたい方は、Project Euler の問題 3からどうぞ。

Project Euler シリーズの3問目です。素因数分解はこのシリーズで何度も出てくる道具で、今回はその最も基本的な形である試し割り法(trial division)をきちんと押さえておきます。コードは十数行ですが、「ループが終わったときに残る値はなぜ素数なのか」「ループは最悪で何回回るのか」をきちんと考えると、見た目以上に学ぶことがあります。

ネタバレ注意: この記事には解法の全体と答えが含まれます。先に自分で解きたい場合は、元の問題から始めてください。このシリーズの記事には広告を載せません。問題文のライセンス(CC BY-NC-SA 4.0)に非営利の条件があるためです。

問題

ある整数を割り切る素数を、その数の素因数といいます。問題に出てくる例では、13195の素因数は5、7、13、29で、そのうち最大のものは29です。同じ意味で 600851475143の最大の素因数を求めるのが問題です。

入力は12桁の整数一つで、求めるのは素数一つです。

前提知識

素数と合成数。 2以上の整数のうち、約数が1と自分自身だけのものが素数、それ以外の2以上の整数が合成数です。1は素数でも合成数でもなく、素因数を一つも持ちません。

素因数分解とユークリッドの補題。 2以上のすべての整数は素数の積で表せ、掛ける順序を除けばその表し方は一通りです(算術の基本定理)。この記事で直接使うのは、その土台にあるユークリッドの補題です。素数 p が積 ab を割り切るなら、p は a か b の少なくとも一方を割り切ります。したがって n = a × b と分けておけば、n の素因数は a の素因数か b の素因数のどちらかです。

合成数は √m 以下の素因数を持ちます。 m が合成数なら m = a × b(1 < a ≤ b < m)と書けます。a ≤ b なので a × a ≤ a × b = m、つまり a ≤ √m です。a の素因数の一つを p とすると p ≤ a ≤ √m で、p は m も割り切ります。対偶をとると、2以上の m が √m 以下のどの素数でも割り切れなければ、m は素数です。これが今回の解法の核心です。

整数の割り算と切り捨て。 正の整数 f、r について r / f を整数の割り算(小数点以下切り捨て、⌊r / f⌋)で計算するとき、f × f ≤ r と f ≤ ⌊r / f⌋ は同じ条件です。f × f ≤ r なら f ≤ r / f で、f は整数なので r / f を切り捨てた値以下でもあります。逆に f ≤ ⌊r / f⌋ ≤ r / f なら、両辺に f を掛けて f × f ≤ r です。この同値関係は、C++ の実装で掛け算のあふれを避けるのに使います。

素朴な方法

定義をそのまま書くなら、2から n までのすべての整数 d について「d が n を割り切るか」「d が素数か」を調べ、両方を満たす d のうち最大のものを選べば済みます。

def is_prime(m: int) -> bool:
    return m >= 2 and all(m % d != 0 for d in range(2, m))

n = 13195
best = 0
for d in range(2, n + 1):
    if n % d == 0 and is_prime(d):
        best = d

小さな例なら問題なく29が出ますが、n が600851475143だと外側のループだけで6000億回近く回ることになります。外側のループだけで O(n) で、約数ごとに素数判定を別に行えばさらにかかります。n の桁が一つ増えるごとに仕事が10倍になるので、n が大きくなるとこの方法は使えません。

上から下へ数えていき、最初に見つかった素数の約数で止める形に変えても、最大の素因数が n よりずっと小さければその間の数をほとんど調べることになるので、本質的には改善しません。

解法の導出

ステップ1: 小さな素因数を割り取ると問題が小さくなる。 13195を小さい数から順に割ってみます。2では割り切れず(奇数)、3でも割り切れません(各桁の和が19)。2で割り切れないので4でも当然割り切れません。5で割ると

13195 = 5 × 2639

です。これで13195の残りの素因数はすべて2639の素因数になります(ユークリッドの補題)。元の数の代わりに、より小さい2639だけを見ればよくなりました。

ステップ2: 同じことを続ける。 2639はもう5では割り切れないので6に進み、6でも割り切れないので7で割ります。

2639 = 7 × 377

377はもう7では割り切れません。8、9、10、11、12はどれも割り切れず、13で

377 = 13 × 29

となります。残った値は29です。

ステップ3: いつ止めてよいかを見る。 次の候補は14ですが、14 × 14 = 196 > 29 です。ここまでで2から13までの数で割れるだけ割ったので、29は13以下のどの素数でも割り切れません。√29 ≈ 5.4 なので、√29 以下の素数はすでにすべて試しています。前提知識の事実から29は素数です。29に対してこれ以上試す必要はありません。

まとめると 13195 = 5 × 7 × 13 × 29 で、割り取った5、7、13より最後に残った29のほうが大きいので、答えは29です。問題の例と一致します。

ステップ4: 合成数の候補は気にしなくてよい。 ステップ1では4や6のような合成数も候補として試しましたが、割り切れませんでした。これは偶然ではありません。候補 f が合成数なら、f の素因数はすべて f より小さく、それらの素数は f の番が来る前に残りの値からすでに割り取られています。ですから残りの値が f で割り切れることはありません。結果として、実際に割り切れる候補は常に素数です。素数だけを試すために素数の一覧を別に作る必要はありません。

最終的なアルゴリズム

  1. n < 2 なら素因数がないので0を返します(後の「境界のケース」)。
  2. remaining = n、largest = 1、factor = 2 から始めます。
  3. factor × factor ≤ remaining の間、remaining が factor で割り切れる限り割り続け、そのたびに largest = factor とします。その後 factor を1増やします。
  4. ループが終わったとき remaining > 1 なら largest = remaining とします。
  5. largest を返します。
  6. 答えを出力する前に小さな入力を確かめます。1 → 0、2 → 2、9 → 3、64 → 2、97 → 97、13195 → 29。一つでも違えば0以外の終了コードで終わります。

C++ の実装では手順3の条件を factor ≤ remaining / factor と書きます。前提知識で見たとおり同じ条件ですが、掛け算がないのであふれません。その理由は後の「整数の範囲」で説明します。

正しさの証明

ループが候補 f で条件を判定する時点(f = 2, 3, …)ごとに、次が成り立ちます(不変条件)。

  • (a) n = (これまでに割り取った素数の積) × remaining で、割り取った素数はすべて f より小さい。
  • (b) remaining は f より小さいどの素数でも割り切れない。
  • (c) 一つでも割り取っていれば largest はその中で最大のもの、そうでなければ1。

最初(f = 2)は何も割り取っておらず、remaining = n、largest = 1 なので成り立ちます。条件が真で f の番を処理するとき、内側のループは remaining が f で割り切れなくなるまで f を割り取ります。f が実際に割り切ったなら、ステップ4で見たとおり f は素数で、割るたびに f を一つ積の側へ移すので (a) が保たれます。内側のループが終わると remaining は f で割り切れないので、f + 1 で (b) が成り立ちます。f はそれまでに割り取ったどの素数よりも大きいので、largest = f とすれば (c) も保たれます。remaining は1回割るごとに少なくとも半分になるので、内側のループは必ず終わります。

外側のループでは f が1ずつ増え remaining は増えないので、いつか f × f > remaining となって終わります。終わったときの f を F とすると、二つの場合があります。

  • remaining = 1: (a) より n は割り取った素数の積です。ユークリッドの補題により n の素因数はそれらの素数のどれかなので、(c) の largest が最大の素因数です。
  • remaining > 1: (b) より remaining は F より小さい素数で割り切れず、F × F > remaining なので √remaining 以下の素数はすべて F より小さいです。前提知識の事実から remaining は素数です。また remaining の最小の素因数(= remaining 自身)は F 以上なので、F より小さい割り取った素数よりも大きいです。n の素因数は割り取った素数か remaining のどちらかなので、最大の素因数は remaining です。

どちらの場合も、返す値は n の最大の素因数です。

計算量

最悪の反復回数。 外側のループの本体は f × f ≤ remaining ≤ n となる f についてしか回らないので、多くても ⌊√n⌋ − 1 回です。内側のループで実際に割る回数は、n を素因数分解したときの素因数の個数(重複を含む)に等しく、素因数はどれも2以上なので log₂ n 以下です。それ以外の演算は、外側のループ1回ごとに余りの判定1回と条件の判定1回です。固定幅の整数ではこれらは定数時間なので、時間計算量は O(√n)、保持する値は3つだけなので追加の空間計算量は O(1) です。

この上限には実際に到達します。n が素数なら remaining は一度も減らず、f は √n 付近まで上がる必要があります。「割るたびに remaining が減るので対数時間」と言うのは誤りです。割る回数は対数で抑えられますが、割り切れない候補を通り過ぎる回数はそうではありません。

入力によって早く終わる場合。 ループは remaining が小さくなるほど早く止まります。素因数を重複を含めて大きさの順に並べ、最大のものを p₁、2番目を p₂ とすると、ループはおおよそ max(p₂, √p₁) まで回ります。p₂ までは割り取る必要があり、その後は remaining に p₁ だけが残るので、f が √p₁ を超えたところで止まるからです。小さな素因数が多い数ほど、√n よりずっと早く終わります。この問題の入力で実際に何回回るかは「答えと検証」で数えます。

整数の範囲。

  • 600851475143は 2^39 ≈ 5.5 × 10^11 より大きく 2^40 ≈ 1.1 × 10^12 より小さい40ビットの数です。32ビット整数(最大約 4.3 × 10^9)には収まらないので、64ビット整数が必要です。
  • C++ は unsigned long long を使います。C++ 標準はこの型が少なくとも64ビットであることを保証し(C 標準の ULLONG_MAX の最小値に従う)、符号なし整数の演算は 2^N を法として定義されているので未定義動作になりません。ただしあふれると黙って誤った値になります。リテラルには ULL 接尾辞を付けて、最初から同じ型にしておきます。
  • ループ内の割り算と余りの演算は値を大きくしないので、あふれることはありません。危ないのは条件 factor * factor <= remaining の掛け算です。remaining が 2^64 に近い素数なら factor は 2^32 まで上がり、そのとき factor × factor = 2^64 は0に回り込んで条件が再び真になります。するとループが止まるべきところで止まりません。そこで C++ の実装は同じ条件を factor <= remaining / factor と書きます。割り算しか使わないのであふれず、前提知識の同値関係により結果も同じです。
  • 条件が真の間、factor は √remaining < 2^32 以下なので、++factor もあふれません。したがって C++ の関数は unsigned long long のすべての値で正確です。ただし 2^64 に近い素数を渡すとループは約 2^32 ≈ 43億回回る必要があるので、正確ではあっても遅くなります。
  • Python の int は精度に上限がなくあふれないので、条件は factor * factor <= remaining のまま書いています。アルゴリズムは同じです。

境界のケース

  • n < 2: 0と1は素因数を持たないので、「最大の素因数」は定義されません。二つの実装はどちらもこの場合に0を返すことにしました。0は素数ではないので、正常な答えと紛れることはありません。Python では負の数も同じ扱いになり、C++ は符号なしの型なので負の値は受け取りません。小さな例の確認に 1 → 0 が入っています。
  • n が素数: ループの中で一度も割れずに終わり、残った n 自身が答えになります(2 → 2、97 → 97)。2は最初から 2 × 2 > 2 なので、ループの本体は一度も回りません。
  • n が素数のべき: 64 = 2^6 は2で6回割ると remaining が1になり、答えは最後に割った2です。最後の remaining > 1 の判定が、1を答えと取り違えるのを防ぎます。
  • n が素数の2乗: 9は f = 3 で 3 × 3 ≤ 9 が真と判定されて初めて割り取れます。条件を誤って < と書くと9を素数とみなして9を返してしまうので、このケースを小さな例の確認に入れています。

C++17 による実装

#include <iostream>

// Largest prime factor of n. Values below 2 have no prime factor; return 0 for them.
// Divide out every factor from 2 upward while factor * factor <= remaining.
// The loop tests factor <= remaining / factor, the same condition without a
// multiplication that could overflow. Whatever is left above 1 at the end is
// prime and larger than every factor divided out before it.
unsigned long long largest_prime_factor(unsigned long long n) {
    if (n < 2) {
        return 0;
    }
    unsigned long long remaining = n;
    unsigned long long largest = 1;
    for (unsigned long long factor = 2; factor <= remaining / factor; ++factor) {
        while (remaining % factor == 0) {
            remaining /= factor;
            largest = factor;
        }
    }
    if (remaining > 1) {
        largest = remaining;
    }
    return largest;
}

int main() {
    if (largest_prime_factor(1) != 0 || largest_prime_factor(2) != 2 ||
        largest_prime_factor(9) != 3 || largest_prime_factor(64) != 2 ||
        largest_prime_factor(97) != 97 || largest_prime_factor(13195) != 29) {
        return 1;
    }
    std::cout << largest_prime_factor(600851475143ULL) << '\n';
    return 0;
}

小さな例の確認には assert ではなく if を使いました。assert は NDEBUG を定義したビルドで消えるからです。確認が一つでも外れれば return 1 で0以外の終了コードを返します。

Python 3 による実装

def largest_prime_factor(n: int) -> int:
    """Largest prime factor of n. Values below 2 have no prime factor; return 0 for them.

    Divide out every factor from 2 upward while factor * factor <= remaining.
    Whatever is left above 1 at the end is prime and larger than every factor
    divided out before it.
    """
    if n < 2:
        return 0
    remaining = n
    largest = 1
    factor = 2
    while factor * factor <= remaining:
        while remaining % factor == 0:
            remaining //= factor
            largest = factor
        factor += 1
    if remaining > 1:
        largest = remaining
    return largest


def main() -> None:
    if (
        largest_prime_factor(1) != 0
        or largest_prime_factor(2) != 2
        or largest_prime_factor(9) != 3
        or largest_prime_factor(64) != 2
        or largest_prime_factor(97) != 97
        or largest_prime_factor(13195) != 29
    ):
        raise SystemExit(1)
    print(largest_prime_factor(600851475143))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ小さな例です。割り算は切り捨て除算 //= で行い、値が float に変わらないようにしています。/= を使うと結果が浮動小数点数になり、大きな数で精度を失うことがあります。小さな例が合わなければ raise SystemExit(1) で終了コード1を返します。assert は Python を -O 付きで実行すると取り除かれるので使っていません。

答えと検証

600851475143の最大の素因数は 6857 です。

アルゴリズムがたどる過程は次のとおりです。

factor = 2 … 70      割り切れない
factor = 71          600851475143 = 71 × 8462696833
factor = 72 … 838    割り切れない
factor = 839         8462696833   = 839 × 10086647
factor = 840 … 1470  割り切れない
factor = 1471        10086647     = 1471 × 6857
factor = 1472        1472 × 1472 > 6857 → 停止、残った 6857 は素数

つまり600851475143は、表に出てくる四つの素数の積です。外側のループの本体は factor = 2 から1471まで1470回回り、実際の割り算は3回です。最初の上限 ⌊√600851475143⌋ = 775146 よりずっと早く終わるのは、「計算量」で見たとおり、2番目に大きい素因数が1471と小さいからです。

検証は三つの方法で行います。

  1. 小さな例。 二つのプログラムはどちらも、答えを出力する前に1、2、9、64、97、13195の結果が0、2、3、2、97、29かどうかを確かめ、違えば0以外の終了コードで終わります。13195 → 29は「解法の導出」で手でもたどりました。
  2. 二つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。両方の出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則として、この検査に通ったコードだけを記事に載せます。
  3. 手計算。 上の表の四つの素因数を掛け直すと 71 × 839 = 59569、後ろの二つの積は10086647で、59569 × 10086647 = 600851475143 と元の数になります。四つの数が素数であることも前提知識の事実で確かめられます。たとえば最大の素因数は平方根が約82.8で、それ以下の素数2、3、5、…、79のどれでも割り切れません。

発展

  • 候補を半分にする。 2をすべて割り取った後は、奇数の候補だけを試せば十分です。正しさの証明のどこを直す必要があるか、反復回数がどれだけ減るかを考えてみてください。2と3を先に処理して 6k ± 1 の形だけを試すとどうなるでしょうか。
  • すべての素因数と指数。 同じループで割るたびに(素数, 回数)を記録すれば、素因数分解の全体が得られます。そこから約数の個数や和を求める公式につながり、このシリーズの1〜100番の中にも約数の個数や和を扱う問題(たとえば12番、21番)があります。
  • たくさんの数をまとめて。 1から N までのすべての数の最小の素因数をあらかじめふるいで求めておけば、各数はその表をたどって数回割るだけで分解できます。エラトステネスのふるいをどう変えればよいか考えてみてください。素数を多く扱うこのシリーズの7番、10番でも、ふるいが役に立ちます。
  • 最悪の入力。 64ビットの範囲で、このコードが最も長くかかる入力はどんな数でしょうか。大きな素数と、同じくらいの大きさの二つの素数の積とで、それぞれいつ止まるかを「計算量」の max(p₂, √p₁) で比べてみてください。
  • あふれの実験。 C++ の条件を factor * factor <= remaining に変え、2^64 未満で最大の素数付近の値を渡すと何が起きるか考えてみてください(実際に実行すると長い時間がかかります)。

出典とライセンス

  • Project Euler, Problem 3 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、原文と同じ CC BY-NC-SA 4.0 で共有します。問題の帰属表示は Project Euler です。
  • Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません。
  • 解法、証明、コードは独自に書いたものです。Project Euler の解説(overview)、フォーラム、会員限定の資料やほかの解答サイトは参照していません。Motitory は Project Euler と提携・スポンサー関係になく、Project Euler がこの記事を確認したこともありません。
  • C++17 最終作業ドラフト N4659 — §6.9.1 基本型(符号なし整数の 2^N を法とする演算)、§21.3.6 <climits>(範囲は C 標準に従う)、§22.3 assert(NDEBUG)
  • Python 標準ライブラリ: Numeric Types — int の無制限の精度、//(切り捨て除算)と /(真の除算)
  • Python 言語リファレンス: The assert statement — 最適化モードでは assert が取り除かれる

ブログへ