· 開発

Project Euler 4: 3桁の数2つの積で最大の回文数を大きい方から探す

3桁の数2つの積のうち最大の回文数を、大きい積から順に探して枝刈りしながら求めます。枝刈りで答えを見落とさない理由の証明、偶数桁の回文数と11の関係、整数の範囲、C++17とPython 3の実装まで扱います。

Project Eulerアルゴリズム数学回文数

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

Project Euler シリーズの4問目です。すべての組み合わせを試してもすぐ終わる大きさなので、今回の関心は「もっと速い公式」よりも探索をどこで打ち切ってよいかです。大きい積から順に見ていき、見る必要のない部分を切り捨てますが、切り捨てた中に答えが隠れていないことをきちんと証明しておきます。この考え方は後の問題でも繰り返し使います。

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

問題

前から読んでも後ろから読んでも同じになる数を回文数(palindrome)と呼びます。問題文の例では、2桁の数2つの積で作れる最大の回文数は 9009 = 91 × 99 です。同じように、3桁の数2つの積で作れる最大の回文数を求めるのが問題です。

2つの因数はそれぞれ100から999までの整数で、求めるのはその積のうち回文数であるものの最大値1つです。

前提知識

回文数の判定。 数 n の10進の各桁を逆順に並べた数を rev(n) と書くと、n が回文数であることは rev(n) = n と同じです。rev(n) は文字列に変換しなくても整数演算だけで作れます。n % 10 で最後の桁を取り出して結果の末尾に付け(結果 × 10 + 桁)、n を10で割った商に置き換える、という操作を n が0になるまで繰り返します。たとえば9009は 9 → 90 → 900 → 9009 となって回文数、9019は 9 → 91 → 910 → 9109 となって回文数ではありません。末尾が0の数(たとえば10)は、逆にすると先頭の0が消えて元の数と違う値になるので、回文数ではありません。先頭が0の数は存在しないので、これが正しい判定です。

積の単調性。 正の整数 a を固定すると、a × b は b について増加します。b を上から1ずつ減らすと、積は a ずつ減ります。ですから b を下げながら見ていき、ある積が基準値以下になったら、その後の積もすべて基準値以下です。同じ理由で、b の上限が hi なら a × b の最大値は a × hi で、a が小さくなるほどこの値も小さくなります。

偶数桁の回文数は11の倍数です。 4桁の回文数 abba と6桁の回文数 abccba を位ごとに展開すると

abba   = 1000a + 100b + 10b + a          = 1001a + 110b          = 11 × (91a + 10b)
abccba = 100001a + 10010b + 1100c                                = 11 × (9091a + 910b + 100c)

となります(1001 = 11 × 91、100001 = 11 × 9091、10010 = 11 × 910、1100 = 11 × 100)。11は素数なので、偶数桁の回文数を x × y と書けば、x と y の少なくとも一方は11の倍数です。最終的なアルゴリズムはこの事実を使いませんが、小さな例を手で検算するときと「発展」で使います。

素朴な方法

100から999までの a と b をすべて組み合わせて積を作り、回文数であるもののうち最大のものを選べば済みます。

best = 0
for a in range(100, 1000):
    for b in range(100, 1000):
        if str(a * b) == str(a * b)[::-1]:
            best = max(best, a * b)

3桁の数は900個なので積は 900 × 900 = 810000 個あり、それぞれに桁数に比例する回文判定を行います。今のコンピュータならすぐ終わるので、間違った方法ではありません。ただ、桁数を d に一般化すると積の個数は (9 × 10^(d−1))² で、d が1増えるごとに100倍になります。また、この方法は答えが出た後も残りの積を最後まで調べます。最大の回文数を探しているのに、大きな積がどこにあるかという、すでにわかっている情報を使っていないのです。

解法の導出

因数の範囲を lo = 10^(d−1)、hi = 10^d − 1 とします(3桁なら100と999)。素朴な方法の中から、省いてよい作業を一つずつ探していきます。

ステップ1: 順序を入れ替えた組は一度だけ見る。 a × b = b × a なので、a ≤ b の組だけを見れば十分です。積の個数がほぼ半分になります。

ステップ2: 大きい積から見る。 a を hi から lo まで下げ、各 a について b を hi から a まで下げます。前提知識の単調性により、1つの a の中(以下「a の行」)では積は減り続けます。ですから a の行で最初に見つかった回文数がその行で最大の回文数で、行の残りは見る必要がありません。

