· 開発

Project Euler 5: 1から20までのすべてで割り切れる最小の数を最小公倍数から組み立てる

1から20までのすべての数で割り切れる最小の正の整数を、1つずつ試すのではなく、ユークリッドの互除法と lcm(a, b) = a / gcd(a, b) × b の繰り返しで組み立てます。この方法がちょうど最小公倍数を与える理由の証明、整数の範囲、C++17とPython 3の実装まで扱います。

Project Eulerアルゴリズム数学最小公倍数

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

Project Euler シリーズの5問目です。「すべての数で割り切れる」という条件はそのまま公倍数の定義で、その中で最小のものが最小公倍数です。そこで今回は、答えを探し回る代わりに答えを直接組み立てます。その途中で、ユークリッドの互除法と2数の最小公倍数の公式がなぜ正しいのかを一から証明しておきます。どちらも後の問題でよく出てきます。

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

問題

問題文の例では、1から10までのすべての整数で余りなく割り切れる最小の数は2520です。同じように、1から20までのすべての整数で割り切れる最小の正の整数を求めるのが問題です。

入力はなく、上限は20に固定されています。求めるのは、その条件を満たす正の整数の最小値1つです。

前提知識

割り切れることと最小公倍数。 整数 m が正の整数 k の倍数であるとき、つまり m = k × q となる整数 q があるとき、「k は m を割り切る」と言い、k | m と書きます。いくつかの数の公倍数とはそのすべての倍数のことで、最小公倍数(lcm、least common multiple)は正の公倍数のうち最小のものです。正の整数をいくつか全部掛ければ公倍数が1つはできるので、最小公倍数は必ず存在します。

ユークリッドの互除法。 2数の最大公約数(gcd)は次の性質で求めます。a を b(> 0)で割った商を q、余りを r(0 ≤ r < b)とすると

gcd(a, b) = gcd(b, r),    r = a − q × b
gcd(a, 0) = a

です。理由は、公約数の集合が同じだからです。d が a と b を割り切れば a − q × b = r も割り切り、逆に d が b と r を割り切れば q × b + r = a も割り切ります。公約数がすべて同じなので、その最大のものも同じです。余りは毎回真に小さくなるのでいずれ0になり、そのとき残っている数が最大公約数です。たとえば gcd(840, 9) は、840 = 93 × 9 + 3、9 = 3 × 3 + 0 なので3です。

最大公約数は2数の整数結合です(ベズーの等式)。 互除法の各余り r = a − q × b は直前の2数の整数結合で、直前の2数も最初の a、b の整数結合なので、順にたどればすべての余りが s × a + t × b(s、t は整数)の形になります。最後の0でない余りが最大公約数なので、gcd(a, b) = s × a + t × b となる整数 s、t があります。上の例では 3 = 1 × 840 − 93 × 9 です。この事実は下の証明で一度使います。

素因数分解で見る最小公倍数(手計算での確認用)。 各素数の指数を比べると、最小公倍数は素数ごとに最大の指数を選んだ数です。たとえば 2520 = 2³ × 3² × 5 × 7 で、8 = 2³ と 9 = 3² は10以下で最大の2のべきと3のべきです。最終的なアルゴリズムは素因数分解を使いません。結果を手で確かめるときだけ使います。

素朴な方法

m = 1, 2, 3, … を順に試し、2から20までのすべての k で割り切れる最初の m を選べば済みます。

m = 1
while not all(m % k == 0 for k in range(2, 21)):
    m += 1

必ず終わることは保証されています。1から20までを全部掛けた数が条件を満たすので、それまでに必ず止まります。問題は、試す候補の個数が答えそのものだという点です。候補1つあたり最大19回の剰余演算なので、時間は O(答え × n) です(n は上限)。しかも上限が素数 p を越えるたびに答えは p 倍になる(p より小さい数は p を因数に持たない)ので、n が大きくなると答えも候補の数も急速に増えます。n = 10 なら候補は2520個だけですが、n = 20 ではずっと多くなります。正確に何個かは「答えと検証」で見ます。

