# Writeups HackTM CTF 2023

## d-phi-enc

In this challenge, there are two polynomials who share `phi` as root.

Based on the following information:

e*d = phi*2 + 1
e = 3

We define the following polynomials:

f = (2*x+1)^3 - c
g = x^3 - enc_phi

Ideally, we would like to compute the gcd based on these values but sage is not able to compute it over a field. Therefore, we used a groebner basis to recover the root of the polynomials.

from sage.rings.polynomial.toy_d_basis import d_basis
from Crypto.Util.number import long_to_bytes
n = 24476383567792760737445809443492789639532562013922247811020136923589010741644222420227206374197451638950771413340924096340837752043249937740661704552394497914758536695641625358888570907798672682231978378863166006326676708689766394246962358644899609302315269836924417613853084331305979037961661767481870702409724154783024602585993523452019004639755830872907936352210725695418551084182173371461071253191795891364697373409661909944972555863676405650352874457152520233049140800885827642997470620526948414532553390007363221770832301261733085022095468538192372251696747049088035108525038449982810535032819511871880097702167
enc_d = 23851971033205169724442925873736356542293022048328010529601922038597156073052741135967263406916098353904000351147783737673489182435902916159670398843992581022424040234578709904403027939686144718982884200573860698818686908312301218022582288691503272265090891919878763225922888973146019154932207221041956907361037238034826284737842344007626825211682868274941550017877866773242511532247005459314727939294024278155232050689062951137001487973659259356715242237299506824804517181218221923331473121877871094364766799442907255801213557820110837044140390668415470724167526835848871056818034641517677763554906855446709546993374
enc_phi = 3988439673093122433640268099760031932750589560901017694612294237734994528445711289776522094320029720250901589476622749396945875113134575148954745649956408698129211447217738399970996146231987508863215840103938468351716403487636203224224211948248426979344488189039912815110421219060901595845157989550626732212856972549465190609710288441075239289727079931558808667820980978069512061297536414547224423337930529183537834934423347408747058506318052591007082711258005394876388007279867425728777595263973387697391413008399180495885227570437439156801767814674612719688588210328293559385199717899996385433488332567823928840559
enc_flag = 24033688910716813631334059349597835978066437874275978149197947048266360284414281504254842680128144566593025304122689062491362078754654845221441355173479792783568043865858117683452266200159044180325485093879621270026569149364489793568633147270150444227384468763682612472279672856584861388549164193349969030657929104643396225271183660397476206979899360949458826408961911095994102002214251057409490674577323972717947269749817048145947578717519514253771112820567828846282185208033831611286468127988373756949337813132960947907670681901742312384117809682232325292812758263309998505244566881893895088185810009313758025764867
# create two polynomial that share a root
P.<x> = PolynomialRing(IntegerRing(), 1, order='lex')
c = int(enc_d) * 27
f = (2*x+1)^3 - c
g = x^3 - enc_phi
I = ideal(f,g)
gb = d_basis(I)
R = P.change_ring(IntegerModRing(n))
gb = [R(f) for f in gb if R(f)]
recovered_phi = (-gb[0][0]) % n
recovered_d = pow(e, -1, recovered_phi)
long_to_bytes(int(pow(enc_flag, recovered_d,n)))

## kaitenzushi

The first step was to recover `e`. We used the fact that the rotation of the vectors does not affect the norm. Therefore, we used the provided equations to retrieve the public exponent:

a1, a2 = x
b1, b2 = y
e = (2*n - a1^2-a2^2) / (b1^2+b2^2)
e = e.integer_part()

Based on the challenge the initial vectors `x` and `y` were rotated by a random angle. We came up with an algorithm to compute the `theta` value and compute the initial vectors and yes it was ugly.

Once the vectors and the public exponent were known we looked at the following checks:

assert x1 ** 2 + e * y1 ** 2 == p * q
assert x2 ** 2 + e * y2 ** 2 == p * q

Dixon's factorization technique came to our minds after hours of trying to recover the private key, since the idea behind this method fits perfectly with the equations above.

