椭圆曲线

咕咕ing

常用算法

ps.下述的算法在不同领域有很多不同的变体但核心思想是一致的

1. 大整数的分解
2. 离散对数的求解
3. 椭圆曲线上的离散对数求解

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 factors

2.离散对数

咕咕ing

3.椭圆曲线

咕咕ing

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$