果然。。。都有解了呀（）

还得是Nowayback和HashTeam（

Quicker_20240602_181752

EZ_RSA#

题目需要我们算 $p1 + p2 + q1 + q2$ 的和的md5值，所以我们需要算出这四个素数才行。

于是这里分成两部分去计算：

gen1#

在加密代码里，出题人直接把这函数写成一行，这人真的坏（）

所以我们需要把这个gen1转成正常的加密代码才好做题：

1
rand = lambda n: bytes(random.getrandbits(1) for _ in range(n))
2

3
def gen1(bits):
4
    p = getPrime(bits // 2, randfunc=rand)
5
    q = getPrime(bits // 2, randfunc=rand)
6
    return p, q

因为x使用的是 getrandbits(1)；因此，x里会有b’\x00’ 和 b’\x01’ 这两种可能的bytes的组合。

而且如果自己尝试用gen1生成素数会发现：对应的bytes类型值里都是b”\x80”开头，且除去开头后的长度都是63。

故我们可以通过剪枝的方法去爆破出 $p_1$ 、 $q_1$ ：

点击展开代码

1
from Crypto.Util.number import bytes_to_long
2

3

4
def findflag(p, q, n):
5
    if len(p) == 63:
6
        pp = bytes_to_long(b"\x80" + p)
7
        if n % pp == 0:
8
            print(pp)
9
    else:
10
        L = len(p)
11
        pp = bytes_to_long(p)
12
        qq = bytes_to_long(q)
13
        if pp * qq % (256**L) == n % (256**L):
14
            findflag(b"\x00" + p, b"\x01" + q, n)
15
            findflag(b"\x01" + p, b"\x01" + q, n)
16
            findflag(b"\x00" + p, b"\x00" + q, n)
17
            findflag(b"\x01" + p, b"\x00" + q, n)
18

19

20
n1 = 44945076854246685060397710825960160082061127479194994041436997195972585701097443198954359213635892234058786065342178389181538153413878118039445271277476379366294977408981257175008890470376094381644530106799352839565803317977637572325347776636285703517680754624094985374606187797141657688145287340444623176193
21
# 因为n1的最后一位是3，所以p和q的bytes里的最后一位都是b"\x01"
22
findflag(b"\x01", b"\x01", n1)
23
print("over")
24
# 6704108555018235126044943757232820606509394092422982568570889193755927961059256373509251540968761510597955349467469602569838971815245810437719445704016129
25
# 6704109351064115455893295955648034742490075153469647652626278891915863408109665863057399099944508436950445216168956861863970048266975514429634433698038017

gen2#

这个其实是RSA的一个后门：A new idea for RSA backdoors (arxiv.org)，不过不看论文也能推出来的，算是比较巧妙吧（）

对应的解法如下：

获得所有可能的 (q2 mod T, p2 mod T)#

还是一样，先将gen2转成正常的加密代码，这人真的坏（）：

1
def gen2(alpha=512, K=500, T=getPrime(506)):
2
    while True:
3
        q = getPrime(alpha)
4
        r = getPrime(alpha)
5
        for k in range(2, K+1):
6
            p = r + (k * q - r) % T
7
            if isPrime(p):
8
                return p, q, T

现在我们便可以看出：

p_2\equiv kq_2~mod~T

于是代入到 $n_2$ 中便有：

n_2\equiv kq_2^2~mod~T

因为 $k\in[2,K]$ 且 $K$ 很小，那么我们可以遍历 $k$ 的值，

得到所有可能的 $(q_2~mod~T,~p_2~mod~T)$ ，然后通过下面的过程，我们来恢复 $p_2,~q_2$

恢复 p2, q2#

通过gen2以及上述得到的式子，我们可以把 $n_2$ 写成：

n_2=[\pi T+(p_2~mod~T)][vT+(q_2~mod~T)]~~~~~~\left(1\right)

如果我们记： $\delta=[n_2-(p_2~mod~T)(q_2~mod~T)]/T$

那么我们展开便有：

\delta=\pi vT+\pi(q~mod~T)+v(p~mod~T)

然后我们记： $x=\pi,~C=\pi+v,~a=q~mod~T,~b=p~mod~T$ ，

那么我们可以得到：

\delta=x(C-x)T+ax+b(C-x)

我们计算一下 $C$ 的上界，会发现：我们是能穷举 $C$ 的

我们把 $\left(1\right)$ 式张开得到
$n_2=\pi vT^2+vT(p_2~mod~T)+\pi T(q_2~mod~T)+(p_2q_2~mod~T)$
那么明显 $~C=\pi+v\le\pi v\le n_2/T^2$

通过计算（或者直接拿sage去算 $n_2/T^2$ ）可以知道： $n_2/T^2<2^{15}~(small~num)$ ，可以看出 $C$ 值确实很小

于是我们就可以穷举 $C$ ，用下面这个方法来计算 $x$ （即 $\pi$ ），进而得到 $v$ ：

简单来说，就是两步：

1，遍历去找可能的C，使得 $\Delta$ 这个判别式不小于0

2，然后尝试求我们这个二次方程的两个解： $D=(CT+a-b\pm\sqrt{\Delta})/2T$

3，然后判断 $D$ 是否为整数；如果有一个解 $D_i$ 成立，则 $\pi=D_i,~v=C-\pi$

然后再通过下式便可算出 $p_2,~q_2$ ：

p_2=\pi T+b

q_2=vT+a

由于我们有了 $a=q_2~mod~T,~b=p_2~mod~T$ 的所有可能值，故我们可以通过遍历 $(a,~b)$ 进行上述操作，恢复出 $p_2,~q_2$

对应的exp:

点击展开代码

1
import gmpy2
2
from Crypto.Util.number import *
3

4
def Legendre(n, p):
5
    return pow(n, (p - 1) // 2, p)
6

7

8
def Tonelli_Shanks(n, p):
9
    assert Legendre(n, p) == 1
10
    if p % 4 == 3:
11
        return pow(n, (p + 1) // 4, p)
12
    q = p - 1
13
    s = 0
14
    while q % 2 == 0:
15
        q = q // 2
16
        s += 1
17
    for z in range(2, p):
18
        if Legendre(z, p) == p - 1:
19
            c = pow(z, q, p)
20
            break
21
    r = pow(n, (q + 1) // 2, p)
22
    t = pow(n, q, p)
23
    m = s
24
    if t % p == 1:
25
        return r
26
    else:
27
        i = 0
28
        while t % p != 1:
29
            temp = pow(t, 2 ** (i + 1), p)
30
            i += 1
31
            if temp % p == 1:
32
                b = pow(c, 2 ** (m - i - 1), p)
33
                r = r * b % p
34
                c = b * b % p
35
                t = t * c % p
36
                m = i
37
                i = 0
38
        return r
39

40
K=500
41
n = 57784854392324291351358704449756491526369373648574288191576366413179694041729248864428194536249209588548791706980878177790271653262097539281383559433402738548851606776347237650302071287124974607439996041713554182180186588308614458904981542909792071322939678815174962366963098166320441995961513884899917498099
42
T = 150514823288951667574011681197229106951781617714873679347685702558528178681176081082658953342482323349796111911103531429615442550000291753989779754337491
43

44
nt = n % T
45
gammas = []
46
for k in range(2, K + 1):
47
    k_ = inverse(k, T)
48
    if Legendre(nt * k_, T) == 1:
49
        gammas.append(nt * k_ % T)
50
pqs = []
51
for gamma in gammas:
52
    qt1 = Tonelli_Shanks(gamma, T)
53
    assert qt1**2 % T == gamma and qt1 < T
54
    qt2 = T - qt1
55
    assert qt2**2 % T == gamma
56
    pt1 = nt * inverse(qt1, T) % T
57
    pt2 = nt * inverse(qt2, T) % T
58
    pqs.append((qt1, pt1))
59
    pqs.append((qt2, pt2))
60
for a, b in pqs:
61
    begin = gmpy2.iroot(2 * (n // T**2) - 1, 2)[0]
62
    end = n // (T**2)
63
    delta = (n - a * b) // T
64
    for C in range(begin, end + 1):
65
        Delta = (b - a - C * T) ** 2 - 4 * T * (delta - b * C)
66
        if Delta < 0:
67
            continue
68
        xx = gmpy2.iroot(Delta, 2)
69
        if xx[1]:
70
            x1 = C * T + a - b + xx[0]
71
            x2 = C * T + a - b - xx[0]
72
            for x in (x1, x2):
73
                if x % (2 * T) == 0:
74
                    x = x // (2 * T)
75
                    pi = x
76
                    v = C - x
77
                    p1 = pi * T + b
78
                    q1 = v * T + a
79
                    if p1 * q1 == n:
80
                        print(p1, q1)
81
                        exit(0)
82
# 8604143985568971357221106163518321547782942525630490158067993880524661927741225574307260111628133976467492901704516592869940382055272648214920231756723373
83
# 6715932984064668444342570644774156271984002289395510283696469320418962556390690901906940107908954287013902962625382543926197508086756331581472941654687263

last#

最后计算flag即可：

1
import hashlib
2

3
p1 = 6704108555018235126044943757232820606509394092422982568570889193755927961059256373509251540968761510597955349467469602569838971815245810437719445704016129
4
q1 = 6704109351064115455893295955648034742490075153469647652626278891915863408109665863057399099944508436950445216168956861863970048266975514429634433698038017
5
p2 = 8604143985568971357221106163518321547782942525630490158067993880524661927741225574307260111628133976467492901704516592869940382055272648214920231756723373
6
q2 = 6715932984064668444342570644774156271984002289395510283696469320418962556390690901906940107908954287013902962625382543926197508086756331581472941654687263
7
flag = "DASCTF{" + hashlib.md5(str(p1 + p2 + q1 + q2).encode()).hexdigest() + "}"
8
print(flag)
9
# DASCTF{354ed97c5a3d9d16f49ad93fc30e1c6f}

CF#

Boneh and Durfee#

对于相对较大的 $e$ ，并且 $e<n^{0.292}$ ，我们可以利用Boneh and Durfee Attack来做

首先我们可以知道： $e\times d≡1~(mod~\phi(n))$

于是可以转化一下，得到 $e\times d = 1 + k\times\phi(n)$

即 $k\times\phi(n)+1≡0\mod e~——~(1)$

又因为 $\phi(n)=(p-1)(q-1)(r-1)(s-1)=n-t+1$ ，其中 $t$ 为一个包含 $p,q,r,s$ 的多项式

故我们代入到式 $(1)$ 中，得到 $k(n-t+1)+1≡0\mod e$

然后我们记： $A=(n+1)$ ，则上式可以写成：

f(x,y)=x(A-y)+1≡0\mod e

假如我们可以解出这个方程的根，我们便可得到 $\phi(n)$

多元n与phi分解因子#

多元 $n$ 、 $\phi(n)$ 分解因子这一步，需要用到这篇论文：On Some Attacks on Multi-prime RSA；具体位置在 “Equivalence of Factoring and Exposing the Private Key” 这节。

参考其中的方法，使用 $n$ 和 $\phi(n)$ 来分解出因子即可。

不过嘛。。。网上其实是有脚本的：known_phi.py，直接拿来用就行。

完整EXP：

点击展开代码

1
from Crypto.Util.number import *
2
from Crypto.Cipher import AES
3
import hashlib
4
from __future__ import print_function
5

6

7
############################################
8
# Config
9
##########################################
10

11
"""
12
Setting debug to true will display more informations
13
about the lattice, the bounds, the vectors...
14
"""
15
debug = True
16

17
"""
18
Setting strict to true will stop the algorithm (and
19
return (-1, -1)) if we don't have a correct
20
upperbound on the determinant. Note that this
21
doesn't necesseraly mean that no solutions
22
will be found since the theoretical upperbound is
23
usualy far away from actual results. That is why
24
you should probably use `strict = False`
25
"""
26
strict = False
27

28
"""
29
This is experimental, but has provided remarkable results
30
so far. It tries to reduce the lattice as much as it can
31
while keeping its efficiency. I see no reason not to use
32
this option, but if things don't work, you should try
33
disabling it
34
"""
35
helpful_only = True
36
dimension_min = 7 # stop removing if lattice reaches that dimension
37

38
############################################
39
# Functions
40
##########################################
41

42
# display stats on helpful vectors
43
def helpful_vectors(BB, modulus):
44
    nothelpful = 0
45
    for ii in range(BB.dimensions()[0]):
46
        if BB[ii,ii] >= modulus:
47
            nothelpful += 1
48

49
    print(nothelpful, "/", BB.dimensions()[0], " vectors are not helpful")
50

51
# display matrix picture with 0 and X
52
def matrix_overview(BB, bound):
53
    for ii in range(BB.dimensions()[0]):
54
        a = ('%02d ' % ii)
55
        for jj in range(BB.dimensions()[1]):
56
            a += '0' if BB[ii,jj] == 0 else 'X'
57
            if BB.dimensions()[0] < 60:
58
                a += ' '
59
        if BB[ii, ii] >= bound:
60
            a += '~'
61
        print(a)
62

63
# tries to remove unhelpful vectors
64
# we start at current = n-1 (last vector)
65
def remove_unhelpful(BB, monomials, bound, current):
66
    # end of our recursive function
67
    if current == -1 or BB.dimensions()[0] <= dimension_min:
68
        return BB
69

70
    # we start by checking from the end
71
    for ii in range(current, -1, -1):
72
        # if it is unhelpful:
73
        if BB[ii, ii] >= bound:
74
            affected_vectors = 0
75
            affected_vector_index = 0
76
            # let's check if it affects other vectors
77
            for jj in range(ii + 1, BB.dimensions()[0]):
78
                # if another vector is affected:
79
                # we increase the count
80
                if BB[jj, ii] != 0:
81
                    affected_vectors += 1
82
                    affected_vector_index = jj
83

84
            # level:0
85
            # if no other vectors end up affected
86
            # we remove it
87
            if affected_vectors == 0:
88
                print("* removing unhelpful vector", ii)
89
                BB = BB.delete_columns([ii])
90
                BB = BB.delete_rows([ii])
91
                monomials.pop(ii)
92
                BB = remove_unhelpful(BB, monomials, bound, ii-1)
93
                return BB
94

95
            # level:1
96
            # if just one was affected we check
97
            # if it is affecting someone else
98
            elif affected_vectors == 1:
99
                affected_deeper = True
100
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):
101
                    # if it is affecting even one vector
102
                    # we give up on this one
103
                    if BB[kk, affected_vector_index] != 0:
104
                        affected_deeper = False
105
                # remove both it if no other vector was affected and
106
                # this helpful vector is not helpful enough
107
                # compared to our unhelpful one
108
                if affected_deeper and abs(bound - BB[affected_vector_index, affected_vector_index]) < abs(bound - BB[ii, ii]):
109
                    print("* removing unhelpful vectors", ii, "and", affected_vector_index)
110
                    BB = BB.delete_columns([affected_vector_index, ii])
111
                    BB = BB.delete_rows([affected_vector_index, ii])
112
                    monomials.pop(affected_vector_index)
113
                    monomials.pop(ii)
114
                    BB = remove_unhelpful(BB, monomials, bound, ii-1)
115
                    return BB
116
    # nothing happened
117
    return BB
118

119
"""
120
Returns:
121
* 0,0   if it fails
122
* -1,-1 if `strict=true`, and determinant doesn't bound
123
* x0,y0 the solutions of `pol`
124
"""
125
def boneh_durfee(pol, modulus, mm, tt, XX, YY):
126
    """
127
    Boneh and Durfee revisited by Herrmann and May
128

129
    finds a solution if:
130
    * d < N^delta
131
    * |x| < e^delta
132
    * |y| < e^0.5
133
    whenever delta < 1 - sqrt(2)/2 ~ 0.292
134
    """
135

136
    # substitution (Herrman and May)
137
    PR.<u, x, y> = PolynomialRing(ZZ)
138
    Q = PR.quotient(x*y + 1 - u) # u = xy + 1
139
    polZ = Q(pol).lift()
140

141
    UU = XX*YY + 1
142

143
    # x-shifts
144
    gg = []
145
    for kk in range(mm + 1):
146
        for ii in range(mm - kk + 1):
147
            xshift = x^ii * modulus^(mm - kk) * polZ(u, x, y)^kk
148
            gg.append(xshift)
149
    gg.sort()
150

151
    # x-shifts list of monomials
152
    monomials = []
153
    for polynomial in gg:
154
        for monomial in polynomial.monomials():
155
            if monomial not in monomials:
156
                monomials.append(monomial)
157
    monomials.sort()
158

159
    # y-shifts (selected by Herrman and May)
160
    for jj in range(1, tt + 1):
161
        for kk in range(floor(mm/tt) * jj, mm + 1):
162
            yshift = y^jj * polZ(u, x, y)^kk * modulus^(mm - kk)
163
            yshift = Q(yshift).lift()
164
            gg.append(yshift) # substitution
165

166
    # y-shifts list of monomials
167
    for jj in range(1, tt + 1):
168
        for kk in range(floor(mm/tt) * jj, mm + 1):
169
            monomials.append(u^kk * y^jj)
170

171
    # construct lattice B
172
    nn = len(monomials)
173
    BB = Matrix(ZZ, nn)
174
    for ii in range(nn):
175
        BB[ii, 0] = gg[ii](0, 0, 0)
176
        for jj in range(1, ii + 1):
177
            if monomials[jj] in gg[ii].monomials():
178
                BB[ii, jj] = gg[ii].monomial_coefficient(monomials[jj]) * monomials[jj](UU,XX,YY)
179

180
    # Prototype to reduce the lattice
181
    if helpful_only:
182
        # automatically remove
183
        BB = remove_unhelpful(BB, monomials, modulus^mm, nn-1)
184
        # reset dimension
185
        nn = BB.dimensions()[0]
186
        if nn == 0:
187
            print("failure")
188
            return 0,0
189

190
    # check if vectors are helpful
191
    if debug:
192
        helpful_vectors(BB, modulus^mm)
193

194
    # check if determinant is correctly bounded
195
    det = BB.det()
196
    bound = modulus^(mm*nn)
197
    if det >= bound:
198
        print("We do not have det < bound. Solutions might not be found.")
199
        print("Try with highers m and t.")
200
        if debug:
201
            diff = (log(det) - log(bound)) / log(2)
202
            print("size det(L) - size e^(m*n) = ", floor(diff))
203
        if strict:
204
            return -1, -1
205
    else:
206
        print("det(L) < e^(m*n) (good! If a solution exists < N^delta, it will be found)")
207

208
    # display the lattice basis
209
    if debug:
210
        matrix_overview(BB, modulus^mm)
211

212
    # LLL
213
    if debug:
214
        print("optimizing basis of the lattice via LLL, this can take a long time")
215

216
    BB = BB.LLL()
217

218
    if debug:
219
        print("LLL is done!")
220

221
    # transform vector i & j -> polynomials 1 & 2
222
    if debug:
223
        print("looking for independent vectors in the lattice")
224
    found_polynomials = False
225

226
    for pol1_idx in range(nn - 1):
227
        for pol2_idx in range(pol1_idx + 1, nn):
228
            # for i and j, create the two polynomials
229
            PR.<w,z> = PolynomialRing(ZZ)
230
            pol1 = pol2 = 0
231
            for jj in range(nn):
232
                pol1 += monomials[jj](w*z+1,w,z) * BB[pol1_idx, jj] / monomials[jj](UU,XX,YY)
233
                pol2 += monomials[jj](w*z+1,w,z) * BB[pol2_idx, jj] / monomials[jj](UU,XX,YY)
234

235
            # resultant
236
            PR.<q> = PolynomialRing(ZZ)
237
            rr = pol1.resultant(pol2)
238

239
            # are these good polynomials?
240
            if rr.is_zero() or rr.monomials() == [1]:
241
                continue
242
            else:
243
                print("found them, using vectors", pol1_idx, "and", pol2_idx)
244
                found_polynomials = True
245
                break
246
        if found_polynomials:
247
            break
248

249
    if not found_polynomials:
250
        print("no independant vectors could be found. This should very rarely happen...")
251
        return 0, 0
252

253
    rr = rr(q, q)
254

255
    # solutions
256
    soly = rr.roots()
257

258
    if len(soly) == 0:
259
        print("Your prediction (delta) is too small")
260
        return 0, 0
261

262
    soly = soly[0][0]
263
    ss = pol1(q, soly)
264
    solx = ss.roots()[0][0]
265

266
    #
267
    return solx, soly
268

269
def factorize_multi_prime(N, phi):
270
    prime_factors = set()
271
    factors = [N]
272
    while len(factors) > 0:
273
        # Element to factorize.
274
        N = factors[0]
275

276
        w = randrange(2, N - 1)
277
        i = 1
278
        while phi % (2**i) == 0:
279
            sqrt_1 = pow(w, phi // (2**i), N)
280
            if sqrt_1 > 1 and sqrt_1 != N - 1:
281
                # We can remove the element to factorize now, because we have a factorization.
282
                factors = factors[1:]
283

284
                p = gcd(N, sqrt_1 + 1)
285
                q = N // p
286

287
                if is_prime(p):
288
                    prime_factors.add(p)
289
                elif p > 1:
290
                    factors.append(p)
291

292
                if is_prime(q):
293
                    prime_factors.add(q)
294
                elif q > 1:
295
                    factors.append(q)
296

297
                # Continue in the outer loop
298
                break
299

300
            i += 1
301

302
    return list(prime_factors)
303

304

305
def attack(N, e, factor_bit_length, factors, delta=0.25, m=1):
306
    x, y = ZZ["x", "y"].gens()
307
    A = N + 1
308
    f = x * (A + y) + 1
309
    X = int(RR(e) ** delta)
310
    Y = int(2 ** ((factors - 1) * factor_bit_length + 1))
311
    t = int((1 - 2 * delta) * m)
312
    x0, y0 = boneh_durfee(f, e, m, t, X, Y)
313
    z = int(f(x0, y0))
314
    if z % e == 0:
315
        phi = N +int(y0) + 1
316
        factors = factorize_multi_prime(N, phi)
317
        if factors:
318
            return factors
319

320
    return None
321

322

323
n = 12778528771742949806245151753869219326103790041631995252034948773711783128776305944498756929732298934720477166855071150429382343090525399073032692529779161146622028051975895639274962265063528372582516292055195313063685656963925420986244801150981084581230336100629998038062420895185391922920881754851005297105551156140379014123294775868179867798105218243424339964238809811837555910593108364135245826360599234594626605066012137694272914693621191616641820375665250179042481908961611154276842449520816511946371478115661488114557201063593848680402471689545509362224765613961509436533468849519328376263878041094637028661183
324
e = 4446726708272678112679273197419446608921686581114971359716086776036464363243920846432708647591026040092182012898303795518854800856792372040517828716881858432476850992893751986128026419654358442725548028288396111453301336112088168230318117251893266136328216825852616643551255183048159254152784384133765153361821713529774101097531224729203104181285902533238977664673240372553695106609481661124179618839909468411817548602076934523684639875632950838463168454592213740967654900802801128243623511466324869786575827161573559009469945330622017702786149269513046331878690768979142927851424854919322854779975658914469657308779
325
c = b'_\xf7\x16\x00S\x11\xd5\xec\x94+>\x98\x91\x8b\xaeC\xadV3\xf8\x07a\x95\xf6rr\x86\xd4\x1e\x1b\xe7\xf4H\xa0\xd9\x9b\xb5\x05.u\x08\x80\x04\x8d\xee\xec\x98\xf5'
326
p, q, r, s = attack(n, e, 512, 4, 0.127, 10)
327
key = hashlib.md5(str(p + q + r + s).encode()).digest()
328
aes = AES.new(key, mode=AES.MODE_ECB)
329
print(aes.encrypt(c))
330
# DASCTF{d4d0b2c4-b41d-4ce1-871a-b08325900b30}

Pell#

题目给了我们模数N、椭圆曲线点Q和密文C，我们先看加密函数：

1
def encrypt(M, N, e):
2
    xm, ym = M
3
    M = (xm, ym, 0)
4
    a = (1 - xm**3) * inverse(ym**3, N) % N
5
    curve = Pell_Curve(int(a), N)
6
    if curve.is_on_curve(M):
7
        return curve.mul(M, e)
8
    return None

这里是实现了一个Cubic Pell Curve的加密，这一步的解法后边会提及；但在这之前，我们需分解一下N（为后面的解法做准备）。

分解N#

通过ECLCG生成了两个质数作为加密系统的私钥

可以发现ECLCG的生成方式为

Q_n=Q_0+nP

于是我们可以得到对应的两个素数的形式：

p=(Q_0)_x+(n*P)_x

q=(Q_0)_x+[(n_1+1)*n*P]_x

且曲线为 $y^2=x^3$ ，是一个奇异曲线；因此该曲线的dlp计算，我们可以转化成数域上的简单运算 （以下的加法，均为椭圆曲线上的加法）

根据参数方程，我们可以设 $X_1=(t^2,t^3),X_i=((at)^2,(at)^3)$

那么 $X_{i+1}=X_i+X_1$

通过一系列运算化简（当然，上面那步也可以直接用sage去算），可以得到
$X_{i+1}=((\frac{a}{a+1}t)^2,(\frac{a}{a+1}t)^3)$
利用此关系可以简单证明（证明一个数列 $a_{i+1}=a_i/(a_i+1), a_1=1$ 的性质）一下，便可得到下面这个式子
$X_n=(n^{-2}t^2,n^{-3}t^3)$
也即： $nX_1=(n^{-2}t^2,n^{-3}t^3)~~~~~~~~~~~~~\left(2\right)$

对于我们这个随机数发生器，我们构造大量的测试数据可以发现两个质数之间差的n很小，一般不会超过1000；而我们可以由前面的p、q的素数形式知道： $q=p+(n_1*P)_x$ ，且点 $Q_0$ 已知。

（这里为了方便，我们记：参数方程下的 $Q_0=(t^2,~t^3)$ ）

于是我们可以设 $nP=(x^2,~x^3)$ ，然后结合ECLCG的加法，我们就可以得到素数p和q的对应方程

（其中，下式中的 $x_1=(1+n_1)^{-1}x$ ，这个可以根据 $\left(2\right)$ 式和前面的式子推导而得）：

p\equiv[(t^3-x^3)^2(t^2-x^2)^{-2}-x^2-t^2]~mod~N

q\equiv[(t^3-x_1^3)^2(t^2-x_1^2)^{-2}-x_1^2-t^2]~mod~N

然后我们再将上式代入 $N$ ，便能得到一个方程（因为发现太长了，就不列出来了，这里可以参考后面的代码）；然后就会发现：只有一个未知数x，因此我们在爆破 $n_1+1$ 的同时解方程，便可分解出 $p、q$

点击展开代码

1
# recover p, q：
2
import tqdm
3
import gmpy2
4
from Crypto.Util.number import *
5

6
n = 142509889408494696639682201799643202268988370577642546783876593347546850250051841172274152716714403313311584670791108601588046986700175746446804470329761265314268119548997548026516318449862727871202339967955587242463610862701184493904376304507029176806166448249192854001854607465457042204258734279909961546441004233711967226919624405968584449147177981949821415107225952390645278348482729250785152039807053641247569456385545220501027102363800108028762768824577077321340577271010321469215228402821463907345773901277193445125640936231772522681574300491883451795804527966948605710874090658775247402867915876744113646170885038891240778364069379164812880482584571673151293322613478565661348746336931021896668941228934951050789999827329748371987279847108342825214485163497943
7
N = 9909641861967580472493256614158113105414778684219844785944662774988084232380069009372420371597872375863508561123648164278317871844235719752735021659264009
8
Q =  (5725664012637594848838084306454804843458550077896287815106012266176452953193402684379119042639063659980463425502946083139850146060755640351348257807890845,7995259612407104192119579242200802136801092493271952329412936709212369500868134058817979488983954214781719018555338511778896087250394604977285067013758829)
9

10
R.<x> = Zmod(N)[]
11
f1 = x ^ 2 - Q[0]
12
f2 = x ^ 3 - Q[1]
13
for i in f1.roots():
14
    for j in f2.roots():
15
        if i[0] == j[0]:
16
            t = i[0]
17

18
R.<x> = Zmod(N)[]
19
for i in tqdm.tqdm(range(1, 1000)):
20
    xi = x * inverse(i, N)
21
    pp = ((t ^ 3 - x ^ 3) ^ 2 - (t ^ 2 + x ^ 2) * (t ^ 2 - x ^ 2) ^ 2)^3
22
    qq = ((t ^ 3 - xi ^ 3) ^ 2 - (t ^ 2 + xi ^ 2) * (t ^ 2 - xi ^ 2) ^ 2)^2
23
    f = n * (t ^ 2 - x ^ 2) ^ 6 * (t ^ 2 - xi ^ 2) ^ 4 - pp * qq
24
    roots = f.roots()
25
    if roots:
26
        for xx in roots:
27
            x1 = xx[0]
28
            if t != x1:
29
                p1 = (t ^ 3 - x1 ^ 3) ^ 2 * inverse(int(t ^ 2 - x1 ^ 2), N) ^ 2 - t ^ 2 - x1 ^ 2
30
                if n % int(p1) == 0:
31
                    p = int(p1)
32
                    print(p)
33
                    break
34
        if n % int(p1) == 0:
35
            break
36

37

38

39
r = 3
40
s = 2
41
q = int(gmpy2.iroot(n//p^3,2)[0])

Cubic Pell Curve解密+模下开根#

然后就是需要找一个Cubic Pell Curve解密系统

思路出自这篇paper：

A New Public Key Cryptosystem Based on the Cubic Pell Curve

以下是思路：

这里就用到了我们刚刚分解出的p和q；简单来说，这个解密就是以下几步：

1，计算出第一步和第二步中的、模 $p^r$ 和 $q^s$ 下的方程里的两组根： $(a_{p,1},a_{p,2})$ 和 $(a_{q,1},a_{q,2})$ （计算参考论文里的Corollary 3）。

2，然后结合CRT与第三步中所列出来的同余关系，计算出 $a_1、a_2、a_3、a_4$ 。

3，计算 $D=(d_1,d_2,d_3,d_4)$ ，当 $a_i$ 符合某一个 $d_i$ 的条件时，则后面进行椭圆曲线数乘的时候使用这个 $d$ ，其中 $R^3(p)$ 指：在 $Z/pZ$ 下 $p$ 的立方剩余空间（这里直接理解成立方剩余就行）。 $(d_1,d_2,d_3,d_4)$ 的计算公式如下：

4，遍历 $D=(d_1,d_2,d_3,d_4)$ ，并根据第六步的 $M_i=(x_i,y_i,z_i)=d_i⊙(x_C,y_C,z_C)$ 以及 $a_i$ 所对应的Cubic Pell Curve，计算出对应的 $M_i$ ；当我们算出 $M=(x_i,y_i,0)$ 时，我们便恢复出明文 $M$ 了。

该解密系统对应的代码：

点击展开代码

1
# https://eprint.iacr.org/2024/385.pdf
2

3
import gmpy2
4
from Crypto.Util.number import *
5
from sympy.ntheory.modular import crt
6

7

8
def Legendre(n, p):
9
    return pow(n, (p - 1) // 2, p)
10

11

12
def Tonelli_Shanks(n, p):
13
    assert Legendre(n, p) == 1
14
    if p % 4 == 3:
15
        return pow(n, (p + 1) // 4, p)
16
    q = p - 1
17
    s = 0
18
    while q % 2 == 0:
19
        q //= 2
20
        s += 1
21
    z = next(z for z in range(2, p) if Legendre(z, p) == p - 1)
22
    c = pow(z, q, p)
23
    r = pow(n, (q + 1) // 2, p)
24
    t = pow(n, q, p)
25
    m = s
26
    if t % p == 1:
27
        return r
28
    else:
29
        i = 0
30
        while t % p != 1:
31
            temp = pow(t, 2 ** (i + 1), p)
32
            i += 1
33
            if temp % p == 1:
34
                b = pow(c, 2 ** (m - i - 1), p)
35
                r = r * b % p
36
                c = b * b % p
37
                t = t * c % p
38
                m = i
39
                i = 0
40
        return r
41

42

43
class Pell_Curve:
44
    def __init__(self, a, N):
45
        self.a = a
46
        self.N = N
47

48
    def is_on_curve(self, point):
49
        if point is None:
50
            return True
51
        x, y, z = point
52
        return (
53
            x**3 + self.a * y**3 + self.a**2 * z**3 - 3 * self.a * x * y * z
54
        ) % self.N == 1
55

56
    def add(self, P, Q):
57
        x1, y1, z1 = P
58
        x2, y2, z2 = Q
59
        x3 = (x1 * x2 + self.a * (y2 * z1 + y1 * z2)) % self.N
60
        y3 = (x2 * y1 + x1 * y2 + self.a * z1 * z2) % self.N
61
        z3 = (y1 * y2 + x2 * z1 + x1 * z2) % self.N
62
        return (x3, y3, z3)
63

64
    def mul(self, P, x):
65
        Q = (1, 0, 0)
66
        while x > 0:
67
            if x & 1:
68
                Q = self.add(Q, P)
69
            P = self.add(P, P)
70
            x >>= 1
71
        return Q
72

73

74
def psi(p, q, r, s):
75
    psi1 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p**2 + p + 1) * (q**2 + q + 1)
76
    psi2 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p - 1) ** 2 * (q - 1) ** 2
77
    psi3 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p**2 + p + 1) * (q - 1) ** 2
78
    psi4 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p - 1) ** 2 * (q**2 + q + 1)
79
    return (psi1, psi2, psi3, psi4)
80

81

82
def gen(E, Q, r, s):
83
    lcg = LCG(E, Q)
84
    while 1:
85
        p = lcg.get_prime()
86
        q = lcg.get_prime()
87
        if p % 3 == 1 and q % 3 == 1:
88
            N = p**r * q**s
89
            e = 0x20002
90
            return (N, e)
91

92

93
def cubic_residue(a, p):
94
    return pow(a, (p - 1) // gmpy2.gcd(3, p - 1), p) == 1
95

96

97
def judge_d(a, p, q, d):
98
    cr_p = cubic_residue(a, p)
99
    cr_q = cubic_residue(a, q)
100
    if cr_p and cr_q:
101
        return d[1]
102
    elif not cr_p and not cr_q:
103
        return d[0]
104
    elif not cr_p and cr_q:
105
        return d[2]
106
    else:  # cr_p and not cr_q
107
        return d[3]
108

109
# Corollary 3
110
def roots(C, p, n):
111
    xc, yc, zc = C
112
    a = zc**3
113
    b = yc**3 - 3 * xc * yc * zc
114
    c = xc**3 - 1
115
    delta = (b**2 - 4 * a * c) % p
116
    _delta = Tonelli_Shanks(delta, p)
117
    y = (-b + _delta) * gmpy2.invert(2 * a, p) % p
118
    z = (-b - _delta) * gmpy2.invert(2 * a, p) % p
119
    for i in range(1, n):
120
        y = y - (a * y**2 + b * y + c) * gmpy2.invert(
121
            int(2 * a * y + b), p ** (i + 1)
122
        ) % p ** (i + 1)
123
        z = z - (a * z**2 + b * z + c) * gmpy2.invert(
124
            int(2 * a * z + b), p ** (i + 1)
125
        ) % p ** (i + 1)
126
    return (y, z)
127

128
def decrypt(C, p, q, d, r, s):
129
    aps = roots(C, p, r)
130
    aqs = roots(C, q, s)
131
    N = p**r * q**s
132
    A = []
133
    for i in aps:
134
        for j in aqs:
135
            a = crt([p**r, q**s], [int(i), int(j)])[0]
136
            A.append(a)
137
    for i in range(len(A)):
138
        curve = Pell_Curve(A[i], N)
139
        M = curve.mul(C, judge_d(A[i], p, q, d))
140
        if M[2] == 0:
141
            return (int(M[0]), int(M[1]))
142
    return None
143

144

145
ps = psi(p, q, r, s)
146
d = []
147
for i in ps:
148
    d.append(inverse(e, i))
149
print(long_to_bytes(decrypt(C, p, q, d, r, s)[0]))

但代入计算会发现：计算不出来 $d$

回去看加密代码会发现： $e=0x20002$ ，这样的话会有：

\gcd(e,pq(p^2 + p + 1)(q^2 + q + 1)(p - 1)(q - 1))=2

于是我们需要换个方法去转化一下。

这里记 $C$ 为密文， $M$ 为原文， $0x10001$ 为 $l$ ，那么

C=2lM

虽然我们无法直接计算出 $M$ ，但是我们可以知道：

\gcd(e/2,pq(p^2 + p + 1)(q^2 + q + 1)(p - 1)(q - 1))=1

所以我们可以先通过解密系统计算出 $2M$ ，然后再恢复 $M$

而我们去推一下 $2M$ ，可以发现：

2M=(x^2,2xy,y^2)

那么只要对解密完得到的 $M_x$ 开二次方根即可 (我这里用了论文里的Corollary 3)

完整exp:

点击展开代码

1
import tqdm
2
import gmpy2
3
from Crypto.Util.number import *
4
from sympy.ntheory.modular import crt
5

6
n = 142509889408494696639682201799643202268988370577642546783876593347546850250051841172274152716714403313311584670791108601588046986700175746446804470329761265314268119548997548026516318449862727871202339967955587242463610862701184493904376304507029176806166448249192854001854607465457042204258734279909961546441004233711967226919624405968584449147177981949821415107225952390645278348482729250785152039807053641247569456385545220501027102363800108028762768824577077321340577271010321469215228402821463907345773901277193445125640936231772522681574300491883451795804527966948605710874090658775247402867915876744113646170885038891240778364069379164812880482584571673151293322613478565661348746336931021896668941228934951050789999827329748371987279847108342825214485163497943
7
N = 9909641861967580472493256614158113105414778684219844785944662774988084232380069009372420371597872375863508561123648164278317871844235719752735021659264009
8
Q =  (5725664012637594848838084306454804843458550077896287815106012266176452953193402684379119042639063659980463425502946083139850146060755640351348257807890845,7995259612407104192119579242200802136801092493271952329412936709212369500868134058817979488983954214781719018555338511778896087250394604977285067013758829)
9

10
R.<x> = Zmod(N)[]
11
f1 = x ^ 2 - Q[0]
12
f2 = x ^ 3 - Q[1]
13
for i in f1.roots():
14
    for j in f2.roots():
15
        if i[0] == j[0]:
16
            t = i[0]
17

18
R.<x> = Zmod(N)[]
19
for i in tqdm.tqdm(range(1, 1000)):
20
    xi = x * inverse(i, N)
21
    pp = ((t ^ 3 - x ^ 3) ^ 2 - (t ^ 2 + x ^ 2) * (t ^ 2 - x ^ 2) ^ 2)^3
22
    qq = ((t ^ 3 - xi ^ 3) ^ 2 - (t ^ 2 + xi ^ 2) * (t ^ 2 - xi ^ 2) ^ 2)^2
23
    f = n * (t ^ 2 - x ^ 2) ^ 6 * (t ^ 2 - xi ^ 2) ^ 4 - pp * qq
24
    roots = f.roots()
25
    if roots:
26
        for xx in roots:
27
            x1 = xx[0]
28
            if t != x1:
29
                p1 = (t ^ 3 - x1 ^ 3) ^ 2 * inverse(int(t ^ 2 - x1 ^ 2), N) ^ 2 - t ^ 2 - x1 ^ 2
30
                if n % int(p1) == 0:
31
                    p = int(p1)
32
                    print(p)
33
                    break
34
        if n % int(p1) == 0:
35
            break
36

37

38

39
r = 3
40
s = 2
41
q = int(gmpy2.iroot(n//p^3,2)[0])
42

43
def Legendre(n, p):
44
    return pow(n, (p - 1) // 2, p)
45

46

47
def Tonelli_Shanks(n, p):
48
    assert Legendre(n, p) == 1
49
    if p % 4 == 3:
50
        return pow(n, (p + 1) // 4, p)
51
    q = p - 1
52
    s = 0
53
    while q % 2 == 0:
54
        q //= 2
55
        s += 1
56
    z = next(z for z in range(2, p) if Legendre(z, p) == p - 1)
57
    c = pow(z, q, p)
58
    r = pow(n, (q + 1) // 2, p)
59
    t = pow(n, q, p)
60
    m = s
61
    if t % p == 1:
62
        return r
63
    else:
64
        i = 0
65
        while t % p != 1:
66
            temp = pow(t, 2 ** (i + 1), p)
67
            i += 1
68
            if temp % p == 1:
69
                b = pow(c, 2 ** (m - i - 1), p)
70
                r = r * b % p
71
                c = b * b % p
72
                t = t * c % p
73
                m = i
74
                i = 0
75
        return r
76

77

78
class Pell_Curve:
79
    def __init__(self, a, N):
80
        self.a = a
81
        self.N = N
82

83
    def is_on_curve(self, point):
84
        if point is None:
85
            return True
86
        x, y, z = point
87
        return (
88
            x**3 + self.a * y**3 + self.a**2 * z**3 - 3 * self.a * x * y * z
89
        ) % self.N == 1
90

91
    def add(self, P, Q):
92
        x1, y1, z1 = P
93
        x2, y2, z2 = Q
94
        x3 = (x1 * x2 + self.a * (y2 * z1 + y1 * z2)) % self.N
95
        y3 = (x2 * y1 + x1 * y2 + self.a * z1 * z2) % self.N
96
        z3 = (y1 * y2 + x2 * z1 + x1 * z2) % self.N
97
        return (x3, y3, z3)
98

99
    def mul(self, P, x):
100
        Q = (1, 0, 0)
101
        while x > 0:
102
            if x & 1:
103
                Q = self.add(Q, P)
104
            P = self.add(P, P)
105
            x >>= 1
106
        return Q
107

108

109
def psi(p, q, r, s):
110
    psi1 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p**2 + p + 1) * (q**2 + q + 1)
111
    psi2 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p - 1) ** 2 * (q - 1) ** 2
112
    psi3 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p**2 + p + 1) * (q - 1) ** 2
113
    psi4 = p ** (2 * (r - 1)) * q ** (2 * (s - 1)) * (p - 1) ** 2 * (q**2 + q + 1)
114
    return (psi1, psi2, psi3, psi4)
115

116
def cubic_residue(a, p):
117
    return pow(a, (p - 1) // gmpy2.gcd(3, p - 1), p) == 1
118

119

120
def judge_d(a, p, q, d):
121
    cr_p = cubic_residue(a, p)
122
    cr_q = cubic_residue(a, q)
123
    if cr_p and cr_q:
124
        return d[1]
125
    elif not cr_p and not cr_q:
126
        return d[0]
127
    elif not cr_p and cr_q:
128
        return d[2]
129
    else:  # cr_p and not cr_q
130
        return d[3]
131

132
# Corollary 3
133
def roots(C, p, n):
134
    xc, yc, zc = C
135
    a = zc**3
136
    b = yc**3 - 3 * xc * yc * zc
137
    c = xc**3 - 1
138
    delta = (b**2 - 4 * a * c) % p
139
    _delta = Tonelli_Shanks(delta, p)
140
    y = (-b + _delta) * inverse(2 * a, p) % p
141
    z = (-b - _delta) * inverse(2 * a, p) % p
142
    for i in range(1, n):
143
        y = y - (a * y**2 + b * y + c) * inverse(
144
            int(2 * a * y + b), p ** (i + 1)
145
        ) % p ** (i + 1)
146
        z = z - (a * z**2 + b * z + c) * inverse(
147
            int(2 * a * z + b), p ** (i + 1)
148
        ) % p ** (i + 1)
149
    return (y, z)
150

151
C = (81768339111299816705544898152771220210336305743364535623542396932097508874478708007356482559951843443716017684599109593939309497876283954739065532068358640123897297735011312421303760220341679952682608376253590454613919282861879034834442483766217227383792409215337347571227544874051744198403805434968528386779039795337990338248171933970791615195263892724675263032559658819135855374073644306381889879990890042223246077362618291952646985683966244920555989982399613765530011499719074486903003792714562373937144871278164758310693947837335237349195046040995477558132367388842506474592468217861986173383953237474756202802360230890862369060851962186244111055545256271117424905591906972255761770741149563674457745615873496579818814035900990579591845004609499494547080458704584, 84621087300399647293777247835306246465300232341486881635357679809773437325943820311329988605594440622251629971586435278844599108015288735134349648420317858374374591896130432582322507215780484530408523427525797210077752785624079848616300884164345285833494971279538396297733797260240933961493604434803064166573528094704954546014575856837921125063112845773099272164228859908533081610458091806418565502108153124283531626701488036466436102247845200341492584130445948027051529476352653110990934770121255651400555911301783360692285788890607740888376040139286200434818197323063848144168033132174931153362170954175707409126745301216651916596489805505649061280397491087997636237767764403484186207472581036806824115157283392062188592165421921369151939109986184806890233258458794, 107470405748787057257826187107093535161311781207158281438762592876162686482566135109505652982571025667746244660986635749326688338471024529029121466041296205925603803529179856346298760611767192411134153152234712303426575150170977692186997733960581208607060982624871524319162170866870037830416559938924612968969966225954744925337757413696488884826990851697771972617146921133799053964257776476473920346656878321177511107743545375181606366722878715116467369115483252574781605976088763248469134730611983687505906661228606502293949130180171550100994569942435781067167383369188511834406179774120708650048333802855942156250759495072298696263590518886055347952253836780124031369144821306654247715239306949355039924862372681097653701186683219165141054051943634109692683916632353)
152
e = 0x20002
153

154
ps = psi(p, q, r, s)
155
d = []
156
for i in ps:
157
    d.append(inverse(e // 2, i))
158
aps = roots(C, p, r)
159
aqs = roots(C, q, s)
160
A = []
161
for i in aps:
162
    for j in aqs:
163
        a = crt([p**r, q**s], [int(i), int(j)])[0]
164
        A.append(a)
165

166
# 根据Corollary 3改改就能算平方根
167
def roots2(M1, p, n):
168
  a = 1
169
  b = 0
170
  c = -M1
171
  delta = (b**2 - 4 * a * c) % p
172
  _delta = Tonelli_Shanks(delta, p)
173
  y = (-b + _delta) * inverse(2 * a, p) % p
174
  z = (-b - _delta) * inverse(2 * a, p) % p
175
  for i in range(1, n):
176
    y = y - (a * y**2 + b * y + c) * inverse(int(2 * a * y + b), p ** (i + 1)) % p ** (i + 1)
177
    z = z - (a * z**2 + b * z + c) * inverse(int(2 * a * z + b), p ** (i + 1)) % p ** (i + 1)
178
  return (y, z)
179

180
for a in A:
181
    curve = Pell_Curve(a, n)
182
    M1 = curve.mul(C, judge_d(a, p, q, d))
183
    x, y, z = M1
184
    x1 = roots2(x, p, r)
185
    x2 = roots2(x, q, s)
186
    for i in x1:
187
        for j in x2:
188
            m=long_to_bytes(crt([p**r, q**s], [int(i), int(j)])[0])
189
            if m.startswith(b'DASCTF'):
190
                print(m)
191
                exit(0)
192
# DASCTF{3dc7844aafe4e0628ba29ef09501089dd4a2adeec924b916be275cfc3953d681}______This_is_pad:adwd3i2j0fj20ef2j0fj20efj9h2j0fj20efj9h2j0fj2asdaedfqe0efj9h2j0fj20efj9qfqdqwdh2j0fj2qfefwfewfqwedfqwdqwdqw0efj9h2j0fj20efj9h2j0fj20efjwfewfwefwefwe9hj9huiehv89h92j0fj20efj9h8hwd893y198e32