在这里零散地记载一些密码学可能用到的算法脚本（当裁缝），方便之后直接拿来用

RSA中的e,phi不互素问题

参考了这篇文章

换模

当e和p-1或q-1互素时，可以转换到模p或模q下求解

假设e与p-1互素
$$
m^e\equiv c(mod\ n)
$$

$$
m^e=c+kpq
$$

$$
m^e\ mod\ p=c\ mod\ p+kpq\ mod\ p
$$

$$
m^e\equiv c(mod\ p)
$$

python实现：

1
2
3
4
5
6
7
8
9
10
11
import gmpy2
e=
p=
q=
n=
c=
assert gmpy2.gcd(e,p-1)==1
c=c%p
phi=p-1
d=gmpy2.invert(e,phi)
m=pow(c,d,p)

e//gcd iroot(m,2)

当gcd(e,phi)较小时，可以先将e//gcd(e,phi)，使得e和phi互素后，再对算出的m开根

python实现：

1
2
3
4
5
6
7
8
9
10
11
import gmpy2
e=
p=
q=
n=
c=
phi=(p-1)*(q-1)
_gcd=gmpy2.gcd(e,phi)
e=e//_gcd
d=gmpy2.invert(e,phi)
m=gmpy2.iroot(pow(c,d,p),_gcd)[0]

有限域内开方

有时，当e较小时，我们仍然无法用上面的方法得到m。这时，我们可以使用有限域内开方的方法。

前面我们已经证明了，当p为素数时，由
$$
m^e\equiv c(mod\ n)
$$
可以推知
$$
m^e\equiv c(mod\ p)
$$
于是我们可以设法找到所有满足
$$
m^e\equiv x^e\equiv c(mod\ p)
$$
的x，以及所有满足
$$
m^e\equiv y^e\equiv c(mod\ q)
$$
的y，从而列出
$$
\begin{cases}m\equiv x(mod\ p) \ m\equiv y(mod\ q)\end{cases}
$$
的方程组，通过中国剩余定理（CRT）求解。

下面是一个利用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
27
from Crypto.Util.number import *
p = 
q = 
c = 
e = 
n = 

P.<a>=PolynomialRing(Zmod(p),implementation='NTL')
f=a^e-c
mps=f.monic().roots()

P.<a>=PolynomialRing(Zmod(q),implementation='NTL')
g=a^e-c
mqs=g.monic().roots()

flag=[]
for mpp in mps:
    x=mpp[0]
    for mqq in mqs:
        y=mqq[0]
        solution = CRT_list([int(x), int(y)], [p, q])
        flag.append(solution)
for i in flag:
    m=long_to_bytes(i)
    if b'flag'in m:
        print(m)

当然也可以通过观察对算法进行一定的优化，详见上面给出的文章，这里不再赘述。

AMM算法

参考了这篇文章

当e比较大的时候，我们可以使用AMM算法，它能够大大提高运算的速度

原算法

先说明一下AMM算法的原算法，此时e=2，p为奇素数，对于
$$
x^2\equiv r(mod\ p)
$$
我们先对两边开根
$$
x\equiv r^{\frac{1}{2}}(mod\ p)
$$
令p-1=2^{s}t，又由欧拉定理得，
$$
r^{\frac{p-1}{2}}\equiv r^{2^{s-1}t}\equiv 1(mod\ p)
$$

当s=1时，

$$
r^t\equiv 1(mod\ p)
$$

两边同时乘r，再开根，即推出公式
$$
r^{\frac{t+1}{2}}\equiv \pm \sqrt r\equiv \pm x(mod\ p)
$$
将m和c代进去，就是这样(t可以根据p算出来)
$$
\pm m\equiv c^{\frac{t+1}{2}}
$$

当s>1时，

如果直接开根我们会得到一正一负两个式子
$$
r^{2^{s-2}t}\equiv 1(mod\ p)
$$

$$
r^{2^{s-2}t}\equiv -1(mod\ p)
$$

由上面得到的
$$
r^{2^{s-1}t}\equiv 1(mod\ p)
$$
两边同时乘上这个式子，并用k来控制是否要乘（k=0,1）
$$
r^{2^{s-2}t}n^{2^{s-1}tk}\equiv 1(mod\ p)
$$
就这样反复对两边进行开方操作，直至回到前面s=1的情况，即
$$
r^tn^{t*(2k_1+2^2k_2+…+2^{s-1}k_{s-1})}\equiv 1(mod\ p)
$$
两边乘上r再开方
$$
r^{\frac{t+1}{2}}n^{t*(k_1+2k_2+…+2^{s-2}k_{s-1})}\equiv\pm \sqrt r\equiv \pm x(mod\ p)
$$
将m和c代进去，得到
$$
c^{\frac{t+1}{2}}n^{t*(k_1+2k_2+…+2^{s-2}k_{s-1})}\equiv \pm m(mod\ p)
$$