ステップ3: これまでの最大値と比べる。 これまでに見つかった最大の回文数を best とします。a の行を下りていく途中で積が best 以下になったら、その行の残りの積はすべて best 以下なので、行を打ち切ります。

ステップ4: 外側のループも打ち切る。 a の行で最大の積は a × hi です。a × hi ≤ best なら a の行に best を超える積はなく、a より小さい行の積はさらに小さいので、やはりありません。ここで探索全体を終えます。

小さな例で手計算による検算。 2桁(lo = 10、hi = 99)でたどってみます。99の行は 99 × 99 = 9801 の1つだけで、回文数ではありません。98の行から92の行までは、積が回文数になるには前提知識により因数の一方が11の倍数でなければなりませんが、92から99の間にある11の倍数は99だけです。したがって調べるべき積は b = 99 の場合だけです。

92 × 99 = 9108    93 × 99 = 9207    94 × 99 = 9306    95 × 99 = 9405
96 × 99 = 9504    97 × 99 = 9603    98 × 99 = 9702    99 × 99 = 9801

どれも回文数ではないので、99の行から92の行まで best は0のままです(アルゴリズムはこれらの行を最後まで下りて調べます)。91の行の最初の積 91 × 99 = 9009 は回文数なので、best = 9009 として行を打ち切ります。次の90の行では 90 × 99 = 8910 ≤ 9009 なので、ステップ4により探索が終わります。結果は 9009 = 91 × 99 で、問題文の例と一致します。

最終的なアルゴリズム

  1. lo = 10^(d−1)、hi = 10^d − 1、best = 0 から始めます。
  2. a = hi, hi − 1, …, lo の順に:
    1. a × hi ≤ best なら全体を打ち切ります。
    2. b = hi, hi − 1, …, a の順に積 p = a × b を作り、p ≤ best ならこの行を打ち切り、p が回文数なら best = p としてこの行を打ち切ります。
  3. best を返します。回文数が1つもなければ0ですが、この記事で使う d = 1, 2, 3 ではそうなりません。
  4. 答えを出力する前に小さな入力を確認します。回文判定は 9009 → 真、9019 → 偽、7 → 真、10 → 偽で、1桁の結果は9、2桁の結果は9009でなければなりません。1つでも違えば0以外の終了コードで終わります。

正しさの証明

a の行を処理する直前には、毎回次が成り立ちます(不変条件)。

best は、a < x ≤ y ≤ hi を満たすすべての組 (x, y) の積のうち回文数であるものの最大値である(なければ0)。

最初(a = hi)はそのような組がなく best = 0 なので成り立ちます。a の行を処理するとき、行を打ち切る場合は2通りです。積が best 以下になって打ち切った場合は、単調性により飛ばした積はすべて best 以下なので、最大値を変えることはできません。回文数に出会って打ち切った場合は、その積が行の中で最大の回文数で、飛ばした積はそれより小さくなります。行の中で best が変わるのは、best より大きな回文数に出会ったときだけです。どちらの場合も行を終えると best は x ≥ a の組までの最大値になり、次の行(a − 1)で不変条件が成り立ちます。

外側のループがステップ4で止まるとき、残りのすべての組 (x, y) は x ≤ a、y ≤ hi を満たすので x × y ≤ a × hi ≤ best です。残りの積は best を超えられません。ループが lo まで回りきって終わった場合は、すべての組を見たことになります。どちらの場合も best は lo ≤ x ≤ y ≤ hi を満たすすべての組の最大値で、ステップ1の対称性により、順序を問わないすべての組の最大値と一致します。

計算量

最悪の場合。 N = hi − lo + 1 = 9 × 10^(d−1) とすると、a ≤ b の組は N(N + 1) / 2 個で、回文判定1回は積の桁数(2d 以下)に比例します。枝刈りがまったく起きないと仮定すれば時間は O(N² · d) で、素朴な方法のおよそ半分です。3桁では N = 900 なので、最悪でも405450個の積です。保持する値は数個だけなので、追加の空間計算量は O(1) です。

実際にはずっと早く止まります。 大きな回文数が早く見つかれば、ステップ3と4によってほとんどの行は数回の積で終わるか、まったく開かれません。何回で終わるかは入力によって変わり、事前に保証できる式はないので、ここでは最悪の場合の上限だけを計算量として示します。実際の反復回数はこの記事では測っていません。

