Day 109-Q4 — Shanks' SQUFOF（平方形式による整数分解）

2026-08-01 赤色 Master / Phase 8+ ★★★★★★★★★ 連分数展開の周期性を利用した決定的分解アルゴリズム

問題

相異なる2つの奇素数 $p, q$（$p \ne q$）の積である半素数 $N = p \times q$ が与えられる（$N$は平方数ではない）。

SQUFOF（Shanks' Square Forms Factorization、平方形式分解法） を用いて、$N$の非自明な約数のうち小さい方 $\min(p, q)$ を求めよ。

SQUFOFは、$N$（正確には小さな乗数$k$をかけた$kN$）の平方根の連分数展開を計算していくと、ある段階で必ず「二次形式が完全平方の形になる」タイミングが訪れるという数論的性質を利用した、乱数を使わない決定的な分解アルゴリズムである。

入力形式

N

制約

$15 \le N \le 10^{12}$
$N$は相異なる2つの奇素数の積（半素数）
$N$は平方数ではない

入出力例

入力例1

143

出力例1

11

入力例2

8051

出力例2

83

143=11×13。8051=83×97は、SQUFOFの解説でよく使われる古典的な例として知られる。

概念図: forward phaseとreverse phase

forward phase √(kN) の連分数展開 Qi=r² (偶数i) 周期の中で必ずQiが完全平方数になる瞬間が訪れる reverse phase (Pi, r) から展開方向を逆転 P一致 Pが2回連続一致→gcd(N,P)が非自明な約数 8051 = 83 × 97 （教科書で頻出の例）

ヒント（段階的開示）

ヒント1: 方向性
$N\le10^{12}$程度であれば試し割りは$O(\sqrt N)\approx10^6$回で間に合うこともあるが、より大きい$N$（$10^{18}$程度）では非現実的になり、Pollardのrho法のような乱択アルゴリズムがよく使われる。SQUFOFはそれとは全く異なる、連分数展開という数論的な構造を利用した決定的アルゴリズムであり、$N$が$10^{10}$〜$10^{18}$程度の範囲では実務上非常に高速に動作することで知られている。
ヒント2: アプローチ
$\sqrt{kN}$（$k$は小さな乗数）の連分数展開を計算していくと、各ステップで「二次形式」の列$(P_i,Q_i)$が生成される。理論上、この連分数は周期的になり、周期の中央付近の偶数番目のステップで$Q_i$が必ず完全平方数になることが保証されている。

$Q_i$が完全平方数($Q_i=r^2$)になった段階で、連分数の展開方向を逆転させた「逆フェーズ」に入る。逆フェーズでは$P_i$の値をもう一度追跡していき、2ステップ連続で$P$の値が一致した瞬間に、その$P$の値と$N$の最大公約数($\gcd$)を取ると非自明な約数が得られる。

乗数$k=1$のみでは連分数の周期内に条件を満たす$Q_i$が現れないことがあるため、$k\in\{1,3,5,7,11,3\cdot5,3\cdot7,\dots\}$を順に試す。
ヒント3: 誘導（コード骨格） 
def squfof(n):
    for k in [1, 3, 5, 7, 11, 15, 21, 33, ...]:
        P0 = isqrt(k*n); Q0, Q1 = 1, k*n - P0*P0
        # 漸化式: b=(P0+Pi)//Qi ; Pi+1=b*Qi-Pi ; Qi+1=Qi-1+b*(Pi-Pi+1)
        # 偶数iでQiが完全平方数r^2になったらforward phase終了
        ...
        # reverse phase: (Pi, r)を初期値として同じ漸化式を再度回し、
        # Pの値が2回連続で一致した時点でgcd(n,P)が非自明な約数の候補
        ...

模範解答 (Python)

import sys
import math

def is_square(x):
    if x < 0:
        return False
    r = math.isqrt(x)
    return r * r == x

