[TCTF/0CTF2021 final]EzHash 复现

tctf final2021唯一没做出来的题，属于是毫无思路，看了paper越发觉得思路很有借鉴意义，复现一下

题目核心代码

class ezhash:
    def __init__(self):
        self.p = 328145541634163431929233386535821861121
        self.i = 96858741233065729363697642301435517926
        assert pow(self.i, 2, self.p) == self.p - 1
        self.l = 5
        self.g = [
            [[1, -2 * self.i], [-2 * self.i, 1]],
            [[1, -2], [2, 1]],
            [[1 - 2 * self.i, 0], [0, 1 + 2 * self.i]],
            [[1 + 2 * self.i, 0], [0, 1 - 2 * self.i]],
            [[1, 2], [-2, 1]],
            [[1, 2 * self.i], [2 * self.i, 1]]
        ]
        self.pi = [
            [1, 0, 5, 4, 3, 2],
            [0, 5, 4, 3, 2, 1],
            [2, 1, 0, 5, 4, 3],
            [4, 3, 2, 1, 0, 5],
            [3, 2, 1, 0, 5, 4]
        ]
        self.x = [[getRandomRange(1, self.p) for _ in range(2)] for _ in range(2)]
    def base(self, n):
        seq = []
        while n:
            seq.append(n % self.l)
            n //= self.l
        return seq
    def digest(self, msg):
        seq = self.base(int(msg.hex(), 16))
        g_last = self.g[1]
        r = [[1, 0], [0, 1]]
        for _ in range(len(seq)):
            g_cur = self.g[self.pi[seq[_]][self.g.index(g_last)]]
            (a, b), (c, d) = r[0], r[1]
            (e, f), (g, h) = g_cur[0], g_cur[1]
            r[0][0] = (a * e + b * g) % self.p
            r[0][1] = (a * f + b * h) % self.p
            r[1][0] = (c * e + d * g) % self.p
            r[1][1] = (c * f + d * h) % self.p
            g_last = g_cur
        return r
    def check(self, msg):
        d = self.digest(msg)
        s = [inverse(d[0][0], self.p) * self.x[0][0] % self.p, \
        inverse(d[0][1], self.p) * self.x[0][1] % self.p, \
        inverse(d[1][0], self.p) * self.x[1][0] % self.p, \
        inverse(d[1][1], self.p) * self.x[1][1] % self.p]
        if len(list(set(s))) == 1:
            return True
        else:
            return False
class Task(socketserver.BaseRequestHandler):
    def handle(self):
        try:
            signal.signal(signal.SIGALRM, self.timeout_handler)
            signal.alarm(30)
            if not self.proof_of_work():
                self.dosend('You must pass the PoW!')
                return
            signal.alarm(10)
            h = ezhash()
            self.dosend('x: ' + str(h.x))
            msg = self.request.recv(1024).strip()
            msg = bytes.fromhex(msg.decode('utf-8'))
            if h.check(msg):
                self.dosend(flag)
            else:
                self.dosend('Good luck!')

删掉了无关紧要的交互函数，大致是要在1024位内，将一个随机矩阵用6个矩阵基相乘表示。

比赛中的尝试

首先，这里有三个常数，p为模，由于后面也可能涉及离散对数，分解一下p-1,发现并不算很光滑。i² = -1 mod p。结合i的符号，很容易让人想到i相当于模p下的虚数单位i，也就是说这个矩阵我们不该将其看作是模p的普通矩阵，而应当将其看成一个复数矩阵。l = 5，也就是这6个矩阵的行列式（毕竟后面要将六个矩阵相乘，乘法中矩阵特征最容易计算的也就是行列式了）。

观察矩阵，可以发现这些矩阵很好的对称性质，对角线上的两个数共轭，反对角线上的两个数实部相反，虚部相同。很容易想到这些矩阵也许成群，将其乘一下果然仍然有这一性质。因此，第一步肯定是很明朗的，需要将要求的随机矩阵化为上述形式，这一步很简单就能够做到。其次，普通矩阵需要四个数字才能够表示，但是这无法表示我们说的对称性质，因此将矩阵用第一列的两个数的实部和虚部的四个数字组成的四元向量来表示，能够表现出上面的对称性质，尽管正常的矩阵仍然需要四个数字，但是我们的基都能够有两个位置是0，能够更方便。但接下来的分解就非常困难，由于是乘法，在模p下，离散对数问题绝对是大问题，没法解决，即使能够算出，往往需要的矩阵个数也远远超过1024了。

比赛后的思考

