Automne's Shadow.

LCG和模p平方根结合的题目

2020/01/26 Share

Crypto

代码task.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from Crypto.Util.number import *
from flag import flag


state = bytes_to_long(flag.encode())

p = getPrime(128)
q = getPrime(128)
m = p * q
a, b = [getRandomRange(0, m) for _ in range(2)]

with open('output', 'w') as f:
f.write(str(m) + '\n')
for i in range(6):
state = (a*state + b) % m
if i & 1 == 0:
f.write(str(state) + '\n')

输出output

1
2
3
4
55584485421406970301254549526604327461056610434113867867423369210521481715933
3693929062344887200388436902170299368192364992317574256526074702481265639773
23589540966217424118104948577678311766405141841319921607964145788046048746955
47328186237445429020074319580472604276783097851624136984720842711863174734878

代码比较简单易懂,这里有篇介绍LCG线性同余生成器的博文,LCG属于流密码的范畴。

https://zeroyu.xyz/2018/11/02/Cracking-LCG/

不同的是,题目中给到的输出状态不是连续的输出,无法直接套用。

假设初始state,也就是flag记为S0,于是有:

1
2
3
4
5
6
7
8
9
10
11
S1 ≡ (S0*a + b) % m

S2 ≡ (S1*a + b) % m

S3 ≡ (S2*a + b) % m

S4 ≡ (S3*a + b) % m

S5 ≡ (S4*a + b) % m

S6 ≡ (S5*a + b) % m

其中,S1,S3和S5是已知的,于是有:

1
S5-S3 ≡ (S3-S1)*a^2 % m

使用模逆运算计算a^2的值,注意模逆运算实际上使用算法是扩展欧几里得算法,关于欧几里得算法和扩展欧几里得算法可以参考下面的链接加深理解:

https://www.bilibili.com/video/av24248451?from=search&seid=8600306360895299893

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#encoding=utf-8

