Project Euler 8: 1000桁の数で隣り合う13桁の最大の積を、0を数えるスライディングウィンドウで求める

Project Euler シリーズの8問目です。窓ごとに数字を13個掛けてみる方法だけでも、答えはすぐに出ます。この記事はもう一歩進んで、窓を1つずらすたびに出ていく数字で割り、入ってくる数字を掛けるスライディングウィンドウを組み立てます。そこで突き当たるのが0です。0を掛ける代わりに個数だけを数えれば0で割ることが起きないことを示し、積が32ビットに収まらない理由も確かめます。

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

問題

元の問題には1000桁の数が1つ与えられます。その数の中で連続して隣り合う4つの数字の積のうち最大のものは 9 × 9 × 8 × 9 = 5832 だと、例として示されています。同じ数の中で連続して隣り合う13個の数字の積の最大値を求めるのが問題です。

入力はありません。数は元の問題のページにあり、この記事のコードにそのまま写してあります。この記事では桁数を n(= 1000)、窓の長さを k(= 13)と書き、左から i 番目の数字(0から数える)を dᵢ と書きます。長さ k の窓は n − k + 1 = 988 個あります。求めるのは整数1つです。

前提知識

スライディングウィンドウ。 長さ k の窓を左端から右端まで1つずつずらすと、1回ずらすたびに変わるのは左から出ていく数字1つと右から入ってくる数字1つだけです。窓について持ち運ぶ値をこの2つの数字だけで更新できれば、窓ごとに最初から計算し直す必要はありません。

積に0があれば積は0。 窓に0が1つでもあれば、その窓の積は0です。そして0では割れません。そのため「出ていく数字で割る」という考えは、出ていく数字が0のときにはそのままでは通用しません(解法の導出のステップ2)。

