2022 D3CTF 的d3factor，顺便用来学习一元（univariate）的Coppersmith方法，然后尽量讲一下二元（bivariate）的方法；
二元的方法大部分跟一元类似的，只是有些细节还没搞懂。。。（菜
然后二元的方法还有个模的版本，拿了Boneh-Durfee的攻击来做栗子。

全文有点长，不过还是建议从前往后阅读，有些前面详细讲过的东西后面都是直接略了的（逃

附博客链接：https://tover.xyz/p/d3factor-coppersmith/
（有更新/纠错的话应该会先在博客上改）

d3factor

题目链接：https://paste.ubuntu.com/p/VkNv57QGkh/

先讲一下原理，当时队友已经找到这篇论文（[NR15]）了（不过当时不会Coppersmith没搞出来，），原理就在论文的第4章。

首先代码的14行以上是没漏洞的（应该），题目的意思大概是要从16-21行里解出p和q，然后14行前的RSA解m。从代码上看，有N下的两组密钥：
e_1d_1 \equiv 1 \pmod {\varphi(N)} \\ e_2d_2 \equiv 1 \pmod {\varphi(N)}
然后imath[/imath]是一个小的数（1000bits，相对于N是小的。。。），直觉上说是要把这个数搞出来的（而且确定了要用格做的话是肯定要提取一些小的未知数的），所以，把第一条式子乘一个e2，第二条式子乘一个e_1然后相减（有点像扩展Wiener里面Guo 's Approach的操作）：
e_1e_2d_1 \equiv e_2 \pmod {\varphi(N)} \\ e_1e_2d_2 \equiv e_1 \pmod {\varphi(N)} \\ e_1e_2(d_2-d_1) \equiv e_1-e_2 \pmod {\varphi(N)}
PS：构造到这一步，其实有点像之前说用格解方程的格式，但后来算了一下发现不在可解范围里，所以就需要更强大的Coppersmith方法（虽然论文里也有说- ）

接着，设x=d_2-d_1、a_0=e_2-e_1、a_1=e_1e_2，然后把上面最后一条式子整理一下，就是：
a_0+a_1x \equiv 0 \pmod {\varphi(N)}
然后，现在这个式子是模\varphi(N)余0的，而\varphi(N)=p^6(p-1)(q-1)，所以也是模p^6余0（展开一下就可以得到），之所以要转成模p^6的，是因为\varphi(N)是未知的，N是已知的，p^6的一个N的一个因子，一会用Coppersmith时某些p^6的操作可以用已知的N代替（至于具体是什么，一会再说-）。然后上式就可以化为：
a_0+a_1x \equiv 0 \pmod {p^6}
符合一元Coppersmith的格式，虽然，一元Coppersmith里是要求模数已知的，但是调函数时把模数设成N还是可以阴差阳错地解出来- -

一元的Coppersmith在sage里面已经可以直接调了（即官方wp的small_roots，虽然我也现在才知道- ），某个在github上找的代码也可以直接用，这里以后者怼个例子：

# sage
load('coppersmith.sage')

c = 2420624631315473673388732074340410215657378096737020976722603529598864338532404224879219059105950005655100728361198499550862405660043591919681568611707967
N = 1476751427633071977599571983301151063258376731102955975364111147037204614220376883752032253407881568290520059515340434632858734689439268479399482315506043425541162646523388437842149125178447800616137044219916586942207838674001004007237861470176454543718752182312318068466051713087927370670177514666860822341380494154077020472814706123209865769048722380888175401791873273850281384147394075054950169002165357490796510950852631287689747360436384163758289159710264469722036320819123313773301072777844457895388797742631541101152819089150281489897683508400098693808473542212963868834485233858128220055727804326451310080791
e1 = 425735006018518321920113858371691046233291394270779139216531379266829453665704656868245884309574741300746121946724344532456337490492263690989727904837374279175606623404025598533405400677329916633307585813849635071097268989906426771864410852556381279117588496262787146588414873723983855041415476840445850171457530977221981125006107741100779529209163446405585696682186452013669643507275620439492021019544922913941472624874102604249376990616323884331293660116156782891935217575308895791623826306100692059131945495084654854521834016181452508329430102813663713333608459898915361745215871305547069325129687311358338082029
e2 = 1004512650658647383814190582513307789549094672255033373245432814519573537648997991452158231923692387604945039180687417026069655569594454408690445879849410118502279459189421806132654131287284719070037134752526923855821229397612868419416851456578505341237256609343187666849045678291935806441844686439591365338539029504178066823886051731466788474438373839803448380498800384597878814991008672054436093542513518012957106825842251155935855375353004898840663429274565622024673235081082222394015174831078190299524112112571718817712276118850981261489528540025810396786605197437842655180663611669918785635193552649262904644919
a = (e2-e1)*(e1*e2).inverse_mod(N) % N
X = 2^1000