候補を20ずつ飛ばす(答えは20の倍数なので)か2520ずつ飛ばす(答えは1から10までの公倍数で、下の補題Aにより2520の倍数なので)と、候補はその分減ります。それでも答えに比例する個数を調べることに変わりはありません。答えの大きさに関係なく終わる方法が必要です。

解法の導出

ステップ1: 問題を読み直す。 「1から n までのすべての整数で割り切れる正の整数」は 1, 2, …, n の正の公倍数で、その最小のものは定義どおり lcm(1, 2, …, n) です。これを L(n) と書きます。探す必要はなく、L(20) を計算すればよいのです。

ステップ2(補題A): すべての公倍数は最小公倍数の倍数である。 いくつかの数の最小公倍数を l、正の公倍数の1つを m とします。m を l で割って m = q × l + r(0 ≤ r < l)と書くと、それらの数はどれも m と l を割り切るので、r = m − q × l も割り切ります。つまり r は公倍数です。r > 0 なら l より小さい正の公倍数になり l の最小性に反するので r = 0、すなわち l | m です。逆に l の倍数はもちろん公倍数です。したがって公倍数の集合 = l の倍数の集合です。

ステップ3: 1つずつ畳み込んでよい。 L(k) = lcm(L(k − 1), k) です。正の整数 m について

  • m が 1, …, k の公倍数である
  • ⇔ m が 1, …, k − 1 の公倍数で、かつ k | m である
  • ⇔(補題A)L(k − 1) | m かつ k | m である
  • ⇔ m が L(k − 1) と k の公倍数である

なので2つの公倍数の集合は同じで、したがって最小の元も同じです。あとは2数の最小公倍数が求められれば十分です。

ステップ4: 2数の最小公倍数は lcm(a, b) = a / gcd(a, b) × b である。 正の整数 a、b について g = gcd(a, b) とおき、a = g × a′、b = g × b′ と書きます。

  • gcd(a′, b′) = 1 です。1より大きい d が a′ と b′ をともに割り切れば、g × d が a と b の公約数になり、g の最大性に反します。
  • l = g × a′ × b′ は公倍数です。l = a × b′ = a′ × b だからです。また l = (a / g) × b で、g | a なので a / g は整数です。
  • l が最小です。正の公倍数を m = a × x と書くと m = g × a′ × x で、b = g × b′ が m を割り切るので b′ | a′ × x です。ベズーの等式で s × a′ + t × b′ = 1 となる整数 s、t をとると x = s × a′ × x + t × b′ × x で、右辺の2項はどちらも b′ で割り切れるので b′ | x です。x ≥ b′ なので m = a × x ≥ a × b′ = l です。

小さな例で手計算による検算。 問題文の例 n = 10 をたどります。L = 1 から始めて、k = 2, …, 10 を順に畳み込みます。

k    g = gcd(L, k)    L / g × k
2    gcd(1, 2)   = 1     1 / 1 × 2  = 2
3    gcd(2, 3)   = 1     2 / 1 × 3  = 6
4    gcd(6, 4)   = 2     6 / 2 × 4  = 12
5    gcd(12, 5)  = 1    12 / 1 × 5  = 60
6    gcd(60, 6)  = 6    60 / 6 × 6  = 60
7    gcd(60, 7)  = 1    60 / 1 × 7  = 420
8    gcd(420, 8) = 4   420 / 4 × 8  = 840
9    gcd(840, 9) = 3   840 / 3 × 9  = 2520
10   gcd(2520, 10) = 10   2520 / 10 × 10 = 2520

結果の2520は問題文の例と一致します。素因数分解で見ると、4は2の指数を1から2へ、8は2から3へ上げ、9は3の指数を1から2へ上げます。6と10はすでに L の約数なので g = k となり、L は変わりません。最終的な 2520 = 2³ × 3² × 5 × 7 と合っています。