割り切れる除算。 整数 a が a = d × b(d ≠ 0)と書けるなら、a / d は余りなくちょうど b です。整数除算(C++ の /、Python の //)が切り捨てをしても、この場合は切り捨てるものがありません。解法では除算をいつもこの形でしか使いません。

文字から数字へ。 数は文字列として持ちます。文字 '0' から '9' までは文字コードが順に1ずつ大きくなるので、数字の文字 c の値は c − '0' です。プログラムは答えを出す前に、文字列が本当に1000個の '0'〜'9' だけでできているかを確かめます。

固定幅の整数。 符号なし32ビット整数が持てる最大の値は 2³² − 1 = 4294967295、64ビットなら 2⁶⁴ − 1 = 18446744073709551615 です。超えると値が黙って回り込み、間違った答えになります。13個の数字の積がどちらに収まるかは「計算量」で確かめます。

素朴な方法

窓の開始位置ごとに k 個の数字を最初から掛け、その中で最大の値を選びます。

best = 0
for start in range(n - k + 1):
    product = 1
    for j in range(start, start + k):
        product *= d[j]
    best = max(best, product)

正しい方法で、掛け算は (n − k + 1) × k = 988 × 13 = 12844 回なので、この問題では十分に速いです。費用は O(nk) で、窓が長くなるほど1つの窓にかかる仕事が k に比例して増えます。ところが隣り合う2つの窓は k − 1 個の数字を共有しているので、同じ掛け算をほぼ k 回ずつ繰り返していることになります。以下では窓ごとの仕事を定数に減らします。この素朴な方法は捨てずに、実装の中で検算用に使います。

解法の導出

ステップ1: 問題の例を手で確かめる。 元の問題で強調されている連続した4つの数字 9, 9, 8, 9 の積は、9 × 9 = 81、81 × 8 = 648、648 × 9 = 5832 で、問題の例と一致します。この4桁が本当に最大の窓かどうかは手ではすべて確かめられないので、プログラムが k = 4 で計算して5832が出ることを確かめます。

アルゴリズムを手でたどるための小さな例も1つ決めます。数字列 3021945、k = 3 なら窓は5つです。

3 0 2 → 0
0 2 1 → 0
2 1 9 → 18
1 9 4 → 36
9 4 5 → 180
最大値: 180

ステップ2: 割って掛けて窓をずらす考えと、その落とし穴。 窓の積 P が分かっていれば、1つずらすときに出ていく数字 dₒᵤₜ と入ってくる数字 dᵢₙ で、新しい積を P / dₒᵤₜ × dᵢₙ として得られそうです。0がなければ正しいです。しかし上の例では最初の窓 3 0 2 の積が0で、その0が窓から出ていく瞬間に「0 / 0」が必要になります。積がいったん0になると、その中に何が掛けられていたかという情報が失われるからです。

ステップ3: 0は掛けずに数える。 窓について2つの値を持ち運びます。

• Z: 窓の中の0の個数
• P: 窓の中の0でない数字の積(そういう数字がなければ1)

すると窓の積は、Z > 0 なら0、Z = 0 なら P です。1つずらすときは

• 出ていく数字が0なら Z を1減らし、そうでなければ P をその数字で割ります。
• 入ってくる数字が0なら Z を1増やし、そうでなければ P にその数字を掛けます。

なぜ0で割らず、除算が割り切れるのか。 除算は出ていく数字 d が0でないときにしか行いません。そして P は d を含む窓の中の0でない数字の積なので、P = d × (残りの0でない数字の積) であり、P / d は余りなくちょうどその残りの積です。0は一度も P に掛けられないので P はいつも1以上で、窓に0が残っている間も P は「0でない数字の積」という意味をそのまま保ちます。「0が残っているのに割る」という状況そのものが起きません。0は Z を通してだけ出入りします。

ステップ4: 先に除いてから加える。 1つずらすとき、出ていく数字を先に除き、入ってくる数字を後で加えます。すると P に掛けられている数字はどの瞬間にも k 個以下なので、P ≤ 9ᵏ です。順序を逆にすると一時的に k + 1 個が掛けられ、上限が 9ᵏ⁺¹ に大きくなります(「計算量」の整数の範囲)。

3021945、k = 3 でたどってみます。i は入ってくる数字の位置で、i ≥ 3 からは先に dᵢ₋₃ を除きます。

i=0: 3 が入る            Z=0, P=3   (窓がまだ埋まっていない)
i=1: 0 が入る            Z=1, P=3
i=2: 2 が入る            Z=1, P=6   窓 302 → Z>0 なので 0, best=0
i=3: 3 が出て 1 が入る   Z=1, P=2   窓 021 → 0
i=4: 0 が出て 9 が入る   Z=0, P=18  窓 219 → 18, best=18
i=5: 2 が出て 4 が入る   Z=0, P=36  窓 194 → 36, best=36
i=6: 1 が出て 5 が入る   Z=0, P=180 窓 945 → 180, best=180

ステップ1で窓ごとに直接掛けた結果と一致します。i = 3 で P = 6 を3で割った 6 / 3 = 2 は割り切れ、i = 4 で0が出ていくときは割らずに Z だけを減らしました。

最終的なアルゴリズム

1. Z = 0、P = 1、best = 0 から始めます。
2. i = 0 から n − 1 まで:
   • i ≥ k なら dᵢ₋ₖ を除きます(0なら Z − 1、そうでなければ P / dᵢ₋ₖ)。
   • dᵢ を加えます(0なら Z + 1、そうでなければ P × dᵢ)。
   • i + 1 ≥ k(窓が埋まっている)かつ Z = 0 かつ P > best なら best = P とします。
3. best を返します。窓の積が0の場合は best の初期値0と同じなので、別に比べる必要はありません。

前提は 1 ≤ k ≤ n です(「境界のケース」)。

正しさの証明

i 番目の数字を処理した直後の窓を W(i) = { dⱼ : max(0, i − k + 1) ≤ j ≤ i } とします。次の不変条件が成り立ちます。「Z は W(i) の中の0の個数であり、P は W(i) の中の0でない数字の積である。」

• 開始時。 処理前の窓は空で、Z = 0、P = 1(空の積)です。
• i < k のとき。 加えるだけで W(i) = W(i − 1) ∪ {dᵢ} なので、加える規則どおりに Z か P を更新すれば不変条件は保たれます。
• i ≥ k のとき。 W(i) は W(i − 1) から dᵢ₋ₖ を除いて dᵢ を加えたものです。dᵢ₋ₖ が0なら Z − 1 が正しく、0でなければステップ3の理由で P / dᵢ₋ₖ はちょうど残りの0でない数字の積です。その後 dᵢ を加えるのは前の場合と同じです。
• したがって i + 1 ≥ k である各 i で、W(i) はちょうど dᵢ₋ₖ₊₁ … dᵢ の k 個で、その積は Z > 0 なら0、そうでなければ P です。そういう i は k − 1 から n − 1 までの n − k + 1 個で、すべての窓を一度ずつ見ます。best はそれらの積の最大値(すべて0なら0)です。

計算量

• 時間。 数字1つにつき比較が数回と、多くても除算1回と掛け算1回なので O(n) です。窓の長さ k に依存しません。素朴な方法の O(nk) と比べて k 倍減りました。
• 空間。 数字列そのものを除けば Z、P、best などの値がいくつかあるだけなので O(1) です。

整数の範囲。

• P と best は k 個以下の0でない数字の積なので(ステップ4)、9¹³ = 2541865828329(約 2.54 × 10¹²)以下です。これは符号なし32ビットの最大値 4294967295 より大きく、2⁴¹ = 2199023255552 < 9¹³ < 2⁴² = 4398046511104 なので42ビットが必要です。つまり32ビット整数では足りません。 実際の窓の積が32ビットを超えるかどうかは数次第ですが、前もっては分からないので、上限から型を決めます。
• C++ では unsigned long long を使います。C++ 規格はこの型が少なくとも64ビットであることを保証しているので、9¹³ を余裕をもって持てます。検算用の素朴な方法の途中の積も k 個以下の数字の積なので、同じ上限に収まります。
• Python の int は精度に制限がなく、あふれることはありません。
• 添字と Z は std::size_t で、n = 1000 以下です。Z は窓の中の0だけを数えるので、0より下にはなりません(除く0はいつも窓の中にあった0です)。

C++17 による実装

#include <cstddef>
#include <iostream>
#include <string>

// The 1000-digit number from the problem, 50 digits per line.
const std::string DIGITS =
    "73167176531330624919225119674426574742355349194934"
    "96983520312774506326239578318016984801869478851843"
    "85861560789112949495459501737958331952853208805511"
    "12540698747158523863050715693290963295227443043557"
    "66896648950445244523161731856403098711121722383113"
    "62229893423380308135336276614282806444486645238749"
    "30358907296290491560440772390713810515859307960866"
    "70172427121883998797908792274921901699720888093776"
    "65727333001053367881220235421809751254540594752243"
    "52584907711670556013604839586446706324415722155397"
    "53697817977846174064955149290862569321978468622482"
    "83972241375657056057490261407972968652414535100474"
    "82166370484403199890008895243450658541227588666881"
    "16427171479924442928230863465674813919123162824586"
    "17866458359124566529476545682848912883142607690042"
    "24219022671055626321111109370544217506941658960408"
    "07198403850962455444362981230987879927244284909188"
    "84580156166097919133875499200524063689912560717606"
    "05886116467109405077541002256983155200055935729725"
    "71636269561882670428252483600823257530420752963450";

// The value of the digit at position i of a string of '0'..'9'.
unsigned long long digit_at(const std::string& digits, std::size_t i) {
    return static_cast<unsigned long long>(digits[i] - '0');
}

// Largest product of k adjacent digits, for 1 <= k <= digits.size().
// Rolling window: zeros counts the zeros in the window and nonzero_product is the
// product of its other digits, so the window's product is 0 if zeros > 0 and
// nonzero_product otherwise. Only nonzero digits are ever divided out.
unsigned long long largest_product(const std::string& digits, std::size_t k) {
    unsigned long long nonzero_product = 1;
    std::size_t zeros = 0;
    unsigned long long best = 0;
    for (std::size_t i = 0; i < digits.size(); ++i) {
        if (i >= k) {
            // Remove the digit leaving the window before adding the new one,
            // so the window never holds more than k digits.
            const unsigned long long outgoing = digit_at(digits, i - k);
            if (outgoing == 0) {
                --zeros;
            } else {
                nonzero_product /= outgoing;
            }
        }
        const unsigned long long incoming = digit_at(digits, i);
        if (incoming == 0) {
            ++zeros;
        } else {
            nonzero_product *= incoming;
        }
        if (i + 1 >= k && zeros == 0 && nonzero_product > best) {
            best = nonzero_product;
        }
    }
    return best;
}

// Multiplies every window from scratch. Used only to cross-check largest_product.
unsigned long long largest_product_brute(const std::string& digits, std::size_t k) {
    unsigned long long best = 0;
    for (std::size_t start = 0; start + k <= digits.size(); ++start) {
        unsigned long long product = 1;
        for (std::size_t j = start; j < start + k; ++j) {
            product *= digit_at(digits, j);
        }
        if (product > best) {
            best = product;
        }
    }
    return best;
}

int main() {
    if (DIGITS.size() != 1000) {
        return 1;
    }
    for (const char c : DIGITS) {
        if (c < '0' || c > '9') {
            return 1;
        }
    }
    if (largest_product("3021945", 3) != 180 || largest_product("1020304", 3) != 0) {
        return 1;
    }
    if (largest_product(DIGITS, 4) != 5832) {
        return 1;
    }
    for (std::size_t k = 1; k <= 13; ++k) {
        if (largest_product(DIGITS, k) != largest_product_brute(DIGITS, k)) {
            return 1;
        }
    }
    std::cout << largest_product(DIGITS, 13) << '\n';
    return 0;
}

1000桁の数は50桁ずつ20行に分けて書きました。隣り合う文字列リテラルはコンパイル時に1つにつながります。行が抜けたり数字を写し間違えたりした場合に備えて、main はまず長さが1000であること、すべての文字が '0'〜'9' であることを確かめます。

小さな例は3つです。「解法の導出」で手でたどった 3021945(k = 3 → 180)、すべての窓に0が入っていて答えが0になるはずの 1020304(k = 3)、そして問題の例(k = 4 → 5832)です。続けて k = 1 から13まで、スライディングウィンドウと largest_product_brute(素朴な方法)が同じ値を出すかを照合します。Z や P を更新する順序、除算の条件、窓が埋まったかを見る条件に誤りがあれば、この照合で引っかかります。確認は assert ではなく if と return 1 で行いました。assert は NDEBUG が定義されたビルドでは何もしないからです。

Python 3 による実装

# The 1000-digit number from the problem, 50 digits per line.
DIGITS = (
    "73167176531330624919225119674426574742355349194934"
    "96983520312774506326239578318016984801869478851843"
    "85861560789112949495459501737958331952853208805511"
    "12540698747158523863050715693290963295227443043557"
    "66896648950445244523161731856403098711121722383113"
    "62229893423380308135336276614282806444486645238749"
    "30358907296290491560440772390713810515859307960866"
    "70172427121883998797908792274921901699720888093776"
    "65727333001053367881220235421809751254540594752243"
    "52584907711670556013604839586446706324415722155397"
    "53697817977846174064955149290862569321978468622482"
    "83972241375657056057490261407972968652414535100474"
    "82166370484403199890008895243450658541227588666881"
    "16427171479924442928230863465674813919123162824586"
    "17866458359124566529476545682848912883142607690042"
    "24219022671055626321111109370544217506941658960408"
    "07198403850962455444362981230987879927244284909188"
    "84580156166097919133875499200524063689912560717606"
    "05886116467109405077541002256983155200055935729725"
    "71636269561882670428252483600823257530420752963450"
)


def digit_at(digits: str, i: int) -> int:
    """The value of the digit at position i of a string of '0'..'9'."""
    return ord(digits[i]) - ord("0")


def largest_product(digits: str, k: int) -> int:
    """Largest product of k adjacent digits, for 1 <= k <= len(digits).

    Rolling window: zeros counts the zeros in the window and nonzero_product is the
    product of its other digits, so the window's product is 0 if zeros > 0 and
    nonzero_product otherwise. Only nonzero digits are ever divided out.
    """
    nonzero_product = 1
    zeros = 0
    best = 0
    for i in range(len(digits)):
        if i >= k:
            # Remove the digit leaving the window before adding the new one,
            # so the window never holds more than k digits.
            outgoing = digit_at(digits, i - k)
            if outgoing == 0:
                zeros -= 1
            else:
                nonzero_product //= outgoing
        incoming = digit_at(digits, i)
        if incoming == 0:
            zeros += 1
        else:
            nonzero_product *= incoming
        if i + 1 >= k and zeros == 0 and nonzero_product > best:
            best = nonzero_product
    return best


def largest_product_brute(digits: str, k: int) -> int:
    """Multiplies every window from scratch. Used only to cross-check largest_product."""
    best = 0
    for start in range(len(digits) - k + 1):
        product = 1
        for j in range(start, start + k):
            product *= digit_at(digits, j)
        if product > best:
            best = product
    return best


def main() -> None:
    if len(DIGITS) != 1000:
        raise SystemExit(1)
    for c in DIGITS:
        if c < "0" or c > "9":
            raise SystemExit(1)
    if largest_product("3021945", 3) != 180 or largest_product("1020304", 3) != 0:
        raise SystemExit(1)
    if largest_product(DIGITS, 4) != 5832:
        raise SystemExit(1)
    for k in range(1, 14):
        if largest_product(DIGITS, k) != largest_product_brute(DIGITS, k):
            raise SystemExit(1)
    print(largest_product(DIGITS, 13))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ確認です。括弧の中に並べて書いた文字列リテラルは、C++ と同じく1つにつながります。数字の値は C++ の c - '0' と同じ意味の ord(c) - ord("0") で求め、除算は割り切れる場合にしか使わないので、整数除算 // が切り捨てるものはありません。小さな例が失敗すれば raise SystemExit(1) で終了コード1を返します。assert は python -O で実行すると取り除かれるので使っていません。

答えと検証

隣り合う13桁の数字の積の最大値は 23514624000 です。

手でたどれる部分は次のとおりです。

• どの窓か。 この値を出す窓は、1000桁の数の左から198番目から210番目(1から数える)の数字 5576689664895 です。コードの4行目の末尾の 557 と5行目の先頭の 6689664895 で、すぐ次の211番目の数字は0です。
• 積。 これらの数字は5が3個、7が1個、6が4個、8が2個、9が2個、4が1個なので、積は 5³ × 7 × 6⁴ × 8² × 9² × 4 = 125 × 7 × 1296 × 64 × 81 × 4 = 2¹² × 3⁸ × 5³ × 7 で、計算すると上の値になります。9¹³ 以下で32ビットの最大値より大きいので、「計算量」の整数の範囲の議論が実際に必要になる場合です。
• これより大きな窓がないことの手による点検(参考)。 0を含む窓の積は0なので、0で区切られた区間のうち長さが13以上のものだけを見れば十分です。この値に届くには、13個の数字の幾何平均がおよそ6.28以上でなければなりません。数字ごとに9に比べてどれだけ足りないかを常用対数で測ると(1はおよそ0.95、2はおよそ0.65、3はおよそ0.48、4はおよそ0.35)、1つの窓でその合計がおよそ2.03以下でなければならないので、1・2・3をいくつも含む窓はすぐに外れます。そうして区間ごとに調べると、上の窓だけが残ります。ただしこれは手計算なので誤りがありえます。決定的な確認は下のプログラムによる検証です。

検証は三重です。

1. 小さな例と素朴な方法との照合。 2つのプログラムはどちらも答えを出力する前に、写した数が1000個の数字であること、3021945(k = 3)が180、1020304(k = 3)が0、問題の例(k = 4)が5832であることを確かめ、k = 1 から13まで、スライディングウィンドウと素朴な方法の結果がすべて一致するかを確かめます。1つでも違えば0以外の終了コードで終わります。
2. 2つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。2つの出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則上、この検査を通ったコードだけを記事に載せます。
3. 証明。 スライディングウィンドウがすべての窓の積を正確に見ること、0で割ったり余りの出る除算をしたりしないことは、「解法の導出」と「正しさの証明」で示しました。

境界のケース

• k = 1: 窓は数字1つで、答えは最大の数字です。出ていく数字を除くと P は1に戻ります(0でなければ d / d = 1、0なら Z だけが変わる)。プログラムの照合は k = 1 を含みます。
• すべての窓に0があるとき: Z > 0 の窓は best を変えないので、best は初期値の0のままで、これが正しい答えです(1020304、k = 3)。
• 0が出ていったばかりの窓: Z が0に戻った瞬間、P はすでに窓の中の0でない数字の積なので、そのまま窓の積です。たどった例の i = 4 がこの場合です。
• k > n: 窓が1つもなく0を返します。しかし「窓がない」と「最大の積が0」を区別できないので、関数は k ≤ n を前提にしています。
• k = 0: 前提の外です。i ≥ 0 が常に真なので、まだ窓に加えていない数字を除くことになり、不変条件が崩れます(C++ では、その数字が0なら Z が0より下になって符号なしの値が回り込みます)。再利用するなら先に入力を確かめる必要があります。
• 除く順序と加える順序: 加えてから除いても答えは同じですが、一時的に k + 1 個が掛けられて上限が 9¹⁴ に大きくなります。この記事の順序は上限を 9ᵏ に抑えます。

発展

• もっと長い窓。 P ≤ 9ᵏ で、9²⁰(約 1.22 × 10¹⁹)< 2⁶⁴ − 1(約 1.84 × 10¹⁹)< 9²¹(約 1.09 × 10²⁰)なので、64ビットの unsigned long long で安全なのは k ≤ 20 までです。それより長い窓なら、C++ では大きな整数を自分で扱うか、積ではない別の比べ方が必要です。
• 対数の和で比べる。 積の大小は対数の和の大小と同じなので、あふれずに窓を比べられます。ただし浮動小数点の誤差でほぼ等しい2つの窓を取り違えることがあるので、候補を選んだ後で正確な整数の積で確かめ直すほうが安全です。
• 除算なしで。 スタック2つで作るキューに、各スタックの累積の積を一緒に持たせると、除算なしで窓の積を保てます。元に戻せない演算(例えば窓の中の最大値)にもそのまま使えるという利点があります。どうすれば数字1つあたり償却定数時間になるか考えてみてください。
• 積の代わりに和。 同じ問題を積ではなく和で問うと、0の問題は消え、引き算だけで窓をずらせます。積では0が「元に戻せない値」であることが、この記事の Z を必要にしたのです。

出典とライセンス

• Project Euler, Problem 8 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、コードに写した1000桁の数も元の問題から取ったものです。これらの部分は原文と同じ CC BY-NC-SA 4.0 で共有します。問題の著作者表示は Project Euler です。
• Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません(ads: false)。
• Python の int の精度は Numeric Types — int, float, complex、assert と -O は The assert statement に従いました。C++ 側の事実(数字の文字 '0'〜'9' の値が順に1ずつ大きくなること、unsigned long long が少なくとも64ビットであること、NDEBUG が定義されると assert が何もしないこと)は C++ 規格(ISO/IEC 14882)の規定です。
• 解法、証明、コードは自分で書いたものです。Project Euler の解説(overview)、フォーラム、会員限定の資料や他の解答サイトは参照していません。Motitory は Project Euler と提携・後援の関係になく、Project Euler がこの記事を確認したこともありません。