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) == 1p = 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)]=outputsA = 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 *= -1assert 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)*Vk1,k0,l,lu11,lu10=Xu11=lu11/lu10=lu10/lP = x**2 + u11 * x + u10 # this is the "u"-part of j1+j2print(f"j1+j2 recovered: {P}")# which is equal to (x-flag1)*(x-flag2)# so we simply factorizeflag = []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 + v0f = f2(x=x)a = o2[0](x=x)T=(f-V**2)-l*(x-s)**5*aI=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 gss = p - (set(gs)-{1}).pop()print(s)flag.append(int(s).to_bytes(64,"big"))    ```