# NOT Mordell primes
## task.sage

from Crypto.Util.number import bytes_to_long
from secrets import k, FLAG


p = 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
a = 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
b = 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177

E = EllipticCurve(GF(p),[a,b])

P = E(11283606203023552880751516189906896934892241360923251780689387054183187410315259518723242477593131979010442607035913952477781391707487688691661703618439980, 12748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)
Q = k*P
R = (k+1)*P

p = int(Q[0])
q = int(R[0])

assert is_prime(p)
assert is_prime(q)

e = 0x10001
N = p*q
m = bytes_to_long(FLAG)
c = pow(m,e,N)

print(f'N = {N}')
print(f'c = {c}')

# Solution

Our goal is to figure out the value of $p$ and $q$ so we can decrypt the ciphertext $c$ and get the flag.

We have the following relations

$$p = Q_x\\q = R_x\\N = p*q\\Q = k*P\\R = (k + 1)*P$$

We don't know the value of $k$, but we know the value of $P$. Notice that

$$R = Q + P$$

This means we can take the ECC point addition formula and generate new equations.

From https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication#Point_addition, we have the following equations

$$\lambda = \frac{y_Q-y_P}{x_Q-x_P}
\\x_R = \lambda^2 - x_P-x_Q
\\y_R = \lambda(x_P-x_R)-y_P$$

Note that the last 2 equations are over the field of the curve, `GF(p)`, or $\bmod p$.

Our goal is to get a univariate polynomial equation over `GF(p)` in either $x_R$ or $x_Q$, which we can then use sage to find the roots of. As long as we recover $x_R$ or $x_Q$, we can get the other value by using $N$.

There are 2 main ways to do this. Pure algebra, or using Groebner basis/resultants.

## Method 1 - Algebra

Let's focus on the equation for $x_R$. We have
$$x_R = \lambda^2 - x_P-x_Q$$

We substitute $\lambda$ and get
$$x_R = \left(\frac{y_Q-y_P}{x_Q-x_P}\right)^2 - x_P-x_Q\\
x_R * (x_Q-x_P)^2 = \left(y_Q-y_P\right)^2 - (x_P+x_Q)*(x_Q-x_P)^2$$

We can get rid of $x_R$ by multiplying both sides by $x_Q$:
$$N * (x_Q-x_P)^2 = x_Q*\left(\left(y_Q-y_P\right)^2 - (x_P+x_Q)*(x_Q-x_P)^2\right)$$

Now in this equation we know the value of $N, x_P,$ and $y_P$, but we still have 2 unknowns $x_Q$ and $y_Q$. Let's get rid of $y_Q$ by using the equation for the elliptic curve and substitute $x_Q$ instead:

$$y_Q^2 = \left(x_Q^3 + ax_Q + b\right) \bmod p$$

Now we have:

$$N * (x_Q-x_P)^2 = x_Q*\left(\left(y_Q^2-2y_Qy_P+y_P^2\right) - \left(x_Q^3-x_Q^2x_P-x_Qx_P^2+x_P^3\right)\right)\\
N * (x_Q-x_P)^2 = x_Q*\left(\left(x_Q^3 + ax_Q + b-2y_Qy_P+y_P^2\right) - \left(x_Q^3-x_Q^2x_P-x_Qx_P^2+x_P^3\right)\right)$$

Cleaning up a removing some like terms:
$$N * (x_Q-x_P)^2 = x_Q*\left(\left(ax_Q + b-2y_Qy_P+y_P^2\right) + \left(x_Q^2x_P+x_Qx_P^2-x_P^3\right)\right)$$

We still have $y_Q$, so we will put it alone on one side of the equation and square both sides to get rid of it.

