crypto-Pell方程（第一类佩尔方程:x 2 −D∗y 2 =1）

定义

若一个方程具有形式:
$$x^2-D\ast y^2=1(D为正整数，且D不是完全平方数）$$

则称次方程为第一类佩尔方程，要求第一类佩尔方程的解都是正整数解（因此限制D）。

（其实有第二类，但太难了，一般用不到就不展开讲了）

性质

显然佩尔方程 x^2 - Dy^2 = 1有无数个解，化简一下就变成了
$$x=\sqrt{1+D\ast y^2}$$
从y=1一个个试过去，如果√1+D*y^2是整数，则此时（x,y)为方程的一个解，当（x,y)最小时，也就是**最小正整数解(x1,y1)，**设解为(x1,y1)，(x2,y2)，……。所有的解(xi,yi)可表示成如下形式：
$$x^i+\surd D\ast Y^i=(x1+y1\ast D{)}^i$$
或者由以下的递推式得到:
$$x_{(}i+1)=x_1\ast x_i+n\ast y_1\ast y_i$$

$$y_{(}i+1)=x_1\ast y_i+y_1\ast x_i$$

由此递推式可以用前几组接就能得到后面的解

求最小正整数解的脚本如下：

import numpy as np
from collections import deque

d =
m = int(np.sqrt(d))
dq = deque()
dq.append(m)
n0 = n1 = d - m * m
m1 = m
while 1:
    q, m2 = divmod(m1 + m, n1)
    dq.appendleft(q)
    m1 = -m2+m
    n1 = (d-m1*m1)//n1
    if m1 == m and n1 == n0:
        break

dq.popleft()
x = 1
y = 0
for i in dq:
    x1 = y + x * i
    y = x
    x = x1
    y1=y
if x*x-d*y*y==-1:
   b=(x**2+d*y**2)
   y1=2*x*y
   x1=b
print('x1=',x1)
print('y1=',y1)

递推式推导过程

先取方程的相邻的两个解(x1,y1),(x2,y2)，显然有：
$$x_1^2-d\ast y_1^2=1$$
$$x_2^2-d\ast y_2^2=1$$

然后将两式相乘可以得到
$$(x_1^2-d\ast y_1^2)(x_2^2-d\ast y_2^2)=1$$
再对上式进行展开
$$x_1^2{x}_2^2+d^2{y}_1^2{y}_2^2+2d{x}_{1}x2y_1{y}_2-d(d{x}_1^2{y}_2^2+y_1^2{x}_2^2+2x_1{x}_2{y}_1{y}_2)=1$$
再合并得到
$$(x_1{x}_2+d{y}_1{y}_2{)}^2-d(x_1{y}_2+x_2{y}_1{)}^2=1$$
合并后发现是佩尔方程的形式，于是得到
$$x_3=x_1{x}_2+d{y}_1{y}_2$$

$$y_3=x_1{y}_2+x_2{y}_1$$

例题

[HZNUCTF 2023 final]一步到喂

from Crypto.Util.number import *
from secret import flag,x,y
p=getPrime(512)
q=getPrime(512)
n=p*q
e=x
assert 1293023064232431070902426583269468463*pow(x,2)==105279230912868770223946474836383391725923*pow(y,2)+1293023064232431070902426583269468463
m=bytes_to_long(flag)
c=pow(m,e,n)
print(p,q)
print(c)
#9850212100620338486478343799784951721581757936621772924971218807300819701941457605399099898870264241518769370682330612103782092302148525012450902351701339
#10749429992014823019923966246896511618886613763258781706004694804949547801668777655988055847885755337127548775758133804022361510427909703124161450470578543
#66847321508502017574023222490247591875197038421108556106531660421662626233521063441647157067220450229816184622038471812597874582613385516838920632450015292570673816423432903604941781889308906893966588610214614726388822851471695742453496232748358301888465563812947038856742838097152549971517159475947566599664

仔细阅读题目条件，会发现

1293023064232431070902426583269468463*pow(x,2)==105279230912868770223946474836383391725923*pow(y,2)+1293023064232431070902426583269468463

很像佩尔方程的形式，于是我们应该会下意识的去约分，经过尝试后会发现是可行的，于是得到
$$x^2-81421y^2=1$$
用上面介绍的代码可以求出整数解(x1,y1),然后求出e，最后按RSA算法的步骤求解

exp

import numpy as np
from collections import deque
import gmpy2
from Crypto.Util.number import *
d =61
m = int(np.sqrt(d))
dq = deque()
dq.append(m)
n0 = n1 = d - m * m
m1 = m
while 1:
    q, m2 = divmod(m1 + m, n1)
    dq.appendleft(q)
    m1 = -m2+m
    n1 = (d-m1*m1)//n1
    if m1 == m and n1 == n0:
        break

dq.popleft()
x = 1
y = 0
for i in dq:
    x1 = y + x * i
    y = x
    x = x1
    y1=y
if x*x-d*y*y==-1:
   b=(x**2+d*y**2)
   y1=2*x*y
   x1=b
#x1=706458746417678962043621845971467865659328419354757470370053791959400836343459667183452696881749102272803495561977766862981211957547107314899837090421877848369768126470488899783907413049018564965473411482320934865623338739674557285559179500448017722749519381238449907463225595410791751764554609076679816077991393296019025846726102143252917723127850421484059625048891687413514675773883846965460650321320194376953023128886680992904874680874592398084000621265936352579085146973168432725924933271736161246147210167219273151533549827720768127049
e=x1
p=9850212100620338486478343799784951721581757936621772924971218807300819701941457605399099898870264241518769370682330612103782092302148525012450902351701339
q=10749429992014823019923966246896511618886613763258781706004694804949547801668777655988055847885755337127548775758133804022361510427909703124161450470578543
c=66847321508502017574023222490247591875197038421108556106531660421662626233521063441647157067220450229816184622038471812597874582613385516838920632450015292570673816423432903604941781889308906893966588610214614726388822851471695742453496232748358301888465563812947038856742838097152549971517159475947566599664
n=p*q
phi_n =(q-1)*(p-1)
d = gmpy2.invert(e, phi_n)
m = gmpy2.powmod(c, d, n)
flag = long_to_bytes(m)
print(flag)
#b'HZNUCTF{D0_y0u_know_p3ll_3qu4tion!!!}'