Project Euler 1: 3または5の倍数の和をO(1)で求める

Project Euler シリーズの最初の問題です。1行のループでも解けますが、同じ答えをループなしで定数回の演算から求めようとすると、包除原理と等差数列の和という二つの道具がそろって登場します。どちらも後の問題で繰り返し使うので、最初にきちんと導いておきます。

問題

1000 未満の自然数のうち、3の倍数または5の倍数であるものをすべて足した値を求めます。問題には小さな例が付いています。10未満で条件を満たす自然数は3、5、6、9で、その和は23です。

境界に注意が必要です。「1000未満」なので1000そのものは含みません。1000は5の倍数なので、境界を「1000以下」と読み違えると答えがちょうど1000だけ大きくなります。

前提知識

商(床関数による割り算)。 ⌊a / b⌋ は a を b で割った商、つまり a/b 以下の最大の整数です。正の整数 n と k について、「n 未満の k の正の倍数」は k, 2k, …, mk で、m は mk ≤ n − 1 を満たす最大の整数なので m = ⌊(n − 1) / k⌋ になります。「未満」なので n ではなく n − 1 を割ります。たとえば n = 10、k = 5 なら ⌊9 / 5⌋ = 1 で、5だけを数えます。⌊10 / 5⌋ = 2 を使うと10まで入ってしまいます。

1から m までの和。 1 + 2 + … + m = m(m + 1) / 2 です。同じ和を逆順に書いて上下に足すと理由がわかります。

  1   +   2   + … +   m
  m   + (m−1) + … +   1
(m+1) + (m+1) + … + (m+1)   ← m 個

2行の合計は m(m + 1) で、もとの和はその半分です。m と m + 1 のどちらかは必ず偶数なので m(m + 1) は常に偶数で、2で割っても整数が残ります。つまりコードで整数の割り算を使っても切り捨てられる部分はありません。

等差数列の和。 k, 2k, …, mk は公差 k の等差数列で、各項から k をくくり出すと k(1 + 2 + … + m) = k · m(m + 1) / 2 になります。

包除原理。 二つの集合 A、B について |A ∪ B| = |A| + |B| − |A ∩ B| です。A と B の両方に属する要素は |A| と |B| で1回ずつ、合わせて2回数えられるので、1回分を引きます。要素の個数ではなく要素の値の和でも同じことが成り立ちます。共通の要素の値が2回足されているので、1回引けばよいのです。

最小公倍数。 3の倍数でもあり5の倍数でもある数は3と5の公倍数で、公倍数は最小公倍数の倍数です。3と5は互いに素なので、最小公倍数は 3 × 5 = 15 です。

素朴な方法

1から999までを順に見て、3か5で割り切れれば足します。

total = 0
for i in range(1, 1000):
    if i % 3 == 0 or i % 5 == 0:
        total += i

上限が n なら n − 1 個の数を1回ずつ調べるので、時間計算量は O(n)、空間は O(1) です。n = 1000 なら一瞬で終わります。しかし n が 10^12 や 10^18 のように大きくなると、ループ回数もそれだけ増え、現実的な時間では終わりません。

よくある間違いもあります。「3の倍数の和」と「5の倍数の和」を別々のループで求めて足すと、15、30、45のように両方に属する数が2回足されます。上のコードは or で1回だけ足すのでこの問題はありませんが、下の公式による解法ではこの重複を明示的に処理する必要があります。

解法の導出

ステップ1: 小さな例を集合に分ける。 n = 10 のとき、10未満の3の倍数は {3, 6, 9}、5の倍数は {5} です。両方に属する数、つまり15の倍数は10未満にはありません。

ステップ2: 包除原理を当てはめる。 「n 未満の k の正の倍数の和」を S(k, n) とすると、

答え(n) = S(3, n) + S(5, n) − S(15, n)

です。n = 10 なら S(3, 10) = 3 + 6 + 9 = 18、S(5, 10) = 5、S(15, 10) = 0 で、18 + 5 − 0 = 23 と問題の例に一致します。

ステップ3: S(k, n) を閉じた式にする。 前提知識の二つの式を合わせると

m = ⌊(n − 1) / k⌋
S(k, n) = k + 2k + … + mk = k · m(m + 1) / 2