$$N * (x_Q-x_P)^2 = x_Q*\left(ax_Q + b-2y_Qy_P+y_P^2 + x_Q^2x_P+x_Qx_P^2-x_P^3\right)\\
N * (x_Q-x_P)^2 - ax_Q^2 - bx_Q -x_Qy_P^2 -x_Q^3x_P-x_Q^2x_P^2+x_Qx_P^3= x_Q*\left(-2y_Qy_P\right)\\
\left(x_Q*\left(-2y_Qy_P\right)\right)^2 = \left(N * (x_Q-x_P)^2 - ax_Q^2 - bx_Q -x_Qy_P^2 -x_Q^3x_P-x_Q^2x_P^2+x_Qx_P^3\right)^2\\
4x_Q^2y_Q^2y_P^2 = \left(N * (x_Q-x_P)^2 - ax_Q^2 - bx_Q -x_Qy_P^2 -x_Q^3x_P-x_Q^2x_P^2+x_Qx_P^3\right)^2\\
4x_Q^2y_P^2\left(x_Q^3 + ax_Q + b\right) = \left(N * (x_Q-x_P)^2 - ax_Q^2 - bx_Q -x_Qy_P^2 -x_Q^3x_P-x_Q^2x_P^2+x_Qx_P^3\right)^2$$

It looks ugly, but we finally only have an equation in just $x_Q$, we can clean it up with SageMath and try to get the roots of this univariate polynomial.

We also need to use the `pycrypto` python module. With sage we can easily install this package by doing

sage --python -m pip install pycrypto

### Solution in SageMath

from Crypto.Util.number import long_to_bytes

p = 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
a = 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
b = 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177

N = 22607234899418506929126001268361871457071114354768385952661316782742548112938224795906631400222949082488044126564531809419277303594848211922000498018284382244900831520857366772119155202621331079644609558409672584261968029536525583401488106146231216232578818115404806474812984250682928141729397248414221861387
e = 0x10001
c = 15850849981973267982600456876579257471708532525108633915715902825196241000151529259632177065183069032967782114646012018721535909022877307131272587379284451827627191021621449090672315265556221217089055578013603281682705976215360078119427612168005716370941190233189775697324558168779779919848728188151630185987

E = EllipticCurve(GF(p),[a,b])
P = E(11283606203023552880751516189906896934892241360923251780689387054183187410315259518723242477593131979010442607035913952477781391707487688691661703618439980, 12748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)

xp = int(P[0])
yp = int(P[1])
P.<xq> = PolynomialRing(Zmod(p))

f = 4*xq^2*yp^2 * (xq^3 + a*xq + b) - (N * (xq - xp)^2 - a*xq^2 - b*xq - xq*yp^2 - xq^3*xp - xq^2*xp^2 + xq*xp^3)^2
print(f.roots())

for root, multiplicity in f.roots():
    q_guess = int(root^multiplicity)
    if N % q_guess == 0:
        print("found it")
        q = q_guess
        p = N // q
        phi = (p - 1) * (q - 1)
        d = inverse_mod(e, phi)
        print(long_to_bytes(pow(c, d, N)))
        break

Running the script above, we get the flag

## Method 2: Resultants/Groebner basis

If you don't want to do some nasty algebra, you can basically ask SageMath to simply the multivariate polynomials you have into univariate ones, and then easily get the roots. This is primarily done through the Groebner basis.

First let's list the multivariate polynomial equations we have:

$$x_R * (x_Q-x_P)^2 = \left(y_Q-y_P\right)^2 - (x_P+x_Q)*(x_Q-x_P)^2\\
y_R * (x_Q-x_P) = (y_Q-y_P) * (x_P-x_R) - y_P* (x_Q-x_P)\\
N= x_Q*x_R\\
y_Q^2 = x_Q^3 + a*x_Q + b\\
y_R^2 = x_R^3 + a*x_R + b$$

Now we ask SageMath to attempt to calculate the Groebner basis for these 5 multivariate polynomial equations. Also make sure to get rid of the fractions, since the equations are under $\bmod p$.

### SageMath code to do Groebner basis

from Crypto.Util.number import long_to_bytes

p = 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
a = 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
b = 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177

N = 22607234899418506929126001268361871457071114354768385952661316782742548112938224795906631400222949082488044126564531809419277303594848211922000498018284382244900831520857366772119155202621331079644609558409672584261968029536525583401488106146231216232578818115404806474812984250682928141729397248414221861387
e = 0x10001
c = 15850849981973267982600456876579257471708532525108633915715902825196241000151529259632177065183069032967782114646012018721535909022877307131272587379284451827627191021621449090672315265556221217089055578013603281682705976215360078119427612168005716370941190233189775697324558168779779919848728188151630185987

