Day 153-Q1 | ペル方程式（連分数展開による基本解の構成）

問題

平方数でない正整数 $D$ が与えられる。ペル方程式

$x^2 - D y^2 = 1$

を満たす正整数の組 $(x, y)$ のうち、$x$ が最小となるもの（基本解）を求めよ。$\sqrt{D}$ の連分数展開 $[a_0; a_1, a_2, \dots]$ は必ず周期的になり、その収束分数（convergent）$h_n/k_n$ の中に基本解が必ず現れることが知られている。この性質を使って構成的に解を求めること（$D$ を大きくしても解が指数的に巨大になりうるため、単純な全探索は使えない）。

入力形式

制約

$2 \le D \le 10^6$

$D$ は平方数ではない

入出力例

入力例1

出力例1

3 2

入力例2

出力例2

1766319049 226153980

概念図：連分数の収束分数が解を生む過程

ヒント（段階的開示）

ヒント1（方向性）

$x, y$ を小さい方から総当たりで探すと、$D=61$ のように基本解が $10^9$ 級になる入力では絶対に終わらない。全探索ではなく、$\sqrt{D}$ の連分数展開そのものを計算し、その収束分数を使う必要がある。

ヒント2（アプローチ）

$\sqrt{D}$ の連分数展開は $m_0=0, d_0=1, a_0=\lfloor\sqrt{D}\rfloor$ から始め、$m_{n+1}=d_na_n-m_n$, $d_{n+1}=(D-m_{n+1}^2)/d_n$, $a_{n+1}=\lfloor(a_0+m_{n+1})/d_{n+1}\rfloor$ の漸化式で計算できる。収束分数 $h_n/k_n$ も $h_n=a_nh_{n-1}+h_{n-2}$（$k_n$も同様）で更新でき、$h_n^2-Dk_n^2=1$ となった時点の $(h_n,k_n)$ が基本解。

ヒント3（誘導）

import math

def solve_pell(D):
    a0 = math.isqrt(D)
    m, d, a = 0, 1, a0
    h_prev2, h_prev1 = 1, a0
    k_prev2, k_prev1 = 0, 1
    if h_prev1 * h_prev1 - D * k_prev1 * k_prev1 == 1:
        return h_prev1, k_prev1
    while True:
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d
        h = a * h_prev1 + h_prev2
        k = a * k_prev1 + k_prev2
        if h * h - D * k * k == 1:
            return h, k
        h_prev2, h_prev1 = h_prev1, h
        k_prev2, k_prev1 = k_prev1, k

模範解答 (Python)

import sys
import math


def solve_pell(D):
    a0 = math.isqrt(D)
    m, d, a = 0, 1, a0
    h_prev2, h_prev1 = 1, a0
    k_prev2, k_prev1 = 0, 1

    if h_prev1 * h_prev1 - D * k_prev1 * k_prev1 == 1:
        return h_prev1, k_prev1

    while True:
        m = d * a - m
        d = (D - m * m) // d
        a = (a0 + m) // d
        h = a * h_prev1 + h_prev2
        k = a * k_prev1 + k_prev2
        if h * h - D * k * k == 1:
            return h, k
        h_prev2, h_prev1 = h_prev1, h
        k_prev2, k_prev1 = k_prev1, k


def main():
    D = int(sys.stdin.readline())
    x, y = solve_pell(D)
    print(x, y)


main()

Step-by-Step 解説

1連分数展開の初期値
$D=2$: $a_0=\lfloor\sqrt2\rfloor=1$。$h_{-1}=1,h_0=1,k_{-1}=0,k_0=1$。$h_0^2-2k_0^2=1-2=-1$ で不成立。

21回反復して収束分数を更新
$m_1=1,d_1=1,a_1=2$。$h_1=2\cdot1+1=3$、$k_1=2\cdot1+0=2$。

3判定して確定
$h_1^2-2k_1^2=9-8=1$ で基本解 $(3,2)$ が確定、入力例1と一致。

4D=61 の「暴れる」ケース
周期長11で基本解が $x=1766319049, y=226153980$ という巨大な数になる。$D$が小さくても解が巨大になりうる古典的反例で、連分数法の必要性を裏付ける。

よくあるミス

ミス	原因	正しい書き方
$(x,y)$ を小さい方から線形探索する	$D=61$型の入力で終わらない	連分数展開の収束分数だけを追う
$n=0$の特殊ケースの判定を忘れる	$a_0$自体が基本解を与える場合がある	ループ前に $n=0$ もチェックする
整数除算の代わりに浮動小数点演算を使う	大きい$D$で誤差が出る	すべて `//`, `math.isqrt` で統一
「1周＝解」と誤解する	停止条件の誤解	各収束分数ごとに $h^2-Dk^2=1$ を直接チェック
出力を`float`にキャストする	大きい解で精度が壊れる	Pythonのint（多倍長整数）のまま出力

次のステップ

発展: $x^2-Dy^2=-1$（負のペル方程式）の解の存在条件（連分数周期の偶奇）を調べ、存在しない$D$の例を確認する。
次回予告: Suffix Array構築（Manber-Myersのダブリング法）

自己評価

完全理解おおむね理解要復習

Day 153-Q1 — ペル方程式（連分数展開による基本解の構成）

問題