#The code takes as input the equation defining a hyperelliptic curve and returns the #corresponding tau function, according to the Sato formulation. We follow the description #in https://arxiv.org/abs/1808.06748

#The following equations correspond in order to curves of genus 2,3,4. Uncomment the one you want to try out. 
hyp := (x-1)*(x-1-e)*(x-2)*(x-2-e)*(x-3)*(x-3-e):
#hyp := (x + 1 + e)*(x + 1 + 2*e)*(e^2 + e + x + 1)*(e^2 + 2*e + x + 1)*(x + 2 + e)*(x + 2 + 2*e)*(e^2 + e + x + 2)*(e^2 + 2*e + x + 2):
#hyp := (x - 1)*(x - 1 - e)*(x - 2)*(x - 2 - e)*(x - 3)*(x - 3 - e)*(x - 4)*(x - 4 - e)*(x - 5)*(x - 5 - e):




M := 3:    #this value determines the highest order of e appearing in the next calculation
N := 40:


d := degree(hyp, x):
x := 1/z:
y := expand(factor(hyp*z^d))^(1/2):
S := mtaylor(mtaylor(y, e, M), z, N + d/2):
n := d/2:
m := n - 1:


#define the array of the coefficients alpha_i
a := Array(0 ..N + n + m):
a[0] := 1:
for i to N + n do
    a[i] := factor(coeff(S, z, i))
end do:


#define the array of the functions g_n,...,g_(N+n+m)
g := Array(1 ..N + n + m):
for i to N + n do
    for j to i + n 
        do g[i + m] := alpha[j - 1]*x^(i + n - j) + g[i + n - 1] 
    end do:
end do:


#define the array of the functions f_n,...,f_(N+n+m)
f := Array(1 ..N + n + m):
for i to N + n do
    f[i + m] := expand(factor(1/2*(x^(i + m)*y + g[i + m])*z^(i + m)))
end do:


#define the array with the series expansions of the functions 1,f_n,...,f_(N+n+m)
Cmatrix := Array(1 ..N + n + 1):
Cmatrix[1] := taylor(z^m, z, N):
for i to N + n do
    Cmatrix[i + 1] := mtaylor(mtaylor(f[i + m], e, M), z, i + 1 + N + m)
end do:


#compute and display the matrix corresponding to the frame computed above
xi := (i, j) -> factor(coeff(Cmatrix[N + n + 2 - j], z, i - j)):
A := Matrix(2*(N + n), N + n + 1, xi);
numRows := 2*(N + n) + 1:


#compute the elementary Schur polynomials in three variables up to partitions of partLength 
partLength := 8:
p = Array(1 ..partLength):
for i to partLength do
    p[i] := 0
end do:
n1 := 0: n2 := 0: n3 := 0:
for i to partLength do
    for n1 from 0 to i do
        for n2 from 0 to i do
            for n3 from 0 to i do
                if n1 + 2*n2 + 3*n3 = i then
                   p[i] := p[i] + X^n1*Y^n2*T^n3/(n1!*n2!*n3!)
                end if: 
            end do:
        end do:
    end do:
end do:


with(combinat, partition): with(ListTools):
partitions := FlattenOnce([seq(partition(j), j = 1 .. partLength)]):
h := numelems(partitions):

#compute the Maya diagrams corresponding to the partitions
Maya = Array(1 ..h):
for i to h do
    r := numelems(partitions[I]):
    Maya[i] := Array(1 ..r):
    for j to r do
        Maya[i][j] := partitions[i][r - j + 1] - j + N + 4
    end do:
end do:
eqns := {e^4 = 0}:
tau := 0:

#compute the coordinates, the Schur polynomial for the non zero coordinates and the tau function
Coordinates = Array(1 ..h):
with(LinearAlgebra): with(ArrayTools): with(linalg):
for i to h do
    r := numelems(Maya[i]):
    increasingMaya := convert(FlipDimension(Maya[i], 2), list):
    if Maya[i][1] < numCol then
        Coordinates[i] := simplify(Determinant(SubMatrix(A, increasingMaya, [N + n + 1 - r + 1 .. N + n + 1])), eqns):
    else
        Coordinates[i] := 0:
    end if:
    if Coordinates[i] <> 0 then
        l := Array(1 ..r):
        v := Vector(1 ..r):
        for k to r do
            l[k] := partitions[i][k] + k - 1:
            v[k] := p[l[k]]:
        end do:
        Wr := wronskian(v, X):
        schur := det(Wr):
        tau := schur*Coordinates[i] + tau:
    end if:
end do:
tau;




#this part of the code plugs tau function into the Hirota bilinear form

t := (X, Y, T) -> tau:
taux := diff(t(X, Y, T), X):
tauxx := diff(taux, X):
tauxxx := diff(tau;xx`, X):
tauxxxx := diff(tau;xxx`, X):
taut :=diff(t(X, Y, T), T):
tauxt := diff(tau;x`, T):
tauy := diff(t(X, Y, T), Y):
tauyy := diff(tau;y`, Y):
Hirota := simplify(tauxxxx*tau - 4*tauxxx*taux + 3*tauxx^2 + 4*(taut*taux-tauxt*tau) + 3*(-tauy^2 +tauyy*tau), eqns);