比赛后才知道是paper题，而且这真的是个在md5被攻破以后提出又被废弃的hash方案LPS，查hash就能查到（鬼知道啊，这一看还以为只是考线代复变的，以为题目乱取的，居然还真有人拿这当hash方案）首先是dbt在discord找到了出题人说参考的《Full Cryptanalysis of LPS and Morgenstern Hash Functions》^[1]但这篇很多都没讲清楚，在参考文献中又找到了一篇《Collisions for the LPS expander graph hash function》^[2]两篇结合加上一点自己的想法就可以做出这道题了。第二篇是注重碰撞，而碰撞只需要分解单位矩阵I即可，而第一个注重找原像，也就是题目中的分解。

哈希方案LPS

题目中选取的G为PSL(2,Fp)（模p的2*2矩阵），S的选取为那6个矩阵。

这六个矩阵是由一条方程组生成的

该方程共有\(\sum_{d|l}d\)个解，证明可见参考文献[3]中的第二章内容。在这里l为5,也就是有6个解，则六个矩阵则由这6个解按照如下方式生成

当我们将kM看作于M相等时，这些矩阵满足上面的条件。

分解思路

不管是寻求该哈希函数的碰撞，还是找一个输出的原像，它们本质上都是寻找一个矩阵的分解。（找碰撞等价于找单位矩阵的分解，找原像则是找任意矩阵的分解也就是题目中我们要做的）但如果在PSL(2,Fp)上，离散对数问题是没法解决的，所以文献[2]中提供的思路是将其转化到Z[i]上，在Z[i]上进行分解，分解完后便自然在PSL(2, Fp)上成立。因此问题转化为两步，即PSL到Z[i]的转化与Z[i]上的分解。

分解前准备的工具函数与变量

题目中的全局变量

p = 328145541634163431929233386535821861121
i = 96858741233065729363697642301435517926
l = 5
g = [
    [ 1, 0, 0,-2],
    [ 1, 0,-2, 0],
    [ 1,-2, 0, 0],
    [ 1, 2, 0, 0],
    [ 1, 0, 2, 0],
    [ 1, 0, 0, 2]]
ZERO = [1 , 0 , 0 , 0]

这里g是将矩阵转化为四元向量后的形式，并且g[i] = g[5-i]^-1，后面求逆很方便。

矩阵两种形式互相转化的函数

def PSLp_to_Zi(M):
    a = (M[0][0] + M[1][1]) * inverse(2 , p) % p
    b = (M[0][0] - M[1][1]) * inverse(2*i , p) % p
    c = (M[1][0] - M[0][1]) * inverse(2 , p) % p
    d = (M[1][0] + M[0][1]) * inverse(2*i , p) % p
    return [a,b,c,d]
def Zi_to_PSLp(M):
    a,b,c,d = M
    res = [[0,0] , [0,0]]
    res[0][0] = (a+b*i)%p
    res[0][1] = (-c+d*i)%p
    res[1][0] = (c+d*i)%p
    res[1][1] = (a-b*i)%p
    return res

乘法（选择用四元向量的形式）

def mul(A ,B):
    a1,b1,c1,d1 = [int(i) for i in A]
    a2,b2,c2,d2 = [int(i) for i in B]
    a3 = a1*a2 - b1*b2 - c1*c2 - d1*d2
    b3 = a1*b2 + a2*b1 - c1*d2 + c2*d1
    c3 = c1*a2 - d1*b2 + a1*c2 + b1*d2
    d3 = c1*b2 + d1*a2 + n - b1*c2
    return [a3,b3,c3,d3]

Z[i]上的分解

首先需要知道解决这个问题的可行性。即为什么在Z[i]上分解会比在PSL(2, Fp)上容易。粗略的看，我们需要分解的是一个随机矩阵，本来就不是由6个矩阵基相乘得到的，有很大概率是无解的，但其实却不然，文献[3]中给出了一个定理，也就是满足如下条件的矩阵，一定能够被l+1个矩阵基分解。（证明方法还没看到，粗略的看应该是把复数域再扩展了一下，扩展出一个四元域，然后在上面证明，过段时间再补充）

如果我们知道有解，设F = g_k0g_k1...g_kn,然后我们爆破六个矩阵基，分别计算Fg_i的值，当 g_i = g_kn^-1时，F = 5g_k0g_k1...g_kn-1，此时F中的四个数均能被5整除，即使存在i使得g_i != g_kn^-1时，F也能被5整除，那此时F的行列式仍然在下降，不断下降一定会下降到1.

代码很简单

