第一届强网杯密码数学专项赛——赛题一解题笔记

第一届（2022）强网杯密码数学专项赛，赛题一RSA 公钥密码体制的攻击的WP（AK）。

赛题涉及多种常见的RSA攻击，还是挺有学习价值的。

首先附上我们答辩的Slide：Slide.pdf

和赛题数据：RSA-data.txt

由于我只拿到纸质版的赛题，所以下面先码一下赛题描述的关键部分（也可以看Slide）

赛题描述·

赛题其实是一个KEM-DEM结构的混合加密，DEM部分使用128比特的密钥$K$，KEM部分是一个RSA2048-KEM，对密钥$K$进行加密。假设用户A给用户B发送$M$，其加密过程如下：

• A随机产生128比特密钥$K$，加密$M$得到$C_M$；

• A选定RSA2048公钥$(N, e)$，使用固定算法$\text{padding}$填充$K$至2048比特然后使用$(N, e)$加密：

  $$C_K \equiv (\text{padding}(K))^e \pmod N$$

• A发送$(C_M, C_K)$给B；

以上$N=PQ$，其中$P$和$Q$的产生方法比较特殊，以$P$为例，首先有一固定比特序列$(x_0, x_1, ...)$，有一起点$s$，从$i=0$开始，设$P = (1, x_{s+i}, x_{s+i+1}, ..., x_{s+i+1021}, 1)_2$，并检测$P$的素性，若$P$为素数则生成结束，否则$i$增加$1$再检测。$Q$的产生方法类似，不过其起点记作$t$。

现知道多个用户的公钥$(N_i, e_i)$、对应的起点$(s_i, t_i)$和密文$C_{K_i}$，求$K_i$，即攻击OW-CPA安全性。（赛题数据的6组数据中，每列分别对应：序号、$s$、$t$、$N$、$e$、$C_K$）

WriteUp·

数据一和数据三——求公因子·

首先$s_3 - s_1 = 1011 - 1010 = 1$，按上面的素数生成方法大概率会有$P_1 = P_3$（Slide中用素数密度估算，实际直接试一下就好了），而$t_3 - t_1 = 8746 - 3927$，即大概率$Q_1 \ne Q_3$，假设相同则可求$N_1$和$N_3$的公因子恢复$P_1$和$P_3$:

$$\begin{aligned} &P_1 = P_3 = \text{GCD}(N_1, N_3) = \text{GCD}(P_1 Q_1, P_3 Q_3) \\ &Q_1 = N_1 / P_1 \\ &Q_3 = N_3 / P_3 \end{aligned}$$

参考代码：

# sage
import libnum

n1 = 0x62d9353e32425edd17be2ce4d98a6c0407d309673b2949b55adcc51a0d63d2793f920533ad7a79ad23dbb2d733eecde2c1694d1255dccb1a90e3035ac7b755993866a32a606dcab916589b4eb93b42fb01091b1acef3fc5b65598bfccd1081ed3b8509691f935bad0b6360baac1c0050255b5de55fbbf7950df9ef11ac95033f38a9a571b0bd351a154cca2d6f6e7b89519e25b46fde8bf6fa8907a32eeefd614a44da3878c8e04881260b9fc2588b5b71770184f0b1d58b35363829cc39b146758e18f933ea30f0f5b80f7274a50b0449d4ce0b052e052062e14474186daea9d30e30afb2599e9ec51af67eaa9bd7655c0e7e31f3c66d29d5a7159304782e63
n3 = 0x6a3402d4aa029ac85eee66c7f3ac3f39d308a369b8f785bee365fea1e2ad009a89878437d95f45916e801489b224ec15b5772f7fc0a46ae0be1620420602cbd399b05933892cb528a1de807540734e1338d42caf08f92f7f0b9e480b4a2218d42f43797d43f8e0d4c44e6352b15b41f923b32941ab5725f4bdfa3f9f26dbce1528ba1a5d0c209ef494a651bbb4ff878dc207a265385d418330cf0b696da8461693bd415a17673c40acd68e2c9de6afacf82efd78bae6e43cdcc2e1c0d2d8fda932cf987f8516900b734f41d7c93e1fa4278a0fb8957c9d7900e2596d944d7c11f34649ee51cfb7dac1fdf0a410a9ebe70827bad4708762776b463cc9302e38c1

e1 = 31
e3 = 131

