# Barak
Le challenge revient à résoudre une instance de ECDLP sur une courbe Hessienne.

La courbe est d'équation
$$
(E): x^3 + y^3 + c - dxy = 0
$$
avec
```python
p = 73997272456239171124655017039956026551127725934222347
d = 68212800478915688445169020404812347140341674954375635
c = 1
```

Note, ce n'est techniquement **pas** une courbe elliptique à cause du $y^3$.

Pour les curieux, voici à quoi ressemble une courbe Hessienne:
C'est très similaire au Folium de Descartes, mais c'est pas exactement la même courbe.

L'addition de points est similaire à l'addition classique pour les courbes elliptiques.
On trace la droite passant par les deux points, on trouve la 3e intersection et
on effectue une symétrie axiale, cette fois-ci autour de $(\Delta): y=x$.

Il nous est donné deux points $P$, $Q$ de sorte que $P = Q\cdot m$ pour un certain $m$.
L'objectif est de retrouver $m$, qui encode la solution du challenge.

```
P = (71451574057642615329217496196104648829170714086074852,
69505051165402823276701818777271117086632959198597714)
Q = (40867727924496334272422180051448163594354522440089644,
56052452825146620306694006054673427761687498088402245)
```

## Trouver l'ordre du point P

Il est possible de construire une équivalence birationnelle entre $(E)$ et une
courbe elliptique[1].
Le script suivant illustre la construction de cette courbe elliptique:

```python=
p = 73997272456239171124655017039956026551127725934222347
G = GF(p)
d = G(68212800478915688445169020404812347140341674954375635)
D = d/3
c = 1
# Hessian Curve:
# x³ + y³ + 1 = 3Dxy

# Weierstrass equivalent
a = -27*D*(D^3 + 8)
b = 54*(D^6 - 20*D^3 - 8)
E2 = EllipticCurve(G, [a, b])
print(E2)
# Elliptic Curve defined by y^2 = x^3 + 44905983883632632311912975168565494049729462391119290*x + 4053170785171018449128623853386306889464200866918538 over Finite Field of size 73997272456239171124655017039956026551127725934222347
```

Les deux courbes étant équivalentes, cela nous assure qu'elles possèdent le même
ordre, que nous pouvons calculer grâce à sage:

```python=
order = E2.order()
# 73997272456239171124655016995459084401465136460086688
```

Par théorème de Lagrange, l'ordre de $P$ doit obligatoirement diviser l'ordre de $E$,
donc l'ordre de $E2$.
Son ordre va donc être le plus petit diviseur $k$ de `order` tel que $P\cdot k = O$.

Le script suivant implémente les divers opérateurs sur la courbe Hessienne,
puis calcule l'ordre $n$ du point $P$:

```python=
class HessianPoint:
def __init__(self, x, y):
self.x = G(x)
self.y = G(y)

assert (x, y) == (0, 0) or x^3 + y^3 + c - d*x*y == 0

def __eq__(P, Q):
return (P.x, P.y) == (Q.x, Q.y)

def __add__(P, Q):
if (P.x, P.y) == (0, 0):
return Q
if (Q.x, Q.y) == (0, 0):
return P
if (P.x, P.y) == (Q.y, Q.x):
return HessianPoint(0, 0)

x1, y1 = P.x, P.y
x2, y2 = Q.x, Q.y
if P == Q:
x3 = y1 * (c - x1^3) / (x1^3 - y1^3)
y3 = x1 * (y1^3 - c) / (x1^3 - y1^3)
else:
x3 = (y1^2*x2 - y2^2*x1) / (x2*y2 - x1*y1)
y3 = (x1^2*y2 - x2^2*y1) / (x2*y2 - x1*y1)
return HessianPoint(x3, y3)

def __mul__(P, m):
if (P.x, P.y) == (0, 0):
return P
R = HessianPoint(0, 0)
while m != 0:
if m & 1:
R = P + R
m = m >> 1
if m != 0:
P = P + P
return R

def __neg__(self):
return HessianPoint(self.y, self.x)

def __rmul__(P, m):
return P*m

def __repr__(self):
return f"({self.x}, {self.y})"

def __hash__(self):
return hash((self.x, self.y))

P = HessianPoint(71451574057642615329217496196104648829170714086074852, 69505051165402823276701818777271117086632959198597714)
Q = HessianPoint(40867727924496334272422180051448163594354522440089644, 56052452825146620306694006054673427761687498088402245)
O = HessianPoint(0, 0)

for k in order.divisors():
if P * k == O:
break
n = k
# 3083219685676632130193959041477461850061047352503612
```

