总结了一下和Pollard’s p-1算法原理相似的几个算法，以及最近出现过的一些题目

Pollard’s p-1·

Pollard’s p-1算法大家应该都很熟悉，我之前在第一届强网杯密码数学专项赛的文章中也介绍过，大概思路是：

令

n=p⋅qn = p \cdot q

如果存在一个整数LL，使得

p−1∣Lq−1∤Lp - 1 \mid L \\ q - 1 \nmid L

那么根据费马小定理，大概率会存在一个aa，使得

aL≡1(modp)aL≢1(modq)a^L \equiv 1 \pmod p \\ a^L \not\equiv 1 \pmod q \\

通过计算

p=gcd(aL−1(modn),n)p = gcd(a^L - 1 \pmod n, n)

就可以实现分解

当然，算法的关键是怎么找到这个LL

传统的Pollard’s p-1的做法是，枚举i=1,2,3,⋯i = 1, 2, 3, \cdots，然后计算

a0=1ai=ai⋅ai−1(modn)\begin{aligned} a_0 &= 1 \\ a_i &= a^i \cdot a_{i-1} \pmod n \end{aligned}

最后结果相当于

L=∏iiaL≡aL(modn)\begin{aligned} L = \prod_i i \\ a_L \equiv a^L \pmod n \end{aligned}

2021羊城杯的Easy RSA·

题目：https://www.nssctf.cn/problem/237

显然，传统的Pollard’s p-1只能对付p−1p-1或者q−1q - 1是平滑数的分解，因为LL至少要枚举到p−1p - 1或者q−1q - 1的最大因子

但是如果p−1p-1或者q−1q - 1的倍数被暴露出来的话，依然可以使用Pollard’s p-1算法的后半部分进行分解，比如2021羊城杯的Easy RSA

看题目，首先（a,b<220a, b < 2^{20}）

p=2ga+1q=2gb+1p = 2 g a + 1 \\ q = 2 g b + 1

那么

n=pq=(2ga+1)(2gb+1)=4g2ab+2g(a+b)+1=2g(2gab+a+b)+1\begin{aligned} n &= p q \\ &= (2 g a + 1) (2 g b + 1) \\ &= 4 g^2 a b + 2 g (a + b) + 1 \\ &= 2 g (2 g a b + a + b) + 1 \end{aligned}

按题目中令

h=2gab+a+bh = 2 g a b + a + b

的话，就是

n−1=2ghn - 1 = 2 g h

观察题目的参数，a,b<220a, b < 2^{20}可以枚举，所以可以枚举aa，然后算

a(n−1)=2gah=(p−1)ha (n - 1) = 2 g a h = (p - 1) h

就可以得到一个p−1p - 1的倍数，不妨令

L=a(n−1)L = a (n - 1)

那么就有

p−1∣Lp - 1 \mid L

所以通过计算

p=gcd(aL−1(modn),n)p = gcd(a^L - 1 \pmod n, n)

即可分解

参考代码

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
#!/usr/bin/env sage
n = 84236796025318186855187782611491334781897277899439717384242559751095347166978304126358295609924321812851255222430530001043539925782811895605398187299748256080526691975084042025794113521587064616352833904856626744098904922117855866813505228134381046907659080078950018430266048447119221001098505107823645953039
e = 58337
c = 13646200911032594651110040891135783560995665642049282201695300382255436792102048169200570930229947213493204600006876822744757042959653203573780257603577712302687497959686258542388622714078571068849217323703865310256200818493894194213812410547780002879351619924848073893321472704218227047519748394961963394668

# a, b < 2^20
# h = pb + a = qa + b
# n = 2gh + 1
# (n-1)/2 = gh
# p - 1 = 2ag | a (n-1)
from tqdm import tqdm
W = pow(2, n-1, n)
w = 1
for a in tqdm(range(1, 2^20)):
  w = w * W % n
  p = Integer(gcd(n, w-1))
  if p == 1:
    continue
  q = n // p
  if p * q == n:
    print(a)
    print('p = %d' % p)
    print('q = %d' % q)
    break

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

William’s p+1·

William’s p+1算法是类似的思路，如果找到一个RR满足

p+1∣Rp + 1 \mid R

则可实现分解

William’s p+1并没有Pollard’s p-1直接使用费马小定理那么直观，证明有点复杂，所以下面理论部分看不懂的话可以先略过

理论部分·

证明部分我就直接抄William的书上的证明了

Part.1: Lucas Function·

William’s p+1算法依赖于Lucas Function这个工具