です。

ステップ4: 小さな例で公式を検算する。 n = 10 で

• k = 3: m = ⌊9 / 3⌋ = 3、S = 3 · 3 · 4 / 2 = 18
• k = 5: m = ⌊9 / 5⌋ = 1、S = 5 · 1 · 2 / 2 = 5
• k = 15: m = ⌊9 / 15⌋ = 0、S = 15 · 0 · 1 / 2 = 0

となり、18 + 5 − 0 = 23 です。ステップ2で書き出して足した値と同じです。あとは n = 1000 に同じ式を使うだけです。

最終的なアルゴリズムと計算量

1. S(k, n) を計算する関数を作ります。n が0なら0を返し、そうでなければ m = ⌊(n − 1) / k⌋ を求めて k · (m(m + 1) / 2) を返します。
2. S(3, n) + S(5, n) − S(15, n) を返します。
3. 出力の前に n = 10 の結果が23かを確かめ、違えば0以外の終了コードで終わります。

正しさ。 n 未満の自然数 x が和に何回入るかを場合分けします。

• x が3の倍数で5の倍数でない: S(3) で1回足され、S(5) と S(15) には現れません。1回。
• x が5の倍数で3の倍数でない: S(5) でだけ1回。1回。
• x が3と5の両方の倍数: 15の倍数なので S(3) と S(5) で1回ずつ足され、S(15) で1回引かれます。1 + 1 − 1 = 1回。
• x がどちらの倍数でもない: どの項にも現れません。0回。

条件を満たす数はちょうど1回、満たさない数は1回も足されないので、結果は求める和と一致します。各 S(k, n) が n 未満の k の倍数を漏れも重複もなく足すことは前提知識で示しました。

時間と空間。 n の大きさに関係なく、割り算3回と数回の掛け算・足し算だけなので、固定幅の整数では時間計算量 O(1)、空間計算量 O(1) です。Python の整数には大きさの上限がないため、非常に大きな n では一つの演算のコストが桁数に応じて増えますが、演算の回数は定数のままです。

整数の範囲。 n = 1000 なら最大の中間値でも数十万程度なので、unsigned long long はもちろん、標準的な32ビット符号なし整数にも余裕で収まります。それでも一般の n を考えて範囲を確認しておきます。

• C++ の実装は unsigned long long を使います。C++ 標準はこの型が少なくとも64ビットであることを保証し(C 標準の ULLONG_MAX の最小値に従う)、符号なし整数の演算は 2^N を法として定義されているので、符号付き整数と違って未定義動作になりません。ただし範囲を超えると黙って誤った値になるので、範囲内で使う必要があります。
• S(k, n) の中では m(m + 1) / 2 を先に計算してから k を掛けます。2で割る前の積は k = 3 のときに最大で、n ≤ 2^32 なら m ≤ (2^32 − 1) / 3、m(m + 1) は約 2.05 × 10^18 と 2^64(約 1.84 × 10^19)より小さくなります。ただし全体で最大の中間値は、引き算の前の S(3) + S(5) です。およそ n²/6 + n²/10、つまり n² の27%程度で、n = 2^32 では約 4.92 × 10^18 ですが、これも 2^64 未満です。したがってこの C++ のコードは n ≤ 2^32 で安全です。
• 引き算 S(3) + S(5) − S(15) は符号なし整数で負になる心配がありません。15の倍数はすべて3の倍数でもあるので、S(15) ≤ S(3) だからです。
• Python の int は精度に上限がないので、n がどれほど大きくてもあふれません。// は床関数による割り算で、ここではオペランドがすべて0以上なので C++ の / と同じ結果になります。

C++17 による実装

#include <iostream>

// Sum of the positive multiples of k that are less than limit.
unsigned long long sum_of_multiples_below(unsigned long long k, unsigned long long limit) {
    if (limit == 0) {
        return 0;
    }
    const unsigned long long m = (limit - 1) / k;
    return k * (m * (m + 1) / 2);
}

// Sum of the natural numbers less than limit that are divisible by 3 or 5.
unsigned long long sum_of_multiples_of_3_or_5_below(unsigned long long limit) {
    return sum_of_multiples_below(3, limit) + sum_of_multiples_below(5, limit) -
           sum_of_multiples_below(15, limit);
}

