This week, the wonderful HackTM CTF was up and running, and I tried my best at the various crypto tasks available. In this post, I’ll go over the first two that I managed to solve, they were a ton of fun and I had the opportunity of learning a lot, so once again thanks y011d4 for writing really cool challenges.

d-phi-enc⌗ 

from Crypto.Util.number import bytes_to_long, getStrongPrime
from secret import flag
assert len(flag) == 255
e = 3
p = getStrongPrime(1024, e=e)
q = getStrongPrime(1024, e=e)
n = p * q
phi = (p - 1) * (q - 1)
d = pow(e, -1, phi)
enc_d = pow(d, e, n)
enc_phi = pow(phi, e, n)
enc_flag = pow(bytes_to_long(flag), e, n)
print(f"{n = }")
print(f"{enc_d = }")
print(f"{enc_phi = }")
print(f"{enc_flag = }")

This challenge looks like standard RSA, (note that $e = 3$), but we have some significant additional information, namely $d^{*} = d^{3} \bmod n$ and $\varphi^{*} = \varphi(n)^{3} \bmod n$.

At the start, I was looking for some way to decrypt without factoring, as I thought there were some tricks with $m^{e^{3}d^{3}} \bmod n$. After rejecting the idea, I came up with a better way to go about it (albeit quite tedious and painful). We can essentially write the following:

$$ \varphi(n)^{3} \equiv (p-1)^{3}(q-1)^{3} \equiv - p^3 - q^3 + 3p^2 + 3q^2 - 3p - 3q + 1 \bmod n $$ $$ f(p, q) = - p^3 - q^3 + 3p^2 + 3q^2 - 3p - 3q + 1 - \varphi^{*} \equiv 0 \bmod n$$

Note that $d^{*} = d^{3} - kn$. We can extract another polynomial which has $(p, q)$ as roots, by looking at the follwing chain of equations:

$$ e^{3}d^{*} = e^{3}(d^{3}-kn) = (ed)^{3} - e^{3}kn \equiv (1 + j\varphi(n)) ^ {3} \bmod n \equiv (1 + j(p-1)(q-1))^3 \bmod n $$

Hence, we let $g(p, q) = (1 + j(p-1)(q-1))^3 - e^3d^{*} \equiv 0 \bmod n$, where $0 \leq j \leq e$, easily bruteforcable.

Our last polynomial $h(p, q) = pq \equiv 0 \bmod n$ makes it so we have 3 different polynomials in $(p, q)$ sharing a common root, namely the two prime factors of $n$.

At this point, we can pray the gods and hope to extract some info from doing resultants and gcds of the polynomials. Doing so, shows that the polynomial we obtain after tinkering with the original 3 has a weird property, which is that it factors nicely into $x^3(x + k)$. Heuristics also show that $k = n -(p+q)$, so we can recover $p, q$ and then decrypt the flag.

from Crypto.Util.number import long_to_bytes

def my_gcd(a, b): 
    return a.monic() if b == 0 else my_gcd(b, a % b)

n = 24476383567792760737445809443492789639532562013922247811020136923589010741644222420227206374197451638950771413340924096340837752043249937740661704552394497914758536695641625358888570907798672682231978378863166006326676708689766394246962358644899609302315269836924417613853084331305979037961661767481870702409724154783024602585993523452019004639755830872907936352210725695418551084182173371461071253191795891364697373409661909944972555863676405650352874457152520233049140800885827642997470620526948414532553390007363221770832301261733085022095468538192372251696747049088035108525038449982810535032819511871880097702167
enc_d = 23851971033205169724442925873736356542293022048328010529601922038597156073052741135967263406916098353904000351147783737673489182435902916159670398843992581022424040234578709904403027939686144718982884200573860698818686908312301218022582288691503272265090891919878763225922888973146019154932207221041956907361037238034826284737842344007626825211682868274941550017877866773242511532247005459314727939294024278155232050689062951137001487973659259356715242237299506824804517181218221923331473121877871094364766799442907255801213557820110837044140390668415470724167526835848871056818034641517677763554906855446709546993374
enc_phi = 3988439673093122433640268099760031932750589560901017694612294237734994528445711289776522094320029720250901589476622749396945875113134575148954745649956408698129211447217738399970996146231987508863215840103938468351716403487636203224224211948248426979344488189039912815110421219060901595845157989550626732212856972549465190609710288441075239289727079931558808667820980978069512061297536414547224423337930529183537834934423347408747058506318052591007082711258005394876388007279867425728777595263973387697391413008399180495885227570437439156801767814674612719688588210328293559385199717899996385433488332567823928840559
enc_flag = 24033688910716813631334059349597835978066437874275978149197947048266360284414281504254842680128144566593025304122689062491362078754654845221441355173479792783568043865858117683452266200159044180325485093879621270026569149364489793568633147270150444227384468763682612472279672856584861388549164193349969030657929104643396225271183660397476206979899360949458826408961911095994102002214251057409490674577323972717947269749817048145947578717519514253771112820567828846282185208033831611286468127988373756949337813132960947907670681901742312384117809682232325292812758263309998505244566881893895088185810009313758025764867
e=3

