# NOT Mordell Primes - zer0pts CTF 2021
###### tags: `zer0pts CTF 2021` `crypto`

## overview

This is the twisted version of Mordell Primes in Union CTF 2021. The difference are parameters and the Curve is on the $\mathbb{F}_p$.

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

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

P = E(..., ...)
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

Since $P + Q = R$, we can make mulrivariate modular congruence about additive law of Elliptic Curve. It could be solvable by using groebner basis. Anyway please see the source code below.

## exploit

p = 1340087758828337670882568669105503854580062990023098389939480481911369583331761787997658244256484001548844170708665954112184084618798862963196496982555437
a = 579520793019813043301962853514018350546879414461429763549854888716118598035188444377014576936707774693812675902636517550724517748438949781821496069882991
b = 836220137554822093644088628952860378557613203386578961376911335316893863894144398100294478276291588679916831302112802755833376978399222338686452519253077


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

P = EC(235503147087035059574346244036877145714699241108084268073582256595578731905840469662408045507147931695613045864587510088815875998657156709252830833800517, 288043725515598578535261835984701538316089236142610727440633081936390449983737273523495496717793593528955017368412573825832528719973368667183359113795046)

import ast
with open("output.txt") as f:
    N = ast.literal_eval(f.readline().strip().split(" = ")[1])
    c = ast.literal_eval(f.readline().strip().split(" = ")[1])

# find l
PR.<x2, y2, l> = PolynomialRing(GF(p))
x1, y1 = P.xy()

x3 = l^2 - (x1 + x2)
y3 = l*(x1 - x3) - y1

I = Ideal([x2*x3 - N, (x2^3 + a*x2 + b) - (y2^2), (x3^3 + a*x3 + b) - (y3^2), l*(x2 - x1) - (y2 - y1)])
B = I.groebner_basis()
PR.<x> = PolynomialRing(GF(p))
f = x^2 + int(B[5][l])*x + int(B[5].constant_coefficient())

ls = [r[0]^r[1] for r in f.roots()]

# find y2
for l in ls:
    print("l = {}".format(l))
    PR.<x2, y2> = PolynomialRing(GF(p))
    x3 = l^2 - (x1 + x2)
    y3 = l*(x1 - x3) - y1

    I = Ideal([x2*x3 - N, (x2^3 + a*x2 + b) - (y2^2), (x3^3 + a*x3 + b) - (y3^2), l*(x2 - x1) - (y2 - y1)])

    B = I.groebner_basis()
    f = x^2 + int(B[0][y2])*x + int(B[0].constant_coefficient())
    ys = [r[0]^r[1] for r in f.roots()]

    for y2 in ys:
        print("y2 = {}".format(y2))
        x2 = -B[1](y2=y2).constant_coefficient()
        if N % int(x2) != 0:
            continue

        p = N // int(x2)
        q = N // p
        d = inverse_mod(0x10001, (p-1)*(q-1))
        m = pow(c, d, N)

        print(bytes.fromhex(hex(m)[2:]))