R = Zmod(N)
P.<x> = PolynomialRing(R, 1);
f = x-a
bounds = [X]
xs = small_roots(f, bounds, m=8)

assert len(xs) > 0
x0 = xs[0][0]
p6 = gcd(f(x=x0), N)
p = iroot(int(p6), 6)
assert p[1] == True
p = p[0]
q = N // p^7

print('p = %d' % p)
print('q = %d' % q)

Univariate Coppersmith

一元的Coppersmith是[Cop96a]里提出的，但是现在常用的是[HG97]里简化的版本，下面说的也是[HG97]的版本：

简单例子

考虑下面方程：
x^2+4x+3 \equiv 0 \pmod {37}
如果有条件|x| \le 4的话，简单计算一下，4^2+4*4+3=35<37，即模数37在上面其实是不起作用的，所以可以把\pmod {37}删掉然后写成：
x^2+4x+3 = 0
然后就变成了整数域上的方程就有很多方法可以解了。

简单原理

形式化一点就是，设f(x)=\sum^{d}_{i=0} f^{(i)} x^i，f(x) \equiv 0 \pmod N，另|f|(x)=\sum^{d}_{i=0} |f^{(i)}| x^i为f(x)的1范数（1-norm），若存在正整数界限X使得|f|(X)<N的话，则对于所有符合|x|\le X的x，都有|f(x)| < N，若刚好x是f(x)在模N下的根的话，x也是f(x)在整数域下的根。（有点像NTRU的某一步，逃

通常做题时都不会给出系数符合|f(x)| < N的（未知数可能会给出|x|\le X，或者需要构造一下），不然就可以直接解了，所以就需要有一个可以把系数减小的方法，而说到减小，自然（吗）就会想到格的SVP最短向量问题了。

格规约的可行性

先考虑一下把上面方程写成大一学的线性方程形式（不需要模数，因为只是为了降低系数）：
Ax \equiv 0 \pmod N
另L=UA的话，又可以写成：
Lx \equiv UAx \equiv U*0 \equiv 0 \pmod N
这里的L=UA可以看成是一个换基的行操作，或者直接一点，U就是LLL的操作，如果把A看作是某个系数矩阵的话，L就是格规约（reduce）后的系数矩阵（有更小的系数），而从上面的式子看来，规约操作后x并没改变，即是L系数下的（模N的）根时，也是A系数下的（模N的）根。

格的构造（例子）

到这还有个问题，f(x)只是一条方程，并不能把系数写成上面说的矩阵A的形式（顶多写成一个行向量）。要构造出一个矩阵（/格）的话，则还需要构造出一堆g_i(x)，符合：
g_i(x) \equiv 0 \pmod N
PS：每个g_i(x)的系数可以“转化”出一个行向量，所以构造的g_i(x)个数至少是这堆多项式的项数减1，即，g_i(x)加上f(x)对应这个矩阵A的行，这些多项式里的全部项（monomials）对应矩阵A的列，A要是一个基的话至少是一个方阵。

先给个直观点的例子描述一下这个矩阵是长啥样的，比如现在我有3个方程：
f(x) = 1 + 2x + 3x^2 \equiv 0 \pmod N \\ g_1(x) = 4x + 5x^2 \equiv 0 \pmod N \\ g_2(x) = 6x^2 \equiv 0 \pmod N
那么这个矩阵（/格）就是（PS：注意这里把X代进去了，原因是现在要求的是|f|(X)<N时的系数，也可参考下一节的推导）：
\begin{bmatrix} 1 & 2X & 3X^2 \\ 0 & 4X & 5X^2 \\ 0 & 0 & 6X^2 \\ \end{bmatrix}
从上到下三行分别对应三个多项式f(xX)、g_1(xX)和g_2(xX)的系数，从左到右三列分别对应三个项1、x和x^2。多项式和项的顺序其实是不重要的（不过记得格规约后的顺序跟规约前的保持一样就好），但是构造时为了美观/效率/容易算行列式通常会选择弄成三角矩阵（特别是论文里的，分析的时候大多是需要算行列式的- -）

使用条件的推导

怎么构造这些g_i(x)的话，其实有很多方法（而且不同的论文好像都有不同的方法- ），所以讲怎么构造之前，先来讲一下这个格（/矩阵）需要符合什么条件。

上面说了，对于f(x)=\sum^{d}_{i=0} f^{(i)} x^i，f(x) \equiv 0 \pmod N，要把模数N去掉的话，需要符合|f|(X)<N，再形式化一点，令向量：
v_f = (f^{(0)}, f^{(1)}X, f^{(2)}X^2, ... , f^{(d)}X^d)
则需要符合|v_f|<N，注意这里的|·|是指1范数（绝对值范数），因为格规约时用的是2范数（欧几里德范数，这里用\lVert v_f \rVert表示），所以需要用以下关系转换以下（维基）上有讲，证明我就略了（逃：

