DASCTF 2025上半年赛Crypto部分Write Up

Excessive Security

from hashlib import sha256  
from random import randint  
from ecdsa import SECP256k1  
from Crypto.Util.number import inverse, bytes_to_long, getPrime  
from secret import flag  
  
curve = SECP256k1  
G, n = curve.generator, curve.order  
  
def H(msg): return int.from_bytes(sha256(msg.encode()).digest(), 'big')  
  
def sign(h, priv, k):  
    r = (k * G).x() % n  
    s = (inverse(k, n) * (h + r * priv)) % n  
    return s, r  
  
def split(data):  
    s = data.decode() if isinstance(data, bytes) else data  
    return [s[i::4] for i in range(4)]  
  
def sign_encrypt(plaintext):  
    blocks = split(plaintext)  
    h = list(map(H, blocks))  
  
    x1 = randint(1, n - 1)  
    x2 = randint(1, n - 1)  
      
    k1 = randint(n // 8, n - 1)  
    k2 = randint(n // 8, n - 1)  
  
    s1, r1 = sign(h[0], x1, k1)  
    s2, _  = sign(h[1], x2, k1)  
    s3, r2 = sign(h[2], x1, k2)  
    s4, _  = sign(h[3], x2, k2)  
  
    m = bytes_to_long(plaintext)  
    p, q = getPrime(512), getPrime(512)  
    N = p * q  
    e = 65537  
    c1 = pow(m, e, N)  
    c2 = pow(x1 * m + x2, e, N)  
  
    print(f"N = {N}")  
    print(f"c1 = {c1}")  
    print(f"c2 = {c2}")  
    print(f"(h1, s1, r1) = ({h[0]}, {s1}, {r1})")  
    print(f"(h2, s2, r1) = ({h[1]}, {s2}, {r1})")  
    print(f"(h3, s3, r2) = ({h[2]}, {s3}, {r2})")  
    print(f"(h4, s4, r2) = ({h[3]}, {s4}, {r2})")  
  
if __name__ == '__main__':  
    sign_encrypt(flag)  
  
    # N = ...  
    # c1 = ...    
    # c2 = ...    
    # (h1, s1, r1) = (..., ..., ...)    
    # (h2, s2, r1) = (..., ..., ...)    
    # (h3, s3, r2) = (..., ..., ...)    
    # (h4, s4, r2) = (..., ..., ...)

对签名步骤，有（四式从上到下分别编号为1、2、3、4）：

$$\begin{cases} s_1\equiv k_1^{-1}(h_1+r_1x_1)\pmod{n}\\ s_2\equiv k_1^{-1}(h_2+r_1x_2)\pmod{n}\\ s_3\equiv k_2^{-1}(h_3+r_2x_1)\pmod{n}\\ s_4\equiv k_2^{-1}(h_4+r_2x_2)\pmod{n} \end{cases}$$

可以得到：

$$\begin{cases} k_1s_1\equiv h_1+r_1x_1\pmod{n}\\ k_1s_2\equiv h_2+r_1x_2\pmod{n}\\ k_2s_3\equiv h_3+r_2x_1\pmod{n}\\ k_2s_4\equiv h_4+r_2x_2\pmod{n}\\ \end{cases}$$

1、2式两边乘上$r’$，3、4式两边乘上$r$，有：

$$\begin{cases} r_2k_1s_1\equiv r_2h_1+r_1r_2x_1\pmod{n}\\ r_2k_1s_2\equiv r_2h_2+r_1r_2x_2\pmod{n}\\ r_1k_2s_3\equiv r_1h_3+r_1r_2x_1\pmod{n}\\ r_1k_2s_4\equiv r_1h_4+r_1r_2x_2\pmod{n}\\ \end{cases}$$

1、3式相减，2、4式相减可以消去$rr’x_1,rr’x_2$：

$$\begin{cases} r_2k_1s_1-r_1k_2s_3\equiv r_2h_1-r_1h_3\pmod{n}\\ r_2k_1s_2-r_1k_2s_4\equiv r_2h_2-r_1h_4\pmod{n}\\ \end{cases}$$

1式乘上$s_4$，2式乘上$s_3$后相减：

$$r_2k_1(s_1s_4-s_2s_3)\equiv (r_2h_1-r_1h_3)s_4-(r_2h_2-r_1h_4)s_3\pmod{n}$$

由于$r_1,r_2,s_1,s_2,s_3,s_4,h_1,h_2,h_3,h_4$均已知，那么可以求出$k_1$：

$$k_1\equiv \frac{(r_2h_1-r_1h_3)s_4-(r_2h_2-r_1h_4)s_3}{r_2(s_1s_4-s_2s_3)}$$

从而可以计算出：

$$\begin{cases} x_1\equiv (k_1s_1 - h_1)r_1^{-1}\pmod{n}\\ x_2\equiv (k_1s_2 - h_2)r_1^{-1}\pmod{n}\\ \end{cases}$$

那么对于加密部分我们就可以得到两个方程：

$$\begin{aligned} &f\equiv m^{e}-c_1\equiv0\pmod{N}\\ &g\equiv (x_1m+x_2)^{e}-c_2\equiv 0\pmod{N} \end{aligned}$$

可以看到$f,g$有一个共同的根$m$，那么只需要求出$f,g$的最大公约式再求根就可以得到$m$，从而得到flag了：

from Crypto.Util.number import *
from ecdsa import SECP256k1

n = SECP256k1.order
N = ...
c1 = ...
c2 = ...
(h1, s1, r1) = (..., ..., ...)
(h2, s2, r1) = (..., ..., ...)
(h3, s3, r2) = (..., ..., ...)
(h4, s4, r2) = (..., ..., ...)
e = 65537

k1 = (s4 * (r2*h1 - r1*h3) - s3 * (r2*h2 - r1*h4)) * inverse(r2 * (s1*s4 - s2*s3), n) % n

x1 = (k1 * s1 - h1) * inverse(r1, n) % n
x2 = (k1 * s2 - h2) * inverse(r1, n) % n

R.<x> = Zmod(N)[]
f = x^e - c1
g = (x1 * x + x2)^e - c2

G = R(f._pari_with_name('x').gcd(g._pari_with_name('x')))

m = -list(G.monic())[0]
print(long_to_bytes(int(m)))

我在比赛前没有存HGCD的板子，在搜板子的时候惊喜地发现原来sage是提供了快速计算两多项式的最大公约式的方法的，这样以后求两多项式最大公约式的时候就不用这么麻烦了.

补充说明

实际上在求$x_1,x_2$一步还有另外一种更直接的解法，已知有关系：

$$\begin{cases} k_1s_1\equiv h_1+r_1x_1\pmod{n}\\ k_1s_2\equiv h_2+r_1x_2\pmod{n}\\ k_2s_3\equiv h_3+r_2x_1\pmod{n}\\ k_2s_4\equiv h_4+r_2x_2\pmod{n}\\ \end{cases}$$

移项可以得到：

$$\begin{cases} s_1k_1&&-r_1x_1&&\equiv h_1\pmod{n}\\ s_2k_1&&&-r_1x_2&\equiv h_2\pmod{n}\\ &s_3k_2&-r_2x_1&&\equiv h_3\pmod{n}\\ &s_4k_2&&-r_2x_2&\equiv h_4\pmod{n} \end{cases}$$

转换为矩阵形式就是：

$$\left(\begin{matrix} s_1&0&-r_1&0\\ s_2&0&0&-r_1\\ 0&s_3&-r_2&0\\ 0&s_4&0&-r_2 \end{matrix}\right) \left(\begin{matrix} k_1\\k_2\\x_1\\x_2 \end{matrix}\right)\equiv \left(\begin{matrix} h_1\\h_2\\h_3\\h_4 \end{matrix}\right)\pmod{n}$$

直接求解即可：

from Crypto.Util.number import *
from ecdsa import SECP256k1

n = SECP256k1.order
N = ...
c1 = ...
c2 = ...
(h1, s1, r1) = (..., ..., ...)
(h2, s2, r1) = (..., ..., ...)
(h3, s3, r2) = (..., ..., ...)
(h4, s4, r2) = (..., ..., ...)
e = 65537

A = matrix(Zmod(n), [[s1, 0, -r1, 0], [s2, 0, 0, -r1], [0, s3, -r2, 0], [0, s4, 0, -r2]])
v = vector(Zmod(n), [h1, h2, h3, h4])
k1, k2, x1, x2 = A.solve_right(v)

R.<x> = Zmod(N)[]
f = x^e - c1
g = (ZZ(x1) * x + ZZ(x2))^e - c2

G = R(f._pari_with_name('x').gcd(g._pari_with_name('x')))

m = -list(G.monic())[0]
print(long_to_bytes(int(m)))

Strange RSA

import random  
from secret import flag, hint  
from Crypto.Util.number import getPrime, GCD, bytes_to_long, inverse  
  
def generate_flag_params(keysize=256):  
    while True:  
        p, q = getPrime(keysize), getPrime(keysize)  
        N = p * q  
        A = (p**4 - 1) * (q**4 - 1)  
  
        while True:  
            u = random.randint(2, 1 << 26)  
            v = random.randint(2, 1 << 26)  
            if GCD(u, v) == 1:  
                break  
  
        w = random.randint(-(v * N) + 1, v * N - 1)  
        numerator = w + A * v  
        if numerator % u != 0:  
            continue  
  
        e = numerator // u  
        if e <= 0 or GCD(A, e) != 1:  
            continue  
  
        return {'N': N, 'u': u, 'v': v, 'e': e}  
  
def generate_hint_params(u, v, keysize=512):  
    p, q = getPrime(keysize), getPrime(keysize)  
    N = p*q  
    phi = (p-1)*(q-1)  
    e1, e2 = 0, 0  
    while True:  
        try:  
            x = getPrime(256)  
            d1 = (v+u)*x  
            d2 = (v-u)*x  
            e1 = inverse(d1, phi)  
            e2 = inverse(d2, phi)  
            break  
        except:  
            continue  
  
    return {'e1': e1, 'e2': e2, 'N': N}  
  
  
def encrypt(flag, hint):  
    params = generate_flag_params()  
    N1, u, v, e1 = params['N'], params['u'], params['v'], params['e']  
    e2, e3, N2 = generate_hint_params(u, v)  
  
    flag_ct = pow(bytes_to_long(flag), e1, N1)  
    hint_ct = pow(bytes_to_long(hint), random.choice([e2, e3]), N2)  
  
    return {  
        'flag_ct': flag_ct,  
        'hint_ct': hint_ct,  
        'N1': N1,  
        'N2': N2,  
        'e1': e1,  
        'e2': e2,  
        'e3': e3  
    }  
  
  
if __name__ == '__main__':  
    res = encrypt(flag, hint)  
    print(f'flag ciphertext: {res["flag_ct"]}')  
    print(f'hint ciphertext: {res["hint_ct"]}')  
    print(f'N1: {res["N1"]}')  
    print(f'N2: {res["N2"]}')  
    print(f'e1: {res["e1"]}')  
    print(f'e2: {res["e2"]}')  
    print(f'e3: {res["e3"]}')  
  
    # flag_ct = ...  
    # hint_ct = ...    
    # N1 = ...    
    # N2 = ...    
    # e1 = ...    
    # e2 = ...    
    # e3 = ...

对flag，其生成了两个互质整数$u,v<2^{26}$，并且在$[-vN+1,vN-1)$（$N$在代码中是$N_1$）的范围内生成一个随机整数$w$，得到$e=\frac{v(p^4-1)(q^4-1)+w}{u}$，并使用这个$e$对flag进行加密。而hint就是一般的RSA，但是有$d_1=(v+u)x,d_2=(v-u)x$，可以知道$d_1,d_2< 2^{256+26}=2^{282}\approx N^{0.275}$（实际上在代码中是$N_2$），之后分别对两个私钥计算出对应公钥$e_1,e_2$，并在两个公钥中随机选出一个对hint加密。

首先我们考虑求出hint，因为在对hint加密时对公钥是二选一的，所以两个私钥我们都要求出来，通过上面求出的界，可以考虑使用Boneh-Durfee来求出两个$d_1,d_2$，从而直接解密求出hint：

from Crypto.Util.number import *

"""  
Setting debug to true will display more informations  
about the lattice, the bounds, the vectors...  
"""  
debug = False  
  
"""  
Setting strict to true will stop the algorithm (and  
return (-1, -1)) if we don't have a correct  
upperbound on the determinant. Note that this  
doesn't necesseraly mean that no solutions  
will be found since the theoretical upperbound is  
usualy far away from actual results. That is why  
you should probably use `strict = False`  
"""  
strict = False  
  
"""  
This is experimental, but has provided remarkable results  
so far. It tries to reduce the lattice as much as it can  
while keeping its efficiency. I see no reason not to use  
this option, but if things don't work, you should try  
disabling it  
"""  
helpful_only = True  
dimension_min = 7  # stop removing if lattice reaches that dimension  
  
  
############################################  
# Functions  
##########################################  
  
# display stats on helpful vectors  
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")  
  
# display matrix picture with 0 and X  
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)  
  
# tries to remove unhelpful vectors  
# we start at current = n-1 (last vector)  
def remove_unhelpful(BB, monomials, bound, current):  
    # end of our recursive function  
    if current == -1 or BB.dimensions()[0] <= dimension_min:  
        return BB  
  
    # we start by checking from the end  
    for ii in range(current, -1, -1):  
        # if it is unhelpful:  
        if BB[ii, ii] >= bound:  
            affected_vectors = 0  
            affected_vector_index = 0  
            # let's check if it affects other vectors  
            for jj in range(ii + 1, BB.dimensions()[0]):  
                # if another vector is affected:  
                # we increase the count                if BB[jj, ii] != 0:  
                    affected_vectors += 1  
                    affected_vector_index = jj  
  
            # level:0  
            # if no other vectors end up affected            # we remove it            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  
  
            # level:1  
            # if just one was affected we check            
        	# if it is affecting someone else            elif affected_vectors == 1:  
                affected_deeper = True  
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):  
                    # if it is affecting even one vector  
                    # we give up on this one                    if BB[kk, affected_vector_index] != 0:  
                        affected_deeper = False  
                # remove both it if no other vector was affected and  
                # this helpful vector is not helpful enough                
                # compared to our unhelpful one                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  
    # nothing happened  
    return BB  
  
