Common Prime RSA 笔记

发表于 更新于 本文字数： 7k 阅读时长 ≈ 6 分钟
Common Prime RSA 相关的笔记

Common Prime RSA

简介

Wiener 提出，如果 \(p\) 和 \(q\) 是素数，使得 \(p-1\) 和 \(q-1\) 有一个大数因子，那么具有这种性质的素数可以作为 Wiener Attack 的反制措施。

根据这种理论，出现了一种 \(p-1\) 和 \(q-1\) 拥有一个共同的大素数因子的 RSA 变体，我们称作共素数RSA（Common Prime RSA）。

对于某个大素数 \(g\)，让 \(p=2ga+1\) 且 \(q=2gb+1\) 作为平衡素数，保证 \(a,b\) 互素且 \(h=2gab+a+b\) 是素数。第一个限制确保 \(\gcd(p-1,q-1)=2g\)，第二个限制确保 \((pq-1)/2=gh\) 与 \(N=pq\) 的大小接近。

因为由上述算法生成的素数 \(p,q\) 满足 \(g=\gcd(p-1,q-1)\) 是一个大素数因子，故称 \(p,q\) 为共素数（common primes）。其中 \(g\) 为这两个素数的共因子（common factor）。

我们需要注意到对于共素数RSA有着以下性质： \[ \begin{align} \lambda(pq)=\mathrm{lcm}(p-1,q-1)=\mathrm{lcm}(2ga,2gb)=2gab\\ \varphi(pq)=(p-1)(q-1)=2ga2gb=2g\lambda(pq) \end{align} \] 此外存在额外定义，RSA加密指数和解密指数需要与 \(\lambda(pq)\) 互素。

根据上述定义，可以推导出 \[ N=pq=(2ga+1)(2gb+1)=2g(2gab+a+b)+1=2gh+1 \] 即 \(N-1\) 为 \[ N-1=2g(2gab+a+b)=2gh \] 定义 \(\gamma\) 表示共因子 \(g\) 的相对于 \(N\) 的大小，即 \(g=N^\gamma\)。考虑 \(g\leq N^{1/2}\)，故 \(0\leq\gamma\leq1/2\)。

生成算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from Crypto.Util.number import *
try:
    gcd
except NameError:
    from math import gcd

def gen_prime(nbits: int, gamma: float):
    g = getPrime(int(nbits * gamma))
    alpha = 0.5 - gamma
    while True:
        a = getRandomNBitInteger(int(alpha * nbits))
        p = 2 * g * a + 1
        if isPrime(p):
            b = getRandomNBitInteger(int(alpha * nbits))
            q = 2 * g * b + 1
            h = 2 * g * a * b + a + b
            while not isPrime(q) or isPrime(h) or gcd(a, b) != 1:
                b = getRandomNBitInteger(int(alpha * nbits))
                h = 2 * g * a * b + a + b
                q = 2 * g * b + 1
            return p, q
print(gen_prime(512, 0.48))

攻击

\(\gamma\) 接近于 \(1/2\)

当 \(\gamma\) 接近于 \(1/2\) 时，\(g=N^\gamma\) 接近于 \(N^{1/2}\)，由于共素数RSA的特殊构造我们可以在 \(O(N^{1/4-\gamma/2})\) 分解 \(N\)，故此时算法时间复杂度接近于 \(O(1)\)。

此时我们只需修改 Pollard's rho method 的 \(x_i\) 函数。

在 Mckee&Pinch 的论文《Further Attacks on Server-Aided RSA Cryptosystems》中指出将 \(f(x)=x^2+1\) 修改为 \(f(x)=x^{N-1}+3\pmod{N}\) 是最优解。

由于 \(N-1=2gh\) 且 \(p-1=2ga\)，故最多只有一个值不在 \(x^{N-1}\mod{p}\) 的环中。

算法实现如下：

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
import random
try:
    gcd
except NameError:
    from math import gcd