Lucas Function形如

f=x2−Px+Qf = x^2 - P x + Q

令

α=P+P2−4Q2β=P−P2−4Q2\alpha = \frac{P + \sqrt{P^2 - 4Q}}{2} \\ \beta = \frac{P - \sqrt{P^2 - 4Q}}{2} \\

分别为方程f=0f = 0的两个根

那么根据二次方程的性质，可以得到判别式为

D=P2−4Q=(α−β)2D = P^2 - 4 Q = (\alpha - \beta)^2

不妨令

δ=α−β=D\delta = \alpha - \beta = \sqrt{D}

那么可以进一步化简为

α=P+δ2β=P−δ2\alpha = \frac{P + \delta}{2} \\ \beta = \frac{P - \delta}{2}

另外由韦达定理（Vieta’s formulas）可得

P=α+βQ=αβ\begin{aligned} P &= \alpha + \beta \\ Q &= \alpha \beta \end{aligned}

实际使用的时候，用到的其实是Lucas Sequence，用Binet form套一下，令

Un=αn−βnα−βVn=αn+βn\begin{aligned} U_n &= \frac{\alpha^n - \beta^n}{\alpha - \beta} \\ V_n &= \alpha^n + \beta^n \end{aligned}

就可以得到两个序列

U0=0U1=1Un=PUn−1−QUn−2\begin{aligned} U_0 &= 0 \\ U_1 &= 1 \\ U_n &= P U_{n-1} - Q U_{n-2} \end{aligned}

和

V0=2V1=PVn=PVn−1−QVn−2\begin{aligned} V_0 &= 2 \\ V_1 &= P \\ V_n &= P V_{n-1} - Q V_{n-2} \end{aligned}

PS：这东西和Fibonacci Sequence有点类似，不了解的可以先去翻一下Fibonacci Sequence的证明

Part.2: U_n和V_n的性质·

暂时先不管Lucas Sequence，先看看UnU_n和VnV_n的一些性质，首先需要用到这个

Un的证明·

首先看UnU_n，对UnU_n乘上2nδ2^n \delta可得（注：δ=α−β\delta = \alpha - \beta）

2nδUn=2n(α−β)αn−βnα−β=2n(αn−βn)=(2α)n−(2β)n\begin{aligned} 2^n \delta U_n &= 2^n (\alpha - \beta) \frac{\alpha^n - \beta^n}{\alpha - \beta} \\ &= 2^n (\alpha^n - \beta^n) \\ &= (2\alpha)^n - (2\beta)^n \end{aligned}

由于

2α=2P+δ2=P+δ2β=2P−δ2=P−β2 \alpha = 2 \frac{P + \delta}{2} = P + \delta \\ 2 \beta = 2 \frac{P - \delta}{2} = P - \beta \\

所以

2nδUn=(P+δ)n−(P−δ)n\begin{aligned} 2^n \delta U_n &= (P + \delta)^n - (P - \delta)^n \end{aligned}

二项式定理展开一下

2nδUn=∑k=0nCnkPn−kδk−∑k=0nCnkPn−k(−δ)k=∑k=0nCnkPn−k(δk−(−δ)k)\begin{aligned} 2^n \delta U_n &= \sum_{k=0}^{n} C_n^k P^{n-k} \delta^k - \sum_{k=0}^{n} C_n^k P^{n-k} (-\delta)^k \\ &= \sum_{k=0}^{n} C_n^k P^{n-k} (\delta^k - (-\delta)^k) \end{aligned}

显然根据kk的奇偶性，有