F = Zmod(n)
P.<p, q> = PolynomialRing(F)
p, q = P.gens()
j=2 #after some testing

f1 = - p^3 - q^3 + 3*p^2 + 3*q^2 - 3*p - 3*q + 1 - ZZ(enc_phi)
f2 = - 8*p^3 - 8*q^3 + 36*p^2 + 36*q^2 - 54*p - 54*q + 27 -ZZ((e^3 * enc_d)%n)
ffin = p*q

fout = det(ffin.sylvester_matrix(f1, p)).univariate_polynomial()
fout2 = det(ffin.sylvester_matrix(f2, p)).univariate_polynomial()
out = my_gcd(fout, fout2)
#out is of the form f(x) = x^3 * (x + k). for q to be a root of this poly mod n, we need q+k= 0 mod p
#experimentally, it appears that k = n-(p+q), so we can rearrange and obtain p and q

k = out.list()[3]
pplusq = n - int(k)

p1 = (pplusq + isqrt(pplusq^2 - 4*n))//2
q1 = n // p1
assert(p1*q1 == n)

phi = (p1-1)*(q1-1)
d = pow(e, -1, phi)
pt = pow(enc_flag, d, n)
print(long_to_bytes(int(pt)))

kaitenzushi⌗ 

from math import gcd
from Crypto.Util.number import bytes_to_long, isPrime

from secret import p, q, x1, y1, x2, y2, e, flag

# properties of secret variables
assert isPrime(p) and p.bit_length() == 768
assert isPrime(q) and q.bit_length() == 768
assert isPrime(e) and e.bit_length() == 256
assert gcd((p - 1) * (q - 1), e) == 1
assert x1.bit_length() <= 768 and x2.bit_length() <= 768
assert y1.bit_length() <= 640 and y2.bit_length() <= 640
assert x1 ** 2 + e * y1 ** 2 == p * q
assert x2 ** 2 + e * y2 ** 2 == p * q

# encrypt flag by RSA, with xor
n = p * q
c = pow(bytes_to_long(flag) ^^ x1 ^^ y1 ^^ x2 ^^ y2, e, n)
print(f"{n = }")
print(f"{c = }")


# hints 🍣
F = RealField(1337)
x = vector(F, [x1, x2])
y = vector(F, [y1, y2])
# rotate
theta = F.random_element(min=-pi, max=pi)
R = matrix(F, [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]])
x = R * x
y = R * y
print(f"{x = }")
print(f"{y = }")

Recovering the hidden values⌗

First of all, we have a bunch of hidden values, satisfying a bunch of properties. We also have the counter-clockwise rotation of two vectors, $x = (x_{1}, x_{2})$ and $y = (y_{1}, y_{2})$ in $R^{2}$.

Recall that rotation doesn’t influence the norm of the vectors, so $ ||x||^{2} = ||x_{rot}||^{2} = x_1^2 + x_2^2 $. We can recover $e$ by summing the two equations at the start of the script, and solving with our values of $x_1^2 + x_2^2$ and $y_1^2 + y_2^2$.

In theory, because rotation is done over the reals, we don’t precisely have the norm of the originals vectors, but we have a really good approximation for both of them, and in particular the recovered $e$ has a lot of useless $0$ in the decimal places, so we can just round the number to the nearest integer, and then use the fact that, because our approximation of $||y||$ is more precise than the one of $||x||$ ($y$ has smaller components), we can use the newly recovered $e$ and $||y||$ to compute the effective norms of $x$.

At this point, we can assume the original $x_i$ and $y_i$ to be integers, and now use binary search to look for the angle of rotation $\theta$.

Time to factor, but how?⌗

Notice, first of all, that $ x_1^2 + ey_1^2 = pq$ implies that $ (x_1 + \sqrt{-e}y_1) (x_1 - \sqrt{-e}y_1) = pq$ over $\mathbb{Z}(\sqrt{-e})$.

If we could find a homomorphism $f : \mathbb{Z}(\sqrt{-e}) \rightarrow \mathbb{Z}/n\mathbb{Z}$, such that $f(a) = a, \space \forall a \in \mathbb{Z}/n\mathbb{Z}$, we would be able to map homorphically the zero divisors of $n$ to themselves, ending up with

$$ f((x_1 + \sqrt{-e}y_1) (x_1 - \sqrt{-e}y_1)) = f(x_1 + \sqrt{-e}y_1) f(x_1 - \sqrt{-e}y_1) =f(pq) = f(p)f(q) \equiv 0 \bmod n $$

This way, we can hope that the distinct terms in LHS side map distinctly to the zero divisors in the RHS (and that is basically the core idea behind NFS). This helps us find said homomorphism on page 5.

Notice that $x_1^2 + ey_1^2 = n$ means that $(\frac{x_1}{y_1})^{2} + e \equiv 0 \bmod n$, so our choice of $m$ will be exactly $(\frac{x_1}{y_1}) \bmod n$.