int main() {
    if (sum_of_multiples_of_3_or_5_below(10) != 23) {
        return 1;
    }
    std::cout << sum_of_multiples_of_3_or_5_below(1000) << '\n';
    return 0;
}

limit == 0 を別に扱っているのは、符号なし整数で 0 − 1 を計算すると最大値に回り込むからです。この問題で0を渡すことはありませんが、この判定があるので、limit == 0 のときは回り込まずに正しい結果の0を返します。小さな例の確認には assert ではなく if を使いました。assert は NDEBUG を定義したビルドで消えるからです。

Python 3 による実装

def sum_of_multiples_below(k: int, limit: int) -> int:
    """Sum of the positive multiples of k that are less than limit."""
    if limit <= 0:
        return 0
    m = (limit - 1) // k
    return k * (m * (m + 1) // 2)


def sum_of_multiples_of_3_or_5_below(limit: int) -> int:
    """Sum of the natural numbers less than limit that are divisible by 3 or 5."""
    return (
        sum_of_multiples_below(3, limit)
        + sum_of_multiples_below(5, limit)
        - sum_of_multiples_below(15, limit)
    )


def main() -> None:
    if sum_of_multiples_of_3_or_5_below(10) != 23:
        raise SystemExit(1)
    print(sum_of_multiples_of_3_or_5_below(1000))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ関数構成です。小さな例が合わなければ raise SystemExit(1) で終了コード1を返します。assert は Python を -O 付きで実行すると取り除かれるので使っていません。

答えと検証

n = 1000 に公式を当てはめると次のとおりです。

• k = 3: m = ⌊999 / 3⌋ = 333、S = 3 · (333 · 334 / 2) = 3 · 55611 = 166833
• k = 5: m = ⌊999 / 5⌋ = 199、S = 5 · (199 · 200 / 2) = 5 · 19900 = 99500
• k = 15: m = ⌊999 / 15⌋ = 66、S = 15 · (66 · 67 / 2) = 15 · 2211 = 33165

したがって答えは 166833 + 99500 − 33165 = 233168 です。

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

1. 小さな例。 二つのプログラムはどちらも、答えを出力する前に n = 10 の結果が問題にある23かどうかを確かめ、違えば0以外の終了コードで終わります。「解法の導出」のステップ4でも同じ値を手で確認しました。
2. 二つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。両方の出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則として、この検査に通ったコードだけを記事に載せます。
3. 手計算。 すぐ上で三つの項を手で計算し、同じ値を得ました。個数でも確かめられます。1000未満の3の倍数は333個、5の倍数は199個、15の倍数は66個なので、条件を満たす数は 333 + 199 − 66 = 466 個です。上の素朴なループで個数を数えて比べてみるのもよい練習です。

発展

• 境界を変える。 「1000以下」なら何が変わるでしょうか。関数に limit = 1001 を渡せば済みます。1000は5の倍数なので、答えは1000だけ大きくなります。
• 約数を三つに増やす。 3、5、7のどれかで割り切れる数の和は、包除原理を三つの集合に広げて S(3) + S(5) + S(7) − S(15) − S(21) − S(35) + S(105) で求めます。約数が互いに素でない場合(たとえば4と6)、共通部分は積ではなく最小公倍数(12)の倍数になる点に注意してください。
• 一般化。 約数が d 個なら包除原理の項は 2^d − 1 個です。約数が多いと項の数が急に増えるので、どこかで別の方法のほうが有利になるかもしれません。
• オーバーフローの実験。 n が 2^32 を大きく超えたとき、C++ のコードがどこで最初にあふれるかを計算してみてください。同じ n を Python で計算して比べると違いが確かめられます。
• このシリーズの後の問題。 1から m までの和のような閉じた式は、1〜100番の中の別の問題でもまた出てきます。まずループで解き、閉じた式にできないかを考える習慣をつけておくと役に立ちます。

出典

• Project Euler, Problem 1 — 元の問題
• Project Euler, Copyright — 問題文のライセンス(CC BY-NC-SA 4.0)
• 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 が取り除かれる