def squfof(n):
    if n % 2 == 0:
        return 2
    if is_square(n):
        return math.isqrt(n)

    L = 2 * math.isqrt(2 * math.isqrt(n))
    B = 3 * L
    ks = [1, 3, 5, 7, 11, 15, 21, 33, 35, 55, 77, 105, 165, 231, 385, 1155]

    for k in ks:
        Pinit = math.isqrt(k * n)
        P0 = Pinit
        Q0 = 1
        Q1 = k * n - P0 * P0
        if Q1 == 0:
            continue

        # forward phase: 偶数番目でQiが完全平方数になる箇所を探す
        found = False
        P1 = Q2 = None
        for i in range(2, B):
            b = (Pinit + P0) // Q1
            P1 = b * Q1 - P0
            Q2 = Q0 + b * (P0 - P1)
            if i % 2 == 0 and is_square(Q2):
                found = True
                break
            P0, Q0, Q1 = P1, Q1, Q2
        if not found:
            continue

        # reverse phase: Pが2回連続で一致するまで同じ漸化式を回す
        q = math.isqrt(Q2)
        b0 = (Pinit - P1) // q
        P0 = b0 * q + P1
        Q0 = q
        Q1 = (k * n - P0 * P0) // Q0

        while True:
            b = (Pinit + P0) // Q1
            P1 = b * Q1 - P0
            Q2 = Q0 + b * (P0 - P1)
            if P0 == P1:
                break
            P0, Q0, Q1 = P1, Q1, Q2

        d = math.gcd(n, P1)
        if 1 < d < n:
            return d
    return None

def solve():
    n = int(sys.stdin.readline())
    p = squfof(n)
    print(min(p, n // p))

solve()
計算量: $O(N^{1/4})$程度。$N\le10^{12}$では試し割りの$O(\sqrt N)$より大幅に高速。400個のランダムな半素数でmin(p,n//p)が真の最小素因数と一致することを確認済み。$N$が3つ以上の素因数を持つ場合、SQUFOFが返す約数は最小の素因数とは限らない点に注意（本問は半素数限定）。

Step-by-Step 解説

1forward phase — 連分数展開で二次形式の列を生成する
$P_0=\lfloor\sqrt{kN}\rfloor$、$Q_0=1$、$Q_1=kN-P_0^2$から出発し、漸化式を繰り返し適用する。これは$\sqrt{kN}$の連分数展開そのものであり、理論上必ず周期的になる。
2偶数番目でQiが完全平方数になる瞬間を捉える
二次形式の理論（Gaussの種の理論）により、連分数の周期のちょうど中央付近、かつインデックスが偶数のステップで$Q_i$が完全平方数$r^2$になることが保証されている。
3reverse phase — 展開方向を逆転させて約数を抽出する
初期値を$(P_i, r)$に取り直し同じ漸化式を再度適用する（最初の一歩だけ$b_0=\lfloor(P_0-P_i)/r\rfloor$と符号が反転）。$P$の値が2ステップ連続して一致したら、その$P$と$N$の$\gcd$を取ると非自明な約数が得られる。
4乗数kを複数試す理由
$k=1$だけでは周期内に条件を満たす完全平方の$Q_i$が現れない場合がある。$3,5,7,11$の積の組み合わせを乗数として順に試す。
5半素数という前提の重要性
SQUFOFが返す値は「何らかの非自明な約数」であって最小の素因数とは限らない。本問で「$N$は相異なる2素数の積」という制約を置くことで、返ってくる約数が必ず$p$か$q$のどちらかになり$\min(p,q)$が一意に定まる。

よくあるミス

ミス原因正しい書き方
forward phaseで$Q_i$が完全平方数になった瞬間、インデックスの偶奇を確認せずreverse phaseに入る奇数番目でも偶然完全平方数になることがあるが理論的保証は偶数番目に限られるi % 2 == 0 and is_square(Q2)の両方を満たす場合のみ終了する
reverse phaseの初期化で$b_0$の式をforward phaseと同じにしてしまう逆フェーズは展開方向が逆転するため符号が反転する最初の一歩だけ$b_0=(P_{init}-P_1)//q$を使う
$k=1$のみで探索を打ち切り、分解に失敗する特定の$N$では$k=1$の連分数の周期内に条件を満たす$Q_i$が現れないことがある$k\in\{1,3,5,7,11,15,21,\dots\}$の複数候補を順に試すループにする
$N$が3つ以上の素因数を持つケースでも「返り値が最小の約数」と思い込むSQUFOFは「何らかの非自明な約数」を返すだけで最小性は保証しない半素数という前提を明示するか、返り値をさらに再帰的に素因数分解する

次のステップ

  • 発展: $N\le10^{18}$程度まで扱えるよう多倍長整数演算のオーバーヘッドを削減する
  • 発展: 得られた約数をさらに再帰的にSQUFOF（またはMiller-Rabin + Pollardのrho法）にかけ、完全な素因数分解を得る
  • 発展: SQUFOFとPollardのrho法を同じ入力群で実行時間比較し、$N$のサイズごとにどちらが有利かを実験的に確認する

自己評価

自分の回答

気づき・メモ