#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.*=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. *