整数の範囲。

  • 3桁で最大の積は 999 × 999 = 998001 で、32ビット整数にも収まります。C++ の実装は d を1から9まで使えるように unsigned long long を使いました。その範囲で積は (10^9 − 1)² < 10^18 で、unsigned long long は規格で最低64ビット(最大値は約 1.8 × 10^19 以上)が保証されているので、あふれません。符号なし型なのであふれても未定義動作にはなりませんが、この範囲ではそもそもあふれが起きません。
  • 回文判定で作る逆順の数は、どの段階でも元の数より桁数が多くならないので、10^18 未満に収まります。ステップ4の条件の a × hi も同じ上限の中にあります。
  • ループ変数も符号なし型ですが、a >= lo(lo ≥ 1)と b >= a(a ≥ 1)で止まるので、0より下に行って巨大な値に回り込むことはありません。
  • Python の int は精度に制限がないので、あふれることはありません。アルゴリズムは C++ と同じです。

C++17 による実装

#include <iostream>

// True when the decimal digits of n read the same in both directions.
bool is_palindrome(unsigned long long n) {
    unsigned long long reversed = 0;
    for (unsigned long long rest = n; rest > 0; rest /= 10) {
        reversed = reversed * 10 + rest % 10;
    }
    return reversed == n;
}

// Largest palindrome that is a product of two numbers with exactly `digits` digits.
// Only pairs a <= b are tried. For a fixed a, b goes down from the top, so the
// products go down too: stop at the first palindrome or once a product is no
// larger than the best so far. Stop the outer loop once even a * hi cannot beat it.
// Meant for 1 <= digits <= 9, where every value fits in 64 bits.
unsigned long long largest_palindrome_product(int digits) {
    unsigned long long lo = 1;
    for (int i = 1; i < digits; ++i) {
        lo *= 10;
    }
    const unsigned long long hi = lo * 10 - 1;
    unsigned long long best = 0;
    for (unsigned long long a = hi; a >= lo; --a) {
        if (a * hi <= best) {
            break;
        }
        for (unsigned long long b = hi; b >= a; --b) {
            const unsigned long long product = a * b;
            if (product <= best) {
                break;
            }
            if (is_palindrome(product)) {
                best = product;
                break;
            }
        }
    }
    return best;
}

int main() {
    if (!is_palindrome(9009) || is_palindrome(9019) || !is_palindrome(7) ||
        is_palindrome(10) || largest_palindrome_product(1) != 9 ||
        largest_palindrome_product(2) != 9009) {
        return 1;
    }
    std::cout << largest_palindrome_product(3) << '\n';
    return 0;
}

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

Python 3 による実装

def is_palindrome(n: int) -> bool:
    """True when the decimal digits of n read the same in both directions."""
    reversed_n = 0
    rest = n
    while rest > 0:
        reversed_n = reversed_n * 10 + rest % 10
        rest //= 10
    return reversed_n == n


def largest_palindrome_product(digits: int) -> int:
    """Largest palindrome that is a product of two numbers with exactly `digits` digits.

    Only pairs a <= b are tried. For a fixed a, b goes down from the top, so the
    products go down too: stop at the first palindrome or once a product is no
    larger than the best so far. Stop the outer loop once even a * hi cannot beat it.
    """
    lo = 10 ** (digits - 1)
    hi = lo * 10 - 1
    best = 0
    for a in range(hi, lo - 1, -1):
        if a * hi <= best:
            break
        for b in range(hi, a - 1, -1):
            product = a * b
            if product <= best:
                break
            if is_palindrome(product):
                best = product
                break
    return best


def main() -> None:
    if (
        not is_palindrome(9009)
        or is_palindrome(9019)
        or not is_palindrome(7)
        or is_palindrome(10)
        or largest_palindrome_product(1) != 9
        or largest_palindrome_product(2) != 9009
    ):
        raise SystemExit(1)
    print(largest_palindrome_product(3))


if __name__ == "__main__":
    main()

C++ と同じアルゴリズム、同じ小さな例です。回文判定も C++ に合わせて str(n)[::-1] ではなく整数演算で行い、桁を取り除くときは整数除算 //= を使います。range は終端の値を含まないので、range(hi, lo - 1, -1) のように終端を1つ下げています。小さな例が失敗すれば raise SystemExit(1) で終了コード1を返します。assert は python -O で実行すると取り除かれるので使っていません。

答えと検証