""" Returns:  
* 0,0   if it fails  
* -1,-1 if `strict=true`, and determinant doesn't bound  
* x0,y0 the solutions of `pol`  
"""  
  
def boneh_durfee(pol, modulus, mm, tt, XX, YY):  
    """  
    Boneh and Durfee revisited by Herrmann and May  
    finds a solution if:    
    * d < N^delta    
    * |x| < e^delta    
    * |y| < e^0.5    
    whenever delta < 1 - sqrt(2)/2 ~ 0.292    
    """  
    # substitution (Herrman and May)  
    PR.<u,x,y> = PolynomialRing(ZZ)  
    Q = PR.quotient(x * y + 1 - u)  # u = xy + 1  
    polZ = Q(pol).lift()  
  
    UU = XX * YY + 1  
  
    # x-shifts  
    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()  
  
    # x-shifts list of monomials  
    monomials = []  
    for polynomial in gg:  
        for monomial in polynomial.monomials():  
            if monomial not in monomials:  
                monomials.append(monomial)  
    monomials.sort()  
  
    # y-shifts (selected by Herrman and May)  
    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)  # substitution  
  
    # y-shifts list of monomials    for jj in range(1, tt + 1):  
        for kk in range(floor(mm / tt) * jj, mm + 1):  
            monomials.append(u ^ kk * y ^ jj)  
  
    # construct lattice B  
    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)  
  
    # Prototype to reduce the lattice  
    if helpful_only:  
        # automatically remove  
        BB = remove_unhelpful(BB, monomials, modulus ^ mm, nn - 1)  
        # reset dimension  
        nn = BB.dimensions()[0]  
        if nn == 0:  
            print("failure")  
            return 0, 0  
  
    # check if vectors are helpful  
    if debug:  
        helpful_vectors(BB, modulus ^ mm)  
  
    # check if determinant is correctly bounded  
    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)")  
  
    # display the lattice basis  
    if debug:  
        matrix_overview(BB, modulus ^ mm)  
  
    # LLL  
    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!")  
  
    # transform vector i & j -> polynomials 1 & 2  
    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):  
            # for i and j, create the two polynomials  
            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)  
  
            # resultant  
            PR.<q> = PolynomialRing(ZZ)  
            rr = pol1.resultant(pol2)  
  
            # are these good polynomials?  
            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)  
  
    # solutions  
    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 sol(N, e, delta = .271):  
    m = 8  # size of the lattice (bigger the better/slower)  
    t = int((1 - 2 * delta) * m)  # optimization from Herrmann and May  
    X = 2 * floor(N ^ delta)  # this _might_ be too much  
    Y = floor(N ^ (1 / 2))  # correct if p, q are ~ same size  
    P.<x,y> = PolynomialRing(ZZ)  
    A = int((N + 1) / 2)  
    pol = 1 + x * (A + y)  
  
    solx, soly = boneh_durfee(pol, e, m, t, X, Y)  
  
    d = int(pol(solx, soly) / e)  
    return d

