# pbctf2023 writeup
- Crypto (blocky4, remake)

## blocky4

part of source (`task.py`)

def main():
    signal.signal(signal.SIGALRM, handler)
    signal.alarm(60)
    key = get_key()

    with open('flag', 'rb') as f:
        flag = f.read()
        flag += b'\x00' * ((-len(flag)) % 9)
    
    flag_cipher = BlockCipher(key, 20)
    enc_flag = b''
    for i in range(0, len(flag), 9):
        enc_flag += flag_cipher.encrypt(flag[i:i+9])
    print(f"enc_flag: {enc_flag.hex()}")

    cipher = BlockCipher(key, 4)
    while True:
        inp = bytes.fromhex(input("> "))
        assert len(inp) < (3 ** 7)
        assert all(b < 243 for b in inp)

        enc = b''
        for i in range(0, len(inp), 9):
            enc += cipher.encrypt(inp[i:i+9])
        
        print(f"Result: {enc.hex()}")

if __name__ == "__main__":
    main()

part of source (`Cipher.py`)

from GF import GF

SBOX, INV_SBOX = dict(), dict()
for i in range(3 ** 5):
    v = GF(23) + (GF(0) if i == 0 else GF(i).inverse())
    SBOX[GF(i)] = v
    INV_SBOX[v] = GF(i)

class BlockCipher:
    def __init__(self, key: bytes, rnd: int):
        assert len(key) == 9
        sks = [GF(b) for b in key]
        for i in range(rnd * 9):
            sks.append(sks[-1] + SBOX[sks[-9]])
        self.subkeys = [sks[i:i+9] for i in range(0, (rnd + 1) * 9, 9)]
        self.rnd = rnd

    def _add_key(self, l1, l2):
        return [x + y for x, y in zip(l1, l2)]

    def _sub_key(self, l1, l2):
        return [x - y for x, y in zip(l1, l2)]
    
    def _sub(self, l):
        return [SBOX[x] for x in l]

    def _sub_inv(self, l):
        return [INV_SBOX[x] for x in l]
    
    def _shift(self, b):
        return [
            b[0], b[1], b[2],
            b[4], b[5], b[3],
            b[8], b[6], b[7]
        ]
    
    def _shift_inv(self, b):
        return [
            b[0], b[1], b[2],
            b[5], b[3], b[4],
            b[7], b[8], b[6]
        ]
    
    def _mix(self, b):
        b = b[:] # Copy
        for i in range(3):
            x = GF(7) * b[i] + GF(2) * b[3 + i] + b[6 + i]
            y = GF(2) * b[i] + b[3 + i] + GF(7) * b[6 + i]
            z = b[i] + GF(7) * b[3 + i] + GF(2) * b[6 + i]
            b[i], b[3 + i], b[6 + i] = x, y, z
        return b
    
    def _mix_inv(self, b):
        b = b[:] # Copy
        for i in range(3):
            x = GF(86) * b[i] + GF(222) * b[3 + i] + GF(148) * b[6 + i]
            y = GF(222) * b[i] + GF(148) * b[3 + i] + GF(86) * b[6 + i]
            z = GF(148) * b[i] + GF(86) * b[3 + i] + GF(222) * b[6 + i]
            b[i], b[3 + i], b[6 + i] = x, y, z
        return b
    
    def encrypt(self, inp: bytes):
        assert len(inp) == 9
        b = [GF(x) for x in inp]
        
        b = self._add_key(b, self.subkeys[0])
        for i in range(self.rnd):
            b = self._sub(b)
            b = self._shift(b)
            if i < self.rnd - 1:
                b = self._mix(b)
            b = self._add_key(b, self.subkeys[i + 1])
        
        return bytes([x.to_int() for x in b])
    
    def decrypt(self, inp: bytes):
        assert len(inp) == 9
        b = [GF(x) for x in inp]

        for i in reversed(range(self.rnd)):
            b = self._sub_key(b, self.subkeys[i + 1])
            if i < self.rnd - 1:
                b = self._mix_inv(b)
            b = self._shift_inv(b)
            b = self._sub_inv(b)
        b = self._sub_key(b, self.subkeys[0])

        return bytes([x.to_int() for x in b])

(`GF.py` is skipped: an implementation of operations over $GF(3^5)$)