E = EllipticCurve(GF(p),[a,b])
P = E(11283606203023552880751516189906896934892241360923251780689387054183187410315259518723242477593131979010442607035913952477781391707487688691661703618439980, 12748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)

xp = int(P[0])
yp = int(P[1])
P.<l, xq, yq, xr, yr> = PolynomialRing(Zmod(p))

f1 = xr * (xq - xp)^2 - (yq - yp)^2 + (xp + xq) * (xq - xp)^2
f2 = yr * (xq - xp) - (yq - yp) * (xp - xr) + yp * (xq - xp)
f3 = N - xq * xr
f4 = yq^2 - xq^3 - a*xq - b
f5 = yr^2 - xr^3 - a*xr - b
I = ideal([f1, f2, f3, f4, f5])
B = I.groebner_basis()

for poly in B:
    print(poly)

for poly in B:
    print(poly.monomials())

### Groebner basis output

yr^3 + 5666880893547136480255841922859400992206656889082093956636278576669736114777428697989942658325415146429860357899712285975825246928664028840249043769127428*yr^2 + 2632440716432087076619220864393394255755157490518721713965710726543835465963482081212183644474370295202492965626732841187037457979174409411984621954779106*xq + 4511479125640753138127959318606337708642809213596589531445870113941324454817644765549294465322299952086957508784669390906563133846246213225143408272039664*yq + 10639832660929164541090955908121894517208112163344703417931340239334092280805447579091340814045621299806960527909993987043026674648403935049318794253233598*xr + 8171152544770799679817087082265627968638669066642180045880275643022347334799593867069759934835527405337003238723076005323325239347345840349682198501435215*yr + 2466519724204860491894653707127033493511805739176129685833961625734531873863189369751393937246570554336731325189576186258626721734447039920866953619750893
xq^2 + 11538897321743466555075964145532585570510503241506054041704174646293099946784255444406477889567578817042410114309828904409446536528290998736447760795591100*xq + 10863866088592178737453716108835010998879807001346631266944590059577566867491016719390508617544204226971600483061513162364642980144831495092085507828169005*yq + 1763282894498093488335953418281655272275522879423943691312771766626221157294832805869600883835305784318187396642888426756906943957420064166607272774539093*xr + 12744677188652596235574621834656759009258995058736558873594026398868214999344249550883590950148766601460853277677647924477255434765046360498510881442258638
xq*yq + 2436284434476330385602943631629821197481289040917823607897185142436682678048776388086041263086096725572963986270897709722355403182724342518314531325712902*yr^2 + 3021174139288789062888963594455389552125618629291202040082945567879799201850397225113848921026114226246226017133162946833033831925643326004396431354867512*xq + 4332508984532358096276024117124698182997617715592662919049610841322434443431939912210454480301565895919301539115099224451871653456545542785353608161617338*yq + 7794688391180577871249132186531342281655311127456628651408829298111600746074623228680343542078866967234732756411184158217339071088830705241544490830082359*xr + 4120569042532314092373969326695036974065543646190811796008239594330650342562554874165566794273760076966398646795633992601385786694099637995658732186866739
yq^2 + 139257580323997605812634111158630435068017689896720524952175724115726214305770950255473577084368347352935025758302454868570136853606632301506372662127044*yr^2 + 11907778822901066106686786263966820099894771372537114296612050103459578513481420522932082569396427305874670797283921242045170527560308080002765343557025226*xq + 9476608477470607409826079205773769577115458656582746591179890919195393627314107333112011752359773792403582109283837613236610005771590584298142618433912761*yq + 12493197365497218485004231811350293122366183292480621942604148492363954980775983125021297004181568625795355760781519106683131222638168934873450287807990968*xr + 2036484082992517080502595097122186890688091533285652826520570022986809148578960714840398655076671736421242015002829617877011951398611198492029785343983595
xq*xr + 7025565589500174699832226767203371757172750521753105773236856247759706592273000015147169142705661127188041109902584007538678387992530274133947946085439788
yq*xr + 10610604663045315983484525976558731009686475199429371864104973678372725889561315936506802098342341037755666017407904669512332932482183410339954445067266171*yr^2 + 8826688045919728341383137300251813637215622518162432576438962177331074438270581625103423024368627095066015099397642976807498548437196223903181324602591065*xq + 657290005737162821377545409561744325277102466670026106363138367573108548879181450629746774336351911563482185676929609292323487404549585278684335968579768*yq + 5550227053285295054306572863923636997670424341766808862942089962231813159810042413597160224352578205301920826835269190014361657855518966771931974017098796*xr + 2851584918529830843941922772237243951560352312999007532880981883893555252212792135356214228213304432712265627280437418824325254982291346315628758717573659
xr^2 + 9778556302085837812285363767693571652268098764204978009721248565891941472907939713994091178432053713428500561887004893097816147525060229267772037640284368*yr^2 + 6334853071620470256950806276701356766351232438159601283382064127687770428135640269717024414980196437886276867412920922041203533731838747800462251292864230*xq + 10058403592379415174003562269380057245984960688605852683224479273139023360256475743077832372036235231520322077823715136169680482768485711094064602941450452*yq + 7531137243325883426427805295214070578017938301158893915686439907973107382661453520139505428216151918944337216688334324617489983492268153606747673730452337*xr + 12065229246544804480764969691476407174043486207134216196820155002622505658701133325004961738765665676100113737119343532161718910324316627311453246040882295
xq*yr + 2436284434476330385602943631629821197481289040917823607897185142436682678048776388086041263086096725572963986270897709722355403182724342518314531325712902*yr^2 + 3922174704657691471236096865670311497794170493308520853214436203944328791064937381804759932056803259054591324713508433430177393948269609799880163336185926*xq + 1105992888760930666958408008719910946998420412753917584949633399053112608415651355239854109498953872754705210965958817464583456552870478887922936805959325*yq + 7794688391180577871249132186531342281655311127456628651408829298111600746074623228680343542078866967234732756411184158217339071088830705241544490830082359*xr + 1763282894498093488335953418281655272275522879423943691312771766626221157294832805869600883835305784318187396642888426756906943957420064166607272774539093*yr + 10195304178991815525145546835951308255607411927348187939121176936915853315397300189236629133215133330616364376398364960410363080682616406542640217675405414
yq*yr + 139257580323997605812634111158630435068017689896720524952175724115726214305770950255473577084368347352935025758302454868570136853606632301506372662127044*yr^2 + 1657958018558198809503162468771782950025326844245251546097505835977017223416582886805545299862204460209407446667324724233374176011689131327299659853469886*xq + 23151904625831252132060179093786308778344530371685397761611074703234586627551827840892875325837310931592523442880109256871181630308546033036488544352179*yq + 12493197365497218485004231811350293122366183292480621942604148492363954980775983125021297004181568625795355760781519106683131222638168934873450287807990968*xr + 298026346944226556468235442266427072157971228876242042348760439534005338274573317684660405003007409208023579567650968997012393279441919155207488454202082*yr + 3842252231153082831793796223432621953293879574332695729026501177227899265241326563252139678626287961543511917783840098369490061129102541853908159570778301
xr*yr + 4126212875713162195214270837872526194291390681899273114035983750804762333515641770206242750985217501871140603388568561975334769461453811586934224108243938*yr^2 + 3193686994050199972731031954738786583988716214598915981789360694283843647261350545535144420972285678851462659182546942540960050244653013751470679210020784*xq + 5949398831896570036697135160209936360621833634462386488119726671524052150725491670980100499782964509653264328517363851910684299418718621899569842771412846*yq + 8227141031680986341672289935554971365045765534614590398453963017407893529552824710654543796148579770785615102608542881420919979983359485728957347682252527*xr + 12775053331508383767963595315118228287620060577374384160901357909400018296141059472700020502289317779998160356897429840072884082573478293480477048945850689*yr + 3057985480360776141370967415111665062316223667298995004442571293585104122117989352653786232828790396722441237081391056672314536131251110419264982142804432