d1 = sol(N2, e1, 0.276)
d2 = sol(N2, e2, 0.276)

hint = pow(hint_ct, d2, N2)
print(long_to_bytes(hint).decode())

可以求出hint：Do you know the contributions of Cortan and Teşeleanu to RSA?，直接搜索Cortan and Teşeleanu可以搜索到一篇文章

链接：A New Generalized Attack on RSA-like Cryptosystems

其中介绍了一种通过关系：

$$eu-(p^4-1)(q^4-1)v=w$$

从$N=pq$中分解出$p,q$的方法，其中：

$$uv<\frac{2N^4-49N^2+2}{4N+170N^2},|w|<vN$$

对上述关系进行变形显然就是：

$$e=\frac{v(p^4-1)(q^4-1)+w}{u}$$

正是本题的情况，那么我们可以通过论文中给出的方式分解出$p,q$：

在这里$u,v$我们可以通过hint那一步的$d_1=(v+u)x,d_2=(v-u)x$来构造方程求解：

$$\begin{cases} x=(d_1,d_2)\\ v+u=\frac{d_1}{x}\\ v-u=\frac{d_2}{x} \end{cases}$$

得到$u,v$之后就可以通过如下代码直接分解出$p,q$，从而得到flag了：

p_and_q = floor(sqrt(abs(2 * N1 + sqrt((N1^2 + 1)^2 - (e3 * u) / v))))
p_sub_q = floor(sqrt(abs(-2 * N1 + sqrt((N1^2 + 1)^2 - (e3 * u) / v))))