c1 = 0x3c3623571637d93a9a3e9787ea8c8eefa001cd93a2cca69f9e68335142cc5e0c7e6edb87c23196c71820de7979598b9b27a5161196ff3741021bb9d197e0d4380182de125fa914fb76f4efc381a3328e762ada4f95298cbe511009cf5d1849a3c59a31bfc73b3e28d579d652bf0d94a31416fb357957910d05466e50346e1690a3b900d17ec927fd22f0e28693d65f086ce55577564cbc7025dcda5de645051e99c5205b40ed7886b4fb6f8d126c33ba2fdd9195eb60571a14ed5c53feabf6857f8f1a35efefcbd994662bd30dbe9b0b2b0ae22f9d201143f1c742f9fe583185b8fcbc7c2b103629fd5219e9902574a8da871514363e97dc0886cd555c856f5d
c3 = 0x085e59951e4297601786b3b2036d9d6ee9f7da1e24e7c729f3d3b3824b7c8f979a1f56a65eef4223b821ab51362a1b0d42f2652b60a185f635c957475f30d30b1eb2b81c08f8d6da5ab07db2a8a1f59060fc2c1d7d2bcff410ac3a1f83ddb3a2c1e6c4c914acc9ac33d3707ff12ce4b85a17965d4cdec3642c69537652c30a4094f650915f1f01970bac9a0a8bfc4e51943404629786208d90ba15dbefc818ce0a3e0369bcd847d46695e60ef40f79b156ae62a30fc819484353d75d5bbab30c7ce963e7363ea151341bd61c42c419ce703653b8bf2f23486af75f4c3be36bc3b8d7acd4afaacc28bab973dca45933517143487c21b04b9b34484778abc7b5c7

p1 = p3 = gcd(n1, n3)
q1 = n1 // p1
q3 = n3 // p3
assert p1 != 1 and p1*q1 == n1 and p3*q3 == n3

print('p1 = %s' % hex(p1))
print('q1 = %s' % hex(q1))
print('p3 = %s' % hex(p3))
print('q3 = %s' % hex(q3))

d1 = e1.inverse_mod((p1-1)*(q1-1))
d3 = e3.inverse_mod((p3-1)*(q3-1))