Let's list the monomials for each one

[yr^3, yr^2, xq, yq, xr, yr, 1]
[xq^2, xq, yq, xr, 1]
[xq*yq, yr^2, xq, yq, xr, 1]
[yq^2, yr^2, xq, yq, xr, 1]
[xq*xr, 1]
[yq*xr, yr^2, xq, yq, xr, 1]
[xr^2, yr^2, xq, yq, xr, 1]
[xq*yr, yr^2, xq, yq, xr, yr, 1]
[yq*yr, yr^2, xq, yq, xr, yr, 1]
[xr*yr, yr^2, xq, yq, xr, yr, 1]

These polynomials are new polynomials created from the basis of the original 5 polynomials we put in.

What we can do now is take some of the polynomials and use the resultant to return a new polynomial where a variable is removed. We continue this process until we get a univariate polynomial in either $x_Q$ or $x_R$ that we can take the roots of.

The first step is to choose polynomials that we want to calculate the resultant of. I chose polynomials that had a small amount of monomials and similar variables.

The second polynomial has only 3 variables in it ($x_Q$, $y_Q$, and $x_R$), and the 5th one only has 2 ($x_Q$ and $x_R$), so we will first take 2 polynomials with all 4 variables (any 2 will do, doesn't matter which) and remove $y_R$ using the resultant. Then we will get a polynomial with only $x_Q$, $y_Q$, and $x_R$ in it. Then we take this polynomial and the 2nd polynomial from the Groebner basis and calculate the resultant to remove $y_Q$. Then we take that polynomial and the 5th polynomial from the Groebner basis and calculate the resultant to remove either $x_Q$ or $x_R$ (my code did both).

That last resultant will be a univariate polynomial that we then can just call `.roots()` and get our flag. Note that in SageMath the way that we calculate the resultant results in a multivariate polynomial being returned, so we have to call `.univariate_polynomial()` to get a univariate polynomial and then call `.roots()`. Solution code below.

### SageMath code to calculate resultant and get flag

from Crypto.Util.number import long_to_bytes

p = 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
a = 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
b = 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177

N = 22607234899418506929126001268361871457071114354768385952661316782742548112938224795906631400222949082488044126564531809419277303594848211922000498018284382244900831520857366772119155202621331079644609558409672584261968029536525583401488106146231216232578818115404806474812984250682928141729397248414221861387
e = 0x10001
c = 15850849981973267982600456876579257471708532525108633915715902825196241000151529259632177065183069032967782114646012018721535909022877307131272587379284451827627191021621449090672315265556221217089055578013603281682705976215360078119427612168005716370941190233189775697324558168779779919848728188151630185987

E = EllipticCurve(GF(p),[a,b])
P = E(11283606203023552880751516189906896934892241360923251780689387054183187410315259518723242477593131979010442607035913952477781391707487688691661703618439980, 12748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)

xp = int(P[0])
yp = int(P[1])
P.<xq, yq, xr, yr> = PolynomialRing(Zmod(p))

f1 = xr * (xq - xp)^2 - (yq - yp)^2 + (xp + xq) * (xq - xp)^2
f2 = yr * (xq - xp) - (yq - yp) * (xp - xr) + yp * (xq - xp)
f3 = N - xq * xr
f4 = yq^2 - xq^3 - a*xq - b
f5 = yr^2 - xr^3 - a*xr - b
I = ideal([f1, f2, f3, f4, f5])
B = I.groebner_basis()

'''
for poly in B:
    print(poly)

for poly in B:
    print(poly.monomials())
'''

def resultant(p1, p2, var):
    p1 = p1.change_ring(QQ)
    p2 = p2.change_ring(QQ)
    var = var.change_ring(QQ)
    r = p1.resultant(p2, var)
    return r.change_ring(Zmod(p))

#B[4] - only xq and xr
#B[1] - only xq, yq, and xr

# Get rid of yr
h1 = resultant(B[2], B[5], yr) # only xq, yq, and xr remains
h2 = resultant(h1, B[1], yq) # only xq, and xr remains
h3 = resultant(h2, B[4], xq) # only xr remains
h4 = resultant(h2, B[4], xr) # only xq remains

roots = h3.univariate_polynomial().roots()
roots2 = h4.univariate_polynomial().roots()

for root, multiplicity in (roots + roots2):
    q_guess = int(root)
    if N % q_guess == 0:
        print("found it")
        q = q_guess
        p = N // q
        phi = (p - 1) * (q - 1)
        d = inverse_mod(e, phi)
        print(long_to_bytes(pow(c, d, N)))
        break

# Flag
`zer0pts{7h4nk_y0u_j4ck_7h4nk_y0u_cr0wn}`