Day 101-Q4 — Lehmer's GCD Algorithm（多倍長整数の高速ユークリッド互除法）

2026-07-24 赤色 Master / Phase 8+ ★★★★★★★★★ 先頭桁近似 + ハイブリッド高速化

問題

2つの巨大な非負整数 $A \ge B > 0$（最大300桁程度）の最大公約数 $\gcd(A,B)$ を求めよ。

普通のユークリッド互除法は正しいが、桁数が非常に大きいとき各ステップの`a % b`が低速になる。Lehmer's GCD Algorithmは、$A,B$の先頭の数十〜数百ビットだけを取り出した近似値$A_h,B_h$（機械語1語サイズ）を使って複数回分の商をまとめて計算し、$2\times2$整数行列として蓄積したうえで最後に1回だけ元の巨大な$A,B$に一括適用することで高速化する。

アルゴリズム概要（各反復）:

  1. $B$のビット長が64ビット以下なら通常のユークリッド互除法で仕上げて終了。
  2. そうでなければ$A$の先頭64ビット相当を$A_h$、同じシフト量の$B$を$B_h$とする。
  3. 真の比$A/B$を挟み込む境界ペア$P=(A_h,B_h+1)$（比が最小）と$Q=(A_h+1,B_h)$（比が最大）で並行に互除法を進め、両方が同じ商を返す限り$2\times2$行列$(x_0,x_1,y_0,y_1)$に蓄積。商が食い違ったら打ち切り。
  4. 1つも安全に確定できなければ通常の互除法を1ステップだけ元の$A,B$に適用。そうでなければ蓄積した変換を一括適用。
  5. $B=0$になったら$A$が答え。

入力形式

A B

制約

$1 \le B \le A < 10^{300}$
10進文字列として入力

入出力例

入力例1

270 192

出力例1

6

入力例2

10946 6765

出力例2

1

10946, 6765は連続するフィボナッチ数。連続フィボナッチ数は常に互いに素であり、ユークリッド互除法が最悪ケースの反復回数になる有名な入力例でもある。

概念図

先頭ビットだけで複数ステップ分の商を先取り 先頭64bit=A_h 残り数百桁（この部分は見ない） a B_h b A_h, B_h だけでユークリッド互除法を進め、商が「下位ビット に依存せず確定する」限り 2x2行列 (x0,x1,y0,y1) に蓄積 → 蓄積した変換を a,b に一括適用（複数ステップを1回の巨大整数演算に）

ヒント（段階的開示）

ヒント1: 方向性
Pythonの整数は多倍長対応だから`math.gcd`を呼べば済むが、この問題の目的はアルゴリズム自体を実装して理解すること。実世界の多倍長演算ライブラリ（GMPなど）が、なぜ単純なユークリッド互除法だけでは不十分で、こうした先頭桁の近似計算による高速化が必要になるのかを体感するのが狙い。
ヒント2: アプローチ
「$A,B$の下位ビットを見なくても、先頭の数ビットだけで商$q=\lfloor A/B\rfloor$が確定できる場合が多い」という点がポイント。ただし先頭ビットだけでは下位ビット次第で真の商とズレる危険がある。「真の$a/b$が取りうる最小の比と最大の比、その両極端で互除法を進めても同じ商になるなら、実際の下位ビットが何であっても安全にその商を採用できる」という判定がLehmer法の核心的な安全性チェック。
ヒント3: 誘導（コード骨格） 
P0, P1 = A, B + 1        # 比が最小になる側の境界ペア
Q0, Q1 = A + 1, B        # 比が最大になる側の境界ペア
x0, x1, y0, y1 = 1, 0, 0, 1

while P1 != 0 and Q1 != 0:
    q1 = P0 // P1
    q2 = Q0 // Q1
    if q1 != q2:
        break
    q = q1
    P0, P1 = P1, P0 - q * P1
    Q0, Q1 = Q1, Q0 - q * Q1
    x0, x1, y0, y1 = y0, y1, x0 - q * y0, x1 - q * y1

