d3factor/Boneh-Durfee与Coppersmith学习笔记

2022 D3CTF 的d3factor，顺便用来学习一元（univariate）的Coppersmith方法，然后尽量讲一下二元（bivariate）的方法；

二元的方法大部分跟一元类似的，只是有些细节还没搞懂。。。（菜

然后二元的方法还有个模的版本，拿了Boneh-Durfee的攻击来做栗子。

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

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)}$$

然后$(d_2-d_1)$是一个小的数（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$里元素个数是$(d+1)$）:

$$|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)$的所有项成上一个$(f^{d})^{-1} \pmod N$然后把$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)$，按照行是$(g_0(xX), g_1(xX), g_2(xX), f(xX))$的系数（注意顺序），列是上面对应的项，可以排出矩阵：

$$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}$；
• 最后把两种方法结合的话，就有$N^{t-i} f^i(x) x^j \equiv 0 \pmod {N^t}$，可以简单验证一下，因为$f(x) \equiv 0 \pmod N$，所以存在$k$使得$f(x) = kN$，于是

  $$N^{t-i} f^i(x) x^j = N^{t-i} k^i N^i x^j = k^i x^j N^t \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)

2022.3.24补充·

大概算了一下怎么卡参数（x），补充下过程。YaoLu的方法在论文里已经有给出$|x_0|$的范围了，就直接引结论了（变量名被我换了，且把 $u$和$v$代进了具体的$r$ ，这里的 $\beta$大概是$p$的规模，在d3factor中即$\frac{256}{2048}=\frac{1}{8}$，$\epsilon$是误差）：

$$|x_0| \le N^{r(r-1)\beta^2 - \epsilon}$$

所以这里主要关注HG97的方法的范围。方法还是关注$det(L) < N^{\omega t}$（类似但不完全相同），首先算$det(L)$，方法类似“更优构造”里讲的，或者直接把$d=1$代进去就好了，算出：

$$det(L) = \prod_{i=0}^t N^{t-i}X^i = N^{\frac{t(t+1)}{2}} X^{\frac{t(t+1)}{2}}$$

$\omega$就是维度了即$\omega=t+1$。然后重点关注$N^{\omega t}$里的$N$，前面说了这个$N$本应是使得$g_{i, j} \equiv 0 \pmod {N^t}$的那个模数里的$N$，现在换到YaoLu的方法里的话是$g_k(y) \equiv 0 \pmod {p^{vt}}$，即原本的$N$现在应该是$p^v=p^{r-1}$。所以最终须符合的关系是：

$$det(L) < p^{\omega t (r-1)}$$

带进去具体值：

$$(p^rq)^{\frac{t(t+1)}{2}} X^{\frac{t(t+1)}{2}} < p^{t(t+1)(r-1)}$$

化简（假设$p$和$q$同规模）：

$$X < \frac{p^{r-2}}{q} \approx p^{r-3}$$

放到d3factor里，如果用YaoLu的方法算出来的上界是$1344-\epsilon$ bits（$\epsilon$是某个误差），用HG97算出来的大概是$1024$ bits。

实际测试的话YaoLu的大概是$1225$ bits，HG97的大概是题目的$1000$ bits（small_roots的跟HG97相近）

PS：原来small_roots用的不是HG97那篇论文。。。

最后附测试代码：

# sage
#
# 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))

  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

# gen
size = 1225
from Crypto.Util.number import getPrime
from sympy import nextprime
from gmpy2 import invert, iroot
p = getPrime(256)
q = getPrime(256)
N = pow(p, 7) * q
phi = pow(p, 6) * (p - 1) * (q - 1)
d1 = getPrime(2000)
d2 = nextprime(d1 + getPrime(size))
e1 = invert(d1, phi)
e2 = invert(d2, phi)
print('Gen:')
print('p = %d' % p)
print('q = %d' % q)

N = Integer(N)
e1 = Integer(e1)
e2 = Integer(e2)
print()

# from Nu1L wp
print('small_roots:')
try:
  P.<x> = PolynomialRing(Zmod(N))
  f = e1*e2*x - (e1-e2)
  f = f.monic()
  res = int(f.small_roots(X=2**size, beta=0.75)[0])
  #p = int(gcd(int(f(res)), N))  # ????????
  p = iroot(int(gcd(int(f(res)), N)), 6)
  assert p[1] == True
  p = p[0]
  q = N // p**7
  print('p = %d' % p)
  print('q = %d' % q)
except:
  print('Fail')
print()


print('YaoLU:')
try:
  P.<x> = PolynomialRing(ZZ)
  a = (e1-e2)*(e1*e2).inverse_mod(N) % N
  f = x-a 
  X = 2^size
  u = 7
  v = u-1

  # YaoLu's params
  m = 16
  t = 12

  # 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)

except:
  print('Fail')
print()


# HG97
print('HG97:')
try:
  P.<x> = PolynomialRing(ZZ)
  a = (e1-e2)*(e1*e2).inverse_mod(N) % N
  f = x-a 
  X = 2^size

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

  # 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))

  # 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])

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

  assert len(res) > 0
  x0 = res[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)

except:
  print('Fail')

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

总结·

菜（：

参考·

1. [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.
2. [Cop96a] D. Coppersmith, Finding a Small Root of a Univariate Modular Equation, proceedings of Eurocrypt ’96, vol. 1070, Lecture Notes in Computer Science.
3. [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.
4. [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.
5. [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.
6. [Gal12] Galbraith, Steven D. Mathematics of public key cryptography. Cambridge University Press, 2012.
7. [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.
8. [HPS14] Hoffstein, Jeffrey, et al. An introduction to mathematical cryptography. Vol. 1. New York: springer, 2008.
9. [Jou09] Joux, Antoine. Algorithmic cryptanalysis. Chapman and Hall/CRC, 2009.
10. [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.
11. [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.
12. https://github.com/defund/coppersmith
13. https://github.com/pwang00/Cryptographic-Attacks/blob/master/Public%20Key/RSA/coron.sage