佩尔方程(pell方程)

第一类佩尔方程

形如： ,(d>1,且d不是完全平方数) , 要求第一类佩尔方程的解都是正整数解,即( x , y ) , x>0 , y>0。

​ 方程存在无数个解，要求最小整数解的话，只要枚举y的值直到x为整数即可。

想要得到所有解，我们就要推导通解公式，过程如下。 首先，设为方程的两个特解两式相乘：展开一顿化简可以得到：所以得到通项公式为

这样的话我们只要知道第一个特解( x1, y1)，和第n-1个特解,就能够求出第n个解,我们这时候可以用矩阵快速幂,以O(log n)的复杂度来就快速求取，第n个解。

第二类佩尔方程

形如： 的方程叫做第二类佩尔方程。

解第二类佩恩方程，需要用到第一类方程的解。过程与第一类大同小异。

我们假设有的一个特解又的正整数解为两式相乘：推导得到：所以通解为 可以转化为矩阵乘法的形式 $ = ^{n-1}

$ ,这样就方便用矩阵快速幂求第k个解.

给出一些求解脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#连分数求解pell方程最小整数解
from math import sqrt, floor

def pell_solver(D):
    a0 = floor(sqrt(D))
    if a0 * a0 == D:
        return None
    p0, p1 = 0, 1
    q0, q1 = 1, a0
    xn, yn = 1, 0
    while True:
        an = floor((sqrt(D) + p1) / q1)
        pn = an * q1 - p1
        qn = (D - pn * pn) // q1
        xn, xn_1 = an * xn + xn_1, xn
        yn, yn_1 = an * yn + yn_1, yn
        if xn * xn - D * yn * yn == 1:
            return xn, yn
       
from math import ceil,floor,sqrt


def pell_minimum_solution(n):
    a = []
    m = floor(sqrt(n))
    sq = sqrt(n)
    a.append(m)
    b = m
    c = 1
    i = 1
    while a[i-1] != 2 * a[0]:
        c = (n - b * b) / c
        tmp = (sq + b) / c
        a.append(floor(tmp))
        i += 1
        b = a[i-1] * c - b
    p = 1
    q = 0
    for j in range(i-2,-1,-1):
        t = p
        p = q + p * a[j]
        q = t
    if (i-1) % 2 == 0:
        x0 = p
        y0 = q
    else:
        x0 = 2 * p ** 2 + 1
        y0 = 2 * p * q
    return x0,y0



from Crypto.Util.number import *
from sympy import *

def solve_pell(N, numTry=1500):
    cf = continued_fraction(sqrt(N))
    for i in range(numTry):
        denom = cf.denominator(i)
        numer = cf.numerator(i)
        if numer^2 - N * denom^2 == 1:
            return numer, denom
    return None, None


from sympy.solvers.diophantine.diophantine import diop_DN
# sympy求解 x^2-Dy^2=N 输入(D,N)
x = diop_DN(105279230912868770223946474836383391725923//1293023064232431070902426583269468463,1)[0][0]