行列の更新は$x$側と$y$側を入れ替えながら行う（`x0,x1 = y0,y1`のように）点に要注意。単純に$x$は$x$だけ、$y$は$y$だけで独立更新すると、複数ステップ後に元の$a,b$への正しい変換行列にならない（Lehmer法で最もバグりやすい箇所）。

模範解答 (Python)

import sys


def lehmer_gcd(a, b):
    WORD_BITS = 64
    while b:
        if b.bit_length() <= WORD_BITS:
            while b:
                a, b = b, a % b
            return a

        shift = max(0, a.bit_length() - WORD_BITS)
        A = a >> shift
        B = b >> shift

        P0, P1 = A, B + 1
        Q0, Q1 = A + 1, B
        x0, x1, y0, y1 = 1, 0, 0, 1

        while P1 != 0 and Q1 != 0:
            q1 = P0 // P1
            q2 = Q0 // Q1
            if q1 != q2:
                break
            q = q1
            P0, P1 = P1, P0 - q * P1
            Q0, Q1 = Q1, Q0 - q * Q1
            x0, x1, y0, y1 = y0, y1, x0 - q * y0, x1 - q * y1

        if x1 == 0:
            a, b = b, a % b
        else:
            a, b = x0 * a + x1 * b, y0 * a + y1 * b

    return a


def solve():
    A, B = map(int, sys.stdin.read().split())
    print(lehmer_gcd(A, B))


solve()
計算量: 巨大整数演算の回数が $O(\log(\max(A,B)))$ から実用上大きく削減される（機械語サイズの近似計算で複数ステップをまとめるため）。

Step-by-Step 解説

1小さくなったら通常の互除法へ
`b.bit_length()<=WORD_BITS`なら普通のユークリッド互除法で仕上げ、最終的な正しさを保証。
2先頭64ビットの抽出
`a.bit_length()-WORD_BITS`だけ右シフトして$A$を得、$b$も同じシフト量で$B$を得る。
3境界ペアを使った安全性チェック
真の比$a/b$を挟み込む2つの境界ペア$P=(A,B+1)$,$Q=(A+1,B)$の商が一致する限り、その商は真の$a,b$に対しても正しい。
4蓄積した変換の一括適用
複数回分のステップを1回の巨大整数演算にまとめるのが高速化の本質。
5フォールバック
`x1==0`（安全な商が1つも見つからない）場合は通常の互除法を1ステップだけ適用。

よくあるミス

ミス原因正しい書き方
行列の更新で$x$側・$y$側を独立に更新してしまう「行と行を入れ替えながら更新する」必要性に気づかず`x0,x1=x1,x0-q*x1`のように自己完結で更新する`x0,x1,y0,y1 = y0,y1, x0-q*y0, x1-q*y1`のように新しい$x$側は古い$y$側から計算する
A,Bの抽出シフト量を別々に計算位取りがずれて商の近似が破綻`shift=a.bit_length()-WORD_BITS`という1つの値を両方に使う
境界ペアの片方（Pのみ）しか見ない商が下位ビット次第で変わるケースを見逃すP側・Q側両方の境界から計算した商が完全一致することを確認
閾値を小さくしすぎて高速化の恩恵が出ないガード桁が不足実用的には32〜64ビット程度の閾値を確保

次のステップ

  • 発展: 100桁・300桁のランダムな大整数ペアで`math.gcd`との一致を大量検証してから、素朴な互除法とLehmer法の`a%b`呼び出し回数を比較する
  • 発展: 拡張ユークリッド互除法版（$ax+by=\gcd(a,b)$）にLehmer法を組み込む
  • 次回予告: 3-NTT素数によるCRT合成（Garnerのアルゴリズム）

自己評価

自分の回答

気づき・メモ