p0 = (p_and_q + p_sub_q) // 2

R.<x> = Zmod(N1)[]
f = x + p0
f = f.monic()
pad = f.small_roots(X = 2^16, beta = 0.4)

p = ZZ(p0 + pad[0])
q = N1 // p

phi = (p - 1) * (q - 1)
d = inverse_mod(e3, phi)

m = pow(flag_ct, d, N1)
print(long_to_bytes(m).decode())

整道题的完整求解代码如下：

from Crypto.Util.number import *  
  
flag_ct = ...  
hint_ct = ...  
N1 = ...  
N2 = ...  
e1 = ...  
e2 = ...  
e3 = ...  
  
"""  
Setting debug to true will display more informations  
about the lattice, the bounds, the vectors...  
"""  
debug = False  
  
"""  
Setting strict to true will stop the algorithm (and  
return (-1, -1)) if we don't have a correct  
upperbound on the determinant. Note that this  
doesn't necesseraly mean that no solutions  
will be found since the theoretical upperbound is  
usualy far away from actual results. That is why  
you should probably use `strict = False`  
"""  
strict = False  
  
"""  
This is experimental, but has provided remarkable results  
so far. It tries to reduce the lattice as much as it can  
while keeping its efficiency. I see no reason not to use  
this option, but if things don't work, you should try  
disabling it  
"""  
helpful_only = True  
dimension_min = 7  # stop removing if lattice reaches that dimension  
  
  
############################################  
# Functions  
##########################################  
  