最終的なアルゴリズム

  1. L = 1 から始めます。
  2. k = 2, 3, …, n の順に、g = gcd(L, k) をユークリッドの互除法で求め、L を L / g × k に置き換えます。
  3. L を返します。n < 2 ならループは回らず1です。
  4. 答えを出力する前に小さな入力を確認します。gcd(12, 18) = 6、gcd(840, 9) = 3、lcm(4, 6) = 12、L(1) = 1、L(10) = 2520 でなければならず、素朴な方法で n = 10 を直接探した値も2520でなければなりません。1つでも違えば0以外の終了コードで終わります。

正しさの証明

k 回目の反復を終えた直後には、毎回次が成り立ちます(不変条件)。

L = lcm(1, 2, …, k) である。

反復の前(k = 1)は L = 1 = lcm(1) なので成り立ちます。k − 1 まで成り立つとすると、k 回目の反復で g = gcd(L, k) は前提知識のユークリッドの互除法が正確に計算し、新しい値 L / g × k はステップ4により lcm(L(k − 1), k) で、ステップ3によりこれが L(k) です。したがって k = n まで反復を終えると L = L(n) で、ステップ1によりこれが1から n までのすべての整数で割り切れる最小の正の整数です。反復はちょうど n − 1 回で、互除法は余りが毎回減るので毎回終わります。

計算量

時間。 gcd を n − 1 回呼びます。gcd(L, k) は最初の割り算で (k, L mod k) になり、その後は k より小さい数だけを扱います。互除法で a ≥ b の (a, b) が (b, r) → (r, ·) と2段階進むと、先頭の数は a の半分未満になります。b ≤ a / 2 なら r < b ≤ a / 2、b > a / 2 なら商が1なので r = a − b < a / 2 だからです。したがって gcd(L, k) 1回は O(log k) 回の割り算で、全体は O(n log n) 回の算術演算です(値が機械語の整数に収まる間)。保持する値は数個だけなので、追加の空間計算量は O(1) です。素朴な方法の O(答え × n) と違い、答えの大きさに依存しません。

整数の範囲。

  • 反復中の値はすべて最終結果以下です。L は反復ごとに増えるか変わらず(L(k − 1) | L(k))、途中の値 L / g は L 以下で、L / g × k はその段階の結果 L(k) です。したがって最終結果が収まる型なら途中であふれません。
  • 掛ける前に割ります。 L × k / g と書くと、L × k は結果より最大 k 倍大きくなります。n = 20 では問題になりませんが、たとえば L(43) = 2⁵ × 3³ × 5² × 7 × 11 × 13 × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 = 9419588158802421600(約 9.42 × 10¹⁸)に k = 44 を畳み込むとき、L × 44 は64ビット符号なし整数の最大値 2⁶⁴ − 1(約 1.84 × 10¹⁹)を超えますが、L / 44 × 44 = L は超えません(4と11が L を割り切るので g = 44)。
  • C++ の実装は unsigned long long を使います。規格で最低64ビットが保証されています。上の計算どおり44、45、46は新しい素数のべきではないので L(46) = L(43) は64ビットに収まり、L(47) = 47 × L(46) はあふれます。そのためコードのコメントに limit <= 46 と書きました。符号なし整数はあふれると未定義動作ではなく 2ᴺ(N は型のビット数)で割った余りに回り込みますが、これはプログラムが止まらずに黙って間違った値を出すということなので、範囲を守ることが大切です。
  • Python の int は精度に制限がないので、あふれることはありません。ただし / は浮動小数点の結果を返すので、整数除算 // を使います。g は L を割り切るので、// は正確な商です。

C++17 による実装

#include <iostream>

