from sage.allimport* d = ... cf = continued_fraction(sqrt(d)) a0 = cf[0] i = 1 whileTrue: if cf[i] == 2*a0: c = cf.convergent(i-1) x, y = c.as_integer_ratio() if x**2 - d*y**2 == 1: print((x,y)) break i = i + 1
在这里扩展地提一下广义佩尔方程:
定义
形如的方程称为广义佩尔方程.
求解
通过连分数求出广义佩尔方程的最小正整数解后,可以知道也是该方程的整数解(其中为方程的整数解)
求解代码
python
from sage.allimport * defpell_roots(D: int, C: int = 1): intervals = 2**10 defget_a_root(D: int, C: int): cf = continued_fraction(sqrt(D)) for i inrange(1, intervals): c = cf.convergent(i - 1) x, y = c.as_integer_ratio() if x**2 - D * y**2 == C: return x, y for y inrange(1, intervals): x2 = C + D * y**2 if x2 <= 0: continue x = isqrt(x2) if x ** 2 == x2: return x, y return0, 0 r, s = get_a_root(D, 1) x, y = get_a_root(D, C) whileTrue: yield x, y x, y = x * r + D * y * s, r * y + s * x D = ... C = ... x, y = next(pell_roots(D, C)) print(x**2 - D * y**2 == C) print((x,y))