ペル方程式の解の列挙方法

Project Euler 66 のネタバレです。見たくない人は見ないでください。本質的には全然理解できてないですが、それなりに有用な情報だと思ったので、解き方だけメモします。

ペル方程式とは

$x^2-Dy^2=\pm1$ (ただし D は平方数でない自然数)
という形の不定方程式をペル方程式というそうです。

これを満たす整数 x と y は無数にあります。ですが、D の値によっては最小解でもかなり大きい値になることがあり、1 から順番に探していくことは事実上不可能です。たとえば D = 166 のときは以下が最小解です。

p(1700902565**2 - 166 * 132015642**2)  #=> 1

最小解の見つけ方

最小解を高速に探し出す方法があります (参考: 二次無理数の連分数展開とペル方程式の解の構成) 。D の平方根の連分数表示を使って、漸化式で解けるとのこと。

例えば、14 の平方根の連分数表示は以下のようになります。

$\sqrt{14} = 3 + \frac{1}{1 + \frac{1}{2 + \frac{1}{1 + \frac{1}{6 + \cdots}}}}$
1, 2, 1, 6 の後は循環して、1, 2, 1, 6, 1, 2, 1, 6, ... が続きます。循環する二歩手前までの数列 [3, 1, 2, 1] を q(n) として (最後の 6 は含まない) 、以下の漸化式で表される X(n) と Y(n) の最後の組が、なんと最小の解になります。

X(0) = 1, X(1) = q(0), X(n+1) = q(n) * X(n) + X(n-1)
Y(0) = 0, Y(1) = 0,    Y(n+1) = q(n) * Y(n) + Y(n-1)

計算するとこうなります。

X(2) = 1 * X(1) + X(0) =  4,  Y(2) = 1 * Y(1) + Y(0) = 1
X(3) = 2 * X(2) + X(1) = 11,  Y(3) = 2 * Y(2) + Y(1) = 3
X(4) = 1 * X(3) + X(2) = 15,  Y(4) = 1 * Y(3) + Y(2) = 4

よって、x = 15, y = 4 の時に式が満たされることになります。実際、満たされます。

p(15**2 - 14 * 4**2)  #=> 1

なんでこれで最小解が得られるのか全然わかってないんですが、Tossy-2 に直感的に教えてくれと言ったところ、 $x^2 - Dy^2 = 1$ は $x^2 \approx D y^2$ だから $x \approx \sqrt{D} y$ で、連分数展開して一番近くなったところが解になる、だそうです。なるほどー。ありがとう。

解の列挙

最小の解を $(x_1, y_1)$ とすると、以下を満たす $(x_n, y_n)$ は全て解になり、かつ、解はこれだけだそうです。へー。

$x_n + \sqrt{D}y_n = (x_1 + \sqrt{D}y_1)^n$

プログラム

以上をプログラムにしました。実行例。

n = 14
PellEquation.each_solution(n, true) do |x, y, v|
  puts "#{x}^2 - #{n} * #{y}^2 = #{v}"
end

15^2 - 14 * 4^2 = 1
449^2 - 14 * 120^2 = 1
13455^2 - 14 * 3596^2 = 1
403201^2 - 14 * 107760^2 = 1
12082575^2 - 14 * 3229204^2 = 1
362074049^2 - 14 * 96768360^2 = 1
...

以下コード。

module PellEquation
  def self.each_solution(n, allow_minus_one = true)
    unless block_given?
      return to_enum(:each_solution, n, allow_minus_one)
    end

    return if n <= 1

    # find a = n^(1/2) by newton method
    a, b = n, 1
    a, b = a / 2, b * 2 until a <= b
    a = b + 1
    a, b = b, (b + n / b) / 2 until a <= b

    # no solution when n is square
    return if a * a == n

    # find fundamental solution
    a0, b, f = a, 1, false
    x2, x = 1, a
    y2, y = 0, 1
    loop do
      # continued fraction expansion
      b = (n - a * a) / b
      q = (a0 + a) / b
      a = q * b - a
      f = !f
      break if b == 1

      # recurrence equation
      x2, x = x, x * q + x2
      y2, y = y, y * q + y2
    end

    # enumerate solutions
    x0, y0 = x, y
    (f ? [-1, 1] : [1]).cycle do |v|
      yield(x, y, v) if allow_minus_one || v == 1
      x, y = x0 * x + n * y0 * y, x * y0 + x0 * y
    end
  end
end