First of all, I needed to recover the prime `p`.
I used the third output constructed from the point with x-coordinate `hint`; I knew that `hint` would therefore have to be a zero for the polynomial that is the x-part of this output.
Therefore, I could just set `x=hint` and factor the result:

```
y = var("y")
hint = int.from_bytes(b"Inspired by theoremoon's SECCON 2022 Finals Challenge - Hell. Thank you!", "big")
outputs = [(x^2 + 14762123602851553604749022996287576858980220577795651427829269858766434621297346961387874961427459051934768224338447011128244905975068497090840444625419470*x + 8519674750729750620690035589812482119785861876353468044895414394332293279114303071755954851101633319350193436546144692795403444364414318973131157246232656, y + 17770738679279057916355557895675090129563269633432826251932824463003364931275912702916209480950481351904761364290424406482997835483807402182326014818733821*x + 12306668467523337827805393760490897581559948654643366727345701375757143864825442910779617850907143245102792529282031529618639723158417652048624567379151171),
(x^3 + 13441279154284544764330805782065565325543470739559917045273482055514440837785754044182874902421009026981197721504820302867945812937528249594953326223176272*x^2 + 3795282115520834934850220740151212731596814319504043674340537364041453624883995759365119899076774262882230308591629439035308527946872182029742910504122735*x + 3726617245981099594981815385059428688276726297460965450328320328460867196111587736356492934195556032891106446058683147130913147722036293641303193921962091, y + 2103349591221335944593862709600493681857281410337020721978302326614691696399677635217262732543672829811190387220058078405239568477387817550236173432744263*x^2 + 4784247634355946154999459446762911004042472267922959302672838559247991353014786987556174410735592161587023899368989617780068662559773261109676326152316907*x + 2640959823121300693709616791657128464111647959613642856293592234010564318329382577397798309822254798484629398268742247779165733848105319417195858443049412),
(x + 540914350057159927735436910109553086959907629357396102386062800609858847095404691934517983640846787552171946520460064397364544720996858815571065665255645, y + 541917331856005964100090629475512429550322452567752818120774876171019476274441296070275457561095853517207532108745504694853066426720092700847788666013730)]

possible_p = [i[0] for i in factor(outputs[2][0](x=hint)) if i[0].bit_length()==513]
assert len(possible_p) == 1
p = possible_p[0]
print(f"p recovered: {p}")
```

The next challenge was to recover $f_1$ and $f_2`$
Denoting the divisors as $[(u_1,v_1), (u_2,v_2), (u_3,v_3)]$, I knew that $f_1 \equiv v_1^2 \mod u_1$, $f_2 \equiv v_2^2 \mod u_2$, $f_2 \equiv v_3^2 \mod u_3$ by the basic properties of the Mumford representation of divisors. Using CRT, I could recover $f_2$ modulo $\text{lcm}(u_2, u_3)=u_2\cdot u_3$.

Now write $f_1 = f_2 + sx^7 + tx^6$.
By CRT again (using that $u_3$ is coprime to $u_1\cdot u_2$), we can find $S,T,U$ such that:

$$
S \equiv 0 \mod u_0 \\
S \equiv x^7 \mod u_1 \cdot u_2 \\
T \equiv 0 \mod u_0 \\
T \equiv x^6 \mod u_1 \cdot u_2 \\
U \equiv f_1 \mod u_0\\
U \equiv f_2 \mod u_1 \cdot u_2
$$

Then consider $P := U - sS - tT$. Then on one hand, $P \equiv f_1 + 0s + 0t= f_1\mod u_0$, but on the other $P \equiv f_2 - sx^7 - tx^6 = f_1 \mod u_1\cdot u_2$, so it follows that $P \equiv f_1 \mod u_0\cdot u_1 \cdot u_2$.

But because $\deg f_1 = 5 < 6 = \deg (u_0\cdot u_1 \cdot u_2)$, we have $P = f_1$ (in $GF(p)[x]$) already.

Consider now the vectors
$$U\|\|(C,0,0), S\|\|(0,1,0), T\|\|(0,0,1)$$
$$(p,0,0,\ldots,0)\|\|(0,0,0), (0,p,0,\ldots,0)\|\|(0,0,0), \ldots, (0,0,0,\ldots,p)\|\|(0,0,0)$$
in $\mathbb{Z}^9$, where we identify polyonimals (of degree $\leq 5$) with their coefficient vector in $\mathbb{Z}^6$ and $C$ is a large-ish constant (I chose $C=10^{50}$).

Then clearly, $f_1 \|\| (C, -s, -t)$ is an integral linear combination of those vectors.

Fortunately, this vector is rather small, as the coeffients of $f_1$ and $f_2$ are only 336 bits large, significantly smaller than the 513-bit $p$.
Thus, we can use LLL to find this vector, and recover $f_1$ and $f_2$:

```
F = GF(p)
PR.<x> = PolynomialRing(F)
outputs = [(PR(i[0]),PR(y-i[1])) for i in outputs]
[(u1,v1),(u2,v2),(u3,v3)]=outputs
A = crt([v2**2,v3**2],[u2,u3])
U = crt([v1**2,A],[u1,u2*u3])
S=crt([0,x**7],[u1,u2*u3])
T=crt([0,x**6],[u1,u2*u3])