// Greatest common divisor by Euclid's algorithm: gcd(a, b) = gcd(b, a mod b).
unsigned long long gcd(unsigned long long a, unsigned long long b) {
    while (b != 0) {
        const unsigned long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

// Least common multiple of two positive numbers. Dividing before multiplying
// keeps every intermediate value no larger than the result.
unsigned long long lcm(unsigned long long a, unsigned long long b) {
    return a / gcd(a, b) * b;
}

// lcm(1, 2, ..., limit), folding in one number at a time. Returns 1 when limit < 2.
// Meant for limit <= 46, where the result fits in 64 bits.
unsigned long long lcm_up_to(unsigned long long limit) {
    unsigned long long result = 1;
    for (unsigned long long k = 2; k <= limit; ++k) {
        result = lcm(result, k);
    }
    return result;
}

// Naive search: the first positive m divisible by every k in 1..limit.
// Used only to cross-check the small example.
unsigned long long smallest_by_search(unsigned long long limit) {
    for (unsigned long long m = 1;; ++m) {
        bool divisible = true;
        for (unsigned long long k = 2; k <= limit; ++k) {
            if (m % k != 0) {
                divisible = false;
                break;
            }
        }
        if (divisible) {
            return m;
        }
    }
}

int main() {
    if (gcd(12, 18) != 6 || gcd(840, 9) != 3 || lcm(4, 6) != 12 ||
        lcm_up_to(1) != 1 || lcm_up_to(10) != 2520 ||
        smallest_by_search(10) != 2520) {
        return 1;
    }
    std::cout << lcm_up_to(20) << '\n';
    return 0;
}

smallest_by_search は素朴な方法をそのまま書いた関数で、公式の例 n = 10 で改良した方法と同じ値になるかを照合するためだけに使います。小さな例の確認は assert ではなく if で行いました。assert は NDEBUG が定義されたビルドでは何もしませんが、if の比較と return 1 はプログラムの観測可能な動作なので、-O2 の最適化でも消えません。確認が1つでも失敗すれば、0以外の終了コードを返します。

Python 3 による実装

def gcd(a: int, b: int) -> int:
    """Greatest common divisor by Euclid's algorithm: gcd(a, b) = gcd(b, a mod b)."""
    while b != 0:
        a, b = b, a % b
    return a


def lcm(a: int, b: int) -> int:
    """Least common multiple of two positive numbers.

    Dividing before multiplying keeps every intermediate value no larger than the result.
    """
    return a // gcd(a, b) * b


def lcm_up_to(limit: int) -> int:
    """lcm(1, 2, ..., limit), folding in one number at a time. Returns 1 when limit < 2."""
    result = 1
    for k in range(2, limit + 1):
        result = lcm(result, k)
    return result


def smallest_by_search(limit: int) -> int:
    """Naive search: the first positive m divisible by every k in 1..limit.

    Used only to cross-check the small example.
    """
    m = 1
    while True:
        if all(m % k == 0 for k in range(2, limit + 1)):
            return m
        m += 1


def main() -> None:
    if (
        gcd(12, 18) != 6
        or gcd(840, 9) != 3
        or lcm(4, 6) != 12
        or lcm_up_to(1) != 1
        or lcm_up_to(10) != 2520
        or smallest_by_search(10) != 2520
    ):
        raise SystemExit(1)
    print(lcm_up_to(20))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ小さな例です。a, b = b, a % b は右辺をすべて評価してから代入するので、C++ の一時変数 r と同じ役割を果たします。range(2, limit + 1) は終端の値を含まないので、終端を1つ上げています。小さな例が失敗すれば raise SystemExit(1) で終了コード1を返します。assert は python -O で実行すると取り除かれるので使っていません。

答えと検証

1から20までのすべての整数で割り切れる最小の正の整数は 232792560 です。

n = 10 の結果2520から続けて k = 11, …, 20 を畳み込むと、次のようになります。

k    g = gcd(L, k)    L
11   1      27720
12   12     27720
13   1      360360
14   14     360360
15   15     360360
16   8      720720
17   1      12252240
18   18     12252240
19   1      上の答え
20   20     上の答え

手で確かめられる部分は次のとおりです。

  • 素因数分解。 20以下の素数は 2, 3, 5, 7, 11, 13, 17, 19 で、20以下で最大のべきは 16 = 2⁴、9 = 3²、そして残りの素数そのものです。16 × 9 = 144、× 5 = 720、× 7 = 5040、× 11 = 55440、× 13 = 720720、× 17 = 12252240、× 19 = 上の答え です。
  • すべてで割り切れること。 20以下の k を素因数分解すると、各素数 p のべき pᵉ は k 以下、したがって20以下なので、上の指数を超えません。だから k は答えを割り切ります。
  • 最小であること。 16, 9, 5, 7, 11, 13, 17, 19 のどれか1つでも指数を1つ減らすと、その数で割り切れなくなります。たとえば 2⁴ を 2³ に減らすと16で割り切れません。
  • 2520との関係。 上の答え = 2520 × 92378 で、92378 = 2 × 11 × 13 × 17 × 19 です。n = 10 から n = 20 に進んで増えたのは、16のための2が1つと新しい素数4つだけです。

ここから、素朴な方法でこの値にたどり着くには候補をその数だけ試す必要があることもわかります。実行時間はこの記事では測っていません。

検証は三重です。

  1. 小さな例。 2つのプログラムはどちらも答えを出力する前に gcd(12, 18) = 6、gcd(840, 9) = 3、lcm(4, 6) = 12、L(1) = 1、L(10) = 2520 を確認し、素朴な方法で n = 10 を直接探して2520と一致するかも確認します。1つでも違えば0以外の終了コードで終わります。n = 10 の過程は「解法の導出」で手でもたどりました。
  2. 2つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。2つの出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則上、この検査を通ったコードだけを記事に載せます。
  3. 証明。 最小公倍数を1つずつ畳み込むとちょうど L(n) になることは「正しさの証明」で示しました。

境界のケース

  • n = 1: ループは一度も回らず1を返します。1は1で割り切れる最小の正の整数なので正しい結果です。小さな例に L(1) = 1 が入っています。n = 0 は問題にはありませんが、同じ理由で1を返し、空集合の最小公倍数を1とする慣習と合っています。
  • すでに L の約数である k(6, 10, 12, 14, 15, 18, 20): g = gcd(L, k) = k となり、L / k × k = L のまま変わりません。
  • 新しい素数 k(7, 11, 13, 17, 19): L の素因数はすべて k より小さいので g = 1 で、L は k 倍になります。
  • 素数のべき k(4, 8, 9, 16): 指数が1つだけ増えます。たとえば16では g = 8 なので、L は2倍になるだけです。
  • 掛けてから割る誤り: 「整数の範囲」のように、n が大きくなると L × k があふれることがあります。先に割れば途中の値は結果を超えません。

発展

  • 素数のべきから直接計算する。 答えの確認で使ったとおり、L(n) は n 以下の各素数 p について、n を超えない最大の pᵉ をすべて掛けたものです。エラトステネスの篩で素数を求めてこの積を作るように変えてみてください。gcd を n − 1 回呼ぶ方法と比べて、どちらが演算が少ないでしょうか。素数の扱いはこのシリーズの7問目と10問目で再び出てきます。
  • 順序を変えてもよいか。 ステップ3の議論は「1から k まで」という順序に頼っていません。数をどんな順序で畳み込んでも同じ最小公倍数になることを、同じ方法で示してみてください。
  • ベズーの係数を求める。 互除法を拡張すると、gcd(a, b) = s × a + t × b の s、t も一緒に求められます。840と9で手で試してみてください。
  • 最大公約数と最小公倍数の積。 ステップ4から gcd(a, b) × lcm(a, b) = a × b がすぐに出ます。3つ以上の数でも同じ関係が成り立つでしょうか。小さな反例を探してみてください。

出典とライセンス

  • Project Euler, Problem 5 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、原文と同じ CC BY-NC-SA 4.0 で共有します。問題の著作者表示は Project Euler です。
  • Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません(ads: false)。
  • C++ の整数型の最小幅と符号なし整数の算術は C++ 規格の作業草案 [basic.fundamental]、assert と NDEBUG は [assertions.assert] に従いました。Python の同時代入は Assignment statements、assert と -O は The assert statement、int の精度と /・// の結果は Numeric Types に従いました。
  • 解法、証明、コードは自分で書いたものです。ユークリッドの互除法と最小公倍数の性質は上で直接証明しました。Project Euler の解説(overview)、フォーラム、会員限定の資料や他の解答サイトは参照していません。Motitory は Project Euler と提携・後援の関係になく、Project Euler がこの記事を確認したこともありません。

ブログへ