#SageMath code for the experiments of the genus 4 JS surface in Section 3.3.3 of the arXiv version.

#Importing RiemannTheta package, see https://github.com/nbruin/RiemannTheta for the installation.
from riemann_theta.riemann_theta import RiemannTheta

#Input: Precision of the underlying complex field.
#Output: The value of Schottky-Igusa modular form when evaluated at JSg4.
#The matrix JSg4 is the discrete Riemann matrix of level 7 in Table 6 of the arXiv version.
#For the modular form, see [Igu82] in the references of the arXiv version. 

def Igusa(p):
    #Complex field with p bits of precision.
    CC.<I>=ComplexField(prec=p)
    #The discrete Riemann matrix of level 7 in Table 6 of the arXiv version.
    JSg4=matrix(CC,[[-0.184082280391336+0.964527352840432*I,0.000339057948694669+0.000310458236828142*I,0.000339057948610197+0.000310458236672205*I,-0.183404164493683-0.0348517306848045*I],[0.000339057948721604+0.000310458236985237*I,-0.184082280391383+0.964527352840458*I,0.000339057948558669+0.00031045823671365*I,-0.183404164493691-0.0348517306850943*I],[0.000339057948725268+0.000310458236989523*I,0.000339057948727943+0.000310458236994094*I,-0.184082280391416+0.964527352840407*I,-0.183404164493851-0.03485173068532*I],[-0.183404164493886-0.0348517306854523*I,-0.183404164493833-0.0348517306853228*I,-0.183404164493574-0.0348517306849086*I,-0.550212493481976+2.89544480794274*I]])
    #Riemann theta function in JSg4.
    RT=RiemannTheta(JSg4)
    #Finite field of size 2.
    G=GF(2)
    #The characteristics that form the modular form.
    #e.g. m1 is the characteristic [[1,0,1,0],[1,0,1,0]].
    m1=vector(G,[1,0,1,0,1,0,1,0])
    m2=vector(G,[0,0,0,1,1,0,0,0])
    m3=vector(G,[0,0,1,1,1,0,1,1])
    n1=vector(G,[0,0,0,1,1,1,1,0])
    n2=vector(G,[0,0,1,1,0,0,0,1])
    n3=vector(G,[0,0,1,0,1,0,1,1])
    N=[a*n1+b*n2+c*n3 for a in G for b in G for c in G]
    m1N=[m1+n for n in N]
    m2N=[m2+n for n in N]
    m3N=[m3+n for n in N]
    Pi1=prod([RT(char=c) for c in m1N])
    Pi2=prod([RT(char=c) for c in m2N])
    Pi3=prod([RT(char=c) for c in m3N])
    return Pi1^2+Pi2^2+Pi3^2-2*Pi1*Pi2-2*Pi1*Pi3-2*Pi2*Pi3

#The list of the values when p ranges over [5,8,10,13,16,19,23,27] to estimate Igusa(p).
[Igusa(p) for p in [5,8,10,13,16,19,23,27]]
#Here we keep the precision low as the level 7 approximation achieves a few digits of the entries
and the discrete Riemann matrix does not necessarily lie on the Schottky locus. 

####################################################################################################
#In this part, we experiment to see the number of even vanishing theta constants of JSg4.
#We start with 5 bits of precision. One may increase the precision gradually to keep track of these values. 
#The number of such values appear to be 10. 
#Here we keep the precision low as the level 7 approximation achieves a few digits of the entries and the discrete Riemann matrix does not necessarily lie on the Schottky locus. 

#Complex field with 5 bits of precision.
CC.<I>=ComplexField(prec=5)
#The discrete Riemann matrix of level 7 in Table 6 of the arXiv version.
JSg4=matrix(CC,[[-0.184082280391336+0.964527352840432*I,0.000339057948694669+0.000310458236828142*I,0.000339057948610197+0.000310458236672205*I,-0.183404164493683-0.0348517306848045*I],[0.000339057948721604+0.000310458236985237*I,-0.184082280391383+0.964527352840458*I,0.000339057948558669+0.00031045823671365*I,-0.183404164493691-0.0348517306850943*I],[0.000339057948725268+0.000310458236989523*I,0.000339057948727943+0.000310458236994094*I,-0.184082280391416+0.964527352840407*I,-0.183404164493851-0.03485173068532*I],[-0.183404164493886-0.0348517306854523*I,-0.183404164493833-0.0348517306853228*I,-0.183404164493574-0.0348517306849086*I,-0.550212493481976+2.89544480794274*I]])
#Riemann theta function in JSg4.
RT=RiemannTheta(JSg4)
#The list of even characteristics.
even=[v for v in GF(2)^8 if v[:4]*v[4:] == 0] 
#Here we start with low precision as the level 7 approximation achieves a few digits and the discrete Riemann matrix does not necessarily lie on the Schottky locus. 
#One may experiment further how the values behave by tuning the precision. 
#The list of even theta constants (with precision 5 bits).
valuesEven=[RT(char=c) for c in even]