#really garbage code incoming
n = 990853953648382437503731888872568785013804329239290721076418541795771569507440261620612308640652961121590348037236708702361580700250705591203587939980126323233833431892076634892318387020242015741789265095380967467201291693288654956012435416445991341222221539511583706970342630678909437274145759598920314784293470918464283814408418704426938549136143925649863711450268227592032494660523680280136089617838412326902639568680941504799777445608524961048789627301462833
c = 312168688094168684887530746663711142224819184527420449851136749248641895825646649162310024737395663075921549510262779965673286770730468773215063305158197748549937395602308558217528064655976647148323981103647078862713773074121667862786737690376212246588956833193632937835958166526006128435536115531865213269197137648990987207140262543956087199861542889002996727146832659889656384027201202873352819689303456895088190857667281342371263570535523695457095802010885279
rot_x = (9.93659400123277470926327676478883140697376509010297766512845139881487348637477791719517951397052010880811619509960535668814993293095146708957649613776125686226138447162258666762024346093786649228730054881453449071976210130217897905782845690384638460560301964009359233596889465133986468021963885911072779457835979983964294586954038412718305000570678333513135467257498071686562749882446495426493483289204e230, -1.20540611958254673086539287012513674064476659427085664430224592760592531301348857885707154893714440960111029743010026152632150988429192286517249118913535366887447596463819555191858702861383725310592687577510708180057642425944345656558038998574368521689142109798891989865473206201635908814994474491537093810680632691594902962470061189337645818851446622588020765058461348047229165216450857822980873846637e230)
rot_y = (9.02899744041999015549480362358897037217795303901085937071039171882835297563545959015336648016772002396355451308252077767567617065937943765701645833054147976124287566465577491039263554806622908070370269238064956822205986576949383035741108310668397305286076364909407660314991847716094610949669608550117248147017329449889127749721988228613503029640191269319151291514601769696635252288607881829734506023770e191, 2.82245306887391321716872765000993510002376761684498801971981175059452895101888694909625866715259620501905532121092041448909218372087306882364769769589919830746245167403566884491547911250261820661981772195356239940907493773024918284094309809964348965190219508641693641202225028173892050377939993484981988687903270349415531065381420872722271855270893103191849754016799925873189392548972340802542077635974e192)

F = RealField(1337)
rot_x = vector(F, rot_x)
rot_y = vector(F, rot_y)
nxsq_old =rot_x * rot_x
nysq_old =rot_y * rot_y
e = (2*n - nxsq_old) / nysq_old
e = ZZ(QQ(e.numerical_approx(prec  = 700, digits = 400)))
nysq = ZZ(nysq_old.round())
nxsq = ZZ(2*n - e*nysq)

lb = -pi/2
up = pi/2

for j in range(20):
    theta = (lb+up) / 2
    R = matrix(F, [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]])
    Rinv = R^-1
    newx = Rinv * vector(rot_x)
    newy = Rinv * vector(rot_y)
    tmpx1 = round(newx[0])
    tmpx2 = round(newx[1])
    tmpy1 = round(newy[0])
    tmpy2 = round(newy[1])
    if (tmpx2 < 0 or tmpy2 < 0):
        up = theta
    elif (tmpx1 < 0 or tmpy1 < 0):
        lb = theta
        
    elif (tmpx1.nbits() >768 or tmpy1.nbits()>640):
        up = theta
        
    elif (tmpx2.nbits() >768 or tmpy2.nbits()>640):
        lb = theta
        
for j in range(8000):
    theta = (lb+up) / 2
    R = matrix(F, [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]])
    Rinv = R^-1
    newx = Rinv * vector(rot_x)
    newy = Rinv * vector(rot_y)
    tmpx1 = round(newx[0])
    tmpx2 = round(newx[1])
    tmpy1 = round(newy[0])
    tmpy2 = round(newy[1])
    lhs1 = tmpx1^2 + e*tmpy1^2
    lhs2 = tmpx2^2 + e*tmpy2^2
    if (lhs1 > n or lhs2 <n):
        up = theta
    else:
        lb = theta
    if (abs(lhs1 - n) == 0 and abs(lhs2 - n)==0):
        print("found")
        x1 = tmpx1
        x2 = tmpx2
        y1 = tmpy1
        y2 = tmpy2
        break
    

assert x1.nbits() <= 768 and x2.nbits() <= 768
assert y1.nbits() <= 640 and y2.nbits() <= 640
assert x1 ** 2 + e * y1 ** 2 == n
assert x2 ** 2 + e * y2 ** 2 == n
assert x1^2 + x2^2 ==nxsq
assert y1^2 + y2^2 ==nysq

tmp = int(x1 * pow(y1, -1, n))
assert((tmp^2 + e) % n == 0)
#map: f(a + var *b) = a + tmp*b
pfake =(x2 + tmp*y2) % n
qfake =(x2 - tmp*y2) % n
p = gcd(p, n)
q = n//p
phi = (p-1)*(q-1)
d = pow(e, -1, phi)
m = pow(c, d, n)
m = int(m) ^^ x1 ^^ x2 ^^ y1 ^^ y2
from Crypto.Util.number import long_to_bytes
print(long_to_bytes(m))