This challenge task is to decrypt encrypted flag with original block cipher on by using arbitrary known plaintext-ciphertext pairs (encryption oracle). This cipher is based on AES, but over $GF(3^5)$ instead over $GF(2^8)$.

The known effective attack for AES 4-round is [Square attack (from davidwong.fr)](https://www.davidwong.fr/blockbreakers/square_2_attack4rounds.html). It uses the property that sum of some set of 3-round (pre last round) encrypted plaintexts equals 0. This can be applied if:

- sum of all elements of base field equals 0 (all finite field satisfied)
- mixcolumn is one row linear multiplication
- sbox is permutation of all elements of base field

Then, we can apply the Square attack directly. Though I experienced I needed a few multiple sets for reducing key space on similar other challenge, this challenge needs only one set ($3^5$ plaintext-ciphertext pairs).

Here is the code.

from typing import List
import itertools

from pwn import *

from GF import GF
from Cipher import SBOX, INV_SBOX, BlockCipher

proc = remote('blocky-4.chal.perfect.blue', 1337)

encflag = bytes.fromhex(proc.recvline().strip().split(b': ')[1].decode())
(proc.recvuntil(b'> '))

ctlst = []
for i in range(3**5):
    pt = bytes([i]) + bytes([0 for _ in range(1, 9)])
    proc.sendline(pt.hex().encode())
    ctlst.append(bytes.fromhex(proc.recvline().strip().split(b': ')[1].decode()))

print('gen ctlst done')

def _sub_key(l1, l2):
    return [x - y for x, y in zip(l1, l2)]

def _shift_inv(b: List[GF]) -> List[GF]:
    return [
        b[0], b[1], b[2],
        b[5], b[3], b[4],
        b[7], b[8], b[6]
    ]

def _sub_inv(l: List[GF]) -> List[GF]:
    return [INV_SBOX[x] for x in l]

def after_shift_idx(i):
    return [0, 1, 2, 5, 3, 4, 7, 8, 6][i]

lastkeycand = [[] for _ in range(9)]
for i in range(3**5):
    keycand = [GF(i) for _ in range(9)]
    after_3round_result = [GF(0) for _ in range(9)]
    for ct in ctlst:
        state = [GF(ele) for ele in ct]
        state = _sub_key(state, keycand)
        state = _shift_inv(state)
        state = _sub_inv(state)
        for j in range(9):
            after_3round_result[j] += state[j]
    for j in range(9):
        if after_3round_result[j] == GF(0):
            lastkeycand[after_shift_idx(j)].append(GF(i))
print(lastkeycand)

rnd = 4
for key in itertools.product(*lastkeycand):
    subkeys = [None for i in range(0, rnd * 9)] + list(key)
    for i in range((rnd+1)*9-1, 9-1, -1):
        subkeys[i-9] = INV_SBOX[subkeys[i] - subkeys[i-1]]
    sks = bytes([x.to_int() for x in subkeys[:9]])
    cipher = BlockCipher(sks, 4)
    flag_cipher = BlockCipher(sks, 20)
    pt = b''
    for chunk in range(0, len(encflag), 9):
        pt += (flag_cipher.decrypt(encflag[chunk:chunk+9]))
    print(pt)

## remake (unsolved during the CTF)

source (`task.py`)

import os 

flag = open("flag", "rb").read()
flag = flag.lstrip(b"pbctf{").rstrip(b"}")
assert len(flag) == 192

while True:
    p = random_prime(1 << 513, lbound = 1 << 512)
    coefs = [int.from_bytes(os.urandom(42), "big") for _ in range(8)]
    PR.<x> = PolynomialRing(GF(p))

    g1, g2 = 2, 3
    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))

    flag1 = GF(p)(int.from_bytes(flag[:64], "big"))
    flag2 = GF(p)(int.from_bytes(flag[64:128], "big"))
    flag3 = GF(p)(int.from_bytes(flag[128:], "big"))
    hint = GF(p)(int.from_bytes(b"Inspired by theoremoon's SECCON 2022 Finals Challenge - Hell. Thank you!", "big"))

    pol1 = x * x - f1(flag1)
    pol2 = x * x - f1(flag2)
    pol3 = x * x - f2(flag3)
    pol4 = x * x - f2(hint)

    if len(pol1.roots()) * len(pol2.roots()) * len(pol3.roots()) * len(pol4.roots()) == 0:
        continue 

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

    HC2 = HyperellipticCurve(f2, 0)
    J2 = HC2.jacobian()(GF(p))

    P1 = HC1((flag1, pol1.roots()[0][0]))
    P2 = HC1((flag2, pol2.roots()[0][0]))
    P3 = HC2((flag3, pol3.roots()[0][0]))
    P4 = HC2((hint, pol4.roots()[0][0]))

    print(2 * J1(P1) + 2 * J1(P2))
    print(5 * J2(P3))
    print(J2(P4))
    break