def egcd(a, b):
if a == 0:
return (b, 0, 1)
else:
g, x, y = egcd(b % a, a)
return (g, y - (b // a) * x, x)

def modinv(b, n):
g, x, _ = egcd(b, n)
if g == 1:
return x % n

def crack_unknown_multiplier(states, modulus):
multiplier = (states[2] - states[1]) * modinv(states[1] - states[0], modulus) % modulus
return multiplier

print crack_unknown_multiplier([3693929062344887200388436902170299368192364992317574256526074702481265639773, 23589540966217424118104948577678311766405141841319921607964145788046048746955, 47328186237445429020074319580472604276783097851624136984720842711863174734878], 55584485421406970301254549526604327461056610434113867867423369210521481715933)

求得:

1
2
3
a^2 ≡ 

29895665113034070831270870121097218673769521820755738885706707818988524697486 % m

直接开方求a的值是不好使的。搜索模p平方根找到了计算模平方根的算法:Tonelli–Shanks算法,从这里找到了这个算法的python实现:https://rosettacode.org/wiki/Tonelli-Shanks_algorithm#Python

但是直接带入

29895665113034070831270870121097218673769521820755738885706707818988524697486和55584485421406970301254549526604327461056610434113867867423369210521481715933

求解是报错的

automne

看一下算法约束条件可以发现要求模数p是奇质数

automne

调用sympy库查看,发现这里的p不是质数

automne

使用yafu分解p得到质因子P和Q:

1
2
3
P = 225198214939471987016766287890315801997

Q = 246824715890162716916904917666448620689

既然p不是质数,那么上面的算法就不适用,根据Rabin加密系统(维基百科搜Rabin cryptosystem)的加解密算法,可以进行求解。该算法要求质因子P和Q满足模4余3的条件。

但是这里的P和Q并不是模4余3的,尽管如此同样可以解,参考:https://xz.aliyun.com/t/5113

a^2 ≡

29895665113034070831270870121097218673769521820755738885706707818988524697486 % m

分解为两个同余式:

a^2 ≡

29895665113034070831270870121097218673769521820755738885706707818988524697486 % 225198214939471987016766287890315801997

a^2 ≡

29895665113034070831270870121097218673769521820755738885706707818988524697486 % 246824715890162716916904917666448620689

这里的两个同余式里的模数都是质数,所以是可解的

automne

二次同余方程解法的通用算法可以参考Cipolla’s algorithm,和上面的Tonelli-Shanks algorithm都是解决x^2≡n mod p问题的,但是对于不同的约束条件求解效率不同。

automne

通过Cipolla’s algorithm算法求出解之后,再使用中国剩余定理求解a的值,得到a的值之后,再根据LCG的特点:

1
S3 ≡ ((S1*a+b)*a + b) % m = (S1*a^2 + (a+1)*b) % m

1
(a+1)*b ≡ (S3-S1*a^2) % m

可求得b的值

然后再通过

1
S1 ≡ (S0*a + b) % m => a*S0 ≡ (S1-b) % m

求得S0,也就是flag。

于是最终的解题脚本:

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
from random import randint
from libnum import *
from gmpy2 import *
from Crypto.Util.number import *


def power(s1, s2, k1, k2, w, p):
return ((s1*k1+s2*k2*w)%p, (s1*k2+s2*k1)%p)

def Cipolla_algorithm(p, n):
a = randint(1, p)
w = a ** 2 - n

while pow(w, (p-1)/2, p) !=p-1 :
a = randint(1, p)
w = a ** 2 - n

times = (p+1)/2

k1 = 1
k2 = 0

first = True

sum1 = 1
sum2 = 0
while times != 0:
if first:
k1, k2=power(k1, k2, a, 1, w, p)
first = False

else:
k1, k2=power(k1, k2, k1, k2, w, p)

if times & 1:
sum1, sum2 = power(sum1, sum2, k1, k2, w, p)
times >>= 1

return sum1

def CRT(c, n):
for i in range(len(n)):
for j in range(i + 1, len(n)):
assert gcd(n[i], n[j]) == 1
assert len(c) == len(n)

N = reduce(lambda a, b: a*b, n)
x = 0
for i, j in zip(c, n):
N_i = N/j
N_i_1 = invert(N_i, j)
x+=i*N_i*N_i_1
return x % N

def egcd(a, b):
if a == 0:
return (b, 0, 1)
else:
g, x, y = egcd(b % a, a)
return (g, y - (b // a) * x, x)

def modinv(b, n):
g, x, _ = egcd(b, n)
if g == 1:
return x % n

if __name__ == '__main__':
p = 225198214939471987016766287890315801997
q = 246824715890162716916904917666448620689
n = p * q
c = 29895665113034070831270870121097218673769521820755738885706707818988524697486
s1 = 3693929062344887200388436902170299368192364992317574256526074702481265639773
s3 = 23589540966217424118104948577678311766405141841319921607964145788046048746955
s5 = 47328186237445429020074319580472604276783097851624136984720842711863174734878

get_x1 = Cipolla_algorithm(p, c)
get_x2 = Cipolla_algorithm(q, c)


assert pow(get_x1, 2, p) == c % p
assert pow(get_x2, 2, q) == c % q

c11 = get_x1
c12 = p-get_x1
c21 = get_x2
c22 = q-get_x2

print 'possible a :' + str(CRT([c11, c21], [p, q]))
print 'possible a :' + str(CRT([c11, c22], [p, q]))
print 'possible a :' + str(CRT([c12, c21], [p, q]))
print 'possible a :' + str(CRT([c12, c22], [p, q]))
a = [int(CRT([c11, c21], [p, q])),int(CRT([c11, c22], [p, q])),int(CRT([c12, c21], [p, q])),int(CRT([c12, c22], [p, q]))]
for i in a:
#print i
b = (s3-s1*c) * modinv(i+1,n) % n
#print b
flag = (s1-b) * modinv(i,n) % n
print long_to_bytes(flag)

automne

CATALOG