3桁の数2つの積で作れる最大の回文数は 906609 = 913 × 993 です。以下ではこの値を単に「結果」と呼びます。

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

  • 積。 913 × 993 = 913 × 1000 − 913 × 7 = 913000 − 6391 で、これがちょうど結果に一致します。913と993はどちらも3桁です。
  • 回文数。 結果は6桁で、上3桁を逆から読むと下3桁になります。
  • 11との関係。 結果は 11 × 82419 で、913 = 11 × 83 なので、因数の一方が11の倍数という前提知識の事実と合っています。さらに分解すると結果は 3 × 11 × 83 × 331 で、993 = 3 × 331 です。

この結果が「最大」であることは手計算ではなく、上の正しさの証明とプログラムの実行で保証します。検証は三重です。

  1. 小さな例。 2つのプログラムはどちらも答えを出力する前に、回文判定4つ(9009、9019、7、10)と1桁の結果9、2桁の結果9009を確認し、違えば0以外の終了コードで終わります。2桁の結果9009は「解法の導出」で手でもたどりました。
  2. 2つの実装の一致。 このシリーズのリポジトリでは、検査スクリプトが C++ のコードを g++ -std=c++17 -O2 -Wall -Wextra -Werror でコンパイルして実行し(警告も失敗扱い)、Python のコードを python3 -I -B で実行します。2つの出力が保存してある answer.txt とバイト単位で一致しなければならず、0以外の終了コードやタイムアウト(C++ 1秒、Python 10秒)は失敗です。シリーズの規則上、この検査を通ったコードだけを記事に載せます。
  3. 証明。 枝刈りで飛ばした積が best を超えられないことは「正しさの証明」で示しました。疑わしければ、「素朴な方法」のコードをそのまま実行して結果を比べてみることもできます。

境界のケース

  • 1桁(d = 1): 1から9までの積のうち1桁のものはすべて回文数で、2桁の回文数 11, 22, …, 99 はすべて11の倍数なので、1桁の数2つの積にはなりません。したがって答えは 9(= 9 × 1 = 3 × 3)です。アルゴリズムは9の行から4の行までで回文数を見つけられず、3の行の最後の積 3 × 3 で9を見つけます。2の行は 2 × 9 = 18 が9より大きいので開かれますが、積が9以下になったところで打ち切られ、1の行は 1 × 9 = 9 ≤ 9 なのでステップ4で探索が終わります。
  • 末尾が0の積: 10や9010のように0で終わる数は回文数ではありません。整数で逆順にすると先頭の0が消えるので、判定は自然に偽になります。小さな例に 10 → 偽 が入っています。
  • 1桁の数の判定: 7のような1桁の数はそれ自体が回文数です。ループが1回回って 7 = 7 になります。
  • 行の中の積が best と等しいとき: ステップ3の条件は <= です。等しい値は best を変えないので、飛ばしてかまいません。

発展

  • 11の倍数だけを試す。 3桁の数2つの積は10000以上なので、その中の回文数は5桁か6桁です。答えが6桁の回文数だとわかれば、a か b の一方は11の倍数でなければなりません。a が11の倍数でないときは、b を11の倍数だけ下りるように変えてみてください。正しさの証明のどこを直す必要があるでしょうか。5桁の回文数は11の倍数とは限らない点もあわせて考える必要があります。
  • 回文数から作る。 逆に、999999、998899、997799、… と6桁の回文数を大きい方から作り、100以上999以下の約数 x を持ち、商も3桁になるかを確かめる方法もあります。上3桁を決めれば回文数が決まることを使います。どちらの方が試す数が少ないか考えてみてください。
  • 桁数を増やす。 d = 4, 5, … と増やすと、枝刈りの効果はどう変わるでしょうか。上の C++ のコードは d ≤ 9 で64ビットに収まりますが、そのときも最悪の上限は O(N² · d) です。
  • ほかの進法。 is_palindrome は10進の桁しか見ていません。進法を引数で受け取るにはどう変えればよいでしょうか。0で終わる数は回文数ではないという議論は、そのまま成り立つでしょうか。

出典とライセンス

  • Project Euler, Problem 4 — 元の問題。この記事の「問題」の節はその内容を独自に要約したもので、原文と同じ CC BY-NC-SA 4.0 で共有します。問題の著作者表示は Project Euler です。
  • Project Euler, Copyright — 問題文のライセンスとリンクの方法。非営利の条件に従い、この記事には広告を載せません(ads: false)。
  • 解法、証明、コードは自分で書いたものです。Project Euler の解説(overview)、フォーラム、会員限定の資料や他の解答サイトは参照していません。Motitory は Project Euler と提携・後援の関係になく、Project Euler がこの記事を確認したこともありません。

ブログへ