e>2

对于
$$
x^e\equiv r(mod\ p)
$$
令p-1=e^{s}t，则有
$$
r^{\frac{p-1}{e}}\equiv r^{e^{s-1}t}\equiv 1(mod\ p)
$$
此时可以找到δ，使得t|(eδ-1)，则
$$
r^{e^{s-1}(eδ-1)}\equiv 1(mod\ p)
$$

当s=1时，

$$
r^{(eδ-1)}\equiv 1(mod\ p)
$$

两边乘r，再开e次方
$$
r^\delta\equiv r^{\frac{1}{e}}\equiv x(mod\ p)
$$

当s>1时，

构造e次非剩余集合
$$
K_i=\rho^{i\frac{p-1}{e}}=\rho^{ie^{s-1}t},0\leq i\leq e-1
$$

$$
K_i^e=\rho^{ie^st}=\rho^{i(p-1)}
$$

所以根据欧拉定理，得
$$
\rho^{i*(p-1)}\equiv \rho^{(p-1)}\equiv1(mod\ p)
$$
由上面的式子又可以推知
$$
\begin{cases}K_i=\rho^{\frac{i*(p-1)}{e}} \ K_{e-i}=\rho^{\frac{(e-i)*(p-1)}{e}}\end{cases}
$$

$$
K_iK_{e-i}=\rho^{p-1}
$$
由欧拉定理又可以得到
$$
K_iK_{e-i}\equiv 1(mod\ p)
$$
所以K_i和K_{e-i}互为逆元

对于前面这个式子
$$
r^{e^{s-1}(eδ-1)}\equiv 1(mod\ p)
$$
两边开e次方得到一个集合K中的数设为K_{e-j}
$$
r^{e^{s-2}(eδ-1)}\equiv K_{e-j}(mod\ p)
$$
两边乘上K_j然后开e次方
$$
r^{e^{s-2}(eδ-1)}K_j\equiv K_{e-j}K_j\equiv 1(mod\ p)
$$

$$
r^{e^{s-2}(eδ-1)}\rho^{j*e^{s-1}t}\equiv 1(mod\ p)
$$

反复进行上述操作，直至回到s=1的情况
$$
r^{(eδ-1)}\rho^{ej_1+e^2j_2+…+e^{s-1}j_{s-1}}\equiv 1(mod\ p)
$$
两边乘r，开e次方
$$
r^δ\rho^{j_1+ej_2+…+e^{s-2}j_{s-1}}\equiv r^{\frac{1}{e}}\equiv x(mod\ p)
$$
将m和c代进去，得到
$$
c^δ\rho^{j_1+ej_2+…+e^{s-2}j_{s-1}}\equiv m(mod\ p)
$$
此时我们便得到了其中一个根，剩余的根可以通过不断乘上集合K得到.

当我们得到了所有的解以后，使用中国剩余定理对下面的方程组求解即可
$$
\begin{cases} m^e\equiv cp(mod\ p) \ m^e\equiv cq(mod\ q)\end{cases}
$$
python代码实现：

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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import random 
import time 
from Crypto.Util.number import * # pip install pycryptodome
import gmpy2 # pip install gmpy2
from tqdm import tqdm # pip install tqdm