δk−(−δ)k={0, k is even2δk, k is odd\delta^k - (-\delta)^k = \begin{cases} 0 &,\ \text{k is even} \\ 2\delta^k &,\ \text{k is odd} \end{cases}

所以

2nδUn=2∑k is oddnCnkPn−kδk\begin{aligned} 2^n \delta U_n &= 2 \sum_{k\text{ is odd}}^{n} C_n^k P^{n-k} \delta^k \end{aligned}

这样有点难看，不妨令

k=2k′+1k = 2 k' + 1

就可以得到

2nδUn=2∑k′=0(n−1)/2Cn2k′+1Pn−2k′−1δ2k′+1=2∑k′=0(n−1)/2Cn2k′+1Pn−2k′−1δ2k′⋅δ=2δ∑k′=0(n−1)/2Cn2k′+1Pn−2k′−1δ2k′=2δ∑k′=0(n−1)/2Cn2k′+1Pn−2k′−1Dk′\begin{aligned} 2^n \delta U_n &= 2 \sum_{k'=0}^{(n-1)/2} C_n^{2k'+1} P^{n-2k'-1} \delta^{2k'+1} \\ &= 2 \sum_{k'=0}^{(n-1)/2} C_n^{2k'+1} P^{n-2k'-1} \delta^{2k'} \cdot \delta \\ &= 2 \delta \sum_{k'=0}^{(n-1)/2} C_n^{2k'+1} P^{n-2k'-1} \delta^{2k'} \\ &= 2 \delta \sum_{k'=0}^{(n-1)/2} C_n^{2k'+1} P^{n-2k'-1} D^{k'} \\ \end{aligned}

左右消去2δ2 \delta，就是

2n−1Un=∑k=0(n−1)/2Cn2k+1Pn−2k−1Dk\begin{aligned} 2^{n-1} U_n &= \sum_{k=0}^{(n-1)/2} C_n^{2k+1} P^{n-2k-1} D^k \end{aligned}

Vn的证明·

VnVn也是类似的，对VnV_n乘上2n2^n，可得

2nVn=2n(αn+βn)=(2α)n+(2β)n\begin{aligned} 2^n V_n &= 2^n (\alpha^n + \beta^n) \\ &= (2\alpha)^n + (2\beta)^n \end{aligned}

二项式定理展开

2nVn=∑k=0nCnkPn−kδk+∑k=0nCnkPn−k(−δ)k=∑k=0nCnkPn−k(δk+(−δ)k)\begin{aligned} 2^n V_n &= \sum_{k=0}^{n} C_n^k P^{n-k} \delta^k + \sum_{k=0}^{n} C_n^k P^{n-k} (-\delta)^k \\ &= \sum_{k=0}^{n} C_n^k P^{n-k} (\delta^k + (-\delta)^k) \end{aligned}

这里反过来，当kk是奇数的时候可以抵消，即

δk+(−δ)k={0, k is odd2δk, k is even\delta^k + (-\delta)^k = \begin{cases} 0 &,\ \text{k is odd} \\ 2\delta^k &,\ \text{k is even} \end{cases}

所以

2nVn=2∑k is evennCnkPn−kδk\begin{aligned} 2^n V_n &= 2 \sum_{k\text{ is even}}^{n} C_n^k P^{n-k} \delta^k \end{aligned}

不妨令

k=2k′k = 2 k'

就可以得到

2nVn=2∑k′=0n/2Cn2k′Pn−2k′δ2k′=2∑k′=0n/2Cn2k′Pn−2k′Dk′\begin{aligned} 2^n V_n &= 2 \sum_{k'=0}^{n/2} C_n^{2k'} P^{n-2k'} \delta^{2k'} \\ &= 2 \sum_{k'=0}^{n/2} C_n^{2k'} P^{n-2k'} D^{k'} \\ \end{aligned}

两边消去22，就是

2n−1Vn=∑k=0n/2Cn2kPn−2kDk\begin{aligned} 2^{n-1} V_n &= \sum_{k=0}^{n/2} C_n^{2k} P^{n-2k} D^k \end{aligned}

Part.3: U_p和V_p的性质·

第二个需要用到的性质是这个

其中ϵp\epsilon_p是Legendre symbol，或者说Jacobi symbol，大概意思是，如果D(modp)D \pmod p是二次剩余，则ϵp=1\epsilon_p = 1，如果D(modp)D \pmod p是非二次剩余，则ϵp=−1\epsilon_p = -1

Up的证明·

首先根据上一Part的证明，有

2p−1Up=∑k=0(p−1)/2Cp2k+1Pp−2k−1Dk\begin{aligned} 2^{p-1} U_p &= \sum_{k=0}^{(p-1)/2} C_p^{2k+1} P^{p-2k-1} D^k \end{aligned}

根据排列组合的性质

Cp2k+1=p!(2k+1)!(p−(2k+1))!C_p^{2k+1} = \frac{p!}{(2k+1)! (p - (2k+1))!}

所以除非p=2k+1p = 2k + 1（即k=(p−1)/2k = (p-1) / 2），或2k+1=02k+1 = 0（显然这个不会成立），否则Cp2k+1C_p^{2k+1}都含因子pp，即

Cp2k+1≡{0(modp), p≠2k+11(modp), p=2k+1C_p^{2k+1} \equiv \begin{cases} 0 \pmod p &,\ p \ne 2k + 1 \\ 1 \pmod p &,\ p = 2k + 1 \end{cases}

而根据费马小定理

2p−1≡1(modp)2^{p-1} \equiv 1 \pmod p

所以综合可得

Up≡1⋅P0⋅D(p−1)/2≡D(p−1)/2≡ϵp(modp)U_p \equiv 1 \cdot P^{0} \cdot D^{(p-1)/2} \equiv D^{(p-1)/2} \equiv \epsilon_p \pmod p

其中最后一个等号根据欧拉准则可得

Vp的证明·

VpV_p的证明类似，首先

2p−1Vp=∑k=0p/2Cp2kPp−2kDk\begin{aligned} 2^{p-1} V_p &= \sum_{k=0}^{p/2} C_p^{2k} P^{p-2k} D^k \end{aligned}

然后

Cp2k=p!(2k)!(p−(2k))!C_p^{2k} = \frac{p!}{(2k)! (p - (2k))!}

所以除非k=0k = 0或p=2kp = 2k（奇素数不成立），否则Cp2kC_p^{2k}都含因子pp，即

Cp2k≡{0(modp), k≠01(modp), k=0C_p^{2k} \equiv \begin{cases} 0 \pmod p &,\ k \ne 0 \\ 1 \pmod p &,\ k = 0 \end{cases}

最后综合可得

Vp≡1⋅Pp⋅D0≡Pp−1⋅P≡P(modp)V_p \equiv 1 \cdot P^{p} \cdot D^0 \equiv P^{p-1} \cdot P \equiv P \pmod p

Part.4: V_{p+1}的性质·

接下来需要使用到这个

因为William’s p+1只用了VV来计算，所以只需要关注下面的2.11

证明也不难，从右边算起，就是

VmVn+DUmUn=(αm+βm)(αn+βn)+(α−β)2αm−βmα−βαn−βnα−β=(αm+βm)(αn+βn)+(α−β)2(αm−βm)(αn−βn)(α−β)2=(αm+βm)(αn+βn)+(αm−βm)(αn−βn)=(αmαn+βmβn)+(αmβn+βmαn)+(αmαn+βmβn)−(αmβn+βmαn)=2(αmαn+βmβn)=2(αm+n+βm+n)=2Vm+n\begin{aligned} V_m V_n + D U_m U_n &= (\alpha^m + \beta^m) (\alpha^n + \beta^n) + (\alpha - \beta)^2 \frac{\alpha^m - \beta^m}{\alpha - \beta} \frac{\alpha^n - \beta^n}{\alpha - \beta} \\ &= (\alpha^m + \beta^m) (\alpha^n + \beta^n) + (\alpha - \beta)^2 \frac{(\alpha^m - \beta^m) (\alpha^n - \beta^n)}{(\alpha - \beta)^2} \\ &= (\alpha^m + \beta^m) (\alpha^n + \beta^n) + (\alpha^m - \beta^m) (\alpha^n - \beta^n) \\ &= (\alpha^m \alpha^n + \beta^m \beta^n) + (\alpha^m \beta^n + \beta^m \alpha^n) + (\alpha^m \alpha^n + \beta^m \beta^n) - (\alpha^m \beta^n + \beta^m \alpha^n) \\ &= 2(\alpha^m \alpha^n + \beta^m \beta^n) \\ &= 2 (\alpha^{m+n} + \beta^{m+n}) \\ &= 2 V_{m+n} \end{aligned}

由前文可知

V1=PU1=1Vp≡P(modp)Up≡ϵp(modp)\begin{aligned} V_1 &= P \\ U_1 &= 1 \\ V_p &\equiv P \pmod p \\ U_p &\equiv \epsilon_p \pmod p \\ \end{aligned}

所以如果把m=pm = p和n=1n = 1代进去的话就是

2Vp+1≡VpV1+DUpU1≡P2+(P2−4Q)ϵp≡(1+ϵp)P2−4ϵpQ(modp)\begin{aligned} 2 V_{p+1} &\equiv V_p V_1 + D U_p U_1 \\ &\equiv P^2 + (P^2 - 4Q) \epsilon_p \\ &\equiv (1+\epsilon_p) P^2 - 4 \epsilon_p Q \pmod p \end{aligned}

如果ϵp=−1\epsilon_p = -1，也即D(modp)D \pmod p是非二次剩余的话，就是

2Vp+1≡4Q(modp)Vp+1≡2Q(modp)\begin{aligned} 2 V_{p+1} &\equiv 4 Q\pmod p \\ V_{p+1} &\equiv 2 Q\pmod p \end{aligned}

这个就是Lemma 15的Case2

PS：根据这个Lemma 15也可以知道，William’s p+1算法其实是可以算p-1的，只不过速度会比Pollard’s p-1慢

Part.5: V_{mn}的性质·

最后需要用到这个

只需要用到VmnV_{mn}的性质，UmnU_{mn}的有兴趣可以自己看书

证明也不难，首先（令α0\alpha_0、β0\beta_0为以PP和QQ为参数的根）

Vmn(P,Q)=α0mn+β0mn=(α0m)n+(β0m)nV_{mn}(P, Q) = \alpha_0^{mn} + \beta_0^{mn} = (\alpha_0^m)^n + (\beta_0^m)^n

令

α1=α0mβ1=β0m\begin{aligned} \alpha_1 &= \alpha_0^m \\ \beta_1 &= \beta_0^m \\ \end{aligned}

稍微算一下可得

α1+β1=α0+β0=Vmα1β1=(α0β0)m=Qm\begin{aligned} \alpha_1 + \beta_1 &= \alpha_0 + \beta_0 = V_m \\ \alpha_1 \beta_1 &= (\alpha_0 \beta_0)^m = Q^m \end{aligned}

由韦达定理反推回去就得到α1\alpha_1、β1\beta_1是以VmV_m和QmQ^m为参数的根，即

Vmn(P,Q)=Vn(Vm,Qm)V_{mn}(P, Q) = V_n(V_m, Q^m)

上一部分中已经推出来了

Vp+1≡2Q(modp)V_{p+1} \equiv 2 Q\pmod p

不妨令m=p+1m = p+1，就得到

Vn(p+1)≡Vn(Vp+1,Qp+1)≡Vn(2Q,Q2)(modp)V_{n(p+1)} \equiv V_n(V_{p+1}, Q^{p+1}) \equiv V_n(2Q, Q^2) \pmod p

令α2\alpha_2、β2\beta_2为以2Q2Q和Q2Q^2为参数的根，由韦达定理可知

{2Q=α2+β2Q2=α2β2\begin{cases} 2Q = \alpha_2 + \beta_2 \\ Q^2 = \alpha_2 \beta_2 \end{cases}

一眼顶针就可以看到

α2=β2=Q\alpha_2 = \beta_2 = Q

于是

Vn(p+1)(P,Q)≡Vn(2Q,Q2)≡α2n+β2n≡Qn+Qn≡2Qn(modp)\begin{aligned} V_{n(p+1)}(P, Q) &\equiv V_n(2Q, Q^2) \\ &\equiv \alpha_2^n + \beta_2^n \\ &\equiv Q^n + Q^n \\ &\equiv 2 Q^n \pmod p \end{aligned}

这里的n(p+1)n(p+1)就是我们想要找的RR，而实际应用中通常会设Q=1Q = 1，于是就

VR≡2(modp)V_R \equiv 2 \pmod p

Part.6: Lucas Sequence·

由上一部分的结论可得，如果知道RR满足p+1∣Rp + 1 \mid R和n=pqn = p q，通过计算

p=gcd(VR−2(modn),n)p = gcd(V_R - 2 \pmod n, n)

即可实现分解

剩下的问题是，如何快速地算得VRV_R

代进VnV_n就是

V0=2V1=PVn=PVn−1−QVn−2\begin{aligned} V_0 &= 2 \\ V_1 &= P \\ V_n &= P V_{n-1} - Q V_{n-2} \end{aligned}

书上有介绍一种快速计算的方法，但按我习惯我会用线性代数的方法去计算，即

[P−Q10][Vn−1Vn−2]=[VnVn−1]\begin{bmatrix} P & -Q \\ 1 & 0 \end{bmatrix} \begin{bmatrix} V_{n-1} \\ V_{n-2} \end{bmatrix} = \begin{bmatrix} V_{n} \\ V_{n-1} \end{bmatrix}

然后就是

[P−Q10]n−1[V1V0]=[VnVn−1]\begin{bmatrix} P & -Q \\ 1 & 0 \end{bmatrix}^{n-1} \begin{bmatrix} V_{1} \\ V_{0} \end{bmatrix} = \begin{bmatrix} V_{n} \\ V_{n-1} \end{bmatrix}

这个和斐波那契的快速计算式方法是一样的

2023 N1CTF的e20k·

题目：https://www.nssctf.cn/problem/4641

普通的William’s p+1的做法是，当p+1p+1是平滑数时，通过枚举

R=1×2×3⋯R = 1 \times 2 \times 3 \cdots

然后计算

p=gcd(VR−2(modn),n)p = gcd(V_R - 2 \pmod n, n)

和Pollard’s p-1类似，如果p+1p + 1或q+1q + 1的倍数被暴露，则也可以把这个倍数看作RR直接计算

比如这题，通过审题可以知道

N=PQ=p(2q2−q)=pq(2q−1)N = P Q = p (2q^2 - q) = p q (2 q - 1)

令

r=2q−1r = 2 q - 1

那么

N=pqrN = p q r

且显然

r+1=2q∣(2q)pr=2Nr + 1 = 2 q \mid (2q) p r = 2N

于是就可以通过计算

r=gcd(V2N−2(modn),n)r = gcd(V_{2N} - 2 \pmod n, n)

分解出rr，然后

q=r+12p=Nqr\begin{aligned} q &= \frac{r+1}{2} \\ p &= \frac{N}{q r} \end{aligned}

实现NN的分解

参考代码

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
#!/usr/bin/env sage
import os
os.environ['TERM'] = 'xterm-256color'
from pwn import remote

def V(k, N, A=5):
  M = matrix(Zmod(N), [
    [A, -1],
    [1, 0]
  ])
  v = vector(Zmod(N), [A, 2])
  vk = M^k * v
  return Integer(vk[1])

def factorN(N):
  A = 5
  while True:
    v2n = V(2*N, N, A)
    g = gcd(v2n-2, N)
    print('[Log] %d - %d' % (A, g))
    if g.nbits() >= 256 and is_prime(g):
      r = g
      q = (r+1) // 2
      p = N // (q * r)
      if p * q * r == N:
        return p, q, r
    A = next_prime(A)

while True:
  #context.log_level = 'debug'
  rm = remote('node4.anna.nssctf.cn', 28738)

  rm.recvuntil(b'N = ')
  N = Integer(rm.recvline()[:-1].decode())
  print('N = %d' % N)
  fN = factorN(N)
  print('fN = %s' % str(fN))

  rm.interactive()
  rm.close()

Bach-Shallit·

令ϕk\phi_k为第kk个分圆多项式，在知道ϕk(p)\phi_k(p)的倍数或者ϕk(p)\phi_k(p)是Smooth的情况下，可以使用Bach-Shallit算法分解

这个算法的资料有点少，主要还是得啃Bach和Shallit的论文，这文章省略的东西有点多，得等我把数论基础补完再回来填坑了（逃

另外可以参考@algellar的代码，这应该算一个模板代码了

ImaginaryCTF 2023的Sus·

题目：https://github.com/maple3142/My-CTF-Challenges/tree/master/ImaginaryCTF%202023/Sus

@Maple的题目，上面的wp其实也挺详细的，可以直接看

这题不用完整地跑完Bach-Shallit算法，所以难度降低了不少

看题，现在知道

n=pqrq=p2+p+1\begin{aligned} n &= p q r \\ q &= p^2 + p + 1 \end{aligned}

其中，因为

ϕ3(x)=x2+x+1\phi_3(x) = x^2 + x + 1

所以

q=ϕ3(p)q = \phi_3(p)

也即nn自身是ϕ3(p)\phi_3(p)的一个倍数，刚好符合Bach-Shallit的条件

首先随机选一个不可约的多项式f(x)f(x)生成环

Rp=(Z/pZ)[x]f(x)R_p = \frac{(\mathbb{Z}/p\mathbb{Z})[x]}{f(x)}

由于f(x)f(x)不可约，所以RpR_p其实是伽罗瓦域，乘法阶为

p3−1=(p2+p+1)(p−1)p^3 - 1 = (p^2 + p + 1) (p - 1)

随机挑选RpR_p中的元素

t←Rpt \leftarrow R_p

则

tn≡tpr⋅(p2+p+1)≡tpr⋅p3p−1t^{n} \equiv t^{p r \cdot (p^2 + p + 1)} \equiv t^{pr \cdot \frac{p^3}{p-1}}

所以tnt^n是乘法阶为p−1p-1的子环的元素

正常的Bach-Shallit算法到这里还要用伽罗瓦群的知识计算一堆东西，但这里因为这个阶为p−1p-1的子环大概率只包含x0x^0的项（即x1x^1和x2x^2的项的系数为0），也即论文中的Om\mathcal{O}_m，所以可以省去这部分的计算

可以搞个代码验证一下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/usr/bin/env sage
k = 3
phik = cyclotomic_polynomial(k)

p = 1997
Rp_ = PolynomialRing(Zmod(p),'xp')
while True:
  fp = Rp_.random_element(degree=k, monic=True)
  if fp.is_irreducible():
    break
Rp = Rp_.quotient(fp, 'yp')
t = Rp.random_element()
tp = t^(phik(p))
print(tp)

由于我们并不知道pp，所以也不知道RpR_p，于是以上运算应该在环

Rn=(Z/nZ)[x]f(x)R_n = \frac{(\mathbb{Z}/n\mathbb{Z})[x]}{f(x)}

上进行，先随机挑选

t←Rnt \leftarrow R_n

然后计算tn (in Rn)t^n \text{ (in } R_n \text{)}，根据上面的分析，x1x^1的项（或x2x^2）的系数在模pp的情况下为00，即在模nn的情况下应该是pp的倍数，令

tn=c0+c1x1+c2x2 (in Rn)t^n = c_0 + c_1 x^1 + c_2 x^2 \text{ (in } R_n \text{)}

则通过计算

p=gcd(c1,n)p = gcd(c_1, n)

即可实现分解

参考代码

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
#!/usr/bin/env sage
n = 1125214074953003550338693571791155006090796212726975350140792193817691133917160305053542782925680862373280169090301712046464465620409850385467397784321453675396878680853302837289474127359729865584385059201707775238870232263306676727868754652536541637937452062469058451096996211856806586253080405693761350527787379604466148473842686716964601958192702845072731564672276539223958840687948377362736246683236421110649264777630992389514349446404208015326249112846962181797559349629761850980006919766121844428696162839329082145670839314341690501211334703611464066066160436143122967781441535203415038656670887399283874866947000313980542931425158634358276922283935422468847585940180566157146439137197351555650475378438394062212134921921469936079889107953354092227029819250669352732749370070996858744765757449980454966317942024199049138183043402199967786003097786359894608611331234652313102498596516590920508269648305903583314189707679
e = 65537
c = 27126515219921985451218320201366564737456358918573497792847882486241545924393718080635287342203823068993908455514036540227753141033250259348250042460824265354495259080135197893797181975792914836927018186710682244471711855070708553557141164725366015684184788037988219652565179002870519189669615988416860357430127767472945833762628172190228773085208896682176968903038031026206396635685564975604545616032008575709303331751883115339943537730056794403071865003610653521994963115230275035006559472462643936356750299150351321395319301955415098283861947785178475071537482868994223452727527403307442567556712365701010481560424826125138571692894677625407372483041209738234171713326474989489802947311918341152810838321622423351767716922856007838074781678539986694993211304216853828611394630793531337147512240710591162375077547224679647786450310708451590432452046103209161611561606066902404405369379357958777252744114825323084960942810

k = 3
phik = cyclotomic_polynomial(k)
Rn_ = PolynomialRing(Zmod(n),'xn')

while True:
  fn = Rn_.random_element(degree=k, monic=True)
  Rn = Rn_.quotient(fn, 'yn')
  t = Rn.random_element()

  tn = t^n
  tn1 = tn.list()[1]
  p = Integer(gcd(tn1, n))
  if p != 1 and n % p == 0:
    q = phik(p)
    r = n // (p * q)
    assert p * q * r == n
    print('p = %d' % p)
    print('q = %d' % q)
    print('r = %d' % r)
    break

PS：或者直接套@algellar的代码也行

Lenstra分解算法·

之前在山石实习的时候也写过一篇文章讲这个算法，可以直接参考

大概思路是，令

n=pqn = p q

那么就可以使用参数(a,b)(a, b)生成三条曲线EnE_n、EpE_p和EqE_q，这三条曲线都形如

Y2=X3+aX+bY^2 = X^3 + aX + b

只是模数不同

假设可以找到一个数ll，这个数能被EpE_p的群阶整除，而不被EqE_q的阶整除的话，就和Pollard’s p-1类似，令PP为曲线EnE_n上的一个随机点、O\mathcal{O}为无穷远点，有

l⋅P=O (in Ep)l⋅P≠O (in Eq)\begin{aligned} l \cdot P = \mathcal{O} \text{ (in } E_p \text{)} \\ l \cdot P \ne \mathcal{O} \text{ (in } E_q \text{)} \\ \end{aligned}

此时，在计算l⋅P=O (in En)l \cdot P = \mathcal{O} \text{ (in } E_n \text{)}时，会有一步需要计算pp模nn的逆元，而因为pp和nn不互素，所以会求逆失败，触发报错，使用报错信息就可以分解

给个测试代码

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
#!/usr/bin/env sage
import re

bits = 64
p = random_prime(2^bits)
q = random_prime(2^bits)
n = p * q
print('p = %d' % p)
print('q = %d' % q)
print('n = %d' % n)

a = 3
b = 7
Ep = EllipticCurve(Zmod(p), [3, 7])
Eq = EllipticCurve(Zmod(q), [3, 7])
En = EllipticCurve(Zmod(n), [3, 7])
P = Ep.random_element()
Q = Eq.random_element()
xn = crt([Integer(P.xy()[0]), Integer(Q.xy()[0])], [p, q])
yn = crt([Integer(P.xy()[1]), Integer(Q.xy()[1])], [p, q])
N = En([xn, yn])
phip = Ep.order()
phiq = Eq.order()
phin = phip * phiq
print('phip = %d' % phip)
print('phiq = %d' % phiq)
print()

try:
  phip * N
except Exception as e:
  E = e
_, n, p, q = [Integer(_) for _ in re.findall(r'\d+', E.args[0])]
assert n == p * q
print('n = %d' % n)
print('p = %d' % p)
print('q = %d' % q)
print()

try:
  phiq * N
except Exception as e:
  E = e
_, n, p, q = [Integer(_) for _ in re.findall(r'\d+', E.args[0])]
assert n == p * q
print('n = %d' % n)
print('p = %d' % p)
print('q = %d' % q)

2023华为杯的2023华为杯的EC_Party-I-chall·

题目：EC_Party-I-chall.zip

几乎就是直接套Lenstra算法，令φp\varphi_p和φq\varphi_q分别为EpE_p和EqE_q的阶，泄露的leak为

leak=φp⋅φqleak = \varphi_p \cdot \varphi_q

注意如果直接用leak的话，会导致

leak⋅C=O (in Ep)leak⋅C=O (in Eq)\begin{aligned} leak \cdot C = \mathcal{O} \text{ (in } E_p \text{)} \\ leak \cdot C = \mathcal{O} \text{ (in } E_q \text{)} \\ \end{aligned}

就不符合Lenstra的条件

一个技巧是，通过Yafu分解leak可以得到φp\varphi_p或φq\varphi_q的小因子，假设φq\varphi_q的一个小因子为aqa_q，则

leakaq⋅C≠O (in Eq)\frac{leak}{a_q} \cdot C \ne \mathcal{O} \text{ (in } E_q \text{)}

满足Lenstra条件

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#!/usr/bin/env sage
import re
a, b, n, C, leak = [138681122158674534796479818810828100269024674330030901179877002756402543027343312824423418859769980312713625658733, 4989541340743108588577899263469059346332852532421276369038720203527706762720292559751463880310075002363945271507040, 762981334990685089884160169295988791471426441106522959345412318178660817286272606245181160960267776171409174142433857335352402619564485470678152764621235882232914864951345067231483720755544188962798600739631026707678945887174897543, (19591102741441427006422487362547101973286873135330241799412389205281057650306427438686318050682578531286702107543065985988634367524715153650482199099194389191525898366546842016339136884277515665890331906261550080128989942048438965, 728465071542637655949094554469510039681717865811604984652385614821789556549826602178972137405550902004858456181137844771163710123158955524137202319902378503104952106036911634918189377295743976966073577013775200078470659428344462772), 762981334990685089884160169295988791471426441106522959345445792076415993922016249232021560266153453470937452118572318136597282436269660557904217923887981072203978473274822142705255987334355747997513083011853917049784914749699536828]
En = EllipticCurve(Zmod(n), [a, b])
C = En(C)

assert leak % 199 == 0
try:
  (leak // 199) * C
except Exception as e:
  E = e
_, n, p, q = [Integer(_) for _ in re.findall(r'\d+', E.args[0])]
assert n == p * q
print('n = %d' % n)
print('p = %d' % p)
print('q = %d' % q)

参考·

  1. Ballot, Christian J-C., and Hugh C. Williams. The Lucas Sequences: Theory and Applications. Vol. 8. Springer Nature, 2023.
  2. Bach, Eric, and Jeffrey Shallit. “Factoring with cyclotomic polynomials.” Mathematics of Computation 52.185 (1989): 201-219.
  3. https://github.com/algellar/Factoring-with-Cyclotomic-Polynomials
  4. Hoffstein J, Pipher J, Silverman J H, et al. An Introduction to Cryptography[M]. Springer New York, 2014.