Python ペル方程式の最小解を求める

ペル方程式とは

$$x^2 - D y^2 = 1$$ はペル方程式と呼ばれており, $\sqrt{D}$ の連分数展開によって最小解を求めることができる。

$$\sqrt{D} = a_0 + \cfrac{1}
{a_1 + \cfrac{1}
{a_2 + \cfrac{1}{\ddots}}
}$$

こちらが非常にわかりやすいです。

こちらの記事はわかりやすく解説されており、連分数を求めるアプリも開発されています。

実装

こちらに置きました。

import math
from decimal import (
    Decimal,
    getcontext
)

def is_square(n):
    is_square = False
    for i in range(2, int(math.sqrt(n)) + 1):
        if n == i ** 2:
            is_square = True
            break
    return is_square

def get_sqrt(D):
    getcontext().prec = 100
    x = int(D)
    for i in range(20):
        x = Decimal(str(x - (x ** 2 - D) / (2 * x)))
    return x

def get_continued_fraction(D):
    continued_fraction = []
    if not is_square(D):
        y = get_sqrt(D) - int(get_sqrt(D))
        continued_fraction.append(int(1 / y))
        z = 1 / y - int(1 / y)
        while True:
            if abs(y - z) < 1e-5:
                break
            else:
                continued_fraction.append(int(1 / z))
                z = 1 / z - int(1 / z)
    return continued_fraction

def get_minimal_solution(D):
    continued_fraction = get_continued_fraction(D)
    length = len(continued_fraction)
    num = length % 2 + 1
    continued_fraction = continued_fraction * num
    an2, bn2 = 0, 0
    an1, bn1 = 0, 0
    an, bn = 1, 0 

    for n in range (1, num * length + 1):
        if n == 1:
            an1, bn1 = int(math.sqrt(D)), 1
        else:
            an2, bn2 = (an + continued_fraction[n - 2] * an1,
                        bn + continued_fraction[n - 2] * bn1)
            an, bn = an1, bn1
            an1, bn1 = an2, bn2
    print_continued_fraction = continued_fraction[:int(len(continued_fraction) / num)]
    print_continued_fraction.insert(0, int(get_sqrt(D)))
    print('連分数展開: {}'.format(print_continued_fraction))
    return an2, bn2

if __name__ == '__main__':
    d = 2
    x, y = get_minimal_solution(d)
    print('最小解\n  x: {}\n  y: {}'.format(x, y))

出力例

$D=2$

連分数展開: [1, 2]
最小解
  x: 3
  y: 2
処理時間: 0.0020055770874023438

$D=7$

連分数展開: [2, 1, 1, 1, 4]
最小解
  x: 8
  y: 3
処理時間: 0.0009946823120117188

$D=14$

連分数展開: [3, 1, 2, 1, 6]
最小解
  x: 15
  y: 4
処理時間: 0.0010349750518798828

$D=61$

連分数展開: [7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14]
最小解
  x: 1766319049
  y: 226153980
処理時間: 0.0010004043579101562

$D=751$

連分数展開: [27, 2, 2, 8, 1, 2, 1, 3, 5, 1, 4, 1, 1, 1, 3, 1, 1, 3, 10, 1, 2, 7, 2, 17, 1, 4, 27, 4, 1, 17, 2, 7, 2, 1, 10, 3, 1, 1, 3, 1, 1, 1, 4, 1, 5, 3, 1, 2, 1, 8, 2, 2, 54]
最小解
  x: 7293318466794882424418960
  y: 266136970677206024456793
処理時間: 0.0009648799896240234

$D=991$

連分数展開: [31, 2, 12, 10, 2, 2, 2, 1, 1, 2, 6, 1, 1, 1, 1, 3, 1, 8, 4, 1, 2, 1, 2, 3, 1, 4, 1, 20, 6, 4, 31, 4, 6, 20, 1, 4, 1, 3, 2, 1, 2, 1, 4, 8, 1, 3, 1, 1, 1, 1, 6, 2, 1, 1, 2, 2, 2, 10, 12, 2, 62]
最小解
  x: 379516400906811930638014896080
  y: 12055735790331359447442538767
処理時間: 0.002000093460083008

参考記事