def onemod(e, q):
    p = random.randint(1, q - 1)
    while ((gmpy2.powmod(p, (q - 1) // e, q)) == 1):
        p = random.randint(1, q)
    return p

def AMM_rth(x, e, N):
    assert ((N - 1) % e == 0)
    p = onemod(e, N)
    
    t = 0
    s = N - 1
    while (s % e == 0):
        s = s // e
        t += 1
    k = 1
    
    while ((s * k + 1) % e != 0):
        k += 1
    delta = (s * k + 1) // e

    a = gmpy2.powmod(p, e ** (t - 1) * s, N)
    b = gmpy2.powmod(x, e * delta - 1, N)
    c = gmpy2.powmod(p, s ,N)
    h = 1
    for i in tqdm(range(1, t)):
        d = gmpy2.powmod(b, e ** (t - 1 -i), N)
        if d == 1:
            j = 0
        else:
            j = (- math.Log(d, a)) % e
        b = b * (c ** (e ** j)) % N
        h = h * (c ** j) % N
        c = c ** e % N
    result = gmpy2.powmod(x, delta, N) * h
    return result

def ALL_ROOT2(r, q):  
    list1 = set()
    while(len(list1) < r):
        p = gmpy2.powmod(random.randint(1, q - 1), (q - 1) // r, q)
        list1.add(p)
    return list1

def ALL_Solution(m, q, rt, cq, e):
    mp = []
    for pr in rt:
        r = (pr * m) % q
        mp.append(r)
    return mp

def calc(mp, mq, e, p, q):
    i = 1
    j = 1
    t1 = gmpy2.invert(q, p)
    t2 = gmpy2.invert(p, q)
    for mp1 in mp:
        for mq1 in mq:
            j += 1
            if j % 100000 == 0:
                print(j)
            ans = (mp1 * t1 * q + mq1 * t2 * p) % (p * q)
            if check(ans):
                return
    return

def check(m):
    try:
        a = long_to_bytes(m)
        if flag_info in a:
            print(a)
            return True
        else:
            return False
    except:
        return False

n = 
c = 
p = 
q = 
e = 
flag_info = b'flag{'

cp = c % p
cq = c % q

print("start")

mp = AMM_rth(cp, e, p)
mq = AMM_rth(cq, e, q)

rt1 = ALL_ROOT2(e, p)  
rt2 = ALL_ROOT2(e, q)

amp = ALL_Solution(mp, p, rt1, cp, e)  
amq = ALL_Solution(mq, q, rt2, cq, e)

calc(amp, amq, e, p, q)
print("run over")

当然，这里还得介绍一下sagemath中有一个很好用的方法.nth_root，可以非常有效地完成域下的高次开根

大概用法如下，可视情况做出修改

1
2
K=Zmod(n)
x=K(c).nth_root(e,all=True)

这样可以返回所有在模n整数环下，满足x^e ≡ c (mod n)的x

共模攻击

用于解决在GF(n)下，明文相同，公钥不同，而模n又很大难以分解的情况（两组公钥和密文已知）
$$
\begin{cases}m^{e_1}\equiv c_1(mod\ n) \ m^{e_2}\equiv c_2(mod\ n)\end{cases}
$$
两边分别同时乘s1,s2次方
$$
\begin{cases}m^{e_1s_1}\equiv c_1^{s_1}(mod\ n) \ m^{e_2s_2}\equiv c_2^{s_2}(mod\ n)\end{cases}
$$
两式相乘
$$
m^{e_1s_1+e_2s_2}\equiv c_1^{s_1}c_2^{s_2}(mod\ n)
$$
这里使用扩展欧几里得算法，我们可以找到能满足e1s1+e2s2=1的s1和s2。又因为一般来说m<n，所以
$$
m=c_1^{s_1}c_2^{s_2}\ mod\ n
$$
代码实现：

1
2
3
4
5
6
7
8
import gmpy2
n=
e1=
e2=
c1=
c2=
_gcd,s1,s2=gmpy2.gcdext(e1,e2)
m=pow(c1,s1,n)*pow(c2,s2,n)%n

维纳攻击

参考这篇文章

用于RSA中e很大的时候

代码实现：

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
import gmpy2
def continuedFra(x, y):
    cf = []
    while y:
        cf.append(x // y)
        x, y = y, x % y
    return cf
def gradualFra(cf):
    numerator = 0
    denominator = 1
    for x in cf[::-1]:
        numerator, denominator = denominator, x * denominator + numerator
    return numerator, denominator
def solve_pq(a, b, c):
    par = gmpy2.isqrt(b * b - 4 * a * c)
    return (-b + par) // (2 * a), (-b - par) // (2 * a)
def getGradualFra(cf):
    gf = []
    for i in range(1, len(cf) + 1):
        gf.append(gradualFra(cf[:i]))
    return gf
def wienerAttack(e, n):
    cf = continuedFra(e, n)
    gf = getGradualFra(cf)
    for d, k in gf:
        if k == 0: continue
        if (e * d - 1) % k != 0:
            continue
        phi = (e * d - 1) // k
        p, q = solve_pq(1, n - phi + 1, n)
        if p * q == n:
            return d
n=
e=
c=
d=wienerAttack(e, n)
m=pow(c, d, n)

Boneh Durfee攻击

类似维纳攻击，但可使d范围更大

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
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
import time
from Crypto.Util.number import *
debug = True
strict = False
helpful_only = True
dimension_min = 7
def helpful_vectors(BB, modulus):
    nothelpful = 0
    for ii in range(BB.dimensions()[0]):
        if BB[ii,ii] >= modulus:
            nothelpful += 1
 
    print (nothelpful, "/", BB.dimensions()[0], " vectors are not helpful")
def matrix_overview(BB, bound):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            a += '0' if BB[ii,jj] == 0 else 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        if BB[ii, ii] >= bound:
            a += '~'
        print (a)
def remove_unhelpful(BB, monomials, bound, current):
    if current == -1 or BB.dimensions()[0] <= dimension_min:
        return BB
    for ii in range(current, -1, -1):
        if BB[ii, ii] >= bound:
            affected_vectors = 0
            affected_vector_index = 0
            for jj in range(ii + 1, BB.dimensions()[0]):
                if BB[jj, ii] != 0:
                    affected_vectors += 1
                    affected_vector_index = jj
            if affected_vectors == 0:
                print ("* removing unhelpful vector", ii)
                BB = BB.delete_columns([ii])
                BB = BB.delete_rows([ii])
                monomials.pop(ii)
                BB = remove_unhelpful(BB, monomials, bound, ii-1)
                return BB
            elif affected_vectors == 1:
                affected_deeper = True
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):
                    if BB[kk, affected_vector_index] != 0:
                        affected_deeper = False
                if affected_deeper and abs(bound - BB[affected_vector_index, affected_vector_index]) < abs(bound - BB[ii, ii]):
                    print ("* removing unhelpful vectors", ii, "and", affected_vector_index)
                    BB = BB.delete_columns([affected_vector_index, ii])
                    BB = BB.delete_rows([affected_vector_index, ii])
                    monomials.pop(affected_vector_index)
                    monomials.pop(ii)
                    BB = remove_unhelpful(BB, monomials, bound, ii-1)
                    return BB
    return BB
def boneh_durfee(pol, modulus, mm, tt, XX, YY):
    PR.<u, x, y> = PolynomialRing(ZZ)
    Q = PR.quotient(x*y + 1 - u) # u = xy + 1
    polZ = Q(pol).lift()
    UU = XX*YY + 1
    gg = []
    for kk in range(mm + 1):
        for ii in range(mm - kk + 1):
            xshift = x^ii * modulus^(mm - kk) * polZ(u, x, y)^kk
            gg.append(xshift)
    gg.sort()
    monomials = []
    for polynomial in gg:
        for monomial in polynomial.monomials():
            if monomial not in monomials:
                monomials.append(monomial)
    monomials.sort()
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            yshift = y^jj * polZ(u, x, y)^kk * modulus^(mm - kk)
            yshift = Q(yshift).lift()
            gg.append(yshift)
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            monomials.append(u^kk * y^jj)
    nn = len(monomials)
    BB = Matrix(ZZ, nn)
    for ii in range(nn):
        BB[ii, 0] = gg[ii](0, 0, 0)
        for jj in range(1, ii + 1):
            if monomials[jj] in gg[ii].monomials():
                BB[ii, jj] = gg[ii].monomial_coefficient(monomials[jj]) * monomials[jj](UU,XX,YY)
    if helpful_only:
        BB = remove_unhelpful(BB, monomials, modulus^mm, nn-1)
        nn = BB.dimensions()[0]
        if nn == 0:
            print ("failure")
            return 0,0
    if debug:
        helpful_vectors(BB, modulus^mm)
    det = BB.det()
    bound = modulus^(mm*nn)
    if det >= bound:
        print ("We do not have det < bound. Solutions might not be found.")
        print ("Try with highers m and t.")
        if debug:
            diff = (log(det) - log(bound)) / log(2)
            print ("size det(L) - size e^(m*n) = ", floor(diff))
        if strict:
            return -1, -1
    else:
        print ("det(L) < e^(m*n) (good! If a solution exists < N^delta, it will be found)")
    if debug:
        matrix_overview(BB, modulus^mm)
    if debug:
        print ("optimizing basis of the lattice via LLL, this can take a long time")
    BB = BB.LLL()
    if debug:
        print ("LLL is done!")
    if debug:
        print ("looking for independent vectors in the lattice")
    found_polynomials = False
    for pol1_idx in range(nn - 1):
        for pol2_idx in range(pol1_idx + 1, nn):
            PR.<w,z> = PolynomialRing(ZZ)
            pol1 = pol2 = 0
            for jj in range(nn):
                pol1 += monomials[jj](w*z+1,w,z) * BB[pol1_idx, jj] / monomials[jj](UU,XX,YY)
                pol2 += monomials[jj](w*z+1,w,z) * BB[pol2_idx, jj] / monomials[jj](UU,XX,YY)
            PR.<q> = PolynomialRing(ZZ)
            rr = pol1.resultant(pol2)
            if rr.is_zero() or rr.monomials() == [1]:
                continue
            else:
                print ("found them, using vectors", pol1_idx, "and", pol2_idx)
                found_polynomials = True
                break
        if found_polynomials:
            break
    if not found_polynomials:
        print ("no independant vectors could be found. This should very rarely happen...")
        return 0, 0
    rr = rr(q, q)
    soly = rr.roots()
    if len(soly) == 0:
        print ("Your prediction (delta) is too small")
        return 0, 0
    soly = soly[0][0]
    ss = pol1(q, soly)
    solx = ss.roots()[0][0]
    return solx, soly
def example(N,e,delta):
    t = int((1-2*delta) * m)
    X = 2*floor(N^delta)
    Y = floor(N^(1/2))
    P.<x,y> = PolynomialRing(ZZ)
    A = int((N+1)/2)
    pol = 1 + x * (A + y)
    if debug:
        print ("=== checking values ===")
        print ("* delta:", delta)
        print ("* delta < 0.292", delta < 0.292)
        print ("* size of e:", int(log(e)/log(2)))
        print ("* size of N:", int(log(N)/log(2)))
        print ("* m:", m, ", t:", t)
    if debug:
        print ("=== running algorithm ===")
        start_time = time.time()
 
    solx, soly = boneh_durfee(pol, e, m, t, X, Y)
    if solx > 0:
        print ("=== solution found ===")
        if False:
            print ("x:", solx)
            print ("y:", soly)
 
        d = int(pol(solx, soly) / e)
        print ("private key found:", d)
    else:
        print ("=== no solution was found ===")
 
    if debug:
        print("=== %s seconds ===" % (time.time() - start_time))
    return d
 
if __name__ == "__main__":
    c =
    n =
    e =
    delta = 0.28 # this means that d < N^delta
    d = example(n,e,delta)
    print(long_to_bytes(int(pow(c,d,n))))

CBC字节翻转攻击

参考这篇文章

用于CBC模式下的AES加密

个人认为上面的文章中这一段讲得已经很清楚了

image-20240721165239569

实现代码：

1
2
3
4
5
6
7
8
def strxor(a1, a2): 
    return bytes([b1 ^ b2 for b1,b2 in zip(a1,a2)])
A=b''
C=b''
CC=b''
AA=strxor(A,C)
AA=strxor(AA,CC)
print(AA)

MT19937伪随机数预测

MT19937，即梅森旋转算法，是一种伪随机数的生成算法，python中的random模块生成“随机数”时使用的就是这种算法。

根据其原理，有一个叫randcrack的python模块可以对其生成的“随机数”进行预测，前提是需要已知已经生成的至少624个32位二进制数，这样才能预测出下一个生成的数会是多少。

实现代码：

1
2
3
4
5
6
7
from randcrack import RandCrack
rc = RandCrack()
Rand=[] #624个32位二进制数
for i in range(624):
    rc.submit(Rand[i])
pre=rc.predict_getrandbits(100) #假设要预测的是下一个100位二进制数
print(pre)

那么前面给出的是312个64位数呢？

其实，random生成64位数的方式就是先生成2个32位数，然后将它们拼起来得到的。因此我们只需要把这个64位数拆开成两个32位数即可。

实现代码：

1
2
3
4
5
6
7
8
9
10
11
from randcrack import RandCrack
rc = RandCrack()
Rand= [] #312个64位数
prng=[]
for i in Rand:
    prng.append(int(i)& (2 ** 32 - 1))
    prng.append(int(i)>> 32)
for i in range(624):
    rc.submit(prng[i])
notrandom=rc.predict_getrandbits(100) #假设要预测的是下一个100位二进制数
print(notrandom)

LCG

参考这篇文章和这篇文章

LCG，线性同余法，是一种生成伪随机数的方法，用一个公式来概括就是
$$
X_{n+1}=(aX_n+b)\ mod\ m
$$
基本围绕以下四个公式，即可解决各类基础的LCG问题

c00bad18-4f1e-4719-81d0-526a8db7030c