把我们的v_f代进去一下就是：
|v_f| \le \sqrt{d+1} \lVert v_f \rVert
为了保证|f|(X)<N一定会成立，联立一下|v_f|<N（注意这里我是在反推，即以下每多一个不等号界可能会更紧，然后最终推出最坏情况下的界）就是（PS：v_f里元素个数是imath[/imath]）:
|v_f| \le \sqrt{d+1} \lVert v_f \rVert < N
根据LLL的性质（取\delta=3/4时），

假设\lVert v_f \rVert是我们用LLL解出来的L中的最短向量，则会满足\lVert v_f \rVert \le 2^{d/4}\ det(L)^{1/(d+1)}，同样为了保证\sqrt{d+1} \lVert v_f \rVert < N一定会成立，再联立一下，就是：
|v_f| \le \sqrt{d+1} \lVert v_f \rVert \le \sqrt{d+1} · 2^{d/4} · det(L)^{1/(d+1) } < N
最后的"<"就是推出的条件，通常在使用时，度d是小的（可以看成常数），N是大的，即夸张点当N \to \infty时，\sqrt{d+1} · 2^{d/4}部分可以忽略，最终就是需要（d \ll N）：
det(L) < N^{d+1}
PS：注意换基不会改变格的行列式（的绝对值），所以也可以看成det(A) < N^{d+1}（A和L是上面例子的A和L）。

从上面结论来看，我们会想要构造出来的格行列式det(L)越小越好，而且从应用来看，会希望X越大越好，但从上节的例子来看，det(L)和X是相关的，即小的det(L)和大的X是矛盾的，所以，要使得构造出来的格在X尽量大的时候符合det(L) < N^{d+1}（应该）是一个挺难的问题。

具体步骤（例子）

这里先不考虑最优的构造，举个简单的例子（[Gal12]19.1.1节），顺便说说Coppersmith/HG97的方法是怎么操作的，而且怎么计算在这一个构造下X的界。

按照上面的说法，我们需要构造一堆g_i(x) \equiv 0 \pmod N，有模N余0的话比较简单且自然（吗）地会想到，只要g_i(x)里面地每一项都含N就好了。这里还有一处细节，因为g_i(x)（加上f(x)）是用来构造格基地，所以需要保证每个行向量v_{g_i}加上v_f是线性无关的。
g_i(x) = Nx^i, 0 \le i < d
然后构造出来的矩阵就是（上面说了，要代X，即矩阵里放的是f(xX)或g_i(xX)的系数，还有一处细节是，f(x)的最高项系数为1，即其实是可以对f(x)的所有项成上一个imath^{-1} \pmod N[/imath]然后把f(x)转化成monic的，这里这个过程略了）：
M = \begin{pmatrix} N & 0 & \cdots & 0 & 0 \\ 0 & NX & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & NX^{d-1} & 0 \\ f^{(0)} & f^{(1)}X & \cdots & f^{(d-1)}X^{d-1} & X^d \\ \end{pmatrix} \\

栗子