# display stats on helpful vectors  
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")  
  
# display matrix picture with 0 and X  
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)  
  
# tries to remove unhelpful vectors  
# we start at current = n-1 (last vector)  
def remove_unhelpful(BB, monomials, bound, current):  
    # end of our recursive function  
    if current == -1 or BB.dimensions()[0] <= dimension_min:  
        return BB  
  
    # we start by checking from the end  
    for ii in range(current, -1, -1):  
        # if it is unhelpful:  
        if BB[ii, ii] >= bound:  
            affected_vectors = 0  
            affected_vector_index = 0  
            # let's check if it affects other vectors  
            for jj in range(ii + 1, BB.dimensions()[0]):  
                # if another vector is affected:  
                # we increase the count                if BB[jj, ii] != 0:  
                    affected_vectors += 1  
                    affected_vector_index = jj  
  
            # level:0  
            # if no other vectors end up affected            
        # we remove it            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  
  
            # level:1  
            # if just one was affected we check            
        # if it is affecting someone else            elif affected_vectors == 1:  
                affected_deeper = True  
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):  
                    # if it is affecting even one vector  
                    # we give up on this one                    if BB[kk, affected_vector_index] != 0:  
                        affected_deeper = False  
                # remove both it if no other vector was affected and  
                # this helpful vector is not helpful enough                
                # compared to our unhelpful one                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  
    # nothing happened  
    return BB  
  
