// this program enables to compute the sign in Lemma 3.6 of "a new proof of a Thomae like-formula
// for plane smooth quartics"


// definition of theta function at z with characteristics M, Riemann matrix tau with a certain
// precision prec

function theta_int(z,M,n,tau)
e1:=M[1];
f1:=M[2];
C:=BaseRing(tau);
I:=C.1;
pi:=Pi(C);
g:=#z;
S:=0;
e:=ChangeRing(e1,C);
f:=ChangeRing(f1,C);
for i:=1 to g do
for j:=1 to g do
S:=S+ (n[i]+e[i]/2)*tau[i,j]*(n[j]+e[j]/2);
end for;
end for;
S:=1/2*S;
for i:=1 to g do
S:=S+(n[i]+e[i]/2)*(z[i]+f[i]/2);
end for;
return Exp(2*I*pi*S);
end function;

function theta(z,M,tau,prec)
g:=#z;
L:=CartesianPower([-prec..prec],g);
S:=0;
for n in L do
S:=S+theta_int(z,M,n,tau);
end for;
return S;
end function;


// definition of the vth derivative of theta function at z
//with characteristics M, Riemann matrix tau with a certain precision prec

function theta_par_int(z,M,n,tau,v)
e1:=M[1];
f1:=M[2];
C:=BaseRing(tau);
I:=C.1;
pi:=Pi(C);
g:=#z;
S:=0;
e:=ChangeRing(e1,C);
f:=ChangeRing(f1,C);
for i:=1 to g do
for j:=1 to g do
S:=S+ (n[i]+e[i]/2)*tau[i,j]*(n[j]+e[j]/2);
end for;
end for;
S:=1/2*S;
for i:=1 to g do
S:=S+(n[i]+e[i]/2)*(z[i]+f[i]/2);
end for;
return 2*I*pi*(n[v]+e[v]/2)*Exp(2*I*pi*S);
end function;


function theta_par(z,M,tau,v,prec)
g:=#z;
L:=CartesianPower([-prec..prec],g);
S:=0;
for n in L do
S:=S+theta_par_int(z,M,n,tau,v);
end for;
return S;
end function;


// definition of Jacobian Nullwert at tau with respect to the g first entries
// of a fundamental system EN

function jacobi_nullwerte(EN,tau,prec)
C:=BaseRing(tau);
I:=C.1;
pi:=Pi(C);
g:=#EN div 2 -1;
z:=[C!0 : i in [1..g]];
M:=ZeroMatrix(C,g,g);
for i:=1 to g do
for v:=1 to g do
M[i,v]:=theta_par(z,EN[i],tau,v,prec);
end for;
end for;
return Determinant(M)/pi^g;
end function;


// definition of the product of Thetanullwerte for a term in the Riemann-Jacobi formula
// related to a fundamental system EN

function theta_product(EN,tau,prec)
C:=BaseRing(tau);
g:=#EN div 2 -1;
z:=[C!0 : i in [1..g]];
P:=1;
for i:=g+1 to 2*g+2 do
P:=P*theta(z,EN[i],tau,prec);
end for;
return P;
end function;

// the normalization of a fundamental system with coefficients 0 or 1

function reduce_mat(e)
g:=NumberOfColumns(e);
return Matrix(2,g,[[e[1][i] mod 2 : i in [1..g]],[e[2][i] mod 2 : i in [1..g]]]);
end function;

function normalized_FS(EN)
return [reduce_mat(e) : e in EN];
end function;

// The operation of translation of a fundamental system EN by v=EN[j]+EN[2g+2]

function translateFS(EN,j)
h:=#EN;
vj:=EN[j]-EN[h];
E:=[M + vj : M in EN];
Eh:=E[h];
Ej:=E[j];
E[j]:=Eh;
E[h]:=Ej;
return E;
end function;

// The specific example in the article

EN:=[Matrix(2,3,[[1,0,0],[1,0,0]]),Matrix(2,3,[[0,1,0],[1,1,0]]),Matrix(2,3,[[0,0,1],[1,1,1]]),
Matrix(2,3,[[1,0,0],[0,0,0]]),Matrix(2,3,[[0,1,0],[1,0,0]]),Matrix(2,3,[[0,0,1],[1,1,0]]),Matrix(2,3,[[0,0,0],[1,1,1]]),Matrix(2,3,[[0,0,0],[0,0,0]])];

ENp:=[Matrix(2,3,[[0,1,1],[0,0,1]]),Matrix(2,3,[[1,0,1],[0,1,1]]),Matrix(2,3,[[1,1,0],[0,1,0]]),
Matrix(2,3,[[1,0,0],[0,0,0]]),Matrix(2,3,[[0,1,0],[1,0,0]]),Matrix(2,3,[[0,0,1],[1,1,0]]),Matrix(2,3,[[0,0,0],[1,1,1]]),Matrix(2,3,[[1,1,1],[1,0,1]])];

EN1:=translateFS(EN,1);
EN1:=normalized_FS(EN1);
EN2:=translateFS(EN,2);
EN2:=normalized_FS(EN2);
EN3:=translateFS(EN,3);
EN3:=normalized_FS(EN3);


ENp1:=translateFS(ENp,1);
ENp1:=normalized_FS(ENp1);
ENp2:=translateFS(ENp,2);
ENp2:=normalized_FS(ENp2);
ENp3:=translateFS(ENp,3);
ENp3:=normalized_FS(ENp3);

N:=50;
prec:=5;
C<I>:=ComplexField(N);

// this is the matrix tau we choose in Lemma 3.6
tau:=Matrix([[5*I,2*I+1,2*I+2],[2*I+1,5*I,3+I],[2*I+2,3+I,4*I]]);

// S([N_0],tau)
jacN:=jacobi_nullwerte(EN,tau,prec);
prodN:=theta_product(EN,tau,prec);
jacN/prodN;

// S([N_1],tau)
jacN1:=jacobi_nullwerte(EN1,tau,prec);
prodN1:=theta_product(EN1,tau,prec);
jacN1/prodN1;

// S([N_2],tau)
jacN2:=jacobi_nullwerte(EN2,tau,prec);
prodN2:=theta_product(EN2,tau,prec);
jacN2/prodN2;

// S([N_3],tau)
jacN3:=jacobi_nullwerte(EN3,tau,prec);
prodN3:=theta_product(EN3,tau,prec);
jacN3/prodN3;


// S([N_0'],tau)
jacNp:=jacobi_nullwerte(ENp,tau,prec);
prodNp:=theta_product(ENp,tau,prec);
jacNp/prodNp;

// S([N_1'],tau)
jacNp1:=jacobi_nullwerte(ENp1,tau,prec);
prodNp1:=theta_product(ENp1,tau,prec);
jacNp1/prodNp1;

// S([N_2'],tau)
jacNp2:=jacobi_nullwerte(ENp2,tau,prec);
prodNp2:=theta_product(ENp2,tau,prec);
jacNp2/prodNp2;

// S([N_3'],tau)
jacNp3:=jacobi_nullwerte(ENp3,tau,prec);
prodNp3:=theta_product(ENp3,tau,prec);
jacNp3/prodNp3;

// \iota(N_0)
(jacN/prodN*jacN1/prodN1*jacN2/prodN2*jacN3/prodN3)/
(jacNp/prodNp*jacNp1/prodNp1*jacNp2/prodNp2*jacNp3/prodNp3);