def rho(N):
    f = lambda x: (pow(x, N-1, N) + 3) % N
    while True:
        # 加快入环速度
        t = random.randint(2, N)
        h = f(t)
        step_times = 0
        step_limit = 2
        while True:
            if not step_times < step_limit:
                step_times = 0
                step_limit *= 2
                t = h
                h = f(h)
            p = gcd(abs(int(t) - int(h)), N)
            if p == N:
                break
            elif p > 1:
                return (p, N // p)
            else:
                h = f(h)
                step_times += 1

print(rho(84236796025318186855187782611491334781897277899439717384242559751095347166978304126358295609924321812851255222430530001043539925782811895605398187299748256080526691975084042025794113521587064616352833904856626744098904922117855866813505228134381046907659080078950018430266048447119221001098505107823645953039))

故事实证明 \(\gamma\) 值的选取不能过于接近 \(1/2\)。

已知 \(a,b\)

已知 \(N=2g(2gab+a+b)+1\)，于是我们构造方程 \[ 4abg^2+2(a+b)g-N+1=0 \] 可以在多项式时间 \(\log(N)\) 内分解 \(N\)。

Sagemath 代码如下：

1
2
3
4
5
6
7
8
9
10
11
a = 1185431345934512
b = 1989628969125971
N = 54692260436051338814890781701826055707958209029414126894070449935683071253184867947357262267840171428710181955973010913204514025135188192484651672240708141692701667242130748316666406528479191422804307020656050201187035928715833163999813216597718706449260040885862566373392398826670863398295350419792842640631
P.<g> = ZZ[]
f = 4 * a * b * g ^ 2 + 2 * (a + b) * g - N + 1
g = f.roots()
if g:
    g = g[0][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

Python 可以借助 sympy 解方程。

已知 \(g\)

当 \(g\geq a+b\) 时，会导致在多项式时间为 \(\log(N)\) 内分解 \(N\)，同时因为素数是平衡的，条件相当于 \(g>N^{1/4}\)。

证明如下：

我们假设 \(g>a+b\)，那么给予 \(N,g\)，令 \(M=(N-1)/(2g)\)，\(c=a+b\)，那么方程 \[ N=2g(2gab+a+b)+1 \] 可以改写成 \[ M=2gab+c \] 因为 \(c=a+b<g\)，根据假设，可以归约为模 \(g\) 域下的方程： \[ c=M\pmod{g} \] 因此，\(c=a+b\) 是已知的。

将 \(b=c-a\) 带回方程 \(N=2g(2gab+a+b)+1\)，整理可得二次方程： \[ 2ga^2-2gca+(N-1)/(2g)-c=0 \] 可以解得 \(a,b\) 为该方程两解。

假设 \(g=a+b\)，我们将方程 \(N=2g(2gab+a+b)+1\) 作等价替换，整理得到方程 \[ \frac{N-1}{4g^2}=ab+\frac{1}{2} \] 我们再进一步替换 \(b=g-a\)，再次整理可得二次方程： \[ a^2-ga+\frac{N-1}{4g^2}-\frac{1}{2}=0 \] 可以解得 \(a,b\)。

代码如下：

\(g>a+b\) 时，

1
2
3
4
5
6
7
8
9
10
11
12
g = 2056971706333850947354991471886113601423457483931388832864204860321308350537317091564919029078296379733989138742162694786565228112885684303
N = 67324909911911622626246005558967775211455024820506932698435813321567574468019013664789401988015894964099052816176029553245881317276340043887466584645914352982274378611180595397686920214079479901514703963131435008906250160656759300390805929849374653321934393399433471228218819498373221757779799476717494079667
M = (N - 1) // (2 * g)
c = M % g
P.<a> = ZZ[]
f = 2 * g * a ^ 2 - 2 * g * c * a + M - c
a = f.roots()
if a:
    a, b = a[0][0], a[1][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

\(g=a+b\) 时，

1
2
3
4
5
6
7
8
9
10
11
g = 2855372645569408464444580237486670388029956719716115953907612135874419892154982850222965560661211729647325085879529571229774148545656169021
N = 159549169988238873893531105042878385551537587717347282632324748268846735710748763722602882823022008548774298858161130258369850715542192739582830583643642436399008902770027668038725347353393047833875066622910131525247842517372845617227325882916166114361718015983671803859502931814932543107911548450229250776542101141849788751722460468073974316977656001286989710480324512919121409123619799426221232443036698458643438020098037548757403
M = (N - 1) // (2 * g)
P.<a> = ZZ[]
f = 2 * a ^ 2 - 2 * g * a + (N - 1) // (2 * g ^ 2) - 1
a = f.roots()
if a:
    a, b = a[0][0], a[1][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

因为 \(p=2ga+1\)，注意到 \(g>N^{1/4}\) 是已知的，故我们也可以用 Coppersmith's factoring 方法分解模数。

而当 \(g<a+b\) 时，我们可以使用时间复杂度为 \(O(N^{1/4-\gamma})\) 的算法分解 \(N\)，每个操作的多项式时间为 \(\log(N)\)。

证明如下：

已知 \(g\)，我们可以通过除法算法计算 \(u\) 和 \(0\leq v\leq2g\)，例如： \[ \frac{N-1}{2g}=2gu+v \] 因为我们知道 \(N=2g(2gab+a+b)+1\)，于是乎 \[ a+b=v+2gc\\ ab=u-c \] 其中 \(c\) 为任意整数。

对于任意与 \(N\) 互素的整数 \(x\)，我们有 \[ x^{u2g}\equiv x^{ab2g+c2g}\equiv x^{c2g}\pmod{N} \] 因为 \(u=ab+c\) 且 \(\lambda(N)=2gab\)，所以 \(x^{2gab}\equiv 1\pmod{N}\)。

让 \(y=x^{2g}\)，我们有 \[ y^u\equiv y^c\pmod{N} \] 根据这个关系，未知的 \(c\) 可以用 Shanks 的小步大步法（baby-step giant-step methodology）求解。对于某些 \(d>\sqrt{c}\)，我们计算得到大步为 \[ y^0,y^d,y^{2d},\cdots,y^{d^2}\mod{N} \] 小步为 \[ y^u,y^{u-1},y^{u-2},\cdots,y^{u-d}\mod{N} \] 在其中搜索碰撞，将会产生一个 \(r\) 和 \(s\) 满足： \[ y^{rd}\equiv y^{u-s}\pmod{N} \] 其中 \(c=rd+s\)。

当 \(c\) 已知，我们可以计算 \(a,b\)。

计算，排序和搜索需要 \(O(d\log(d))\) 操作，其中 \(d>\sqrt{c}\)。

故使用这种算法需求 \(\gamma\) 接近于 \(1/4\)。

本人多次跑这个算法，最终验证可得 \(c\) 的大致范围大概率会落在 \(N^{(1/2-2\gamma)}\) 的附近，取上下浮动两位大小最佳。

Sagemath 如下：

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
from sage.groups.generic import bsgs

g = 28838314918840273611038952856252141148393903858815521699377328440901497
N = 23895816238623712091906549349650291003358466897973970818205381234024572750472476312894286233088847723906345049342176974080465940396626170377635736786780099297743807402105462746828310247313063710700453371777283064626446124365885063073122303271609231832129823571694756234973129694217982548355078107610764879363
nbits = 1024
gamma = 0.23
cbits = ceil(nbits * (0.5 - 2 * gamma))

M = (N - 1) // (2 * g)
u = M // (2 * g)
v = M - 2 * g * u
GF = Zmod(N)
x = GF.random_element()
y = x ^ (2 * g)
# c的范围大概与N^(0.5-2*gamma)很接近
c = bsgs(y, y ^ u, (int(2**(cbits-1)), int(2**(cbits+1))))
ab = u - c
apb = v + 2 * g * c
P.<x> = ZZ[]
f = x ^ 2 - apb * x + ab
a = f.roots()
if a:
    a, b = a[0][0], a[1][0]
    p = 2 * g * a + 1
    q = 2 * g * b + 1
    assert p * q == N

分解 \(N-1\)

当 \(\gamma\) 过小，即 \(g=N^\gamma\) 过小时，因为已知 \(N-1=2gh\)，故分解出 \(g\) 是较为容易的。

可以使用 yafu 分解 \(\frac{N-1}{2}\)，其中当 \(\gamma\) 值约为 0.10 左右时分解迅速。

Mumtaz-Luo攻击

由已知，得 \[ \begin{cases} ed=(p-1)bk+1\\ ed=(q-1)ak+1 \end{cases} \] 可以变换为 \[ \begin{cases} ed-1+bk=pbk\\ ed-1+ak=qak \end{cases} \] 两式相乘，可得 \[ (ed-1+bk)(ed-1+ak)=pbk\cdot qak=abk^2N \] 展开得到完整式 \[ e^2d^2+ed(ak+bk-2)-abk^2(N-1)-ak-bk+1=0 \] 不妨构建多项式函数 \[ f(x,y,z)=e^2x^2+ex(y+z-2)-(y+z-1)-(N-1)yz \] 其中特解 \((x^*,y^*,z^*)=(d,ak,bk)\)，边界为 \[ \begin{align} &|x^*|\leq N^\delta\\ &|y^*|\leq N^{\delta-\gamma+0.5}\\ &|z^*|\leq N^{\delta-\gamma+0.5} \end{align} \] 其中，\(\delta\) 是私钥 \(d\) 的大小。