#SageMath code for the experiments of the genus 3 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 for the underlying complex field.
#Output: The list of even theta constants of JSg3, and the 33th element in the list.
#Here JSg3 is the discrete Riemann matrix of level 7 in Table 5 of the arXiv version.
def computeEven(p):
    #Complex field with p bits of precision.
    CC.<I>=ComplexField(p)
    #The discrete Riemann matrix of level 7 in Table 5 of the arXiv version.
    JSg3=matrix(CC,[[-0.184080315481225+0.964526509924181*I,0.000676150987319488+0.000621759390095499*I,-0.183404164493888-0.0348517306851504*I],[0.000676150987324502+0.000621759390253077*I,-0.184080315481242+0.964526509924194*I,-0.183404164493867-0.0348517306853311*I],[-0.183404164493892-0.0348517306854673*I,-0.183404164493862-0.0348517306853276*I,-0.366808328987806+1.93029653862877*I]])
    #Riemann theta function in JSg3.
    RT=RiemannTheta(JSg3)
    #The list of even characteristics.
    even=[v for v in GF(2)^6 if v[:3]*v[3:] == 0]
    #The list of even theta constants.
    valuesEven=[vector(RT(char=c)) for c in even]
    return valuesEven,valuesEven[33]

#computeEven(p)[1] is the 33th element of the list of even theta constants in precision p.
#The list of these values when p ranges over [5,10,15,20,25,30,35]. 
[computeEven(p)[1] for p in [5,10,15,20,25,30,35]]
#This value is zero up to numerical round-off.