SECCON 2020 OnlineCTF Crypto Write-Ups

I participated in SECCON 2020 Online CTF as a member of HangulSarang. We got 1st place :)

Hangul Day is also my birthday, so that's something I guess :D

The competition collided with ACM-ICPC quals, so I had to start solving at like 7PM. (4 hours late)

During that 4 hours, my teammate solved "This is RSA", which is solved using that each byte can be only 0x30 to 0x39. After that, we can compute solutions in $\pmod{2^k}$ and move upwards using recursion. It was the easiest crypto problem in this CTF, (and the solver wasn't me) so I think I don't need to explain further. The other four crypto problems were quite good, challenging, and very enjoyable, so I will describe my thought process as well. Huge thanks to the problem authors for creating these challenges. :D

koharu

while True:
    p = random_prime(1<<64)
    if is_prime((p+1) // 2):
        break
 
with open("flag.txt", "rb") as f:
    flag = f.read()
flag = int.from_bytes(flag, "big")
 
 
PR.<x> = PolynomialRing(GF(p))
while True:
    P = PR.random_element(degree=64)
    if P.is_irreducible():
        break
 
while True:
    Q = PR.random_element(degree=64)
    if Q.is_irreducible():
        break
 
NP = p**P.degree()
NQ = p**Q.degree()
 
while True:
    R = PR.random_element(degree=64)
    if power_mod(R, (NP-1)//2, P) != 1 and power_mod(R, (NQ-1)//2, Q) != 1:
        break
 
PQ = P*Q
c = []
while flag:
    S = PR.random_element(degree=64)
    if flag & 1:
        c.append((S * S) % PQ)
    else:
        c.append((S * S * R) % PQ)
    flag = flag >> 1
 
print("p =", p)
print("PQ =", PQ)
print("R =", R)
print("c =", c)
 
 

The code screams "quadratic residue", and it's similar to bitcrypto in InterKosenCTF. (write-up in this blog)

The only difference is that we are not using large primes, but polynomials in $GF(p)[x]$. This is already weak. 

The reason is that we can easily factorize polynomials in $GF(p)$. We can ask Sage to do it for us :)

Now we can determine whether a given polynomial is a quadratic residue or not. Combine this to get the flag.

p = 4832823609987476353
F.<x> = PolynomialRing(GF(p))
PQ = 2475361839038406994*x^128 + 1816580044636445865*x^127 + 771106714052997910*x^126 + 2532248969060743840*x^125 + 157159147928168793*x^124 + 1165294508775017303*x^123 + 54498477947855453*x^122 + 564670281176250610*x^121 + 4686383084102262935*x^120 + 4798143559496813901*x^119 + 2373759188753852032*x^118 + 3458843219210551923*x^117 + 3389173528515223367*x^116 + 3175114023644661971*x^115 + 2668820643276713526*x^114 + 1644657084961816584*x^113 + 1949973045428555331*x^112 + 2314884799372359978*x^111 + 1614909032209480656*x^110 + 3706101120120959039*x^109 + 1443476119293487220*x^108 + 507539962924420368*x^107 + 2851578707595377440*x^106 + 2660707099322090529*x^105 + 2275120831055073492*x^104 + 4642644673121099806*x^103 + 780741129747777966*x^102 + 3824963851609159359*x^101 + 1445016816241934269*x^100 + 4706494165496469049*x^99 + 91460120231848540*x^98 + 2033361932245472629*x^97 + 4657205830657809352*x^96 + 627579987075662316*x^95 + 2638155163726745709*x^94 + 773248040814209977*x^93 + 4426134463977473378*x^92 + 1748835523159978170*x^91 + 2545886874835388035*x^90 + 4318027045196127783*x^89 + 529092995613843935*x^88 + 37621695756851259*x^87 + 724317479549357114*x^86 + 235872728824864204*x^85 + 1409136599403563059*x^84 + 984842291673572708*x^83 + 1000642979551429427*x^82 + 2599952022893048437*x^81 + 33489199855748196*x^80 + 2138571356326295553*x^79 + 357904099457660261*x^78 + 1388605866466399741*x^77 + 2123614714168365349*x^76 + 1296407111118101425*x^75 + 3175149128196009486*x^74 + 4407671566428651830*x^73 + 3653949472018283742*x^72 + 2150666969917189331*x^71 + 2425834809198809729*x^70 + 202017664024051124*x^69 + 4656859267960293209*x^68 + 95544718007904685*x^67 + 551963924883187932*x^66 + 1220133766833256737*x^65 + 418789913385574936*x^64 + 3140425594489130574*x^63 + 653426727346469624*x^62 + 2168508737790275670*x^61 + 1350675684196344669*x^60 + 86970043713584944*x^59 + 3125122442296761190*x^58 + 1691082709013935740*x^57 + 14954357710735056*x^56 + 1951640599446313225*x^55 + 3057759244385615044*x^54 + 2842299299534580663*x^53 + 60118912044101305*x^52 + 3791459205438092561*x^51 + 3961025931327708139*x^50 + 3352223936735193809*x^49 + 458087980170556413*x^48 + 303065746752057039*x^47 + 270269323703788403*x^46 + 3435561048914221019*x^45 + 244980776425782882*x^44 + 1756735569264346021*x^43 + 1049402079460555244*x^42 + 1181023304135761892*x^41 + 2480814159047994100*x^40 + 3359295278584507081*x^39 + 1031815312165038169*x^38 + 2284789340145013050*x^37 + 2507227047920435897*x^36 + 4212274843760760739*x^35 + 1874163516348469998*x^34 + 4184876619139253979*x^33 + 2454055493008310058*x^32 + 4810631595605704078*x^31 + 2705618732956794205*x^30 + 4588422028499215564*x^29 + 1362947071518584749*x^28 + 200625668549982104*x^27 + 4162225127389871946*x^26 + 3671964574429446847*x^25 + 497776717675475749*x^24 + 3171362364421276926*x^23 + 4040585504650270495*x^22 + 55143980688943936*x^21 + 1680279432641096886*x^20 + 1141249890787830167*x^19 + 1632171956841566025*x^18 + 4489792289887403690*x^17 + 72863318133800422*x^16 + 3512973315964270180*x^15 + 1880837549990432714*x^14 + 629108155937185931*x^13 + 605563550674482475*x^12 + 3125052390516629852*x^11 + 3434353753938817079*x^10 + 2199180089161294937*x^9 + 4128993677150612079*x^8 + 875038461592559534*x^7 + 1344699457303227348*x^6 + 3605318452000064928*x^5 + 1825112182884559504*x^4 + 4214849563830404245*x^3 + 3018789469914511583*x^2 + 4256870332540451928*x + 3478109193918270445
R = 10529800129354981*x^64 + 4658846300069202283*x^63 + 1343603688498785880*x^62 + 77535778799313918*x^61 + 3909004297055292936*x^60 + 1574062357470841720*x^59 + 2255026177942473610*x^58 + 2913895405335010190*x^57 + 910153010204378491*x^56 + 4823161627331431259*x^55 + 4314926186108070132*x^54 + 3776194104903441585*x^53 + 4218241384907734159*x^52 + 2928099962473177675*x^51 + 3620663369166129209*x^50 + 4671199329340054093*x^49 + 2953252709684913819*x^48 + 1470028746745533363*x^47 + 393509208258687360*x^46 + 2631641671658679748*x^45 + 4823463900549231672*x^44 + 22025139085889956*x^43 + 3905072220448754367*x^42 + 3525611426409694274*x^41 + 1087703571442464513*x^40 + 983613039355879671*x^39 + 2292836760450398296*x^38 + 2429042383184252432*x^37 + 4241866215562144008*x^36 + 3567456235250802214*x^35 + 289826756486726727*x^34 + 3070079221437908111*x^33 + 3164478508626375897*x^32 + 4028195041942471423*x^31 + 1611744044712776226*x^30 + 682031605725048858*x^29 + 2334009162012075842*x^28 + 1056698946696323305*x^27 + 1193918408929283326*x^26 + 1546583097398597126*x^25 + 632624061599387394*x^24 + 3924194912006864689*x^23 + 836241738980292724*x^22 + 2019639656826418643*x^21 + 646182266409329495*x^20 + 3568811299250961381*x^19 + 4024124722170180214*x^18 + 2765626713849083593*x^17 + 830125243533734584*x^16 + 3773807917205041413*x^15 + 4579071273569219071*x^14 + 4169012455774239610*x^13 + 2779202281389813792*x^12 + 1668767138196611027*x^11 + 3668902156196312613*x^10 + 2118966174503976203*x^9 + 2876683474352545557*x^8 + 4749450906737437136*x^7 + 2048549559963146669*x^6 + 2337906091414592304*x^5 + 3234395871197583532*x^4 + 624006023034932764*x^3 + 1020142386943254010*x^2 + 4346889740151908150*x + 2337193413394346074
c = []
P = (x^64 + 2705838326093066801*x^63 + 1861763125820805142*x^62 + 1919270169024731361*x^61 + 728192979251886197*x^60 + 3703504742135431297*x^59 + 608310330267197202*x^58 + 677522369546315305*x^57 + 45111914222503868*x^56 + 3231090245423531905*x^55 + 4439626063971680541*x^54 + 264779255326565930*x^53 + 943573327092647824*x^52 + 3642035360519473519*x^51 + 4624797912514728904*x^50 + 815168423497123035*x^49 + 2058290770523809000*x^48 + 4368972367338353614*x^47 + 1102710837251449034*x^46 + 1838631000574578462*x^45 + 1550208773716319692*x^44 + 4479635398032603580*x^43 + 2547505501081696879*x^42 + 4733577241261296757*x^41 + 1459044726889718801*x^40 + 4736670792998507780*x^39 + 3481084975759672453*x^38 + 4491590348438475003*x^37 + 4286960290474469508*x^36 + 2519824328645383346*x^35 + 722570560813334776*x^34 + 3203376079187925593*x^33 + 2137713042365333594*x^32 + 2529680584881125743*x^31 + 881878615185959251*x^30 + 2648895700342509353*x^29 + 3093613170934869890*x^28 + 1839149659686122740*x^27 + 901352037355979824*x^26 + 3079388294575162468*x^25 + 4316897640303347156*x^24 + 3768144827267250554*x^23 + 1585476600468626452*x^22 + 2408180731465025131*x^21 + 2754322334879778466*x^20 + 1965864600205111832*x^19 + 3016989393277154199*x^18 + 993850365653028982*x^17 + 1661221355151932055*x^16 + 2141520480611688809*x^15 + 636670112723307258*x^14 + 1200822100799196786*x^13 + 2223563845526420680*x^12 + 3134534498746508642*x^11 + 1820327632682349699*x^10 + 4628418849122802568*x^9 + 3731553570235638636*x^8 + 1636534607043587796*x^7 + 1007966122754856335*x^6 + 3571611463638839115*x^5 + 4733247188903796455*x^4 + 3512981852602786831*x^3 + 1560667366459827025*x^2 + 1113839338158290233*x + 4011393002849553527)
Q = (x^64 + 2776622066009961678*x^63 + 1020248994272362724*x^62 + 1812889731002797017*x^61 + 3946133096396475132*x^60 + 1064362775780462120*x^59 + 4267166204741846229*x^58 + 4461168980925876722*x^57 + 3701193932757315736*x^56 + 4004259984657770019*x^55 + 2566923830139634808*x^54 + 2958380329303059106*x^53 + 4642913814072279374*x^52 + 713990683265973444*x^51 + 2282781718594249732*x^50 + 1691679008052617295*x^49 + 4723620313305465430*x^48 + 4052669689859242595*x^47 + 4607757741831461143*x^46 + 3048879536065529044*x^45 + 2012013680568798151*x^44 + 2125237235418450484*x^43 + 2622384625077739224*x^42 + 710661875195936255*x^41 + 375897308743404378*x^40 + 3253268532586707019*x^39 + 3759767504239334681*x^38 + 2945194932005180334*x^37 + 1316716161821289054*x^36 + 2210075866459201344*x^35 + 3421886933443572088*x^34 + 2124192011313760002*x^33 + 3183242335232871177*x^32 + 4722704310996441203*x^31 + 1640862527872462873*x^30 + 292078618156889354*x^29 + 3970255331239899451*x^28 + 290424178543927660*x^27 + 3979382049081683506*x^26 + 3341058157535181184*x^25 + 1891458780676141416*x^24 + 4585931142037966308*x^23 + 2621586816910493860*x^22 + 4526407296014662985*x^21 + 3345825075365423903*x^20 + 433595205227433076*x^19 + 3510443356995660854*x^18 + 1469161865274264871*x^17 + 1968552305256496645*x^16 + 1902262417167822976*x^15 + 3211385257470450715*x^14 + 259183745852362935*x^13 + 1368548986536267534*x^12 + 3726482530039832086*x^11 + 1196244075361051439*x^10 + 3346319329141804238*x^9 + 2362535635162047034*x^8 + 2131037938034625812*x^7 + 3970887869581347678*x^6 + 4428522899784697485*x^5 + 2482987898184812388*x^4 + 3180131420672415636*x^3 + 4690602932003451909*x^2 + 2572790493146370264*x + 802891458181310745)
cv = (p ** 64 - 1) // 2
fin = 0
add = 1
for ply in c:
    Pv = power_mod(ply, cv, P)
    Qv = power_mod(ply, cv, Q)
    if Pv == 1 and Qv == 1:
        fin += add
    add = add * 2
print(fin) ## long to bytes here

urara

from flag import flag
 
p = random_prime(1 << 1024)
q = random_prime(1 << 1024)
n = p * q
 
print("n =", n)
 
# ---
 
x = int.from_bytes(flag, "big")
y = randint(0, n-1)
 
a = randint(0, n-1)
b = (y^2 - (x^3 + a*x)) % n
 
EC = EllipticCurve(Zmod(n), [a, b])
 
P = EC((x, y))
Q = 2 * P
 
print("a, b =", [a, b])
print("Q =", Q.xy())
 
# ---
 
m = int.from_bytes(flag, "big")
t = randint(0, n-1)
 
c = power_mod(m + t, 65537, n)
print("t =", t)
print("c =", c)
 

Very concise problem. We are given $t, e, n$ and $(m+t)^e \pmod{n}$, and a point is double the point with [$x$-coordinate equal to $m$] [which lies on the elliptic curve $y^2 = x^3 + ax + b$]. I broke that sentence up into pieces because it's pretty long and confusing. So anyways, doing actually elliptic curve stuff seems nearly impossible with the lack of hints we have, and we have a huge hint in $(m+t)^e \pmod{n}$. If we have another polynomial equation about $m$ in $\mathbb{Z}_n[x]$, we can solve this problem using polynomial GCD. (Franklin-Reiter related message attack, like padrsa in InterKosen CTF) Of course, we do have another polynomial equation! Directly using the doubling formula, we have $$ (3m^2 + a)^2 - (2m + Q_x)(4(m^3 + am+b)) \equiv 0 \pmod{n}$$ which is the type of equation we want. Write polynomial GCD, and we can easily find $m$. Again, note that if we encounter a problem while performing standard Euclidean Algorithm, it's because the leading coefficient has a non-trivial GCD with $n$. In this case, we can just factorize $n$ to solve the problem. I was stuck in sharsable for so long before solving this one, so I got some energy from it :D 

For the implementation of polynomial GCD, refer to the write-up of the problem padrsa. 

K = Zmod(n)
P.<x> = PolynomialRing(K, implementation='NTL')
f = (3 * x^2 + a)^2 - (2*x + Qx) *(4*(x^3 + a*x+b))
g = power_mod(x + t, 65537, f) - c
print(GCD(f, g, n))
## the remaining details are trivial

sharsable

from Crypto.Util.number import getPrime, GCD
from flag import FLAG
import random
 
def egcd(a, b):
    r0, r1 = a, b
    s0, s1 = 1, 0
    t0, t1 = 0, 1
    while r1 > 0:
        q = r0 // r1
        r0, r1 = r1, r0 % r1
        s0, s1 = s1, s0 - q * s1
        t0, t1 = t1, t0 - q * t1
    return s0, t0
 
def generateKey():
    p = getPrime(512)
    q = getPrime(512)
    n = p * q
    phi = (p-1)*(q-1)
 
    while True:
        d1 = getPrime(int(n.bit_length()*0.16))
        e1 = random.randint(1, phi)
        ed1 = e1 * d1 % phi
 
        d2 = getPrime(int(n.bit_length()*0.16))
        e2, k = egcd(d2, phi)
        e2 = e2 * (phi + 1 - ed1) % phi
        ed2 = e2 * d2 % phi
 
        if GCD(e1, e2) > 10:
            break
 
    assert((ed1 + ed2) % phi == 1)
 
    return (n, (e1, d1), (e2, d2))
 
n, A, B = generateKey()
M = int.from_bytes(FLAG, 'big')
C1 = pow(M, A[0], n)
C2 = pow(M, B[0], n)
assert(pow(C1, A[1], n) * pow(C2, B[1], n) % n == M)
 
import json
print(json.dumps({
    "n": n,
    "A": (A[0], C1),
    "B": (B[0], C2),
    #"d": (A[1], B[1]), # for debug
    }))
 

The code looks convoluted, but it can be compressed (without information loss) into the following results:

• $d_1, d_2$ are quite small, around $n^{0.16}$. 
• $e_1d_1 + e_2d_2 \equiv 1 \pmod{\phi(n)}$.
• We know $e_1, e_2, n$, but we can't use common modulus attack due to $\text{gcd}(e_1, e_2) = 11$.

I didn't even realize we could use common modulus attack here, but that was okay since it doesn't give us useful information anyway. The key clearly lies in the small $d_1, d_2$. At first, my thought was headed to the Wiener's Attack, which I think is justifiable because that attack works with small $d$ as well. Of course, it's hard (if not impossible) to use the continued fraction idea directly. After reading up on some generalizations of Wiener's Attacks, I thought this problem was related to the coppersmith attack. I tried weird stuff like 3-variable coppersmith but it all failed pretty badly. That led to me solving this problem the last out of the four I solved. The idea of the last problem helped me :)

We begin by writing (note that this kind of manipulation is done in Wiener's Attack as well) $$e_1 d_1 + e_2 d_2 = k \phi(n) + 1 = k(n-p-q+1) + 1 \equiv -k(p+q-1) + 1 \pmod{n}$$ We may now note that $k$ is probably around $n^{0.16}$ as well, so the RHS is around $-n^{0.66}$ multiplied by some constant. 

Okay, so we want $d_1, d_2$ to be around $n^{0.16}$, and $e_1d_1 + e_2d_2 \pmod{n}$ to be around $-n^{0.66} \cdot c$ for some "reasonable" constant $c$. Let's first decide the range of $c$ we will look for. We can expect $k$ to be $0.4 \cdot n^{0.16}$ to $2 \cdot n^{0.16}$. We can expect $p+q$ to be $2 \cdot n^{0.5}$ to $3/\sqrt{2} \cdot n^{0.5}$. Combining, we can get something like $c \in [0.8, 4]$. Note that I'm using a lot of handwaving. Now we will fix $c$. 

The result? We have goal values for $d_1, d_2, e_1d_1 + e_2d_2 \pmod{n}$. Sounds like CVP to me now. 

Indeed, we can work this as a CVP with the lattice generated by $$ [e_1, 1, 0], \quad [e_2, 0, 1], \quad [n, 0, 0] $$ and making the result vector as close as possible to $$ [-c \cdot n^{0.66}, n^{0.16}, n^{0.16} ] $$ Here's a problem. The scale between the first column and the next two columns are off by $n^{0.5}$. If we run CVP now, it'll probably ignore the $n^{0.16}$ condition in order to get similar values for the first column. At least, that's what I thought :D 

Therefore, we rescale. We will work with the lattice generated by $$  [e_1, n^{0.5}, 0], \quad [e_2, 0, n^{0.5}], \quad [n, 0, 0] $$ and making our result vector as close as possible to $$ [-c \cdot n^{0.66}, n^{0.66}, n^{0.66} ] $$ which sounds more reasonable. If we run the CVP algorithm, (Babai's) we will have our candidate values for $d_1, d_2$. 

If correct, we will have that $e_1d_1 + e_2d_2 - 1$ is a product of $\phi(n)$, so we can finish easily. So, how do we find the right $c$ to use?

Simple. Just try a whole bunch of them. Refer to the code below for details.  

def Babai_closest_vector(M, G, target):
        small = target
        for _ in range(1):
            for i in reversed(range(M.nrows())):
                c = ((small * G[i]) / (G[i] * G[i])).round()
                small -=  M[i] * c
        return target - small 
    
n  = 142793817321992828777925840162504083304079023834001118099549928854335392622287928254035247188624975743042449746066633491912316354241339908190889792327014012472372654378644158878787350693992259970146885854641856991605625756536504266728483088687985429310233421251081614258665472164668993082471923690196082829593
e1 = 82815162880874815458042429141267540989513396527359063805652845923737062346339641683097075730151688566721221542188377672708478777831586255213972947470222613130635483227797717393291856129771004300757155687587305350059401683671715424063527610425941387424425367153041852997937972925839362190900175155479532582934
C1 = 108072697038795075732704334514926058617161875495016327352871122917196026504758904760148391499245235850616838765611460630089577948665981247735905622903872682862860306107704253287284051312867625831877418240290183661755993649928399992531008191618616452091127799880839665225093055618092869662205901927957599941568
e2 = 84856171747859965508406237198459622554468224770252249975158471902036102010991476445962577679301719179079633469099994226630172251817358960347828156301869905575867853640850107406452911333646573296923235424617864473580743418995994067645338437540627399276292679100115018844287273293945121023787594592185295794983
C2 = 101960082023987498941061751761131381167414505957511290567652602520714324823481487410890478130601013005035303795327512367595187718926017321227779179404306882163521882309833982882201152721855538832465833869251505131262098978117904455226014402089126682222497271578420753565370375178303927777655414023662528363360
 
rat = int(n ** 0.5)
M = Matrix(ZZ, 3, 3)
M[0, 0] = e1
M[1, 0] = e2
M[2, 0] = n
M[0, 1] = rat
M[1, 2] = rat
 
L = []
for i in range(800, 2450):
    Target = vector([-int(n ** 0.66 * i / 1000) , int(n ** 0.16) * rat, int(n ** 0.16) * rat])
    M = M.LLL()
    GG = M.gram_schmidt()[0]
    TT = Babai_closest_vector(M, GG, Target)
    d1 = TT[1] // rat
    d2 = TT[2] // rat
    L.append(e1 * d1 + e2 * d2 - 1)
print(list(set(L))) ## this is c
 
for phi in c:
    d1 = inverse(e1, phi)
    print(long_to_bytes(pow(C1, d1, n)))
    d2 = inverse(e2, phi)
    print(long_to_bytes(pow(C2, d2, n)))

  

crypto01

from sage.all import *
from flag import flag
from functools import reduce
 
def encrypt(m, e, n):
    n = int(n)
    size = n.bit_length() // 2
    m_low = m & ((1 << size) - 1)
    m_high = (m >> size)
 
    b = (m_low**2 - m_high**3) % n
    EC = EllipticCurve(Zmod(n), [0, b])
 
    return (EC((m_high, m_low)) * e).xy()
 
def decrypt(c, d, n):
    n = int(n)
    size = n.bit_length() // 2
 
    c_high, c_low = c
    b = (c_low**2 - c_high**3) % n
    EC = EllipticCurve(Zmod(n), [0, b])
    m_high, m_low = (EC((c_high, c_low)) * d).xy()
    m_high, m_low = int(m_high), int(m_low)
 
    return (m_high << size) | m_low
 
def gen_prime(size):
    p = random_prime(1 << size)
    while p % 3 != 2:
        p = random_prime(1 << size)
 
    q = random_prime(1 << size)
    while q % 3 != 2:
        q = random_prime(1 << size)
 
    if q > p:
        p, q = q, p
 
    return int(p), int(q)
 
 
 
SIZE = 512
HINTSIZE = 96
N = 3
 
flag = int.from_bytes(flag, "big")
assert flag < (1 << SIZE)
 
masks = [randint(1 << (SIZE-1), 1 << SIZE) for _ in range(N)]
masked_flag = reduce(lambda a, b: a ^ b, masks, flag)
 
 
count = 0
ciphertexts = []
while count < N:
    try:
        p, q = gen_prime(SIZE)
        n = p * q
 
        x = random_prime(int(n ** 0.40))
        y = random_prime(int(sqrt(2 * n // (144 * x*x))))
        zbound = -1 * int(round(((p-q) * (n ** 0.25) * y) / (3 * (p + q))))
 
        z_ = zbound + ((p + 1)*(q + 1)*y - zbound) % x
        e = ((p + 1) * (q + 1) * y - z_) // x
        d = inverse_mod(e, (p + 1)*(q + 1))
 
        assert (x*y*x*y < (2 * n // 144))
        assert (gcd(x, y) == 1)
 
        d = inverse_mod(e, (p+1)*(q+1))
        c = encrypt(masks[count], e, n)
        assert decrypt(c, d, n) == masks[count]
 
        ciphertexts.append({
            "n": n,
            "e": e,
            "c": c,
            "hint": p & ((1<<HINTSIZE)-1)
        })
        count += 1
    except KeyboardInterrupt:
        break
    except (ZeroDivisionError, OverflowError):
        pass
 
 
print("masked_flag = " ,masked_flag)
print("ciphertexts = ", ciphertexts)
 

@diff (pcw) noted that the order of the elliptic curve must be $(p+1)(q+1)$ since $p \equiv q \equiv 2 \pmod{3}$.

This is a cool fact that I didn't know, but it could be easily guessed by the definition of $d$ in the script. 

So, first things first : is this really a elliptic curve challenge? The answer : probably not?

I thought this for a few reasons - I thought the generation of $e$ was quite convoluted so there should be a weakness.

The elliptic curve part looked too clean to have a vulnerability, and it's not even over a field. Hard to do anything, really. 

Therefore, I decided to take a look at how the parameter generation behaves.

First, I look at the sizes of the values. It 's clear that $y$ is small, so $zbound$ and $z$ are also quite small compared to $x$.

This leads to the following suspicion, is $e = \lfloor (p+1)(q+1)y / x \rfloor$? Generating parameters on our own, we see that this is true with a very high probability. This means that we can ignore $zbound$ and $z$ completely now! This is a good progress indeed :) 

This simplifies a lot of things, and we have a good, clean equation to work with. Start with $$ xe  = (p+1)(q+1)y - r \equiv (p+q+1)y - r\pmod{n}$$ with $r < x \approx n^{0.4}$. With $y \approx n^{0.1}$, we see that $(p+q+1)y - r \approx n^{0.6}$. So $x \approx n^{0.4}$ satisfies $xe \pmod{n} < \approx n^{0.6}$. 

Heuristically speaking, there should be a fairly small number of $x$ that satisfy that condition! Can we find them though?

Turns out the answer is yes. I guess you can do this with lattices, but here's a method I learned from competitive programming. 

In fact, given $L, R, M, A$, we can find the minimum nonnegative integer $x$ such that $L \le Ax \pmod{M} \le R$. 

The method is not long, but I don't want to type it all out, so here's a link (scroll down to the end of the editorial)

Also, a problem using this algorithm appeared in NWRRC, (problem G) set by tourist. 

Since we know $e, n$, we can select some $CUT \approx n^{0.6}$ and find the minimum $x$ such that $1 \le ex \pmod{n} \le CUT$. Of course, we have to be cautious with constants when selecting $CUT$. It's not too hard, since we can always verify our value of $x$ by checking its primality and size. This concludes the recovery of $x$, which is a huge progress again! Now we try to recover $y$. 

To recover $y$, we use basic inequalities. Note the following two inequalities $$ y = \frac{xe+r}{(p+1)(q+1)} \ge \frac{xe}{n + 3 \sqrt{n/2} + 1}$$ $$ y = \frac{xe+r}{(p+1)(q+1)} \le \frac{xe + x-1}{n + 2\sqrt{n} + 1}$$ which gives us decent lower/upper bounds. Turns out that this is also perfect, as the two bounds are actually equal! This recovers $y$. 

With this information, we can now bound $p+q$. Simply use $$ p+q = \frac{xe+r}{y} - 1 \ge \frac{xe}{y} - 1 $$ $$ p+q  = \frac{xe+r}{y} -1 \le \frac{xe+x-1}{y} - 1$$ which gives us a good bound for $p+q$. Now we proceed as we did in l33tcrypt in DownUnderCTF and use quadratic formula to get a good bound for $p$. Note that $p > q$. If we analyze the two bounds, we see that this gives us about 200 MSBs of $p$. Now we are practically done. 

We already know 96 LSBs of $p$ via the hint. This adds up to 296 bits of information, so we only need to determine 216 bits of $p$. This can be easily done with coppersmith method, i.e. Sage's small_roots. This part is very straightforward. Now we calculate $p, q, d$ and eventually the masks. The conversion of masks to the actual flag was done by @ironore15. This concludes the problem :D

Also, maybe you can recover $x, y$ with continued fractions on $e/n$? Maybe. I didn't try it :P

UPD : Okay I tried it and it works easily. Why didn't I think about this lmao... This is a much easier way to solve. :)

Basically, from $e = \lfloor (p+1)(q+1)y / x \rfloor$, we can guess that $$ \frac{e}{n} \approx \frac{e}{(p+1)(q+1)} \approx \frac{y}{x}$$ so $y/x$ is a good approximation of $e/n$. Find all convergents of $e/n$, and find those that have primes as its numerator and denominator. 

This gives us $x, y$, so we can proceed with the coppersmith attack as above. Much easier :)

UPD : The continued fraction solution is indeed the intended sol. I guess it could be proved using the given bounds :)

UPD : So apparently this thing had a name, KMOV cryptosystem. didn't know that :P

def kthp(n, k):
    if n == 0:
        return 0
    lef = 1
    rig = 2
    while rig ** k < n:
        rig = rig << 1
    while lef <= rig:
        mid = (lef + rig) // 2
        if mid ** k <= n:
            best = mid
            lef = mid + 1
        else:
            rig = mid - 1
    return best
 
def ceil(n, m):
    return (n + m - 1) // m
 
def optf(A, M, L, R):
    if L == 0:
        return 0
    if 2 * A > M:
        L, R = R, L
        A = M - A
        L = M - L
        R = M - R
    cc_1 = ceil(L, A)
    if A * cc_1 <= R:
        return cc_1
    cc_2 = optf(A - M % A, A, L % A, R % A)
    return ceil(L + M * cc_2, A)
 
def decrypt(c, d, n):
    n = int(n)
    size = n.bit_length() // 2
 
    c_high, c_low = c
    b = (c_low**2 - c_high**3) % n
    EC = EllipticCurve(Zmod(n), [0, b])
    m_high, m_low = (EC((c_high, c_low)) * d).xy()
    m_high, m_low = int(m_high), int(m_low)
 
    return (m_high << size) | m_low
 
ciphertexts =  []
for C in ciphertexts:
    n = C['n']
    e = C['e']
    c = C['c']
    hint = C['hint']
    CUT = kthp( (int)(n ** 0.2) // 72, 2) * kthp(n, 2) * 8
    x = optf(e, n, 1, CUT)
    R = ((e * x + x - 1) // (n + 2 * kthp(n, 2) + 1))
    L = ((e * x) // (n + 3 * kthp(n//2, 2) + 1))
    y = (L + R) // 2
    assert L == R and x in Primes() and y in Primes()
    sum_L = (x * e) // y - 1 - n
    sum_R = (x * e + x - 1) // y - 1 - n
    lr = (sum_R + kthp(sum_R * sum_R - 4 * n, 2)) // 2
    sm = (sum_L + kthp(sum_L * sum_L - 4 * n, 2)) // 2
    assert sm <= lr
    assert (sm >> 312) == (lr >> 312)
    p_hint = hint
    K = Zmod(n)
    P.<t> = PolynomialRing(K, implementation='NTL')
    f = (p_hint * inverse_mod(2 ** 96, n)) % n + t + (2 ** (312-96)) * (lr >> 312)
    x0 = f.small_roots(X = (2 ** 220), beta = 0.5, epsilon = 0.03)
    ## print(x0)
    p = p_hint + x0[0] * (2 ** 96) + (2 ** 312) * (lr >> 312)
    p = (int)(p)
    q = n // p
    d = inverse_mod(e, (p+1) * (q+1))
    print(decrypt(c, d, n))
 
## thanks, ironore!
def recover_flag(masks, masked_flag):
    flag = reduce(lambda a, b: a ^ b, masks, masked_flag)
    return flag.to_bytes(512 // 8, 'big')

def get_red(e, n):
    cur_num, cur_den = e, n
    num_1, den_1 = 0, 1
    num_2, den_2 = 1, 0
    while True:
        val = cur_num // cur_den
        nxt_num = cur_den
        nxt_den = cur_num - val * cur_den
        # calculate new convergent
        num_3 = val * num_2 + num_1
        den_3 = val * den_2 + den_1
        if isPrime(num_3) and isPrime(den_3):
            return num_3, den_3
        if den_3 > int(n ** 0.4):
            return -1
        # update convergents
        num_1, den_1 = num_2, den_2
        num_2, den_2 = num_3, den_3
        # update continued fractions
        cur_num, cur_den = nxt_num, nxt_den
 
ciphertexts =  []
for C in ciphertexts:
    n = C['n']
    e = C['e']
    c = C['c']
    hint = C['hint']
    y, x = get_red(e, n)
    R = ((e * x + x - 1) // (n + 2 * kthp(n, 2) + 1))
    L = ((e * x) // (n + 3 * kthp(n//2, 2) + 1))
    assert y == (L + R) // 2
    assert L == R and isPrime(x) and isPrime(y)
    ## continue as the initial solution