椭圆曲线
咕咕ing
常用算法
ps.下述的算法在不同领域有很多不同的变体但核心思想是一致的
- 大整数的分解
- 离散对数的求解
- 椭圆曲线上的离散对数求解
BSGS
1.离散对数
对于这样的一个式子$g^x=h\ mod\ m$对于一般暴力算法的复杂度是$O(\varphi(m))$，对于较大m会计算困难。于是有BSGS算法将复杂度减小到$O(\sqrt{m})$。个人认为其算法思想和MITM攻击类似。
主要操作如下，先将等式变化
$$g^{A\left \lceil \sqrt{m}\right \rceil-B}=h\ mod\ m$$
$$g^{A\left \lceil \sqrt{m}\right \rceil}=hg^B\ mod\ m$$
对等式右边遍历$B(B<\sqrt{m})$并打表，然后遍历等式左侧的A遍历并且找到相同的值，返回解$A\left \lceil \sqrt{m}\right \rceil-B$
def bsgs(g,h,p):
tmp = ceil(sqrt(p))
bs = {}
for B in range(tmp):
bs[h*pow(g,B,p)%p] = B
tmp1 = pow(g,tmp,p)
for A in range(tmp):
if pow(tmp1,A,p) in bs:
return A*tmp - bs[pow(g,tmp*A,p)]修改了一下如果x的大小是已知的可以添加参数x最大值pi来减小复杂度
def bsgs(g,h,p,pi=None):
if pi is None:
pi = p
tmp = ceil(sqrt(pi))
bs = {}
for B in range(tmp):
bs[h*pow(g,B,p)%p] = B
tmp1 = pow(g,tmp,p)
for A in range(tmp):
if pow(tmp1,A,p) in bs:
return A*tmp - bs[pow(g,tmp*A,p)]2.椭圆曲线
对于$Q=xP$，相似的可以得到
$$P=(a\left \lceil \sqrt{m}\right \rceil+b)G$$
$$P-bG=a\left \lceil \sqrt{m}\right \rceil G$$
同理
def bsgs(G,P,p,pi=None):
if pi == None:
pi = p
tmp = ceil(sqrt(pi))
bs = {}
for b in range(tmp):
bs[P-b*G] = b
tmp1 = G*tmp
for a in range(tmp):
if tmp1*a in bs:
return a*tmp+bs[tmp1*a]Pollard rho
这个算法其实有点运用了🐒算法即随机游走
1.大整数分解
🐒算法体现在构造了一个环上的伪随机生成函数例如$f(x)=x^2+c\ mod\ n$，将这个函数用作递推数列产生不断的随机数。由于是在有限域上所以必定最终随机数会进入循环。易得当$x\equiv y$时，$f(x)\equiv f(y)$。
可以使用 Floted 算法判环。在一个ρ形图中，可以使用两个指针x、y，每次让x前进一步，y 前进两步，如果x和y重合了，就说明有环。同时可以对x，y计算gcd判断来猜测是否存在因数
def pollard_rho(n):
c = random.randint(1, n-1)
f = lambda x: (x*x + c) % n
x, y, d = 0, 0, 1
while d == 1 or d == n:
x = f(x)
y = f(f(y))
if x == y:
return pollard_rho(n)
d = GCD(abs(x - y), n)
return d拓展一下有这样的分解因数算法
def factor(n):
factors = []
while n != 1:
p = pollard_rho(n)
times = 0
while n % p == 0:
n //= p
times += 1
factors.append((p,times))
return factors2.离散对数
对于$g^x\equiv h \mod p$，类似的使用相似的思想，先令$x\equiv g^ah^b\mod n$，我们会试图构造一个关于$x$生成函数，同时使用分区映射（为了快速游走，且满足良好的随机性），将整个群$G$分成数量大致相同的三部分$S$，最终得到这样的生成函数
$$ f(x)=f(\varphi(a,b))=\left\{\begin{matrix} \varphi(2a,2b) & x\in S_1 \\ \varphi(a+1,b) & x\in S_2 \\ \varphi(a,b+1) & x\in S_3 \\ \end{matrix}\right. $$然后同样的使用 Floted 算法判环，寻找到$\varphi(a_1,b_1)=\varphi(a_2,b_2)$，基于对$x$的构造，可以推导出
$$ \begin{matrix} & \varphi(a_1,b_1) & \equiv & \varphi(a_2,b_2) & \mod p \\ \Rightarrow & g^{a_1}h^{b_1} & \equiv & g^{a_2}h^{b_2} & \mod p \\ \Rightarrow & g^{a_1-a_2} & \equiv & g^{x(b_2-b_1)} & \mod p \\ \Rightarrow & x & \equiv & \frac{a_1-a_2}{b_2-b_1} & \mod ord(g) \end{matrix} $$下面是示范用的代码实现，不保证一定有解。
def pollard_rho(g, h, p, o=None):
if o is None:
o = p - 1
f = lambda x, a, b: (x**2%p,2*a,2*b) if x%3 == 0 else ((x*g%p,a+1,b) if x%3 == 1 else (x*h%p,a,b+1))
phi = lambda a,b: pow(g, a, p) * pow(h, b, p) % p
a1,b1 = random.randint(0, o),random.randint(0, o)
x1 = phi(a1,b1)
x2,a2,b2 = f(x1,a1,b1)
while x1 != x2:
x1,a1,b1 = f(x1,a1,b1)
x2,a2,b2 = f(*f(x2,a2,b2))
gcd = lambda a,b: gcd(b,a%b) if b else a
tmp = gcd(o, abs(b2-b1))
if (a1-a2) % tmp != 0:
return pollard_rho(g, h, p, o)
res = ((a1-a2)//tmp)*pow((b2-b1)//tmp,-1,o//tmp) % (o//tmp)
return res⚠ 该代码在指定底数的阶的时候一定有解，但是当阶未知的时候，会使用群阶作为元素的阶，此时可能出现无解，所以使用随机初始状态来尝试多个不同的环，来保证能输出解。此外阶未知时，输出的解只能保证$res\equiv x \mod ord(g)$
3.椭圆曲线
同样的使用类似的思想，对于$P=x\cdot G$，先令$X=a\cdot G+b\cdot P$，然后类似的构造相同的函数$f(x)$
$$ f(X)=f(\varphi(a,b))=\left\{\begin{matrix} \varphi(2a,2b) & x\in S_1 \\ \varphi(a+1,b) & x\in S_2 \\ \varphi(a,b+1) & x\in S_3 \\ \end{matrix}\right. $$并且使用 Floted 算法判环，寻找到$X_1=X_2$，基于对$X$的构造，可以推导出
$$ \begin{matrix} & \varphi(a_1,b_1) & \equiv & \varphi(a_2,b_2) & \mod p \\ \Rightarrow & {a_1}\cdot G+{b_1}\cdot P & \equiv & {a_2}\cdot G+{b_2}\cdot P & \mod p \\ \Rightarrow & (a_1-a_2)\cdot G & \equiv & (b_2-b_1)x\cdot G & \mod p \\ \Rightarrow & x & \equiv & \frac{a_1-a_2}{b_2-b_1} & \mod ord(G) \end{matrix} $$下面是示范用的代码实现
def pollard_rho(E, G, P, o=None):
if o == None:
o = E.order()
f = lambda X, a, b: (2*X, 2*a, 2*b) if int(X[0])%3 == 0 else ((X+G, a+1, b) if int(X[0])%3 == 1 else (X+P, a, b+1))
phi = lambda a,b: a*G + b*P
a1,b1 = random.randint(0, o), random.randint(0, o)
X1 = phi(a1, b1)
X2,a2,b2 = f(X1, a1, b1)
while X1 != X2:
X1,a1,b1 = f(X1, a1, b1)
X2,a2,b2 = f(*f(X2, a2, b2))
R.<x> = PolynomialRing(Zmod(o))
f = (b2-b1)*x - (a1-a2)
res = f.roots(multiplicities=False)
if len(res) > 1e6: # 限制增根数量
return pollard_rho(E, G, P, o)
for i in res:
if i*G == P:
return i
else:
return pollard_rho(E, G, P, o)
这里把反求x的过程使用了Sage的求根运算，因为我们只能保证x是该多项式的一个根，而不是唯一的根，同时防止大素数p导致增根数量过多，这里限制为1e6。
Pohlig-Hellman
1.离散对数
条件：$p$为质数，且$\varphi(p)$是光滑数
对于$a^x=b\ mod\ p$，确定原根g，可以转化为$g^{ux}=g^v\ mod\ p$
也就是$x=v\cdot u^{-1}\ mod\ p-1$，因此只需要计算出u，v即可解的x
下面推导如何计算$g^u=a\ mod\ p$中的u
对于$p-1$可以分解因数得到$p-1=p_1^{q_1}p_2^{q_2}\cdots p_n^{q_n}$
可以构造$u$使得满足$u=c_0p_1^0+c_1p_1^1+c_2p_1^2+\cdots+c_{q_1-1}p_1^{q_1-1}\ mod\ p_1^{q_1}$
于是就有
$$ \left\{\begin{matrix} u=c_0+c_1p_1+c_2p_1^2+&\cdots&+c_{q_1-1}p_1^{q_1-1}\ mod\ p_1^{q_1}\\ u=c_0’+c_1’p_2+c_2’p_2^2+&\cdots&+c_{q_2-1}’p_2^{q_2-1}\ mod\ p_2^{q_2}\\ &\cdots&\\ u=c_0^{”}+c_1^{”}p_3+c_2^{”}p_3^2+&\cdots&+c_{q_3-1}^{”}p_3^{q_3-1}\ mod\ p_3^{q_3}\\ \end{matrix}\right. $$求解每个同余方程的c再用crt合并即可得到
下面解c
构造此式子$(g^u)^{\frac{p-1}{p_1^r}}=a^{\frac{p-1}{p_1^r}}\ mod\ p$
于是令$r=1$时有$g^{c_0\frac{p-1}{p_1}}=a^{\frac{p-1}{p_1}}\ mod\ p$
此时只需要爆破即可解出$c_0$，且p是光滑数，显然爆破是可行的
接着调整$r=2$，同理爆破就可以得到$c_1$
所有c解出来之后用crt合并即可得到u
然后回代上式，得到x
from Crypto.Util.number import *
from sage.all import *
def Pohlig_Hellman(a, b, p):
g = int(GF(p).primitive_element())
factors = list(factor(p-1))
mods = [a**b for a,b in factors]
def cal(a):
residues = []
for i in range(len(factors)):
tmp = 0
tmp1 = pow(g,(p-1)//factors[i][0],p)
cs = []
for r in range(1,factors[i][1]+1):
tmp2 = pow(g,tmp,p)
tmp3 = pow(a,(p-1)//factors[i][0]**r,p)
for c in range(factors[i][0]):
if tmp2*pow(tmp1,c,p) == tmp3:
cs.append(c)
tmp += c*(p-1)//factors[i][0]
tmp //= factors[i][0]
break
residues.append(sum(c*factors[i][0]**k for k,c in enumerate(cs))%mods[i])
u = crt(residues,mods)
return u
u = cal(a)
v = cal(b)
if GCD(u,p-1) == 1:
return v*pow(u,-1,p-1)
else:
tmp = 1
while True:
gcd = GCD(u,p-1)
if gcd == 1:
break
tmp *= gcd
u //= gcd
return ZZ(v*pow(u,-1,p-1))//tmp
if __name__ == '__main__':
p = 7863166752583943287208453249445887802885958578827520225154826621191353388988908983484279021978114049838254701703424499688950361788140197906625796305008451719
a = getRandomNBitInteger(256)
x = getRandomNBitInteger(256)
b = pow(a,x,p)
print(x)
print(Pohlig_Hellman(a, b, p))ps.在暴破c的时候还有一说可以用BSGS来优化算法，如下
# Sample
# BSGS & Pohlig_Hellman for DLP
from Crypto.Util.number import *
from sage.all import *
def bsgs(g,h,p,pi):
if pi == None:
pi = p
tmp = ceil(sqrt(pi))
bs = {}
for B in range(tmp):
bs[h*pow(g,B,p)%p] = B
tmp1 = pow(g,tmp,p)
for A in range(tmp):
if pow(tmp1,A,p) in bs:
return A*tmp - bs[pow(g,tmp*A,p)]
def Pohlig_Hellman(a, b, p):
g = int(GF(p).primitive_element())
factors = list(factor(p-1))
mods = [a**b for a,b in factors]
def cal(a):
residues = []
for i in range(len(factors)):
tmp = 0
tmp1 = pow(g,(p-1)//factors[i][0],p)
cs = []
for r in range(1,factors[i][1]+1):
tmp2 = pow(a,(p-1)//factors[i][0]**r,p)*pow(g,-tmp,p)
c = bsgs(tmp1,tmp2,p,factors[i][0])
cs.append(c)
tmp += c*(p-1)//factors[i][0]
tmp //= factors[i][0]
residues.append(sum(c*factors[i][0]**k for k,c in enumerate(cs))%mods[i])
u = crt(residues,mods)
return u
u = cal(a)
v = cal(b)
if GCD(u,p-1) == 1:
return v*pow(u,-1,p-1)
else:
tmp = 1
while True:
gcd = GCD(u,p-1)
if gcd == 1:
break
tmp *= gcd
u //= gcd
return ZZ(v*pow(u,-1,p-1))//tmp
if __name__ == '__main__':
p = 7863166752583943287208453249445887802885958578827520225154826621191353388988908983484279021978114049838254701703424499688950361788140197906625796305008451719
a = getRandomNBitInteger(256)
x = getRandomNBitInteger(512)
b = pow(a,x,p)
print(x)
print(Pohlig_Hellman(a, b, p))2.椭圆曲线
大致思想和DLP相同
对于椭圆曲线上$P=kG$求解$k$时
先确定G基点的阶n，即nG=0，然后同样的对进行因式分解，得到$n=p_1^{q_1}p_2^{q_2}\cdots p_n^{q_n}$
相似的构造如下的同余方程组
$$ \left\{\begin{matrix} k=c_0+c_1p_1+c_2p_1^2+&\cdots&+c_{q_1-1}p_1^{q_1-1}\ mod\ p_1^{q_1}\\ k=c_0’+c_1’p_2+c_2’p_2^2+&\cdots&+c_{q_2-1}’p_2^{q_2-1}\ mod\ p_2^{q_2}\\ &\cdots&\\ k=c_0^{”}+c_1^{”}p_3+c_2^{”}p_3^2+&\cdots&+c_{q_3-1}^{”}p_3^{q_3-1}\ mod\ p_3^{q_3}\\ \end{matrix}\right. $$同理，构造此式子$\frac{n}{p_1^r}P=\frac{n}{p_1^r}nG\ mod\ p$
当$r=1$时有$\frac{n}{p_1}P=c_0\frac{n}{p_1}G$
爆破即可解出c然后依次遍历r
crt合并恢复k
# Sample
# BSGS & Pohlig_Hellman for ECDLP
from Crypto.Util.number import *
def bsgs(G,P,p):
tmp = ceil(sqrt(p))
bs = {}
for b in range(tmp):
bs[P-b*G] = b
tmp1 = G*tmp
for a in range(tmp):
if tmp1*a in bs:
return a*tmp+bs[tmp1*a]
def Pohlig_Hellman(G,P):
n = G.order()
factors = list(factor(n))
mods = []
residues = []
for i in range(len(factors)):
if factors[i][0] > 10**9:
break
tmp = 0
tmp1 = G*(n//factors[i][0])
cs = []
for r in range(1,factors[i][1]+1):
tmp2 = P*(n//factors[i][0]**r)-G*tmp
c = bsgs(tmp1,tmp2,factors[i][0])
cs.append(c)
tmp += c*n//factors[i][0]
tmp //= factors[i][0]
log = sum(c*factors[i][0]**k for k,c in enumerate(cs))%(factors[i][0]**factors[i][1])
print(log)
residues.append(log)
mods.append(factors[i][0]**factors[i][1])
k = crt(residues,mods)
return k感觉对于较大的数e7~e9计算有点慢，似乎是BSGS本身复杂度的问题，但我尝试用sage内置的对数运算却报错显示无解，真复杂，最后还是用自己的BSGS算了
Lazzaro大佬的板子，将BSGS的部分替换成sage内置的对数运算了
这个板子很快不知道为什么
E = EllipticCurve(GF(p), [a, b])
P =
Q =
n = E.order()
factors = list(factor(n))
m = 1
moduli = []
remainders = []
print(f"[+] Running Pohlig Hellman")
print(factors)
for i, j in factors:
if i > 10**9:
print(i)
break
mod = i**j
g2 = P*(n//mod)
q2 = Q*(n//mod)
r = discrete_log(q2, g2, operation='+')
remainders.append(r)
moduli.append(mod)
m *= mod
r = crt(remainders, moduli)
print(r)这个也很慢
E = EllipticCurve(GF(p), [a, b])
P =
Q =
factors, exponents = zip(*factor(E.order()))
primes = [factors[i] ^ exponents[i] for i in range(len(factors))]
print(primes)
dlogs = []
for fac in primes:
t = int(int(P.order()) // int(fac))
dlog = discrete_log(t*Q,t*P,operation="+")
dlogs += [dlog]
print("factor: "+str(fac)+", Discrete Log: "+str(dlog)) #calculates discrete logarithm for each prime order
l = crt(dlogs,primes)
print(l)Shor 量子算法
咕咕ing
常见攻击
Smart’s attack
条件：$E.order()=p$
MOV attack
条件：$E.order()=p+1$
sean神(๑•̀ㅁ•́ฅ)
我来喽!