def factor_Zi(M):
    flist = []
    testM = [i for i in M]
    while testM.count(0) < 3:
        for i in range(6):
            temp = mul(g[i], testM)
            for j in range(4):
                if temp[j] % 5 != 0:
                    break
            else:
                flist.append(5 - i)
                testM = [i // 5 for i in temp]
                break
        else:
            return -1
    return flist

这里的testM最后应当是[a , 0, 0, 0]的形式。

PSL到Z[i]的转化

首先上面的分解提供了我们需要转化的最终结果，所以这一步的需求是，先将随机矩阵转化为四元向量 (A, B, C ,D),然后求解方程

只要找到满足条件的w，x，y，z和k就行（2k其实无所谓）

显然这么多未知数相当难求解，主要是因为wzxy都有一次项与二次项，所以我们先考虑简单情况，即四元向量中，有两个位置是0的情况。

不妨设C = D = 0

原式化为 \[ (A^2+B^2)\lambda^2-l^{2k}+2\lambda(Aw+Bx)p + (w^2+x^2+y^2+z^2)p^2=0 \] 由于我们需要其他项全是正数，因此l^2k需要较大，但k大了会导致分解后矩阵数增多，因此取k = log_l(8p²)左右即可。

可以看到对p的三个次数分别有三个括号，模p和p²依次解决。

模p \[ \lambda = (\frac{l^{2k}}{A^2+B^2})^{\frac{1}{2}} mod\ p \] 由于l=5是模p的二次剩余，所以这里必须要求A² + B²是模p的二次剩余才能有解。

接着模p²有两个未知数。可以随便选取一个w，这里必须要让\(A\lambda+wp\)是奇数，然后计算出x \[ x = \frac{-\frac{(A^2+B^2)\lambda^2-l^{2k}}{p}-2\lambda Aw}{2\lambda B} mod\ p \] 将这些全部代回，只剩最后两个未知数y和z，它们需要满足(n太长了懒得写了，反正就是剩下的) \[ y^2 + z^2 = n \] 可以证明，只要n模4余1,这条式子一定有解，（模4余3一定无解）^[3]。

当n模4余1时，这是一个经典问题，可以用连分数来求解y与z

由于需要计算-1开根，素数会相当容易，对于一般的N需要分解，如果数据过大可以撞出素数N即可。这里我是通过随机的x直接撞素数N

class ContinuedFraction():
    def __init__(self , numerator, denumerator):
        self.numberlist = []    #number in continued fraction
        self.fractionlist = []  #the near fraction list
        self.GenerateNumberList(numerator , denumerator)
        self.GenerateFractionList()
    def GenerateNumberList(self,numerator,denumerator):
        while numerator != 1:
            quotient = numerator//denumerator
            remainder = numerator % denumerator
            self.numberlist.append(quotient)
            numerator = denumerator
            denumerator = remainder
    def GenerateFractionList(self):
        self.fractionlist.append([self.numberlist[0] ,1 ])
        for i in range(1,len(self.numberlist)):
            numerator = self.numberlist[i]
            denumerator = 1
            for j in range(i):
                temp = numerator
                numerator = denumerator + numerator * self.numberlist[i-j-1]
                denumerator = temp
            self.fractionlist.append([numerator,denumerator])
def compute_square_modp(a , p):
    assert is_prime(p)
    if pow(a , (p-1)//2 , p) != 1:
        return None
    R = GF(p)
    return int(R(a).nth_root(2))
def solve_yz(n):
    if not is_prime(n) or n % 4 != 1:
        return -1
    else:
        r = compute_square_modp(-1 , n)
        print(1)
        cf = ContinuedFraction(r , n)
        square = iroot(n , 2)[0]

        for i in range(len(cf.fractionlist)):
            if cf.fractionlist[i][1] <= square and cf.fractionlist[i+1][1] >= square:
                p , q = cf.fractionlist[i]
                a = q
                b = q*r - p*n
                assert a **2 + b ** 2 == n
                return a , abs(b)
def convert_to_Zi(A):
    assert A.count(0) == 2
    nonzero = []
    zero = []
    temp = []
    for i in range(4):
        if A[i] == 0:
            zero.append(i)
        else:
            nonzero.append(i)
            temp.append(A[i])
    a ,b = temp
    k = 224
    if pow(a**2 + b**2 , (p-1)//2 , p) != 1:
        return -1
    lamda = compute_square_modp((l**k) * inverse(a**2+b**2 , p) , p)
    m = (l**k - (a**2+b**2)*lamda**2) // p
    while 1:
        w = randint(1 , p-1)
        if (a*lamda + w) % 2 == 0:
            continue
        x = (m * inverse(2*lamda*b , p) - a * w * inverse(b , p)) % p
        n = (l ** k - (a*lamda +  w*p)**2 - (b*lamda + x*p)**2)// p**2
        if n%4 != 0 or n < 0:
            continue
        n //= 4
        temp = solve_yz(n)
        if temp != -1:
            y , z = temp
            ab = [a*lamda+w*p , b*lamda + x*p]
            res = [0] * 4
            for i in range(4):
                if i in zero:
                    res[i] = temp[zero.index(i)]*2*p
                else:
                    res[i] = ab[nonzero.index(i)]
            return res

这里同时考虑了0的位置不同的问题，开头和结尾是用于B = D = 0或C = D = 0的情况的。但由于这里有很多随机处（二次剩余与n的素数），因此这个函数可能跑的很慢也有可能很快。

有了这个以后，接下来如果我们能够用对角矩阵将任意矩阵变成对角矩阵，就能够使用上面的方法将这些矩阵一一分解。文献[1]使用的是这种方法，使用了7个矩阵来表示任意矩阵(4个对角矩阵与三个基矩阵)但我尝试以后，总是找不到能够同时满足所有的对角矩阵A²+B²为二次剩余的情况，不知道是方法问题还是写出bug了。但我尝试了另一条路。上面的方法并不是一定需要对角矩阵的，只要两个位置是0即可，而对角线相等，反对角线相反的矩阵刚好是B = D = 0的情况，而这种矩阵是很容易构造的。

当\(\alpha,w\)满足 \[ \alpha w = \frac{M_4}{M_1}\\ -\frac{w}{\alpha}=\frac{M_2}{M_3}\\ =>-w^2=\frac{M_4M_2}{M_1M_3} \] 如果右边的是一个二次剩余，则我们能够算出\(\alpha\)与w，从而用三个矩阵来表示出来。如果它不是二次剩余，我们可以通过左乘一个基矩阵，来打散M，由于有加法，因此结果是否为二次剩余是不确定的，粗略看大致为1/2的概率同时，w和a也必须是二次剩余，这里的随机性也得通过左乘来提供。(在尝试中，有些矩阵使用这种方式似乎无效，而其他的矩阵大约在3-4次左右能够成功。)

def diagonalize(M):
    left = 0
    while 1:
        print(left)
        m1 , m2 = M[0]
        m3 , m4 = M[1]
        temp = -m4 * m2 * inverse(m1 * m3 , p) % p
        if pow(temp , (p-1)//2 , p) == 1:
            w = compute_square_modp(temp , p)
            ws = (w,-w)
            assert pow(w,2,p) == temp
            for w in ws:
                a = -w * m3 * inverse(m2 ,p ) % p
                f1 = m1
                f2 = m3 * inverse(a , p) % p
                if pow(a , (p-1)//2,p) ==1 and pow(w , (p-1)//2,p) ==1 and pow(f1**2+f2**2 , (p-1)//2,p) ==1:
                    return left , [[1 , 0] , [0 , a]] , [[f1 , -f2] , [f2 , f1]] , [[1 , 0] , [0 , w]]
        M = Zi_to_PSLp(mul(g[1] , PSLp_to_Zi(M)))
        left += 1

至此，分解完毕。

def preimages(M):
    left , a, b,c = diagonalize(M)
    print(left)
    flist = [4]*left
    for x in (a,b,c):
        temp = PSLp_to_Zi(x)
        tempm = convert_to_Zi(temp)
        print(tempm)
        flist += factor_Zi(tempm)
        print(flist)
    return flist

实验中，仍然有些bug可能，diagonalize函数有概率永远无解，暂未研究原因，其他的也还不知道是否有更好的方法能够提高成功概率。实验时，大多数的矩阵在30s-40s左右跑出结果，但也有少部分矩阵能在3s内秒出，能够满足题目10s的要求。剩下的就没再试着本地起题目试着打了(懒得写交互)

参考文献

[1] Petit C., Lauter K., Quisquater JJ. (2008) Full Cryptanalysis of LPS and Morgenstern Hash Functions. In: Ostrovsky R., De Prisco R., Visconti I. (eds) Security and Cryptography for Networks. SCN 2008. Lecture Notes in Computer Science, vol 5229. Springer, Berlin, Heidelberg.

[2] G. Zémor and J.-P. Tillich. Collisions for the LPS expander graph hash function. To appear in the proceedings of Advances in Cryptology - EUROCRYPT 2008, 2008.

[3] G. Davidoff, P. Sarnak and A. Valette. Elementary Number Theory, Group Theory, and Ramanujan Graphs. Cambridge University Press, 2003.