On remarque que l'ordre de $P$ se factorise simplement:
$$
n = 2^2 × 3^2 × 17 × 2341 × 23497 × 500369 × 5 867327 × 33 510311 × 13824 276503 × 67342 255597
$$

Le plus grand facteur premier de $n$ étant relativement petit, on peut se permettre de résoudre le logarithme discret avec une variante simplifiée de l'algorithme de Pohlig-Hellman.

## Résoudre le Logarithme Discret

Pohlig-Hellman nous permet de simplifier ECDLP en plusieurs plus petites instances pour chaque facteur de $n$ que l'on peut par la suite combiner grâce à un CRT.
Voici une variante simplifiée de l'algorithme de Pohlig-Hellman qui nous est suffisant au vu des contraintes:

Pour chaque facteur $p_i^{e_i}$ de $n$, on cherche à trouver $k_i \equiv m \mod p_i^{e_i}$.
Pour ce faire, on considère les deux points $P' = P \cdot \dfrac n {p_i^{e_i}}$ et $Q' = Q \cdot \dfrac n {p_i^{e_i}}$.
Étant donné que $P' \cdot p_i^{e_i} = P \cdot n$, il est clair que l'ordre de $P'$ est de $p_i^{e_i}$.

Ainsi, on peut résoudre $P' \cdot k_i = Q'$ et obtenir la valeur de $m \mod p_i^{e_i}$.

Afin de résoudre cette plus petite instance, j'utilise l'algorithme baby-step giant-step,
de complexité temporelle $O(\sqrt{p_i^{e_i}})$.

```python=
def baby_step_giant_step(P, Q, order):
"""
Résoud P·x = Q pour x entre 0 et order-1
"""

m = isqrt(order)
# x = am + b
# ⇒ Q = P*(am + b)
# ⇒ P*(am) = Q+P*(-b)

memory = {}
# On y stocke toutes les valeurs possibles de P*(am),
# associées à leur valeur de `a` respective.
M = O
S = P*m # giant step

for a in range(m):
# M = P*(am)
memory[M] = a
M += S

b = 0
M = Q
S = -P # baby step

while M not in memory:
# M = Q+P*(-b)
b += 1
M += S

a = memory[M]
x = a*m + b
assert P * x == Q
return x
```

Et ici se trouve l'algorithme complet en action:

```python=
factors = n.factor()
p = [factor[0] for factor in factors]
e = [factor[1] for factor in factors]

crt_a = []
crt_b = []
for i in range(len(p)):
pi = p[i]
ei = e[i]
P0 = (n // (pi^ei)) * P
Q0 = (n // (pi^ei)) * Q
k = baby_step_giant_step(P0, Q0, pi^ei)
crt_a.append(k)
crt_b.append(pi^ei)
print(f"m = {k} mod {pi^ei}")
# m = 1 mod 4
# m = 1 mod 9
# m = 13 mod 17
# m = 1489 mod 2341
# m = 18879 mod 23497
# m = 339241 mod 500369
# m = 4831116 mod 5867327
# m = 31953114 mod 33510311
# m = 10804270204 mod 13824276503
# m = 57315765050 mod 67342255597

m = CRT(crt_a, crt_b)
# 1780694557271320552511299360138314441283923223949197
assert P * m == Q
```

## Retrouver le flag

Dernière subtilité, il existe plusieurs valeurs de $m$ valides tel que $P \cdot m = Q$.
En effet $P \cdot (m + \lambda n) = P \cdot m + P \cdot n \cdot \lambda = Q + O \cdot \lambda = Q$ pour toute valeur de $\lambda$ entière.
La seule restriction donnée dans le fichier d'entrée est que $m < p$, alors on teste
toutes les valeurs de $m$ valides jusqu'à trouver le flag:

```python=
p = 73997272456239171124655017039956026551127725934222347
while m < p:
flag = long_to_bytes(m)
try:
print('CCTF{' + flag.decode() + '}')
except ValueError:
pass
m += n

# CCTF{_hE5S!4n_f0rM_0F_3CC!!}
```

## Références

[1] https://link.springer.com/content/pdf/10.1007/3-540-44709-1_33.pdf