以[Gal12]的Example 19.1.6作个栗子来讲一下操作步骤（不过变量名被我换了）。模数N=10001，f(x)=x^3+10x^2+5000x-222，可以看出d=3（f(x)的最高次数），取X=10。首先构造g_i(x) = Nx^i, 0 \le i < d：
g_0(x) = Nx^0 = 10001 \\ g_1(x) = Nx^1 = 10001x \\ g_2(x) = Nx^2 = 10001x^2
这里的g_i(x)并没增加更多的项，所以可以简单看出这里的所有项是(1, x, x^2, x^3)，按照行是imath, g_1(xX), g_2(xX), f(xX))[/imath]的系数（注意顺序），列是上面对应的项，可以排出矩阵：
M = \begin{pmatrix} N & 0 & 0 & 0 \\ 0 & NX & 0 & 0 \\ 0 & 0 & NX^{2} & 0 \\ -222 & 5000X & 10X^{2} & X^3 \\ \end{pmatrix} \\
这里的N和X是已知的，把具体值代进去就是：
M = \begin{pmatrix} 10001 & 0 & 0 & 0 \\ 0 & 100010 & 0 & 0 \\ 0 & 0 & 1000100 & 0 \\ -222 & 50000 & 1000 & 1000 \\ \end{pmatrix} \\
对M做LLL（比如取参数\delta=3/4），得到L：
L= \begin{pmatrix} 444 & 10 & -2000 & -2000 \\ 9557 & -10 & 2000 & 2000 \\ -222 & 50000 & 1000 & 1000 \\ 989 & 2500 & 500100 & -500000 \\ \end{pmatrix} \\
根据LLL的性质，L的第一行（应该？）就是格L的最短向量，可以使用这个（行）向量恢复我们想要的小系数，即把这一行中每个元素和(1, X, X^2, X^3)中的元素分别相除，恢复出系数：
(\frac{444}{10^0}, \frac{10}{10^1}, \frac{-2000}{10^2}, \frac{-2000}{10^3}) = (444, 1, -20, -2)
可以简单验证一下这里的除法都是整除的：先把M中每一列的(1, X, X^2, X^3)分离出来，即拆成M=DM'，其中D是(1, X, X^2, X^3)的对角矩阵：
M = DM' = M = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & X & 0 & 0 \\ 0 & 0 & X^2 & 0 \\ 0 & 0 & 0 & X^3 \\ \end{pmatrix} · \begin{pmatrix} N & 0 & 0 & 0 \\ 0 & N & 0 & 0 \\ 0 & 0 & N & 0 \\ -222 & 5000 & 10 & 1 \\ \end{pmatrix} \\
假设U是LLL的操作（上面有说过），则L=UM=UDM'，由于D是对角矩阵，所以D在矩阵乘法中是可交换的，可得L=UDM'=D(UM')，即(1, X, X^2, X^3)可分别整除L中的对应列。

然后把系数(444, 1, -20, -2)和(1, x, x^2, x^3)分别对应回去，即得：
h(x) = -2x^3 - 20x^2 + x + 444
可以简单验证一下|h|(X)<N，把X带进去一下就得到（注意取绝对值）|h|(X)=4454 < 10001 = N。按照最开头说的，这时候的\mod N可以去掉，解h(x) \equiv 0 \pmod N即解h(x)=0\ \ on\ \ \mathbb{Z}。通过解-2x^3 - 20x^2 + x + 444可以解出x_0=4。

根据上面讲“可行性“时说的，h(x)的根与f(x)在模N下是一样的，即x_0=4也是f(x) \equiv 0 \pmod N的根（可以简单代进去验算一下）。

计算条件

在”条件推导“里说道，若d \ll N，则这个方法可以成功使用时，需要det(L) < N^{d+1}，但没有说X需要符合的条件，因为使用不同构造的格算出来的X条件是不一样的。

拿这节构造的格L作栗子，首先根据”条件推导“里的结论，我们只需要算出det(L)就好了。然后，这节构造出来的格基M是个三角矩阵，det(L)即其对角线上的乘积，即：
det(L) = X^d * \Pi_{i=0}^{d-1} (NX^i) = N^{d} \Pi_{i=0}^{d} X^i = N^dX^{\frac{d(d+1)}{2}}
联立，即：
det(L) = N^dX^{\frac{d(d+1)}{2}} < N^{d+1} \\
解一下，最终得到X需要符合的条件是（d \ll N时）：
X < N^{\frac{2}{d(d+1)}}
PS：感觉这个只是x_0的范围，因为把上面栗子的X代进去是不符合的，但x_0符合。

更优的构造

上一节讲的栗子只是个直观的栗子，肯定不是最优的，怎么可以继续优化呢。看回det(L) < N^{d+1}，前面也讲到，det(L)中是含X的，想要同时得到小的det(L)和大的X（似乎）是难的，但是如果可以增大"<"右边呢。有两种方法：增大格的维度（项数）和增大模数：

如果有f(x) \equiv 0 \pmod N的话，那么也会有x^if(x) \equiv 0 \pmod N，通过构造一堆x^if(x)就可以增大多项式的数量和项的数量，从而增大\omega（虽然det(L)也会增大就是了），这种方法又叫“x-shift”；
如果有f(x) \equiv 0 \pmod N的话，那么也会有f^t(x) \equiv 0 \pmod {N^t}，这样就可把原来的模数N增大到N^t（同时项数也会增多，在项数增多后为了弄出个方阵还要增加多项式的数量，所以通常会和上一个方法混着用），其实更广义一点是：N^{t-i}f^i(x) \equiv 0 \pmod {N^t}。

到这里，之前推导的关系就要重写一下了：
det(L) < N^{\omega t}
其中\omega是我们构造的格的维度（原本的d+1），t就是原本的N扩到N^t。

[Cop96]/[HG97]的构造就是混合了上面两种方法，构造了（当然变量名应该和我不一样）：
g_{i, j}(x) = N^{t-i}f^{i}(x)·x^j,\ \ 0 \le i \le t,\ 0 \le j < d
然后（根据上面两个方法说的）就会有g_{i, j} \equiv 0 \pmod {N^t}；多项式数是d·(t+1)，项数是d·t+t=d·(t+1)，即构造出来的格基维度是\omega=d·(t+1)。稍微算一下行列式：
det(L) = \prod_{i=0}^{t} \prod_{j=0}^{d-1} (N^{t-i}·X^{di}·X^j) = N^{dt(t+1)/2} X^{d(t+1)(dt+d-1)/2}

（PS：不会有人手算吧-）

代进det(L) < N^{\omega t}，就有（d \ll N时）：
N^{t} X^{(dt+d-1)} < N^{\frac{2d·(t+1)·t}{d(t+1)}} =N^{2t} \\ X < N^{\frac{t}{dt+d-1}}
即这个构造的成功条件是|x_0| < N^{\frac{t}{dt+d-1}}。

PS：显而易见的是，参数t（注意t是我们选的参数）越大，|x_0|就可以越大：\frac{t}{dt+d-1} = \frac{1}{d+\frac{d-1}{t}}，t越大时整个N的指数部分也越大，当t \to \infin时，"<"右边趋近N^{1/d}，即：
|x_0| < N^\frac{1}{d},\ t \to \infin
当然，t越大构造出来的格也越大。

最后拿[Gal12]的Example 19.1.11做个栗子：

# sage
# stolen from 2022 TQLCTF wp
def matrix_overview(BB):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            if BB[ii, jj] == 0:
                a += ' '
            else:
                a += 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        print(a)

# input
p = 2^30 + 3
q = 2^32 + 15
N = p*q
a0, a1, a2, a3 = (1942528644709637042, 1234567890123456789, 987654321987654321, 1)
P.<x> = PolynomialRing(ZZ)
f = a0 + a1*x + a2*x^2 + a3*x^3
X = 2^14

d = f.degree()
# choose yourself, bigger better and slower
t = 3

# coustruct gi
g_ij = []
for i in range(t+1):
  for j in range(d):
    g_ij.append( N^(t-i) * f^i * x^j )
g_ij = sorted(g_ij)
#print(g_ij)

monomials = []
for g in g_ij:
  monomials += g.monomials()
monomials = sorted(set(monomials))
print('\nmonomials:')
print(monomials)

# construct matrix
w = d * (t+1)
assert w == len(monomials)
M = matrix(ZZ, w, w)
for i in range(w):
  for j in range(w):
    if monomials[j] in g_ij[i].monomials():
      M[i, j] = g_ij[i](x=x*X).monomial_coefficient(monomials[j])
print('\nConstructed matrix:')
matrix_overview(M)

# reduce coeffcients
L = M.LLL()
h = sum([(L[0][i]/monomials[i](x=X)) * monomials[i] for i in range(w)])
print('\nReduced poly:')
print(h)
print('\nRoots of h:')
print(h.roots())
  
# simple check
res = []
for r in h.roots():
  xx = r[0]
  if h(x=xx) < N^t:
    res.append(xx) 

print('\nGot x0:')
print(res)

另

这章的内容主要参考了[Gal12]的19章和[Jou09]的13章，还有一些格的基础内容参考[HPS14]的第7章，还有一些论文（略）。

d3factor again

讲了这么多再回来看一下d3factor，如果上一章有看完的话会发现一个问题，那个多项式方程的模式是要已知的，而d3factor的模数是p^6，p^6肯定是未知的（不然就可以直接解了），然后官方wp用small_roots的时候模数的参数传的是N，还居然可以直接解出来- -

首先要提到的是，[NR15]的其实是用了[LZP15]里3.1的方法，而[LZP15]里构造的g_k是长这样的：

最后还提到，g_k(y) \equiv 0 \pmod {p^{vt}}，即模数是p^{vt}（给没看过这篇论文的观众标注一下，这里g^v=gcd(N, \varphi(N))），仍然是未知的，但是因为g^v是N的一个因子，所以仍然可以用N来扩大模数（参考“更优构造”的第二个方法），但是要注意N中含的p^v的“份量”要大于t-k，即要确保\pmod {p^{vt}}下是余0的。所以就有了上图中的max\{·\}.

而直接用N作模数传进small_roots也可以解出来的原因是，按照上一章的构造，p^v的“份量”只会比[LZP15]的更多，所以还是能保证g_k(y) \equiv 0 \pmod {p^{vt}}，但是显然（吗？）算出来的|x_0|范围会更小，如果题目的输入卡一下的话估计small_roots会算不出来。

直观上来看，对于$det(L) < N^{\omega t}$，两种方法都不会影响"$<$"右边的$N^{\omega t}$，而对于"$<$"左边的$det(L)$，如果使用small_roots的话，$N$的“含量”会增加（没看过源码，但估计是上一章的算法），则$X$的“含量”会相应减小，则$|x_0|$的取值会更小。（具体小多少就懒得算了，逃）

最后附一个[LZP15]解法的代码（和上面的大同小异）：

# sage
# from 2022 tql wp
def matrix_overview(BB):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            if BB[ii, jj] == 0:
                a += ' '
            else:
                a += 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        print(a)

# https://eprint.iacr.org/2014/343.pdf (3.1)
def YaoLu(f, N, m, t, X):
  P = f.parent()
  dim = m+1
  
  # as g_k(x) in YaoLu's paper
  gks = []
  for k in range(dim):
    gks.append( N^max(ceil(v*(t-k)/u), 0) * f^k )
  gks = sorted(gks)

  # order is not important
  # but monomials would be convenient (:
  monomials = []
  for gk in gks:
    monomials += gk.monomials()
  monomials = sorted(set(monomials))
  print(monomials)

  assert dim == len(monomials)
  M = matrix(QQ, dim, dim)
  for i in range(dim):
    for j in range(dim):
      if monomials[j] in gks[i].monomials():
        M[i, j] = gks[i](x=x*X).monomial_coefficient(monomials[j])
  # just debug
  matrix_overview(M)

  # LLL reduce coefficients
  L = M.LLL()
  h = sum([(L[0][i]/monomials[i](x=X)) * monomials[i] for i in range(dim)])

  # simple check
  res = []
  for r in h.roots():
    xx = r[0]
    if h(x=xx) < N^t:  # p^{vt} < N^t, loose
      res.append(xx) 
  return res

if __name__ == '__main__':
  from gmpy2 import iroot
  import libnum

  # known
  c = 2420624631315473673388732074340410215657378096737020976722603529598864338532404224879219059105950005655100728361198499550862405660043591919681568611707967
  N = 1476751427633071977599571983301151063258376731102955975364111147037204614220376883752032253407881568290520059515340434632858734689439268479399482315506043425541162646523388437842149125178447800616137044219916586942207838674001004007237861470176454543718752182312318068466051713087927370670177514666860822341380494154077020472814706123209865769048722380888175401791873273850281384147394075054950169002165357490796510950852631287689747360436384163758289159710264469722036320819123313773301072777844457895388797742631541101152819089150281489897683508400098693808473542212963868834485233858128220055727804326451310080791
  e1 = 425735006018518321920113858371691046233291394270779139216531379266829453665704656868245884309574741300746121946724344532456337490492263690989727904837374279175606623404025598533405400677329916633307585813849635071097268989906426771864410852556381279117588496262787146588414873723983855041415476840445850171457530977221981125006107741100779529209163446405585696682186452013669643507275620439492021019544922913941472624874102604249376990616323884331293660116156782891935217575308895791623826306100692059131945495084654854521834016181452508329430102813663713333608459898915361745215871305547069325129687311358338082029
  e2 = 1004512650658647383814190582513307789549094672255033373245432814519573537648997991452158231923692387604945039180687417026069655569594454408690445879849410118502279459189421806132654131287284719070037134752526923855821229397612868419416851456578505341237256609343187666849045678291935806441844686439591365338539029504178066823886051731466788474438373839803448380498800384597878814991008672054436093542513518012957106825842251155935855375353004898840663429274565622024673235081082222394015174831078190299524112112571718817712276118850981261489528540025810396786605197437842655180663611669918785635193552649262904644919
  a = (e1-e2)*(e1*e2).inverse_mod(N) % N

  P.<x> = PolynomialRing(ZZ)
  f = x-a 
  X = 2^1000
  u = 7
  v = u-1

  # YaoLu's params
  m = 8
  t = 6

  # main
  xs = YaoLu(f, N, m, t, X)
  assert len(xs) > 0
  x0 = xs[0]

  # factor p, q
  pv = gcd(f(x=x0), N)  # gcd(k*p^vt, p^u*q) == p^v
  p = iroot(int(pv), v)
  assert p[1] == True
  p = p[0]
  q = N // p^7
  print('p = %d' % p)
  print('q = %d' % q)

  # hack
  n = p*q
  phi = (p-1)*(q-1)
  e = 65537
  d = e.inverse_mod(phi)
  m = pow(c, d, n)
  flag = libnum.n2s(int(m))
  print(flag)

Bivariate Coppersmith

二元的Coppersmith在[Cop96b]里提出，然后Coron在[Cor04]和[Cor07]里提出了更简化的版本。

上图截自[Cop96b]，注意这里的多项式p是在\mathbb{Z}上的，而不是\pmod N，然后大概就是构造一个q(x, y) \equiv 0 \pmod N，（下图截自[Cop96b]的完整版本-）

通过类似一元Coppersmith的方法规约出h(x, y) \equiv 0 \pmod N，如果h(x, y)的系数足够小的话，模数N可以去掉，就是h(x, y) = 0。最后如果p(x, y)和h(x, y)不相关的话（上面写定理时假设了p(x, y)是不可约的，所以也不会相关），就是两条方程两个未知数，联立解出x、y。

PS：具体的做法其实是下面这样的，先构造Resultant解出x_0，然后x_0代进其中一条方程解出y_0。

模N下的二元Coppersmith

在实际应用（或做题）中，通常需要解的是p(x, y) \equiv 0 \pmod N，比如著名的Boneh-Durfee攻击（[BD00]）。

这种二元模N下的Coppersmith解法其实更像三元的Coppersmith（想一想把二元的同余方程展开不就是整数上的三元方程了- ），而且还少了构造q(x, y) \equiv 0 \pmod N这一步：

上面大概意思是说，我们对p(x, y) \equiv 0 \pmod N使用类似一元Coppersmith的方法进行LLL规约后，如果只拿格基中的第一行构造h_1(x, y)是不够的（因为一个方程两个未知数- ），所以还可以拿格基的第二行（或者其他的行）来构造h_2(x, y)，这样才有两个方程。当X和Y足够小的时候，可以同时保证整数上有h_1(x, y)和h_2(x, y)，联立这两条方程即可解x_0和y_0。

下面大概说说[BD00]的构造（详细的话看[BD00]的Ⅳ章就好了，略- ）：

[BD00]造格时用了两种“偏移”：x-shift的g_{i, k}(x, y)和y-shift的h_{j, k}(x, y)，其中k=0, ..., m、i=0, ..., m-k、j=0, ..., t，这里的m和t都是参数。首先注意[BD00]里的模数是e（即RSA的解密指数），k+(m-k)=m，即g_{i, k}(x, y)和h_{j, k}(x, y)都是模e^m余0的。然后i和j的选择估计是为了把格基弄成方阵（表示不知道怎么想出来的，大概可能的想法可以参考下面维度的计算）。

接下来算一下这个格的维度，首先是多项式的数量（行数），即g_{ik}(x, y)的个数加上h_{jk}(x, y)的个数，减去g_{ik}(x, y)和h_{jk}(x, y)重叠的个数（在i=0和j=0时，g_{0k}(x, y)=h_{0k}(x, y)=f^k(x, y)e^{m-k}）：
\sum_{k=0}^{m}(\sum_{i=0}^{m-k} + \sum_{j=0}^{t} -1) = \frac{ (m+1) (2t+m+2) }{2}

然后是项的数量，直接算项数好像有点困难，但可以取巧一下，比如证明除了在i=0和j=0时，g_{ik}(x, y)和h_{jk}(x, y)的最高项都是不同的。设d_x和d_y分别是f(x, y)里x和y的度，则g_{ik}(x, y)和h_{jk}(x, y)的最高项可以表示如下：
g_{ik}:\ x^{kd_x+i}y^{kd_y} \\ h_{jk}:\ x^{kd_x}y^{kd_y+j}
以上最高项中，g_{ik}的x指数部分总会比h_{jk}的多了i，h_{jk}的y指数部分总会比g_{ik}的多了j，所以可以知道只有在i=0和j=0时上面才会有相同的项出现。而i=0和j=0刚好就是上面算多项式个数时重叠的情况，所以除去这些重叠的部分后即最高项的个数（然后应该还要证一下这些最高项包含所有项的，不会，逃）。随后算出的维度就是：

然后算行列式的话，把g_{ik}(x, y)和h_{jk}(x, y)排列一下的话，格基可以排成上面列出来的“最高项”为对角线的三角矩阵，所以行列式就是上面的“最高项”（删去重叠部分）的乘积了，和[BD00]一样，我把行列式分成det=det_x·det_y来计算。首先是det_x，即g_{ik}的最高项乘积（不会有人手算吧*2）：

然后是det_y，即h_{jk}的最高项乘积，注意这里我删去了和g_{ik}重叠的部分：

最后两个乘起来就是总的行列式：

最最后把行列式代进det<e^{\omega m}就是X和Y需要符合的条件（复鬼杂咯- -）：

附一个从2022 TQLCTF里偷来的代码，应该可以直接用（~~清华改的代码就是不一样~~），另，这里最后对H的两个循环其实就是在试规约后的格的全部行中哪两行是不相关的，其他的应该上面都有说了（逃：

# sage
def matrix_overview(BB):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            if BB[ii, jj] == 0:
                a += ' '
            else:
                a += 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        print(a)

def lattice_attack(PR, pol, e, mm, tt, X, Y):
    x, y = PR.gens()
    polys = []

    for ii in range(mm + 1):
        for jj in range(0, mm - ii + 1):
            poly = e ^ (mm - ii) * x ^ jj * pol ^ ii
            polys.append(poly)

    for ii in range(mm + 1):
        for jj in range(1, tt + 1):
            poly = e ^ (mm - ii) * y ^ jj * pol ^ ii
            polys.append(poly)

    polys = sorted(polys)
    monomials = []
    for poly in polys:
        monomials += poly.monomials()
    monomials = sorted(set(monomials))
    dims1 = len(polys)
    dims2 = len(monomials)
    M = matrix(QQ, dims1, dims2)

    for ii in range(dims1):
        M[ii, 0] = polys[ii](0, 0)
        for jj in range(dims2):
            if monomials[jj] in polys[ii].monomials():
                M[ii, jj] = polys[ii](x * X, y * Y).monomial_coefficient(monomials[jj])
    
    matrix_overview(M)
    print('=' * 128)
    print(len(M.rows()), len(M.columns()))
    #M = M.hermite_form()
    B = M.LLL()
    print('LLL done')

    det = B.det()
    print(f"monomials: {monomials}")
    nn = len(monomials)
    matrix_overview(B)
    H = [(i, 0) for i in range(dims1)]
    H = dict(H)
    for j in range(dims2):
        for i in range(dims1):
            H[i] += (monomials[j] * B[i, j]) / monomials[j](X, Y)
    
    PQ.<q> = PolynomialRing(ZZ)
    H = list(H.values())
    solutions = []
    print(len(H))
    for i in range(len(H)):
        for j in range(i + 1, len(H)):
            pol1 = PR(H[i])
            pol2 = PR(H[j])
            rr = pol1.resultant(pol2, y)
            if rr.is_zero() or rr.monomials() == [1]:
                continue
            sols = rr(q, q).roots()
            for sol in sols:
                solx = sol[0]
                if solx == -1:
                    continue
                try:
                    soly = pol1(solx, q).roots()[0][0]
                    solutions.append((solx, soly))
                    print('=' * 128)
                except:
                    pass
                if len(solutions) > 0:
                    break
            if len(solutions) > 0:
                break
        if len(solutions) > 0:
            break
    return solutions

另外，[BD00]的Ⅴ章里有继续优化的方案，即构造出来的格基据说可以删掉某些行和列的（具体略）；另另外[Cor07]的方法（好像）在构造时不用考虑构造成方阵的，留个坑有空撸个代码（菜，
https://github.com/pwang00/Cryptographic-Attacks/blob/master/Public%20Key/RSA/coron.sage

总结

菜（：

参考

[BD00] Boneh D, Durfee G. Cryptanalysis of RSA with private key d less than N/sup 0.292[J]. IEEE transactions on Information Theory, 2000, 46(4): 1339-1349.
[Cop96a] D. Coppersmith, Finding a Small Root of a Univariate Modular Equation, proceedings of Eurocrypt ’96, vol. 1070, Lecture Notes in Computer Science.
[Cop96b] D. Coppersmith, Finding a Small Root of a Bivariate Integer Equation; Factoring with High Bits Known, proceedings of Eurocrypt’ 96, vol. 1070, Lecture Notes in Computer Science.
[Cor04] Coron J S. Finding small roots of bivariate integer polynomial equations revisited[C]//International Conference on the Theory and Applications of Cryptographic Techniques. Springer, Berlin, Heidelberg, 2004: 492-505.
[Cor07] Coron J S. Finding small roots of bivariate integer polynomial equations: A direct approach[C]//Annual International Cryptology Conference. Springer, Berlin, Heidelberg, 2007: 379-394.
[Gal12] Galbraith, Steven D. Mathematics of public key cryptography. Cambridge University Press, 2012.
[HG97] N. A. Howgrave-Graham, Finding small roots of univariate modular equations revisited. In Cryptography and Coding, volume 1355 of LNCS, pp. 131-142. Springer Verlag, 1997.
[HPS14] Hoffstein, Jeffrey, et al. An introduction to mathematical cryptography. Vol. 1. New York: springer, 2008.
[Jou09] Joux, Antoine. Algorithmic cryptanalysis. Chapman and Hall/CRC, 2009.
[LZP15] Lu Y, Zhang R, Peng L, et al. Solving linear equations modulo unknown divisors: revisited[C]//International Conference on the Theory and Application of Cryptology and Information Security. Springer, Berlin, Heidelberg, 2015: 189-213.
[NR15] Nitaj A, Rachidi T. New Attacks on RSA with Moduli N= p r q[C]//International Conference on Codes, Cryptology, and Information Security. Springer, Cham, 2015: 352-360.
https://github.com/defund/coppersmith
https://github.com/pwang00/Cryptographic-Attacks/blob/master/Public%20Key/RSA/coron.sage