m1pad = libnum.n2s(int(pow(c1, d1, n1)))
m3pad = libnum.n2s(int(pow(c3, d3, n3)))
m1 = m1pad[-(128//8):]
m3 = m3pad[-(128//8):]

print('m1 = %s' % hex(libnum.s2n(m1)))
print(m1)
print('m3 = %s' % hex(libnum.s2n(m3)))
print(m3)

'''
p1 = 0x81652a17874b7b8f19a2f3e192b4a1a8303073435e23423310b1a5a6c44aea72dcfcb61226e36a22b08188413d217c186912e9552d4c9803e78d846443b57c5d23173f383bfa0d418d4b8d0fbb73a12116d04231b3413a1b20f794d70f88ea7ca149c95635b147fb487e0d3f30a90692fc917a7313934979e48e55089ea64561
q1 = 0xc390b88a439fa5392c9b3adc8ff8482cde86781e7edf8652078abae7fe2a081284f7614d535d3217a7af3087689328bd460ddc69ded3e3ef14889dd6e709da4c85539e59914384b57451f614791ed976a370523dcf265f3fb292c219beeb5921fc05c47882919839f73aca9df0f250b7f8350b20f43bdcf71287c079fe2bc643
p3 = 0x81652a17874b7b8f19a2f3e192b4a1a8303073435e23423310b1a5a6c44aea72dcfcb61226e36a22b08188413d217c186912e9552d4c9803e78d846443b57c5d23173f383bfa0d418d4b8d0fbb73a12116d04231b3413a1b20f794d70f88ea7ca149c95635b147fb487e0d3f30a90692fc917a7313934979e48e55089ea64561
q3 = 0xd21db9a921f269df2e7e1313f43055c27795621c555a82f729d0d76e5c2324411127c07cf85f8a5773535fc81511545732241d58be5ca9b471937577d49cc67ef4745649e96de30d2000a495039f52b470ca9a5669dc6284d923d62645736ef1a8ac52310578bba12f3aebe324d45f6aeb035f3991997bfc54e3683c34e54f61
m1 = 0x56eebe53d69efebd4505540305375f07
b'V\xee\xbeS\xd6\x9e\xfe\xbdE\x05T\x03\x057_\x07'
m3 = 0x7cd11086cad330d2cbbe4fc7e59ee2e5
b'|\xd1\x10\x86\xca\xd30\xd2\xcb\xbeO\xc7\xe5\x9e\xe2\xe5'
'''

填充去除·

以上恢复出$P_1$和$P_3$后可以解密出$\text{padding}(K_1)$和$\text{padding}(K_3)$，观察得出算法$\text{padding}(K)$实际为$15$个$K$相接：

$$\text{padding}(K) := 0^{128} \| \underbrace{K \| K \| \ldots \| K}_{\text{15 $K$s fused}}$$

这样的表述不容易利用，所以转个数学表达（可以类比进制转换）：

$$\text{padding}(K) := K \cdot \sum_{i=0}^{15-1} 2^{128 \cdot i}$$

然后，令$C_K \equiv \text{padding}(K)^e \pmod N$，定义以下$\text{unpad}_{N, e}$操作在不解密的情况下去除密文的padding：

$$\begin{aligned} \text{unpad}_{N, e}(C_K) :&= C_K \cdot (\sum_{i=0}^{15-1} 2^{128 \cdot i})^{-e} \\ &\equiv K^e \cdot (\sum_{i=0}^{15-1} 2^{128 \cdot i})^{e} \cdot (\sum_{i=0}^{15-1} 2^{128 \cdot i})^{-e} \\ &\equiv K^e \pmod N \end{aligned}$$

数据5——小加密指数·

首先若$K_i^{e_i} < N_i$，则$K_i^{e_i} \pmod {N_i} = K_i^{e_i}$，取余操作可忽略。

观察数据5有：

$$K_5^{e_5} < 2^{9 \cdot 128} < 2^{2047} < N_5$$

所以可以对$C_{K_5}$$\text{unpad}$后开$e_5$次方解密：

$$\begin{aligned} \sqrt[e_5]{\text{unpad}_{N_5, e_5}(C_{K_5})} &= \sqrt[e_5]{K_5^{e_5} \pmod {N_5}} \\ &= \sqrt[e_5]{K_5^{e_5}} \\ &= K_5 \end{aligned}$$

参考代码：

# sage
import libnum

n5 = 0x87a97539530147fecfc9ca4d2ffe847d72796d81e995378bc607ca6daa98eb24112bb86c64fc8982424e055dd6e7b9b88994be1f60da77fd4fc8b2b7ec065437f864532d92c3a303a79ed3db69dda2d94e626188750b6dc819b567810074de31254dd3883d0ef7593b84ec60677dd666953dff23f4b1c88b9d909a2fbda94fc13c51644efcc19d2f581b828e68d6f53c323896b2dbf3c04768f909098bf20a91e41024d869ad3027e2a7ba7f1e6288de4a581852db75988897aaaa6a7dd17fbe85a2cd6efa659d0efa5e1072df457611f8faff7ee34725e747542059338ee5daba0a2934b1cecfc499c3ae3e2e4ae17d989cd2d1f1e5f99d284f74168c2a45df
c5 = 0x68cdd0dea3b194e90557e482d6af1a1a978d7baebd66c9ba308fdc411271fed18202b0d03e4bd73dfe4f905004aad7df11c55b6ef465988edef2e8ecbdb6703be7f4e0517c84faa79c160d6306ac66cc72f75e45f70c3fc717dec6d4658d48a74c385eb0a6a3b650651840171408a3cbc3bd1ce86226ba9aff86af17a4eb137ae1223f15abda51bb832cb872db55b0f594f39a9ee11d4b49f583bd917f585f6ffbdf5920439403824338f2e1be02e019451109a8a559b408ef927ea8f5374d4e026daa228deb687e18ea2aaee9f4acfd77cd60aeae25a1d9334711af70a686a29bf89a46a461f1307659a57177a593a13f8465bbb4f6c05ec728aa5a2ea0f125
e5 = 9

padding = sum([2^(128*i) for i in range(15)])
c5np = int(c5 * pow(padding.inverse_mod(n5), e5, n5) % n5)

m5 = c5np^(1/e5)
print('m5 = %s' % hex(m5))
m5 = libnum.n2s(int(m5))
print(m5)

'''
m5 = 0x160360acc4078a0b5d46e3860255d30a
b'\x16\x03`\xac\xc4\x07\x8a\x0b]F\xe3\x86\x02U\xd3\n'
'''

数据2——小加密指数扩展·

$K_i^{e_i} > N_i$时虽然不可以用上面方法解密， 但若$K_i$可被分解为$K_i = K_{i1} \cdot K_{i2}$，其中$K_{i1}$为可枚举的小因子、$K_{i2}^{e_i} < N_i$，则可以通过枚举$K_{i1}$把$K_i$降为$K_{i2}$再用上面方法解密：

$$\begin{aligned} \sqrt[e_i]{\text{unpad}_{N_i, e_i}(C_{K_i}) \cdot K_{i1}^{-e_i}} &= \sqrt[e_i]{(K_i \cdot K_{i1}^{-1})^{e_i} \pmod {N_i}} \\ &= \sqrt[e_i]{(K_{i2})^{e_i}} \\ &= K_{i2} \end{aligned}$$

稍微计算一下，$e_2 = 17$，设$K_{i1}$为$x$比特，即只要：

$$(128-x)·17 < 2048 \\ x \ge 8 > 128 - \frac{2048}{17}$$

即可。由于赛题的$K$随机选取，所以$K_2$大概率含有8比特以上的小因子，最终枚举出最小符合条件因子为$K_{21} = 477$。

参考代码：

# sage
import libnum

n2 = 0xcb4539eb70dd8966c1adeb86826e5607a2a98cf138646b3a054e57d4d28c2a992499fe87cb1af2a7906777bd84fc3dbd22fddeb61567891376c4f436015cf333b97cde79becb9d89610cc2e545bb354ea08479e7debb0b7a002c23e833cadae227feb6b75de74e74b87bfd2f7f20849e94dabe687f663c66dc17e737ed493af862ab40e58ae40265b97edb5a1d326b7b407f31691f79a72ed2238ef2219525ead1ba0b0a9ecde6f818f105c743258eec9fdb51cd6a16b99cd2d1ec2dd2bc6ed99ebbaca2fc3aabe709a7d82619cf524637bac7c9933e50870b60ce402b4404d2b4414387d4d857ab2f528d7c96e11dbc27ba6f531b748655f36f375280261537
c2 = 0x6e2a4e7fcb266e585e2f4f31a83440c52a1ba8dbb091c9776891b7e37b3bae4b684a932a4cc69c747a7ef67de6cef09297844fa4d3ea82ab641233dc4e33fdaac096369df91c49bd4917d33fc84d32ae6e0d4f31612a82612a81d43b723dc570da950500ce294cf336d7a27ace2f3427db3b74e8a98415e0375503f1bfca219fa60a0a533241d75620f71e6b51d756e3763ceea7cdec8fe3a4bcdefd1364a6b63bd97595cfcee254768518323b593ac3559af80351c8b5e5d4841115548a6069b4812aa3d718a7d852eafdd96fb6e62d4c05640c177f9ea52cd3cbfa7fa1755dba44e649d6135f552c549f573addfc4fa23d61909e53d3cf4083fe7ee315fc0c
e2 = 17

padding = sum([2^(128*i) for i in range(15)])
c2np = c2 * pow(padding.inverse_mod(n2), e2, n2) % n2

# enum, let c2np = (ab)^e % n
for b in range(2^8, 2^10):
  try:
    b = Integer(b)
    c2a = int(c2np * pow(b.inverse_mod(n2), e2, n2) % n2)

    m2a = c2a^(1/e2)
    
    m2 = m2a*b
    if pow(m2, e2) % n2 == c2np:
      print('b = %s' % b)

      print('m2 = %s' % hex(m2))
      m2 = libnum.n2s(int(m2))
      print(m2)
      break
  except:
    pass

数据4——素因子高位泄露·

由于$s_4 − t_3 = 511$，所以存在$Q_3$低位泄露$P_4$一半以上比特的可能性，如果成立，则可以用“P高位泄露攻击”分解$N_4$，下图摘抄自[Gal12]：

大概思路可以参考之前的Coppersmith学习笔记。

由于该攻击需要知道具体的比特数，所以还需要枚举$Q_3$低位，最终枚举得知$Q_3$低位600比特为$P_4$高位。

参考代码：

# sage
# beta: https://doc.sagemath.org/html/en/reference/polynomial_rings/sage/rings/polynomial/polynomial_modn_dense_ntl.html
# compute a root mod b where b is a factor of n and b >= n^beta (b = p or q, so set 0.49)

def solve(n, ph, pl=1, pbits=1024):
  hbits = ph.nbits()
  lbits = pl.nbits()
  PR.<x> = PolynomialRing(Zmod(n))
  f = ph * 2^(pbits-hbits) + x * 2^lbits + pl
  f = f.monic()
  roots = f.small_roots(X=2^(pbits-hbits-lbits+1), beta=0.4)
  if roots:        
    pm = Integer(roots[0])
    p = ph * 2^(pbits-hbits) + pm * 2^lbits + pl
    if n % p == 0:
      q = n // p
      return p, q
  return None

n4 = 0x7a52ce9f4e62bc8992606cdfd4c085d8e92ea0b7f15f5f2dbae0a23a9781773541db723d53e0039d759c4485a8c6cda817211efa294ca422247e663269eaef819e92a49a7812c62d80ff0e21ce3e43cbd7ea14b92f0d426ce45effb6a8d13c746acb7246605a0a8750e79934aa13a1cf0778ce27cc2ffefdb5b745df4b80de310c78185164efd9703e8790c50897117cd6a47ec213f2eb0e6008e750c0e89d3e24361d425abe304549b31e47da184fc8a3294e1f70fb225fc329bdc7d0758f6c5c672d1596b8fbffb7fa0b485d24178e3303d82b3e9a5e8c359a65cffa17f481bef03154301ee0cdcc2f2674cc9fcfaad5b0e6c1363e9adced2241e6b95f23f1
q3 = 0xd21db9a921f269df2e7e1313f43055c27795621c555a82f729d0d76e5c2324411127c07cf85f8a5773535fc81511545732241d58be5ca9b471937577d49cc67ef4745649e96de30d2000a495039f52b470ca9a5669dc6284d923d62645736ef1a8ac52310578bba12f3aebe324d45f6aeb035f3991997bfc54e3683c34e54f61
q3r = bin(q3)[3:-1]
for i in range(550, 1024):
  ph = Integer(int('1'+q3r[-i:], 2))
  res = solve(n4, ph)
  if res == None:
    pass
    #print('%d: Failed.' % i)
  else:
    print('%d: Yes!' % i)
    p4, q4 = res
    assert p4 * q4 == n4
    print('p4 = %s' % hex(p4))
    print('q4 = %s' % hex(q4))
    break

e4 = 253
c4 = 0x51b09ca05d1bc7ecd616642031c1cbf0ed9e49ec3c605e2c2c1fef37ead042c50e7099fbaeedb2c79240842ca21c624d904df9e24a7c28313fa40af51071e1b9aad8039a3973d0c5a68f36438886c07e011c21937a596644c2df3de76da4fe29c1fde91a2af64b966928d0c59065c7df615e9a1362901cb265be02441d9501e511f1c4d77a2047d5551b51d0f696c371afb0a9241f5634baafdbae9bae81bb5a39c2fb9a20377ee7b0bc0ccd135872d2b93cc6f0ce1bc156074627f354edf6c484cabc026e03acd7c6d0fc41e70beca6aa3972c6d489aa76c09506a2002ec73c51e0d4eab97a74cd737126c500a85b39fd05034a7e83e09337233cd78c32d71f

import libnum
phi4 = (p4-1) * (q4-1)
d4 = e4.inverse_mod(phi4)
m4pad = libnum.n2s(int(pow(c4, d4, n4)))
m4 = m4pad[-(128//8):]
print('m4 = %s' % hex(libnum.s2n(m4)))
print(m4)

'''
600: Yes!
p4 = 0x972a6d1c64dd5df527319fbd1d15927a5b78c34800292540e7d4ad1c32a6959a7718a13648f589915cdbbc6a2b148c415e2ee84bcebaf8c93517dabac0d7ce64665eff1538da0f0d3953d8006e062371889743977352014b735c8d9d8589739f7c66178eac11f3685623566a76a2a6fd0bb5d75854d7a0736a1a8d6474952e89
q4 = 0xcf27ccccc5b57655dc22fd0c7ba6962c1272c557edfa81a71ba95b4dcb65b6e0a9398c876eb8318a17423eac62fbb3f9986da6dab3cfec8954453a6e2f179cd0789d3bd64f21645a1a7d83c8af32fdfb2cc5d6e8f0c45df37cfe390a628bdfb79ce35d769285f1ca28a95e805d21e1fcac1eac3ca6f3eb0035a5a712e6793029
m4 = 0x76114187a6afac7b315847ee4736d545
b'v\x11A\x87\xa6\xaf\xac{1XG\xeeG6\xd5E'
'''

数据6——大整数分解·

数据6其实并不存在以上说的所有漏洞，而且感觉上若随机序列真的随机，该数据还是安全的。

不过，赛题其实有一大章介绍了大整数分解算法，所以最后摆烂了直接上分解，然后用Pollard’s p-1算法分解出来了，实际上对于2048比特可选择的分解方法并不多，排除了一下基本就Pollard’s p-1或Williams’s p+1。（从之前山石冬令营的smooth_rsa也获得了点灵感）

下面大概说一下Pollard’s p-1算法（内容摘抄自[HPS14]），假设存在一个整数$L$，使得

$$P - 1 \mid L \\ Q - 1 \not\mid L$$

则根据费马小定理，大概率（取决于所选的$a$，最好$a$是生成元？）有

$$a^L \equiv 1 \pmod P \\ a^L \not\equiv 1 \pmod Q$$

所以即大概率有

$$P \mid a^L - 1 \\ Q \not\mid a^L - 1$$

于是通过计算$P = \text{GCD}(a^L-1, N)$便可分解$N$。

然后问题就变成，如何找到这样的$L$，一种方法是可以通过枚举$n = 2, 3, ...$，然后设$L = n!$，枚举约到$P-1$（或$Q-1$）的最大素因子时即可成功分解。

以上伪代码是教科书算法，实际代码时可以稍作优化，比如上面提到的$k$次迭代才作$\text{GCD}$，又比如另开线程做步骤4和5。

参考代码（gmpy版本）：

import _thread
from math import inf
from gmpy2 import gcd, powmod
import time

DONE = None
def check(n, a, j, lfn):
  global DONE
  d = gcd(a-1, n)
  if d == n:
    print('Bias so large!')
    DONE = False
    return
  if d > 1 and d < n:
    print('Done:')
    DONE = d   # assert BIAS is large enough, so omit lock
    log = '[Done] d = %d' % d
  else:
    log = '[Debug - %s] j = %d\ta = %d\t d = %d' % (time.asctime(time.localtime(time.time())), j, a, d)
  print(log)
  if lfn != None:
    with open(lfn, 'a') as lf:
      lf.write(log + '\n')
  return

def pollard(n, a=2, j=2, B=inf, lfn=None, BIAS=100000):
  global DONE
  counter = j % BIAS - 1  # assert j > 1
  while j < B:
    counter += 1
    if counter == BIAS:
      if DONE != None:
        break
      counter = 0
      _thread.start_new_thread(check, (n, a, j, lfn, ) )
    a = powmod(a, j, n)
    j += 1
  return DONE

n = 0xdd329082ffbd1dbac0ef27b781686a6e8eaccc2a59d37883c5b380acc68b79c9c36fdb381f453d18dd543bbc5a015b6ac82c311260e80a94d784fe54a45c4dbd0534ac7fa5ebca591cb50c0f12f51e478fd2e6d2d948d78185486b2ee10d78c903ec185d8fd7750d04132a99431d9516a73e5d1a3508f428819b487d812392b317a9cdafea3bc1ec7aa91b870f533cc356f3185e05c6a8f052bd4e48800b4b662494fa40c7f810afc3c7293fcf8ffb552c24ea957ce34b8a14ba96e76c88eeed0bb7b2a14f501c7303fd8734fc179f8f56f6c7e8a282321fdcd51d39b11459d2547081f9ed27d918696e7c2ce9f0041976df9b648299d757ecae99e0b2fb4333
print(pollard(n, lfn='3.log', BIAS=100000))

'''
[Debug - Tue Apr 18 14:37:07 2023] j = 100000	a = 7496900376453995838444422927099975267896510388395470026643515771447895895545762691521162629337730874532379835365388788709040203864203032948244411197319102089332745318709697871483609402894941329667560189480016401908029264593724906134741309837153701754212123601252682432056682038235171581815979551104976857548946825991866581044606825975554440351612070207005169020789701718816531407147657081043472602096002347598764981919599269743440681248167507921125024182521464204696113262284083964135241979918178597287892302860485232065608621611745676294835369995815607124326162696696133126596035260288607530664412090415414657959489	 d = 1
  ... ...
[Debug - Wed Apr 19 02:52:07 2023] j = 1010200000	a = 18621235032600903646632361168505542876588771234092696924524378970540517802247635378529236154187296822009448979894979135909760675673864821398535537329550049434085040432515322206857563806712398629745664773908274163494016954960211867505928515570766269453374014727302069278575833073873835631019266109113578649672214620399771733898811234281514094243555913822109512643485354076996945684586465330799483935022562991620806883703981349570389299595310776968506987598333740591397342141210782987058069599256998729983189894163097432700255490225731404563154832441446069765106665540722542290105513805013210660820652627595580441713666	 d = 1
[Debug - Wed Apr 19 02:52:12 2023] j = 1010300000	a = 22960846295748386771473173898353066638167836943473058803358415403729789219893140746703152729552281832109669257672115479672325420225491260604301324196840724480072367125072579967848103767090835320615043939816759850841271408904991024445695894494032241768772903537776261441888600724408146751340464275276127288046834820175958902987686842702594172834176187657499992508845690619086011597393437715939846411638871130497631581285515411716221638059094294020483949261053515850926891651744376536957946380943441291105784915151955305992014749653298120799150887961734360118559172600485088668984574934814741516912838862486208381425911	 d = 1
[Debug - Wed Apr 19 02:52:16 2023] j = 1010400000	a = 5556981995531096297233758890805669574865961211699630646887648193803243412326419934922504688282573434915261545134000082083979362558004758148547668031752021628246906349615648649928265656056355078645895172941239055403440222079027995406027117310310109245473270971269488662940378761736533375879010240559780270053753167583276700726677204815929567782597185523322071111916165080791324336931834767054863076964593789743772327563858596480545024833748198741099697137679017962957794007704163052256117310123990687313471374424865038097367076195143376080953453353259355390036319361100877585014146502680485302883397328065561324041671	 d = 1
[Done] d = 173392280438546792683818136702600967448366743521806189327229178974568108822276784275748977768555394587467983783851310057943432666398430604814797445730219422150371757816773423848241811463624899203965492680109656320607902960451073955585153506896911302285757276495392556781174371601361441446529516356977237812613
'''

另外提供C语言版本的代码，实际测试和gmpy2的速度差不多，由于gmpy2基于Cython，所以底层也是C。

参考代码（C版本）：

// gcc pollard2.c -lgmp -lm -lpthread -o pollard2.exe
// gcc pollard2.c -lgmp -lm -lpthread -O3 -o pollard2 && ./pollard2
#include <time.h>
#include <stdio.h>
#include <string.h>
#include <pthread.h>
#include <gmp.h>
#define BIAS 500000         // smaller than B
#define LFN  "c2.log"
#define B    1099511627776  // 2^40
mpz_t z1, z2, zbias;
mpz_t n, a, res;
unsigned long int j;
char buffer[1000];
int DONE = 0;

void* check(void* arg){
    mpz_t d;
    mpz_init(d);
    // gcd
    mpz_sub(d, a, z1);
    mpz_gcd(d, d, n);

    // check
    if (mpz_cmp(d, n) == 0) {
        printf("Bias so large!\n");
        mpz_clear(d);
        pthread_exit(0);
    }
    if (mpz_cmp(d, z1) > 0 && mpz_cmp(d, n) < 0) {
        DONE = 1;
        printf("Done:\n");
        gmp_printf("%Zd\n", d);
        mpz_set(res, d);
        mpz_clear(d);

        FILE *fp = NULL;
        fp = fopen(LFN, "a+");
        gmp_sprintf(buffer, "[DONE] - d = %Zd\n", d);
        fputs(buffer, fp);
        fclose(fp); 
        pthread_exit(0);
    }
    //gmp_printf("[+++]%lu\t%Zd\n", j, d);
    // log
    time_t t;
    struct tm *timeinfo;
    time(&t);
    timeinfo = localtime(&t);
    sprintf(buffer, "[DEBUG] - %s ", asctime(timeinfo));
    gmp_sprintf(buffer + strlen(buffer), "j = %lu\ta = %Zd\n", j, a);
    printf("%s", buffer);

    // log file
    FILE *fp = NULL;
    fp = fopen(LFN, "a+");
    fputs(buffer, fp);
    fclose(fp); 

    // clear
    mpz_clear(d);
    pthread_exit(0);
}

void pollard() {
    unsigned long int counter = BIAS - j + 1;
    pthread_t tid2;
    while (j < B) {
        counter -= 1;
        if (counter == 0) {
            if (DONE != 0) {
                break;
            }
            counter = BIAS;
            //new thread
            pthread_create(&tid2, NULL, check, NULL);
            //gmp_printf("[---]%lu\t%Zd\n", j, res);
        }
        mpz_powm_ui(a, a, j, n);
        j += 1;
    }
    return;
}

int main() {
    mpz_init_set_str(z1, "1", 10);
    mpz_init_set_str(z2, "2", 10);
    mpz_init_set_ui(zbias, BIAS);

    mpz_init_set_str(res, "0", 10);
    mpz_init_set_str(n, "27923599681213021407034926611637566469874523487438156297829021871854267185336835555834417699638124062081862065978435358275789755449922797880942372540343149490351082513378216121657620510652337065296119029212679608293116863721614934226652937927332211856317133563520857800844724599974454798616386551458666742419594273729122176671137165710858043854085715611469987477104346783638814358091242759009090423077469924079120414624179345633165061041336969051478949273930390780624527965481748666910510642670099841373434666890288335574645730875375871875051442416833323104303435247592133718723309813270278494345959197697574411387699", 10);
    //gmp_printf("%Zd\n", n);
    mpz_init_set_str(a, "2", 10);
    //mpz_init_set_str(j, "2", 10);
    j = 2;
    pollard();
    gmp_printf("%Zd\n", res);

    mpz_clears(z1, z2, zbias, n, a, res, NULL);
    return 0;
}

最终RSA解密：

# sage

 
n6 = 0xdd329082ffbd1dbac0ef27b781686a6e8eaccc2a59d37883c5b380acc68b79c9c36fdb381f453d18dd543bbc5a015b6ac82c311260e80a94d784fe54a45c4dbd0534ac7fa5ebca591cb50c0f12f51e478fd2e6d2d948d78185486b2ee10d78c903ec185d8fd7750d04132a99431d9516a73e5d1a3508f428819b487d812392b317a9cdafea3bc1ec7aa91b870f533cc356f3185e05c6a8f052bd4e48800b4b662494fa40c7f810afc3c7293fcf8ffb552c24ea957ce34b8a14ba96e76c88eeed0bb7b2a14f501c7303fd8734fc179f8f56f6c7e8a282321fdcd51d39b11459d2547081f9ed27d918696e7c2ce9f0041976df9b648299d757ecae99e0b2fb4333
p6 = 173392280438546792683818136702600967448366743521806189327229178974568108822276784275748977768555394587467983783851310057943432666398430604814797445730219422150371757816773423848241811463624899203965492680109656320607902960451073955585153506896911302285757276495392556781174371601361441446529516356977237812613
q6 = n6 // p6
assert p6 * q6 == n6
print('p6 = %s' % hex(p6))
print('q6 = %s' % hex(q6))


import libnum
e6 = 89
c6 = 0x39807efece486e643e434cfd26f38398465aaf3bb2c361782de9ef2c9d34d83721bf3eb2be08273b69f3d1c57e3acb79e002ebe4477b834e5b46aed77dffcab6704aa8b5aa33e7d6657a2d1a399abb46252757734077527987183ced85ad95510e04fb75c137b407426e5eca9e40ab22a3ce61eb52ee9172b68604f108b60d41ad622505923aea5e06e4c8a4904a45200ef5b89c3f243a475198ac666709c32fc6b0f016c896a77905e8d4c1c97fab1b4c4915741e77ab0f8181bc1ff8f5fb1779318d00f94934dcbe7e841875954be45414452b11ae8eabfea75dc3bd3bf1a15347d6bb7fc475ca171138cb3daf4d9b9a14159e65213a875bbbc7e709469728
phi6 = (p6-1) * (q6-1)
d6 = e6.inverse_mod(phi6)
m6pad = libnum.n2s(int(pow(c6, d6, n6)))
m6 = m6pad[-(128//8):]
print('m6 = %s' % hex(libnum.s2n(m6)))
print(m6)

'''
m6 = 0xa9ad7e791d2e9d7ee3c11330851f5897
b'\xa9\xad~y\x1d.\x9d~\xe3\xc1\x130\x85\x1fX\x97'
'''

效率·

贴一下上面算法的复杂度即在我渣机上的运行时间：

$$\begin{array}{|l|l|r|l|} \hline \text{数据} & \text{理论复杂度} & \text{实际运行时间} & \text{备注} \\ \hline \text{1} & \text{O}(\text{log}(N_1)) & 1.98s & \text{$N_1$与$N_3$同规模} \\ \text{2} & \text{O}(k_{21}\text{log}(e_2)) & 3.67s & \text{$k_{21}$为$K_2$中$8$比特以上的最小因子} \\ \text{3} & \text{O}(\text{log}(N_3)) & 1.98s & \text{$N_1$与$N_3$同规模} \\ \text{4} & \text{O}(\text{log}^2(N_4)) & 4.40s & \\ \text{5} & \text{O}(\text{log}(e_5)) & 1.95s & \\ \text{6} & \text{O}(B_6\text{log}(B_6)) & 12h15m09s & \text{$B_6$为$P_6-1$（或$Q_6-1$）最大素因子的上界} \\ \hline \end{array}$$

总结·

结果AK了也只拿了第四，差两分就季了（