from math import gcd
from Crypto.Util.number import bytes_to_long, isPrime
n = 990853953648382437503731888872568785013804329239290721076418541795771569507440261620612308640652961121590348037236708702361580700250705591203587939980126323233833431892076634892318387020242015741789265095380967467201291693288654956012435416445991341222221539511583706970342630678909437274145759598920314784293470918464283814408418704426938549136143925649863711450268227592032494660523680280136089617838412326902639568680941504799777445608524961048789627301462833
c = 312168688094168684887530746663711142224819184527420449851136749248641895825646649162310024737395663075921549510262779965673286770730468773215063305158197748549937395602308558217528064655976647148323981103647078862713773074121667862786737690376212246588956833193632937835958166526006128435536115531865213269197137648990987207140262543956087199861542889002996727146832659889656384027201202873352819689303456895088190857667281342371263570535523695457095802010885279
x = (9.93659400123277470926327676478883140697376509010297766512845139881487348637477791719517951397052010880811619509960535668814993293095146708957649613776125686226138447162258666762024346093786649228730054881453449071976210130217897905782845690384638460560301964009359233596889465133986468021963885911072779457835979983964294586954038412718305000570678333513135467257498071686562749882446495426493483289204e230, -1.20540611958254673086539287012513674064476659427085664430224592760592531301348857885707154893714440960111029743010026152632150988429192286517249118913535366887447596463819555191858702861383725310592687577510708180057642425944345656558038998574368521689142109798891989865473206201635908814994474491537093810680632691594902962470061189337645818851446622588020765058461348047229165216450857822980873846637e230)
y = (9.02899744041999015549480362358897037217795303901085937071039171882835297563545959015336648016772002396355451308252077767567617065937943765701645833054147976124287566465577491039263554806622908070370269238064956822205986576949383035741108310668397305286076364909407660314991847716094610949669608550117248147017329449889127749721988228613503029640191269319151291514601769696635252288607881829734506023770e191, 2.82245306887391321716872765000993510002376761684498801971981175059452895101888694909625866715259620501905532121092041448909218372087306882364769769589919830746245167403566884491547911250261820661981772195356239940907493773024918284094309809964348965190219508641693641202225028173892050377939993484981988687903270349415531065381420872722271855270893103191849754016799925873189392548972340802542077635974e192)

a1, a2 = x
b1, b2 = y

# rotation does not change the norm
# otherwise you can derive the same equations manually
e = (2*n - a1^2-a2^2) / (b1^2+b2^2)
e = e.integer_part()
assert isPrime(e) and int(e).bit_length() == 256
#e = 111578009802636409437123757591617048189760145423552421418627338749835916561801


F = RealField(1337)

a1, a2 = x
b1, b2 = y

t = var("t")
ff = a1^2*cos(t)^2 + a1*a2*sin(2*t) + a2^2*sin(t)^2 + e*(b1^2*cos(t)^2 + b1*b2*sin(2*t)+b2^2*sin(t)^2) - n

f = lambda param : abs(F(ff.subs({t:param})))

# worst approx algorithm pls dont copy
def approx(delta, start, N):
    best_v = 10e1000
    best_x0 = None

    for x0 in range(N):
        x0 = start + delta / N * x0
        v = F(f(x0))
        if abs(v) < abs(best_v):
            best_v = v
            best_x0 = x0

    return best_v, best_x0
    
N = 1000
delta = 2*pi
start = -pi

for j in range(1000):
    v_n, x_n = approx(delta, start, N)
    delta /= 2
    start = x_n - delta / 2

    if j % 50 == 0:
        print(j)
    if v_n == 0:
        break

# iterate possible solutions for theta
L = 4
for offset in range(L):
    theta_r = x_n + (2*pi) / L * offset
    R_ = matrix(F, [[cos(theta_r), sin(theta_r)], [-sin(theta_r), cos(theta_r)]])
    x = vector(F, x)
    y = vector(F, y)
    x_recovered = R_ * x
    y_recovered = R_ * y

    x1 = x_recovered[0].round()
    x2 = x_recovered[1].round()
    y1 = y_recovered[0].round()
    y2 = y_recovered[1].round()

    if x1 ** 2 + e * y1 ** 2 - n == 0 and x2 ** 2 + e * y2 ** 2 - n == 0:
        if x1 > 0 and x2 > 0 and y1 > 0 and y2 > 0:
            print("\n\n", theta_r.n())
            print("x1 =", x1)
            print("x2 =", x2)
            print("y1 =", y1)
            print("y2 =", y2)

#e = 111578009802636409437123757591617048189760145423552421418627338749835916561801
#x1 = 993315378106395196440156892634615357425859001976376351903878161126954317590016249318316631584063366449446002974804447367756266228508159317926113473123770241598131922105753478630709094061327843793983555725542453353312556415777678937
#x2 = 123343431936894440973263647479974540141395074556779828339916509613682879668610901423506961118285523166037774054833601787794419590891163752205158573276826154790166536984681500991748749778629881670438838666011425669518792357094873553
#y1 = 193518098174342694414424160720807163740044134017573004218248685165604434384710484681124817651698709818703976889508767807895216618103609127904817977547152172876909535027087606807328610207963608
#y2 = 2957028917590401838272414886210261099554152128524012256631787151968768935090286908219944634008304129914083074684507666539700290047827545862670465906725813971398170535104589598065683927537059268

# use dixon factorization
from Crypto.Util.number import long_to_bytes
p = gcd(y1*x2 - y2*x1, n)
print(p)
# 957509848415776008506125961998120495161250346184055094697245571121876444575553394581756735245207167681344755095903616730328731358607257251854603846193989936802222147961302618645021044609662945352893811478461448918625795339911124621
q = n // p

d = pow(e, -1, (p - 1) * (q - 1))
m = int(pow(c, d, n)) ^^ x1 ^^ y1 ^^ x2 ^^ y2
print(long_to_bytes(m))