""" Returns:  
* 0,0   if it fails  
* -1,-1 if `strict=true`, and determinant doesn't bound  
* x0,y0 the solutions of `pol`  
"""  
  
def boneh_durfee(pol, modulus, mm, tt, XX, YY):  
    """  
    Boneh and Durfee revisited by Herrmann and May  
    finds a solution if:    
    * d < N^delta    
    * |x| < e^delta    
    * |y| < e^0.5    whenever delta < 1 - sqrt(2)/2 ~ 0.292    
    """  
    # substitution (Herrman and May)  
    PR.<u,x,y> = PolynomialRing(ZZ)  
    Q = PR.quotient(x * y + 1 - u)  # u = xy + 1  
    polZ = Q(pol).lift()  
  
    UU = XX * YY + 1  
  
    # x-shifts  
    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()  
  
    # x-shifts list of monomials  
    monomials = []  
    for polynomial in gg:  
        for monomial in polynomial.monomials():  
            if monomial not in monomials:  
                monomials.append(monomial)  
    monomials.sort()  
  
    # y-shifts (selected by Herrman and May)  
    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)  # substitution  
  
    # y-shifts list of monomials    for jj in range(1, tt + 1):  
        for kk in range(floor(mm / tt) * jj, mm + 1):  
            monomials.append(u ^ kk * y ^ jj)  
  
    # construct lattice B  
    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)  
  
    # Prototype to reduce the lattice  
    if helpful_only:  
        # automatically remove  
        BB = remove_unhelpful(BB, monomials, modulus ^ mm, nn - 1)  
        # reset dimension  
        nn = BB.dimensions()[0]  
        if nn == 0:  
            print("failure")  
            return 0, 0  
  
    # check if vectors are helpful  
    if debug:  
        helpful_vectors(BB, modulus ^ mm)  
  
    # check if determinant is correctly bounded  
    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)")  
  
    # display the lattice basis  
    if debug:  
        matrix_overview(BB, modulus ^ mm)  
  
    # LLL  
    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!")  
  
    # transform vector i & j -> polynomials 1 & 2  
    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):  
            # for i and j, create the two polynomials  
            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)  
  
            # resultant  
            PR.<q> = PolynomialRing(ZZ)  
            rr = pol1.resultant(pol2)  
  
            # are these good polynomials?  
            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)  
  
    # solutions  
    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 sol(N, e, delta = .271):  
    m = 8  # size of the lattice (bigger the better/slower)  
    t = int((1 - 2 * delta) * m)  # optimization from Herrmann and May  
    X = 2 * floor(N ^ delta)  # this _might_ be too much  
    Y = floor(N ^ (1 / 2))  # correct if p, q are ~ same size  
    P.<x,y> = PolynomialRing(ZZ)  
    A = int((N + 1) / 2)  
    pol = 1 + x * (A + y)  
  
    solx, soly = boneh_durfee(pol, e, m, t, X, Y)  
  
    d = int(pol(solx, soly) / e)  
    return d  
  
d1 = sol(N2, e1, 0.276)  
d2 = sol(N2, e2, 0.276)  
  
x = gcd(d1, d2)  
  
v_and_u = d1 // x  
v_sub_u = d2 // x  
v, u = (v_and_u + v_sub_u) // 2, (v_and_u - v_sub_u) // 2  
  
hint = pow(hint_ct, d2, N2)  
print(long_to_bytes(hint).decode())  
  
p_and_q = floor(sqrt(abs(2 * N1 + sqrt((N1^2 + 1)^2 - (e3 * u) / v))))  
p_sub_q = floor(sqrt(abs(-2 * N1 + sqrt((N1^2 + 1)^2 - (e3 * u) / v))))  
  
p0 = (p_and_q + p_sub_q) // 2  
  
R.<x> = Zmod(N1)[]  
f = x + p0  
f = f.monic()  
pad = f.small_roots(X = 2^16, beta = 0.4)  
  
p = ZZ(p0 + pad[0])  
q = N1 // p  
  
phi = (p - 1) * (q - 1)  
d = inverse_mod(e3, phi)  
  
m = pow(flag_ct, d, N1)  
print(long_to_bytes(m).decode())

这题的$e_1,e_2,e_3$的顺序在题目代码中实际上应该为$e_3,e_1,e_2$