This challenge is related to hyperelliptic curves. We need to recover defined polynomial and find flag which is encoded to a small multiple of jacobian of point representation.

I was not so familiar with hyperelliptic curve, so I looked at [Imaginary hyperelliptic curve (from wikipedia)](https://en.wikipedia.org/wiki/Imaginary_hyperelliptic_curve). Then, I remembered the jacobian group is the quotient of divisor group. And I knew the divisor group on elliptic curve is normal elliptic curve group with addition formula.

Then, I tried to understand how the jacobian group is generalized from addition formula on elliptic curve. The wikipedia says reduced divisors (that is, elements of the jacobian group) are represented as Mumford representation.

A divisor $D$ is uniquely represented as $D=(u(x), y-v(x))$ with $u(x), v(x) \in K[x]$ such as

- $\deg{v} < \deg{u} \le g$
- $u(x)$ is monic
- $((v(x)^2 + h(x)*v(x)) - f(x)) \mod{u(x)} = 0$
where the hyperelliptic curve is defined as $y^2+h(x)*y=f(x)$ over $K$ for $\deg{f}=2*g+1, \deg{h}\le g$.

(Note that some article uses $D=(u(x), v(x))$ cause it is determined only $u,v$. On Sagemath, a divisor is printed as our representation, but $Jac(P)[1]$ returns $v(x)$.)
From now, assume $h(x)=0$.

On elliptic curve, each non-infinity point $P=(x_0, y_0)$ corresponding to the divisor $(x-x_0, y-y_0)$. (It is written as $[P]-[\infty]$.) This seems same on hyperelliptic curve. In fact, all non-infinity points on hyperelliptic curve has similar representation.

Then, I experimented for $[i]*Jac(P)$ on $i\le g+1$. I realized

- For $i\le g$, $[i]*Jac(P)[0] = (x-x_0)^i$
- For $i=g+1$, not equal to $(x-x_0)^i$ anymore

For elliptic curve, doubling formula is based on cross point of the tangent line of $P$. So I assume that it is related to some Taylor expansion of $y=\sqrt{f(x)}$. So I put as
$\bar{v}(x)=-(y_0 + \frac{1}{1!}*\frac{dy}{dx}(P)*(x-x_0) + \frac{1}{2!}*\frac{d^2y}{dx^2}(P)*(x-x_0)^2 + \ldots + \frac{1}{g!}*\frac{d^gy}{dx^g}(P)*(x-x_0)^g)$ and confirmed

$(x-x_0)^{(g+1)}*u(x) = (f(x) - {\bar{v}(x)}^2)/\mathbb{lc}(f)$. ($\mathbb{lc}(f)$ is leading coefficient of $f$)
(I divides $\mathbb{lc}(f)$ for $u(x)$ as monic.)
I also confirm that similar formula is satisfied on $i=g+2$. (but we need to divide $-\mathbb{lc}(\bar{v})^2$ instead of $\mathbb{lc}(f)$.)

During the CTF, I applied this formula **without knowing $f(x)$**, but we did not find $f(x)$.

After the CTF, we saw the [writeup](https://rkm0959.tistory.com/285) by @rkm0959. We really depressed that finding $f(x)$ is easy by using LLL and Mumford representation rule. I forgot the check for bit length of coefficients of $f(x)$ are fairly small compared to $p$. Maybe the reason is that I too focused on hyperelliptic curve thing...

However, I notice that @rkm0959 approach uses equation similar to us. But direct substituting $x_0, y_0$ is really wasted time. I rethinks how to solve equation.

Firstly considering $(U2(x), y-V2(x))=5*J2(P3)$ ($P3=(a, \sqrt{f(a)})$, $a=\mathbb{flag3}$). Since this is $5=g_3+2$ case, I have the following formula.

$(x-a)^5 * U2(x) = (f(x)-\bar{V2}(x))/(-\mathbb{lc}(\bar{V2})^2)$.

$\deg{\bar{V2}}=5-1, \deg{U2}=3(=g_3), \deg{f}=7(=2*g_3+1)$. By Mumford represention rule, $f(x)={\bar{V2}(x)}^2 = {V2(x)}^2 \mod{U2(x)}$. Then, I write as

$\bar{V2}(x) = -V2(x) + (A+B*x)*U2(x)$
with unknown $A,B$ and $\mathbb{lc}(\bar{V2})=B$. (I uses $-V2(x)$ cause it is just similar to elliptic curve addition formula. Which one uses does not affect the system.)
By substituting $\bar{V2}(x)$ to above equation and comparing each terms for $1,x,x^2,\ldots,x^8$, we have 9 equations with unknown variable $A, B, a$. (The degree 8 equation is totally useless.) It can be solved by Sagemath.

Then, switch to $(U1(x), y-V1(x))=2*J1(P1)+2*J1(P2)$, ($P1_x=b=\mathbb{flag1}, P2_x=c=\mathbb{flag2}$).
Just thinking similar to above formula,

${((x-b)*(x-c))}^2 * U1(x) = (f(x)-{\bar{V1}(x)}^2)/(-\mathbb{lc}(\bar{V1})^2)$

By using symmetricity, write $\alpha=b+c, \beta=b*c$.

$(x^2-\alpha*x+\beta)^2*U1(x)=(f(x)-{\bar{V1}(x)}^2)/(-\mathbb{lc}(\bar{V1})^2)$

By substituting $\bar{V1}(x) = -V1(x) + (A+B*x)*U1(x)$ as just above, we have 7 equations with unknown variable $A, B, \alpha, \beta$. Unfortunatelly, Sagemath cannot solve the system (maybe Singular issue over the finitefield with large prime characterstic), so I tossed the system to [Magma online calculator (from magma.maths.usyd.edu.au)](http://magma.maths.usyd.edu.au/calc/). It solved the system instantly. Then, I find the solution $b,c$ with quadratic univariate equation from $\alpha,\beta$. so I obtained all flag.

## precompute by J2_P4[0].subs({x:hint, y:y}) = 0 mod p
p = 20206629497318640302613717177493021659164910368667008120702312168658959729889280963268869447056316491892128508524975162223724948508181049113205403468161303

PR = PolynomialRing(GF(p), 2, ["x", "y"])
x,y = PR.gens()

## from output.txt
two_J1_P1_two_J1_P2 = (x^2 + 14762123602851553604749022996287576858980220577795651427829269858766434621297346961387874961427459051934768224338447011128244905975068497090840444625419470*x + 8519674750729750620690035589812482119785861876353468044895414394332293279114303071755954851101633319350193436546144692795403444364414318973131157246232656, y + 17770738679279057916355557895675090129563269633432826251932824463003364931275912702916209480950481351904761364290424406482997835483807402182326014818733821*x + 12306668467523337827805393760490897581559948654643366727345701375757143864825442910779617850907143245102792529282031529618639723158417652048624567379151171)
five_J2_P3 = (x^3 + 13441279154284544764330805782065565325543470739559917045273482055514440837785754044182874902421009026981197721504820302867945812937528249594953326223176272*x^2 + 3795282115520834934850220740151212731596814319504043674340537364041453624883995759365119899076774262882230308591629439035308527946872182029742910504122735*x + 3726617245981099594981815385059428688276726297460965450328320328460867196111587736356492934195556032891106446058683147130913147722036293641303193921962091, y + 2103349591221335944593862709600493681857281410337020721978302326614691696399677635217262732543672829811190387220058078405239568477387817550236173432744263*x^2 + 4784247634355946154999459446762911004042472267922959302672838559247991353014786987556174410735592161587023899368989617780068662559773261109676326152316907*x + 2640959823121300693709616791657128464111647959613642856293592234010564318329382577397798309822254798484629398268742247779165733848105319417195858443049412)
J2_P4 = (x + 540914350057159927735436910109553086959907629357396102386062800609858847095404691934517983640846787552171946520460064397364544720996858815571065665255645, y + 541917331856005964100090629475512429550322452567752818120774876171019476274441296070275457561095853517207532108745504694853066426720092700847788666013730)

U1 = two_J1_P1_two_J1_P2[0]
V1 = -two_J1_P1_two_J1_P2[1].subs({x:x, y:0})
U2 = five_J2_P3[0]
V2 = -five_J2_P3[1].subs({x:x, y:0})
U3 = J2_P4[0]
V3 = -J2_P4[1].subs({x:x, y:0})

hint = GF(p)(int.from_bytes(b"Inspired by theoremoon's SECCON 2022 Finals Challenge - Hell. Thank you!", "big"))

assert J2_P4[0].subs({x:hint, y:y}) == GF(p)(0)

g1, g2 = 2, 3


### step1: compute defining polynomial for hyperelliptic curve
hc_def_pol_coef_rng = PolynomialRing(GF(p), 2*g2+2, ["c%d" % i for i in range(2*g2+2)])
ci = hc_def_pol_coef_rng.gens()

hc_def_pol_rng = PolynomialRing(hc_def_pol_coef_rng, "X")
X = hc_def_pol_rng.gens()[0]

f1symbpol = sum(ci[i] * (X ** i) for i in range(2 * g1 + 2))
f2symbpol = sum(ci[i] * (X ** i) for i in range(2 * g2 + 2))

## Mumford representation relation: f = V^2 mod U
mumford_rel_1 = (f1symbpol - V1.subs({x:X, y:0})^2).quo_rem(U1.subs({x:X, y:0}))[1].list() # g1-equations
mumford_rel_2 = (f2symbpol - V2.subs({x:X, y:0})^2).quo_rem(U2.subs({x:X, y:0}))[1].list() # g2-equations
mumford_rel_3 = (f2symbpol - V3.subs({x:X, y:0})^2).quo_rem(U3.subs({x:X, y:0}))[1].list() # 1-equations

mumford_rel = mumford_rel_1 + mumford_rel_2 + mumford_rel_3
## linear equations, but 2-dim shortage... but each coefficients are (42*8) bits << 512bits -> LLL!

## from rkm Inequality_Solving_with_CVP
load("/home/sage/inequ_cvp_solve.sage")

coef, monos = Sequence(mumford_rel).coefficient_matrix()

mat = []
linequlen = coef.nrows()
linequvars = 2 * g2 + 2
assert coef.ncols() == linequvars + 1

for i in range(linequlen):
    mat.append([int(coef[i][j].lift()) for j in range(linequvars+1)] + [0]*i + [p] + [0]*(linequlen-i-1))
for i in range(linequvars):
    mat.append([0]*i + [1] + [0]*(linequvars-i-1) + [0] + [0]*linequlen)
mat.append([0]*linequvars + [1] + [0]*linequlen)

mat = matrix(ZZ, mat).transpose()
lb = [0]*linequlen + [0]*linequvars + [1]
ub = [0]*linequlen + [2**(8*42)]*linequvars + [1]

print("computing def poly by LLL")
sols = solve(mat, lb, ub)
sol = sols[2] # var solution (sols[0] is equ solution)

f1 = sum(GF(p)(sol[i]) * (X ** i) for i in range(2 * g1 + 2))
f2 = sum(GF(p)(sol[i]) * (X ** i) for i in range(2 * g2 + 2))

assert (f1 - V1.subs({x:X})^2) % U1.subs({x:X}) == 0
assert (f2 - V2.subs({x:X})^2) % U2.subs({x:X}) == 0
assert (f2 - V3.subs({x:X})^2) % U3.subs({x:X}) == 0


# step2: find flag3
polsymbrng_flag3 = PolynomialRing(GF(p), 3, ["A", "B", "a"])
A, B, a = polsymbrng_flag3.gens()
hc_def_pol_flag3 = PolynomialRing(polsymbrng_flag3, "X")
X_flag3 = hc_def_pol_flag3.gens()[0]

## v_bar_flag3: deg=4 (for reducing (X-a)^5), f=v_bar_flag3^2 = V2^2 mod U2
v_bar_flag3 = -V2.subs({x:X_flag3}) + (A + B*X_flag3) * U2.subs({x:X_flag3})
f2symbpol_flag3 = sum([GF(p)(ele) * X_flag3^i for i, ele in enumerate(f2.list())])
lefthand = f2symbpol_flag3 - v_bar_flag3^2
righthand = -B^2 * (X_flag3-a)^5 * U2.subs({x:X_flag3})
lefthand_list = lefthand.list()
righthand_list = righthand.list()

assert len(lefthand_list) == (1+3)*2 + 1 # deg8
assert len(righthand_list) == 5 + 3 + 1 # deg8

equ = []
for i in range(len(lefthand_list)):
    equele = lefthand_list[i] - righthand_list[i]
    equ.append(equele)

I_flag3 = polsymbrng_flag3.ideal(equ)
print("computing flag3")
V_flag3 = I_flag3.variety()
print(V_flag3)
flag3_root = V_flag3[0][a]

flag3 = int.to_bytes(int(flag3_root.lift()), 64, 'big')
print(flag3)


# step3: find flag1, flag2
polsymbrng_flag12 = PolynomialRing(GF(p), 4, ["A", "B", "a", "b"])
A, B, a, b = polsymbrng_flag12.gens()
hc_def_pol_flag12 = PolynomialRing(polsymbrng_flag12, "X")
X_flag12 = hc_def_pol_flag12.gens()[0]

## v_bar_flag12: deg=3 (for reducing (X^2-a*X+b)^2), f=v_bar_flag12^2 = V1^2 mod U1
v_bar_flag12 = -V1.subs({x:X_flag12}) + (A + B*X_flag12) * U1.subs({x:X_flag12})
f1symbpol_flag12 = sum([GF(p)(ele) * X_flag12^i for i, ele in enumerate(f1.list())])
lefthand = f1symbpol_flag12 - v_bar_flag12^2
righthand = -B^2 * (X_flag12^2 - a*X_flag12 + b)^2 * U1.subs({x:X_flag12})
lefthand_list = lefthand.list()
righthand_list = righthand.list()

assert len(lefthand_list) == (1+2)*2 + 1 # deg6
assert len(righthand_list) == 2*2 + 2 + 1 # deg6

equ = []
for i in range(len(lefthand_list)):
    equele = lefthand_list[i] - righthand_list[i]
    equ.append(equele)

I_flag12 = polsymbrng_flag12.ideal(equ)
print("computing flag12")
G_flag12 = I_flag12.groebner_basis()


## for Magma program
G12out = open('G12out.txt', 'w')
G12out.write(f"p:={p};\n")
G12out.write(f"P<A, B, a, b>:=PolynomialRing(GF(p), 4);\n")
G12out.write(f"I := ideal<P |\n")
G12out.write(f"{str(G_flag12[:5])[1:-1]}>;\n")
G12out.write(f"Variety(I);\n")

## Sage9.7 variety does not compute solution... instead use Magma Online calc
#V_flag12 = I_flag12.variety()
#print(V_flag12)
sol = [(13251201369422734297366526607962697431123459669915815849542141065457700206210679127885890919842739624682727419516349360728543108726009316584250248278486350, 16729053782777399746765327111114721998547116314370943300162668871471785975926363509665997824645808370479979293664335945059595374048602182572228181517703089, 5224076322553305114437830665526577139615879548772477449710363729574538835829096345899844681709029008011430318958967087007237927563994860535590585262112322, 1508986926034823363823763034485115530088825680173293957835708459247418883522363409626628316500047427774631513161686709517586272538191205520278401195650306)]


flag12_add_root = GF(p)(sol[0][2])
flag12_mul_root = GF(p)(sol[0][3])

flag12_pol = PolynomialRing(GF(p), "flag12")
flag12 = flag12_pol.gens()[0]
flag12_roots = (flag12^2 - flag12_add_root * flag12 + flag12_mul_root).roots()
print(flag12_roots)

flag1 = int.to_bytes(int(flag12_roots[0][0].lift()), 64, 'big')
print(flag1)
flag2 = int.to_bytes(int(flag12_roots[1][0].lift()), 64, 'big')
print(flag2)

Above method can be applied repeatedly for more general challenge. The method would be justified by [Group Law Computations on Jacobians of Hyperelliptic Curves](https://eprint.iacr.org/2011/306.pdf). The jacobian addition is related to tangent curve and crossing curve (interpolating polynomial).
@rkm0959 points out this challenge is releated to divide-2 for genus 3 problem. It requires 2 steps above method. (and divide-3 for genus 3, 3 steps.) I am thinking rkm method can be worked with reducing variables. (instead of vanishing out $T(x)$, but positively writes $T(x)$ with $u(x), v(x)$.)

Anyway, thank you for nice challenges to the author @rbtree, @rkm0959.