L=[]
L.append([int(i) for i in S.coefficients()]+[1,0,0])
L.append([int(i) for i in T.coefficients()]+[0,1,0])
L.append([int(i) for i in U.coefficients()]+[0,0,10**50])

for i in range(6):
L.append(i*[0]+[p]+(5-i+3)*[0])
L = matrix(L,ring=ZZ).LLL()

coefs = L[0]
if coefs[8] < 0:
coefs *= -1
assert coefs[8] == 10**50
s = -coefs[6]
t = -coefs[7]
coefs = list(coefs[:6]) + [t,s]
assert all(0<=i<2**336 for i in coefs)
print(f"Coefficients recovered: {coefs}")

f1 = sum(coefs[i] * (x ** i) for i in range(2 * g1 + 2))
f2 = sum(coefs[i] * (x ** i) for i in range(2 * g2 + 2))

HC1 = HyperellipticCurve(f1, 0)
J1 = HC1.jacobian()(GF(p))

HC2 = HyperellipticCurve(f2, 0)
J2 = HC2.jacobian()(GF(p))
o1 = J1(outputs[0])
o2 = J2(outputs[1])
o3 = J2(outputs[2])
```

For the first two parts of the flag, I implemented [this paper](https://www.sciencedirect.com/science/article/pii/S1071579715000787) to compute a "bisection" of the first output, recovering $j_1+j_2$ from $2(j_1+j_2)$:

```
F_split = f1(x=x**2).splitting_field("t")
f1_split=f1.change_ring(F_split)
roots=[i[0] for i in f1_split.roots()]
w=[u1(x=i).sqrt() for i in roots]
M=matrix([[roots[i]*u1(x=roots[i]),u1(x=roots[i]),-roots[i]**2*w[i],-roots[i]*w[i],-w[i]] for i in range(5)])
V=vector([v1(x=roots[i]) for i in range(5)])
X=M**(-1)*V
k1,k0,l,lu11,lu10=X
u11=lu11/l
u10=lu10/l
P = x**2 + u11 * x + u10 # this is the "u"-part of j1+j2
print(f"j1+j2 recovered: {P}")
# which is equal to (x-flag1)*(x-flag2)
# so we simply factorize
flag = []
for (root,_) in P.roots():
flag.append(int(root).to_bytes(64,"big").decode())
```

Now all that is left is to recover $j_2$ from $5\cdot j_2$.

This is probably hard in general (or at least, there is no easily google-able way to do this), but it helps that $j_2$ is of the form $[P]-[O]$. Denote its Mumford representation as $(u, v)$; then $u$ is just some linear polynomial $x + s$, and $v$ a constant.

It turns out, by the algorithm for adding divisors, that $5\cdot j_2$ has an *unreduced* representation of the form $(U,V)$ with $U=u^5$, and $V$ some polyonimal of degree at most $4$.

Now consider the algorithm for reducing some representation $(u,v)$:
- While $\deg u > \deg f$ do:
- Set $u \leftarrow \frac{f - v²}{u}$
- Set $v \leftarrow -v \mod u$
- Output $(u',v)$, where $u'$ is $u$ divided by its leading coefficient

By experimenting with generating my own values, I observed that when running on $(U,V)$, the loop of this algorithm will only run once before giving the output $(u,v)$. Denote by $l$ the leading coefficient of $u$ in the final step.

Because $(u_2,v_2)$ is the reduced representation of $5\cdot j_2$, we have:
$$u_2 = \frac{1}{l}\frac{f_2 - V^2}{U} = \frac{f_2 - V^2}{l(x+s)^5}$$
or $l \cdot (x+s)^5 \cdot u_2 = f_2 - V^2$.

Here, we know $f_2$ and $u_2$, but not $l, s$ or $V=v_4x^4+v_3x^3+v_2x^2+v_1x+v_0$.
As this is an equation of polyomials in $x$, we can look at each coefficient (with regards to $x$) seperately, and get a system of $\deg (f_2 - V^2) = 8$ equations.
But as there are only 7 unknowns, this is (likely) solvable.
I did so by computing the Groebner basis for this system of equations:

```
PRcrazy.<s,l,v4,v3,v2,v1,v0> = PolynomialRing(GF(p))
PRcrazyx = PRcrazy[x]
x = PRcrazyx(x)
V = x**4 * v4 + x**3 * v3 + x**2 * v2 + x*v1 + v0
f = f2(x=x)
a = o2[0](x=x)
T=(f-V**2)-l*(x-s)**5*a
I=PRcrazy.ideal(T.coefficients())
G=I.groebner_basis()
gs = [i for i in G if i.degree(s)==1][0]
# hopefully, we have something like "s + {constant} in there"
gs = gs.coefficients()
assert len(gs)==2 and 1 in gs
s = p - (set(gs)-{1}).pop()
print(s)
flag.append(int(s).to_bytes(64,"big"))
```