# Reconstruction of curves from their theta hyperplanes in genera $6$ and $7$

Türkü Özlüm Çelİk  Koç University, Istanbul<EMAIL_ADDRESS>and  David
Lehavi  Tel Aviv<EMAIL_ADDRESS>
###### Abstract.

We derive a formula for reconstructing a generic complex canonical curve $C$
of genus 6 and 7 in terms of the theta hyperplanes of $C$. Hence, we get a
generic inverse to the Torelli map, as well as a complete description of the
Schottky locus in these genera. The computational part of the proof relies on
a certified numerical argument.

###### Key words and phrases:

Curves, Theta hyperplanes, Reconstruction, Inverse Torelli, Schottky,
Certification

###### 2020 Mathematics Subject Classification:

Primary: 14H40, 14H42, 65G40; Secondary: 14H45, 14Q05, 14-04

##  1\. Introduction

A generic complex algebraic curve of genus $g$ possesses $2^{g-1}(2^{g}-1)$
odd theta characteristics, in other words, half canonical divisors of
dimension 1. These theta characteristics uniquely determine the curve, as
shown in [CS03]. By identifying the curve with its canonical image in
$|K_{C}|^{*}\cong\mathbb{P}^{g-1}$, one can consider the _theta hyperplanes_
as the corresponding geometric objects on the canonical model, which meet the
curve in the projective space at $g-1$ points with multiplicity 2 at each
point. This classical topic in algebraic geometry can be studied from various
perspectives. One interesting approach is to reconstruct the canonical image
of the curve given all the theta hyperplanes. This problem has been
extensively investigated for non-hyperelliptic curves of genus 3, with
pioneering work by Aronhold [Aro64] and Coble [Cob29]. For recent expositions
see [Dol12, 6.1.2] and Remark 2.4 respectively. For curves of higher genera,
explicit constructions have been studied for genera 4 and 5, as discussed in
[Leh15] and later work [ÇKRS19, Leh22]. In this article, we extend these
results to curves of genera 6 and 7, providing new insights into the general
reconstruction problem.

In Section 2 we outline the proof of our main result, namely Theorem 2.7. It
explicitly describes generic complex canonical curves of genera 6 and 7 in
terms of their theta hyperplanes. Corollary 2.10 provides a solution to the
Schottky problem in genera 6 and 7, which asks for a characterization of the
locus of Jacobian varieties among principally polarized abelian varieties. As
in previous works [Leh15, Leh22] the proof of Theorem 2.7 starts by proving
the claim for a specific curve, which relies on both symbolic and numerical
computational techniques, and concludes by using a classical degeneration
argument. When $g=6$ and $7$, we use Wiman’s sextic and the Fricke-Macbeath
curve as the specific curves respectively, which are presented in Section 3.
In the same section, we also review the connection between the theta
hyperplanes and the Riemann theta function, emphasizing their analytical
representations in terms of the theta functions, which come with a complete
level 2 structure.

Using computational tools, which rely on numerical techniques, we apply our
formula to reconstruct these specific curves from non-certified analytical
representations of the hyperplanes. Hence, to conclude our proof, we must
certify the computations i.e., bound the accuracy of the hyperplanes and prove
the correctness of the level 2 structure. In Section 4, we build a
certification tool for the theta hyperplanes, which may be of independent
interest in both computing and certifying intersection points acquiring
multiplicities. In Section 5 we show how to certify the level 2 structure, and
in Section 6 we show how to certify the reconstruction formula itself. We
provide detailed outlines of these certifications at the beginning of each of
these sections. Finally, in Section 7 we conclude the proof of the main
theorem by presenting the actual computed bounds needed from previous
sections.

##  2\. Notations, statement of the result, and outline of the proof

We start by establishing our notation, which will be used throughout this
paper unless otherwise specified. Assume that $C$ is a generic complex curve
of genus $g$, and let $K_{C}$ denote the canonical divisor of $C$. As $C$ is
generic, an _odd theta characteristic_ $\theta$ of $C$ is a half-canonical
divisor, $2\theta=K_{C}$, with ${\dim}H^{0}(C,\theta)=1$. By considering the
canonical model of the curve in $|K_{C}|^{*}\cong\mathbb{P}^{g-1}$, we can
identify $\theta$ with a unique hyperplane $H_{\theta}\subset|K_{C}|^{*}$ such
that the intersection $H_{\theta}\cdot C=2\theta$, where $C$ also represents
its canonical model. We call $H_{\theta}$ the _theta hyperplane_ arising from
$\theta$. We denote the homogeneous linear form in $g$ variables that defines
$H_{\theta}$ by $l_{\theta}$.

An important object encapsulating relations between theta hyperplanes is the
_Steiner set_. Denote the group of 2-torsion points of the Jacobian of $C$ by
$\mathrm{Jac}(C)[2]$. Let $0\neq\alpha\in\mathrm{Jac}(C)[2]$. The Steiner set
$\Sigma_{\alpha}$ is composed of pairs ${\theta,\theta+\alpha}$ with both
being odd. Symbolically, we have

(1)
$\Sigma_{\alpha}=\\{\\{\theta,\theta+\alpha\\}:{\dim}H^{0}(C,\theta)={\dim}H^{0}(C,\theta+\alpha)=1\\}.$

See [Dol12, 5.4.2] for further discussion. It is worth noting that if we take
two pairs of odd theta characteristics from a Steiner set and add the four
theta characteristics, we obtain $2K_{C}$. Such quadruples which sum to
$2K_{C}$, are called the _syzygetic_ ones, and the non-syzygetic ones are
called _azygetic_. For such a pair $\theta,\theta+\alpha$, we let the
quadratic form
$q_{\\{\theta,\theta+\alpha\\}}\in|{\mathcal{O}}_{|K_{C}|^{*}}(2)|$ be the
image of $\\{[H_{\theta}],[H_{\theta+\alpha}]\\}\in\mathrm{Sym}^{2}|K_{C}|$
under the map $\mathrm{Sym}^{2}|K_{C}|\to|{\mathcal{O}}_{|K_{C}|^{*}}(2)|$,
namely the projectivization of the product of the linear forms
$l_{\theta}\cdot l_{\theta+\alpha}$. We recall some notation that will furnish
our objects of study. We denote by $\Delta$ the diagonal map:

(2) $\Delta:H^{0}(K_{C}+\alpha)\to\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha).$

Let $\nu_{n}$ be the Veronese map from $\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)$
to $\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)$. We consider also the
following canonically determined map

(3) $\pi_{K_{C}+\alpha}:\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)\to H^{0}(2K_{C}),$

where the subscript signifies the domain of the map. Last but not least, we
will use the $n$-th symmetric power of the map $\pi_{K_{C}+\alpha}$, which we
denote by $\pi^{n}_{K_{C}+\alpha}$. Namely, we have

(4)
$\pi^{n}_{K_{C}+\alpha}:\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)\to\mathrm{Sym}^{n}H^{0}(2K_{C}).$

We display all these maps, as well as some objects to be defined below, in the
following diagram

(5)
${{\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)}}$${{\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)}}$${{V_{n,\alpha}\subset\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C})}}$${{\mathrm{Sym}^{n}H^{0}(2K_{C})\supset
U_{n,\alpha}}}$${{H^{0}(K_{C}+\alpha).}}$$\scriptstyle{{\pi^{n}_{K_{C}+\alpha}}}$$\scriptstyle{\nu_{n}}$$\scriptstyle{\Delta}$$\scriptstyle{{\pi^{n}_{K_{C}}}}$

Finally, let $I_{n}(\Delta(V))$ be the ideal of degree $n$ forms containing
the image $\Delta(V)\subset\mathrm{Sym}^{2}V$ of $V$ under the diagonal map
inside $\mathrm{Sym}^{2}V$. This is the ideal corresponding to the Veronese
image of this object, namely $\nu_{n}(\Delta(V))$. We are interested in the
case when $V=H^{0}(K_{C}+\alpha)$. We now introduce our main object of study.

######  Definition 2.1.

Let
$U_{n,\alpha}:=\pi^{n}_{K_{C}+\alpha}\circ\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha))$.
Denote by $V_{n,\alpha}$ the linear span of the inverse image of
$U_{n,\alpha}$ in ${\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C})}$ under the
map $\pi^{n}_{K_{C}}$ as shown in diagram (5).

######  Definition 2.2.

We say that the curve $C$ has properties A, B, C with respect to $n$ if the
following statements hold respectively.

  * (A)

The codimension of $V_{n,\alpha}$ is

(6) $d_{\alpha}:=\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha))-\dim\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)+\dim\mathrm{Sym}^{n}H^{0}(2K_{C}).$

  * (B)

The scheme $V_{n,\alpha}$ is spanned by the images of
$q_{\\{\theta,\theta+\alpha\\}}$.

  * (C)

The intersection
$\cap_{\alpha\in\mathrm{Jac}(C)[2]\smallsetminus\\{0\\}}V_{n,\alpha}$ is the
kernel of $\pi^{n}_{K_{C}}$.

Loosely speaking, these statements determine whether the pairs of theta
hyperplanes are in _general linear position_ , up to known dimension
constraints. This makes these claims amenable to both specialization and
computational techniques. We make the first of these two observations precise
below.

######  Lemma 2.3.

The codimension of $V_{n,\alpha}$ is at least (6).

###### Proof.

We have the following equality by definition

$\dim\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha))=\dim\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)-\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha)).$

This implies that

$\dim\pi_{K_{C}+\alpha}\circ\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha))\leq\dim\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)-\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha)).$

Then, by flipping the inequality, we have

$\displaystyle\mathrm{codim}\pi_{K_{C}+\alpha}\circ$
$\displaystyle\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha))\geq$
$\displaystyle\dim\mathrm{Sym}^{n}H^{0}(2K_{C})-\left(\dim\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)-\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha))\right).$

As the pull back preserves the codimension, the following codimensions are the
same:

$\mathrm{codim}{(\pi^{n}_{K_{C}}})^{-1}\pi_{K_{C}+\alpha}\circ\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha))=\mathrm{codim}\pi_{K_{C}+\alpha}\circ\nu_{n}\circ\Delta(H^{0}(K_{C}+\alpha)),$

which concludes the proof. ∎

Property B is a general position property since it describes when a finite
number of points in a space spans this space. Property C can also be seen as a
general position property: Given some spaces $V_{n,\alpha}$ that are strictly
contained in $\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C})$ and are all
containing some subspace $V$, the property asserts that intersection
$\cap_{\alpha\in\mathrm{Jac}(C)[2]\smallsetminus\\{0\\}}V_{n,\alpha}$ is
merely $V$.

######  Remark 2.4 ($n=1$).

In this case, we have $\dim I_{n}(\Delta(H^{0}(K_{C}+\alpha))=0$, and so the
codimension (6) in Property A becomes

$\dim H^{0}(2K_{C})-\dim\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)=3g-3-g(g-1)/2.$

This quantity is non-positive when $g>5$. This agrees with the map

(7) $\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)\to H^{0}(2K_{C})$

being surjective, which has been known for genus $g>6$ since [LS96], and the
case of genus 6 is the first case of the Prym-Green conjecture. This case was
settled in [CEFS13].

For $g=3$ the content of Properties A, B for a generic curve has been known to
Coble [Cob29, Chapter IV Section 47 (6)], [ACGH85, Chapter VI exercises
F15-F23]. If one takes a smooth plane quartic $C$ and a Steiner set
corresponding to some non trivial $\alpha\in\mathrm{Jac}(C)[2]$, which is
formed by six pairs of bitangents to $C$ in this case, then the six
intersection points of the pairs of bitangents sit on a conic in the same
plane as $C$, and the double cover of the conic ramified along these six
points is the Prym curve corresponding to $C,\alpha$. Note that in this case
Property C is trivial as the quartic does not sit on any conic.

For $g=4$, the generic case was proved in [Leh15] over the complex numbers,
and later [ÇKRS19] extended the construction that has been carried out in
[Leh15] to space sextics on singular quadrics and over more general fields.

The generic $g=5$ case was proven over the complex numbers in [Leh22].

######  Proposition 2.5.

For $n=2$ the following equalities hold:

$\displaystyle\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha))=\binom{g-1}{2}+(g-1)\binom{g-2}{2}+2\binom{g-1}{4},$
$\displaystyle\dim\mathrm{Sym}^{2}H^{0}(2K_{C})=\binom{3g-2}{2},\quad\dim\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)=\binom{\binom{g}{2}+1}{2}.$

###### Proof.

Only the first equality is not immediate. Suppose $X_{1},\dots,X_{g-1}$ be a
basis of sections for the space $H^{0}(K_{C}+\alpha)$, then the monomials
$(X_{i}X_{j})^{2}$, $X_{i}^{2}X_{j}X_{k}$, $X_{i}X_{j}X_{k}X_{l}$ form a basis
of sections of the linear system that corresponds to the ideal $I_{n}$. This
follows also from [Dol03, Remark 5.1]. ∎

######  Remark 2.6.

Property C fails for $g<6$ as the tacit surjectivity assumption of the map (7)
fails in these cases. Computing the dimensions above for genus $6,7,8,9$, the
codimensions $d_{\alpha}$ from Equation 6 in Property A become 50, 45, 21, -30
respectively. Thus, our codimension computation is absurd already for genus 9.
We also compute the increasing sequence of the dimensions of the spaces ${\rm
Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})$ for genus $6,7,8$:

$\binom{\binom{6+1}{2}+1}{2}=231,\quad\binom{\binom{7+1}{2}+1}{2}=406,\quad\binom{36+1}{2}=666.$

On the other hand, the numbers of pairs of theta hyperplanes arising from a
Steiner set in genera $6,7,8$ are $496,2016,8128$ respectively. Hence, in
principle, Properties A, B, C may all hold for genera $6,7,8$.

We now present the main result of the current study:

######  Theorem 2.7.

Properties A, B, C hold for $n=2,g=6,7$.

We now describe our approach to prove this theorem using ideas from [Leh15,
Leh22], [ÇKRS19]. Indeed, we will make use of the following lemma. Later, in
Remark 2.11, we discuss key differences between the current work and these
previous ones beyond the obvious technical difficulty of moving from $n=1$ to
$n=2$ in property C.

######  Lemma 2.8.

Let ${\mathcal{V}}/X$ be a vector bundle over a base $X$, and let
${\mathcal{V}}_{1},\dots\mathcal{V}_{n}$ be sub-bundles of ${\mathcal{V}}$.
Then the function
$\dim\langle{\mathcal{V}}_{1}|_{x},\ldots,{\mathcal{V}}_{n}|_{x}\rangle$ is
lower semi-continuous on $X$, and the function
$\dim(\cap_{i=1}^{n}{\mathcal{V}}_{i}|_{x})$ is upper semi-continuous on $X$.

###### Proof.

Classical – see e.g. [Leh15, Corrolary 7]. ∎

### The ingredients of the proof of Theorem 2.7

We start by arguing that Theorem 2.7 follows from the statement in the theorem
for a _specific_ curve thanks to Lemma 2.8. Indeed, each
$q_{\\{\theta,\theta+\alpha\\}}$ generates a rank 1 sub-bundle of
$V_{2,\alpha}$ over the family of genus $g$ curves $C$. So it suffices to show
that the dimension of the span of the image of
$q_{\\{\theta,\theta+\alpha\\}}$ is the same as the dimension of
$V_{2,\alpha}$ for a particular curve $C$, and similarly, it is enough to show
that given some subset
$\\{\alpha_{i}\\}\subset\mathrm{Jac}(C)[2]\smallsetminus\\{0\\}$ the
intersection $\cap_{i}V_{2,\alpha_{i}}C$ has codimension at least
$\dim\mathrm{Sym}^{2}H^{0}(2K_{C})$ for a particular curve. For $g=6$ and $7$,
we use Wiman’s sextic and the Fricke-Macbeath curve respectively to carry out
the proof. The examples are presented in Example 3.2 and Example 3.3.

In our computations, we utilize numerical computational tools, which require
certification to conclude the proof. We first computed a Riemann matrix for
each of our examples via SageMath [NBZ19]. This machinery is based on
certified methods. Using the Riemann matrix, we estimated the theta
hyperplanes in terms of the theta functions, as in (11), via the SageMath
library [GB21]. These linear forms come with their labels, which identify the
level 2 structure. Namely, they form the syzygetic quadruples, which are
defined in Equation 1. However, both the linear forms and level 2 structure
lack certification. Our certification task thus breaks into the following
parts:

  1. (1)

In Section 4 we build a certification tool for the theta hyperplanes. This
relies on having a non-certified theta hyperplane approximation, together with
non-certified approximations of the tangency points, which we computed with
`Julia homotopy continuation` [BT17]. As already mentioned in Section 1, this
certification tool may be of independent interest in both computing and
certifying intersection points acquiring multiplicities – see Lemma 4.22,
Remark 4.34 Item 4 and Section 7.6 Item 3.

  2. (2)

In Section 5, given a non-certified level 2 structure and the certified
intersection points of the certified hyperplanes, we build a certification
tool for subsets of some of the Steiner sets.

  3. (3)

In Section 6 we build a certification tool to certify the singular value
decomposition (SVD) computation on the 2nd powers of the quadratic forms
$q_{\theta,\theta+\alpha}$ i.e., the products of the linear forms representing
the certified hyperplanes. This is to certify that the dimension $d_{\alpha}$
in (6) is a lower bound of $\dim V_{2,\alpha}$, which is needed to prove
Property B. The upper bound, which coincides with this lower bound, is proved
in Lemma 2.3. We then discuss the certification of the inequality
$\mathrm{codim}\cap V_{2,\alpha}\geq\mathrm{codim}\ker(\pi^{2}_{K_{C}})$,
which is needed to prove Property C using the SVD on the result of the
previous SVD. The lower bound is provided by computing
$\mathrm{codim}\ker\pi^{2}_{K_{C}}$, which is the difference in the dimensions
of the domain and target spaces.

  4. (4)

In Section 7 we discuss the numerical errors appearing in practice in the
input, and the ones incurred by performing numerical linear algebra
computations. We then apply the certification tools we have built in the
previous sections to conclude the proof of Theorem 2.7. Finally, we also
discuss in this section the testing of the certification code.

######  Corollary 2.9.

For a generic curve of genus 6 or 7, Theorem 2.7 effectively determines the
curve.

###### Proof.

By Theorem 2.7, we can recover the kernel of the map $\pi^{2}_{K_{C}}$ defined
from $\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})$ to
$\mathrm{Sym}^{2}H^{0}(2K_{C})$. Since we work over the complex curves, we can
choose a basis of the kernel, and complete it to a basis of
$\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})$. We thus identify the
complement of the kernel with $\mathrm{Sym}^{2}H^{0}(K_{C})$. Note that we
prefer to work with unitary bases of these spaces, and unitary complements
instead of general ones for the sake of numerical stability. So, in essence,
we have a decomposition

$\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})=\mathrm{Sym}^{2}I_{2}(C)\oplus
I_{2}(C)\otimes H^{0}(2K_{C})\oplus\mathrm{Sym}^{2}H(2K_{C}),$

and we attempt to recover the first summand when we know the last summand. To
this end, we observe that

$\displaystyle\mathrm{Sym}^{2}(\mathrm{Sym}^{2}I_{2}(C)\oplus I_{2}(C)\otimes
H^{0}(2K_{C}))\cap(\mathrm{Sym}^{2}H^{2}(2K_{C})\otimes(\mathrm{Sym}^{2}I_{2}(C)\oplus
I_{2}(C)\otimes H^{0}(2K_{C})))$
$\displaystyle=\mathrm{Sym}^{2}I_{2}(C)\otimes\mathrm{Sym}^{2}H^{0}(2K_{C}).$

Having recovered
$\mathrm{Sym}^{2}I_{2}(C)\otimes\mathrm{Sym}^{2}H^{0}(2K_{C})$ one recovers
$\mathrm{Sym}^{2}I_{2}(C)$ by decomposing in a tensor basis of the spaces
$\mathrm{Sym}^{2}(I_{2}(C)\oplus H^{0}(2K_{C}))$ and
$\mathrm{Sym}^{2}(H^{0}(2K_{C}))$. Finally, note that a generic canonical
curve $C$ is cut out by quadrics unless it is trigonal or a plane quintic. The
quadrics generate the ideal $I_{2}(C)$, and we have recovered
$\mathrm{Sym}^{2}I_{2}(C)$. ∎

Given a principally polarized abelian variety, one can define the theta
hyperplanes via the Gauss map. In explicit terms, given a period matrix that
parametrizes the abelian variety in the Siegel upper half space, the theta
hyperplanes are expressed in terms of the Riemann theta function.

Using the notation introduced at the beginning of Section 2, we state a
description of the Schottky locus, which is the locus of Jacobians in the
moduli space of abelian varieties.

######  Corollary 2.10.

Let $g$ be $6$ or $7$. The Schottky locus in genus $g$ is the closure of the
principally polarized Abelian varieties $(A,\Theta)$ such that the
intersection

$\bigcap_{\alpha\in
A[2]\smallsetminus\\{0\\}}\mathrm{span}\\{\nu_{2}(q_{\theta,\theta+\alpha})\\}_{\\{\theta,\theta+\alpha\\}\in\Sigma_{\alpha}}$

is a non-reduced scheme of multiplicity $2^{g-1}$, whose reduction is smooth
genus $g$ curve. The associated Jacobian variety is $(A,\Theta)$, and the
theta hyperplanes are the theta hyperplanes of $A$ with the same level $2$
structure.

######  Remark 2.11.

The construction of curves in genus 4 [Leh15] was completely algebraic,
moreover, unlike [Leh22, ÇKRS19], the proofs were not computer-aided. In
[Leh22], although the proof was computer-aided, a witness was provided that
furnishes the proof of the main result, which could be verified by rank
computations on non-mechanized scale. Moreover, as the theta hyperplanes were
found algebraically, the proof and also the witness verification could have
been done over any field, with the exception of some bad characteristics.

In the current study, we deviate from this pleasant scheme in two manners:
First, we find the theta hyperplanes by approximate analytic methods via theta
functions, which is presented in Section 3. In particular, we do not know over
which number field they are defined. Second, the techniques we use to certify
our numerical computations involve large-scale linear algebra computations.
Nevertheless, we stress that Theorem 2.7 could be stated over any field of
characteristic other than 2.

Considering applicative aspects, our work is relevant to constructing
algebraic curves from their period matrices. Recently, polygons have been used
as a source of the periods for the construction [ÇFM23].

We close this section by outlining some directions to pursue in future.

######  Question 2.12.

As pointed in Remark 2.6, the dimension counts break down for genus strictly
greater than 8. It would be pointless to consider $I_{n}(\Delta(V))$ for
$n>2$, as this ideal is generated in degree $2$ for any $n$. One might attempt
instead to use syzygies for the ideal, namely relations among the $2$-forms.

######  Question 2.13.

Assuming the claims are correct in the case of genus 8, our certification
procedure would work. There are two challenges in tackling this problem: The
extensive computations involved, and the lack of an obvious candidate curve to
apply our procedure on. To tackle the second issue one may follow [Muk88] to
explicitly construct a general curve of genus 8. The Macalauy 2 package [Sch]
features a function for the construction.

######  Question 2.14.

For numerical reasons we want our examples to have widely spread theta
hyperplanes in their ambient space, as well as distinct enough tangency
points. We also suspect that our curves should be Brill-Noether general to
satisfy the main theorem. Vaguely speaking, non-hyperelliptic curves with many
automorphisms tend to have all these properties. In addition, the
automorphisms provide a useful perspective for the part of our computation of
approximating the tangency points. Therefore, it would be beneficial to have
such an example for the next genera.

######  Question 2.15.

Our proof relies extensively on numerical and symbolic computations. The
numerical part necessitates certification, which adds more to this scheme. As
the underlying genus increases, the complexity of these computations grows
exponentially, making the process even more intricate for higher genera. In
short, is there some conceptual argument which proves our main result?

##  3\. Computing theta hyperplanes numerically

We begin this section by recalling the Riemann theta function and its
connections to theta hyperplanes. Our primary reference for the background
material is [Mum07]. We also mention their numerical computations with current
trends, which provide tools for the computational part of this work. Last but
not least, we will present examples of the curves in genus 6 and 7, which are
cruxes of the proof of the main theorem, namely Wiman’s sextic and the Fricke-
Macbeath curve respectively.

Let $C$ be a complex smooth algebraic curve. Let $\omega_{1},\dots,\omega_{g}$
be a basis of $H^{0}(K_{C})$, i.e., the space of holomorphic differential one-
forms. Let $\alpha_{1},\dots,\alpha_{g},\beta_{1},\dots,\beta_{g}$ be a
symplectic basis of $H_{1}(C,\mathbb{Z}$). The following $g\times 2g$ matrix
is called the _period matrix_ of the curve (with the given choice of basis):

(8)
$(\tau_{1}|\tau_{2}):=\left(\left(\int_{\alpha_{i}}\omega_{j}\right)\Big{|}\left(\int_{\beta_{i}}\omega_{j}\right)\right),$

whereas $\tau:=\tau_{1}^{-1}\tau_{2}$ is called a _Riemann matrix_ of the
algebraic curve $C$. The Riemann matrix lies in the set of $g\times g$
symmetric matrices with positive definite imaginary parts with complex
entries, the so-called Siegel upper half space $\mathbb{H}_{g}$. The theta
function with characteristic $\varepsilon,\delta\in\\{0,1\\}^{g}$ is a
complex-valued function defined on $\mathbb{C}^{g}\times\mathbb{H}_{g}$:

(9) $\theta\begin{bmatrix}\varepsilon\\\ \delta\end{bmatrix}({\bf
z}\,|\,\tau)\,\,\,=\,\,\,\sum_{{\bf n}\in\mathbb{Z}^{g}}{\rm exp}\left(\pi
i\left({\bf n}+\frac{\varepsilon}{2}\right)^{T}\tau\left({\bf
n}+\frac{\varepsilon}{2}\right)+\left({\bf
n}+\frac{\varepsilon}{2}\right)^{T}\left({\bf
z}+\frac{\delta}{2}\right)\right).$

When $\varepsilon=\delta=0$, this is nothing but the _Riemann theta function_
, and differs by an exponential factor from the latter. The characteristic is
called even, odd if $\varepsilon\cdot\delta\equiv 0,1\pmod{2}$ respectively.
So there are $2^{g-1}(2^{g}+1)$ odd and $2^{g-1}(2^{g}+1)$ even
characteristics. For fixed $\tau$, the values
$\theta\begin{bmatrix}\varepsilon\\\ \delta\end{bmatrix}({\bf 0}\,|\,\tau)$ at
${\bf z}={\bf 0}$ are known as _theta constants_. We will also call theta
constant the evaluation of the first-order derivatives of the theta function
at ${\bf 0}$:

(10) $\theta^{\varepsilon,\delta}_{i}:=\frac{\partial}{\partial
z_{i}}\theta\begin{bmatrix}\varepsilon\\\ \delta\end{bmatrix}({\bf
z}\,|\,\tau)_{\big{|}_{{\bf z}={\bf 0}}}.$

The following statement is classical, which follows from the Riemann
singularity theorem. It formulates the theta characteristic divisors in terms
of the canonical model [ACGH85, Page 228]. See for its trim formulation
[Guà02, Theorem 2.2].

######  Theorem 3.1.

Let $\tau_{1}$ be the first $g\times g$ part of the period matrix (8). Let
$\theta$ be an effective theta characteristic divisor of degree $g-1$ with
$\dim H^{0}(\theta)=1$. The equations of the theta hyperplanes $H_{\theta}$
arising from $\theta$ on the canonical model are given as follows:

(11)
$\begin{pmatrix}\theta^{\varepsilon,\delta}_{1}&,\dots,&\theta^{\varepsilon,\delta}_{g}\end{pmatrix}\cdot\tau_{1}^{-1}\cdot\begin{pmatrix}X_{1}\\\
\vdots\\\ X_{g}\end{pmatrix}=0,$

where the characteristic $\begin{bmatrix}\varepsilon\\\ \delta\end{bmatrix}$
ranges over the odd ones.

Some mathematical software packages to numerically compute the Riemann
matrices and the theta functions are available as [NBZ19, Dv01, FK11, MN19,
vW06] and [GB21, AC21, DHB+04, FJK19] respectively. We used the SageMath
packages [NBZ19, GB21] to carry out our computations which involve theta
functions, in particular evaluating the theta hyperplanes in (11).

We close this section by presenting the examples of genus 6 and 7 curves that
lie in the heart of the proof of the main result.

######  Example 3.2 (Wiman’s sextic).

We choose Wiman’s sextic [Wim96] for our proof in genus 6. Its plane equation
is given as:

(12)
$x^{6}+y^{6}+z^{6}+(x^{2}+y^{2}+z^{2})(x^{4}+y^{4}+z^{4})-12x^{2}y^{2}z^{2}=0.$

Its normalization defines a smooth model in $\mathbb{P}^{5}$, whose
automorphism group becomes the symmetric group $\mathbf{S}_{5}$. Via the
SageMath packages [NBZ19, GB21], we computed a period matrix for the curve
numerically using the plane equation in (12), and then computed
$2^{5}(2^{6}-1)$ many linear forms that shall approximate the linear forms
defining the theta hyperplanes as in (11) in finite-digit precision. We used
Julia homotopy continuation [BT17] to approximate the tangency points of the
hyperplanes, which served us to certify these computations. For more on the
certification matters, see Section 4 and Section 5.

######  Example 3.3 (Fricke–Macbeath curve).

For genus 7, our sample is the Fricke-Macbeath [Fri99, Mac65] curve. We learnt
its plane equation from [Hid18], who attributes it to Bradley Brock:

$z^{8}+7xyz^{6}+21x^{2}y^{2}z^{4}+35x^{3}y^{3}z^{2}+28x^{4}y^{4}+2x^{7}z+2y^{7}z=0.$

This is the unique Hurwitz curve of genus 7, namely, its automorphism group is
of order $84(7-1)$. Similar to what we did in Example 3.2, we computed the
8128 linear forms in terms of the estimated theta constants. With these
approximated theta hyperplanes, we ran our algorithm to construct their curve
after certifying underlying numerical computations.

##  4\. Numerically certifying theta hyperplanes

Given a complex algebraic curve, we estimate the theta hyperplanes numerically
using theta functions via the SageMath library [GB21], and determine the
tangency points through homotopy continuation using Julia [BT17]. Both of
these estimation methods are not certified inherently. In this section, we are
given a curve which is an intersection of quadrics in some ambient space. In
practice, this space is given by the dual canonical system, which is not
required from the point of view of this section. Under this setting, we
present a certification process for both the hyperplanes and the points. The
process deals with the two different types of numerical inaccuracies discussed
as follows. The hyperplanes in question are merely putative approximations of
the theta hyperplanes. While denoting the theta hyperplanes by $H$, we denote
$H^{\prime}$ as the approximations of them. The input points are only
approximations of the intersection points of the input hyperplanes and the
curve. We denote the intersection points of $H^{\prime}\cap C$ by $P_{i}$ and
the approximations by $P^{\prime}_{i}$. It is noteworthy to mention that we do
not introduce a notation for the exact multitagency points $H\cap C$. These
points are not used in our certification scheme. The fact that we do not have
exact hyperplanes is the more challenging one to deal with. We argue this
briefly as follows. Under the assumption that the tangency points are not
flexes, which we verify, we locally approximate the curve about each
approximate tangency point $P_{i}$ with a conic. We refer to this with _conic
approximation_. The analysis of this process enables us to establish bounds on
the accuracy of the hyperplanes; thus ensuring their numerical stability. The
plan above presents three problems; the first two originate from numerical
inaccuracies in the input hyperplanes $H^{\prime}$, and the last from
numerical inaccuracies in the input points $P^{\prime}_{i}$. First, we need to
find a conic approximation and prove that it is adequate. Second, we need to
prove that having a hyperplane “close to being _multitangent_ ” to the conic
approximation means that the “candidate” of a theta hyperplane is indeed close
to the theta hyperplane. Finally, we need to address numerical issues arising
from the inaccuracies of the points.

The section is organized as follows. We start by recalling classical numerical
linear algebraic results. We then construct the local conic approximations of
the curve in Lemma 4.4, and bound its accuracy in Corollary 4.10 and Corollary
4.20. We use these bounds to build a certification tool to certify the theta
hyperplanes under the assumption of no numerical errors in the intersection
points of the putative hyperplane and the curve. That is to say, we assume
that we are given the points $P_{i}$ and not the points $P^{\prime}_{i}$,
which form the content of Lemma 4.22. Only then, in Lemma 4.32, we build the
certifying tool for the theta hyperplanes with computational errors in the
points $P^{\prime}_{i}$.

Our analysis requires heavy usage of the properties of the SVD computations.
We begin by recalling one classical result on singular values as well as two
elementary propositions that serve us as preliminaries. Unless otherwise
specified, we will denote the singular values of a matrix by $\sigma_{i}$
which will be assumed to be in decreasing order. In what follows, we use the
operator norm and the Frobenius norm for the matrices, denoted by $|\cdot|$
and $||\cdot||_{F}$ respectively. The same notation $|\cdot|$ will be used
also for the Euclidean norm, which will be distinguishable from the objects.

######  Theorem 4.1.

For any two matrices of the same dimensions $A,E$, we have
$|\sigma_{i}(A+E)-\sigma_{i}(A)|\leq\sigma_{1}(E)=|E|$.

######  Proposition 4.2 (Pseudo kernels are continuous).

Let $A,E$ be matrices of the same dimensions. Suppose that there is some $k$
such that $\sigma_{i}>a$ for $i<k$, and $\sigma_{k}<\epsilon$ for some
$\epsilon$. Assume that the projection of a vector $v$ on the span of the
first $k-1$ singular vectors of $A$ is trivial, then its projection on the
first $k-1$ singular vectors of $A+E$ has norm at most
$\frac{\epsilon+|E|}{a-|E|}|v|$.

###### Proof.

We have $|(A+E)v|\leq(\epsilon+|E|)|v|$. Since the $i$th singular value of
$A+E$ is at least $\sigma_{i}-|E|>a-|E|$, the statement follows. ∎

######  Proposition 4.3 (Solutions of linear equations are continuous).

Suppose that $A$ is an invertible matrix and $Ax=y$, and
$(A+E)x^{\prime}=y+\delta_{y}$, where the smallest singular value of $A$
satisfies $\sigma_{r}(A)>a$. Then
$|x^{\prime}-x|<(|E|/a+\epsilon^{\prime})/(a-|E|)$.

###### Proof.

Setting $B=A+E$, we have
$A^{-1}-B^{-1}=B^{-1}BA^{-1}-B^{-1}AA^{-1}=B^{-1}(B-A)A^{-1}$. Whence

$A^{-1}y-(A+E)^{-1}(y+\delta_{y})=(A+E)^{-1}EA^{-1}y-(A+E)^{-1}\delta_{y}.$

This implies that
$|A^{-1}y-(A+E)^{-1}(y+\delta_{y})|\leq\frac{|E|\cdot|y|}{a(a-|E|)}+\frac{|\delta_{y}|}{a-|E|}.$
∎

### Construction the conic approximation to $C$ at a point $P$

We assume that $C$ is an intersection of quadrics $\\{Q_{i}\\}_{i=1}^{\dim
I_{2}(C)}$ in $|K_{C}|^{*}\cong\mathbb{P}^{g-1}$, where $I_{2}(C)$ is the
linear system of quadrics on $|K_{C}|^{*}$ that contain $C$. We will abuse
notations and denote the $g\times g$ symmetric matrices corresponding to the
quadrics by $Q_{i}$ as well. Then $C$ is given by the null set of
$P\mapsto\\{P^{t}Q_{i}P\\}_{i}$. In what follows, we work with projective
coordinates and represent points as column vectors of size $g$.

######  Lemma 4.4 (The conic approximation of an intersection of quadrics).

Let $P$ be a non flex point on $C$, and let $T,R$ be so that $P+xT+x^{2}R$ is
a conic approximation of $C$ at $P$ then $T,R$ satisfy

(13) $0=2P^{t}Q_{i}T,\qquad 0=2T^{t}Q_{i}T+P^{t}Q_{i}R,\qquad\text{ such that
}T,R\neq 0.$

###### Proof.

The point $P$ lies on each quadric $Q_{i}$. Hence $T$ and $R$ are solutions of
the following equations for any $i$:

$\displaystyle 0=\frac{d}{dx}(P+xT)^{t}Q_{i}(P+xT)|_{x=0},\quad
0=\frac{d^{2}}{dx^{2}}(P+xT+x^{2}R)^{t}Q_{i}(P+xT+x^{2}R)|_{x=0}.$

These equations imply the equations in (13). Finally, $R\neq 0$ since $P$ is
not a flex point. ∎

######  Definition 4.5.

Since our points reside in $\mathbb{P}^{g-1}$, we may assume $|P|=1$.
Moreover, we will take $T$ to be orthogonal to $P$ with $|T|=1$, and consider
$R$ to be the smallest norm solution of the linear system arising from (13).
These last requirements may be formulated as follows: First, denote by $M$ the
$\dim I_{2}(C)\times g$ matrix whose rows are ${\\{P^{t}Q_{i}\\}}_{i}$. We
note that the co-rank of $M$ is 2. Indeed, if there were other vectors in the
kernel than $P,T$ then $P$ would not be a smooth point of $C$. To solve $T$,
we look at the kernel of the $(\dim I_{2}(C)+1)\times g$ matrix
$M_{T}:=\left(\begin{smallmatrix}M\\\ P^{\dagger}\end{smallmatrix}\right)$,
where $\dagger$ denotes the complex transpose operator. To solve $R$, denote
$M_{R}:=\left(\begin{smallmatrix}M_{T}\\\
T^{\dagger}\end{smallmatrix}\right)$, and solve the equation

(14) $\quad M_{R}R=-\frac{1}{2}\left(\begin{smallmatrix}T^{t}Q_{1}T\\\
\vdots\\\ T^{t}Q_{\dim I_{2}(C)}T\\\ 0\\\ 0\end{smallmatrix}\right).$

We now introduce a couple of approximation results for points close to our
curve $C$. The need for these results is twofold. The first is that we are
about to express the curve locally as a “conic plus something small”. The
second is that we will have to bound the effect of computational errors to
certify our computations, as in practice, these are our input points. To this
end, we will assume that
$\sqrt{\sum_{i}|{\tilde{P}}^{t}Q_{i}\tilde{P}|^{2}}<\epsilon$ for some
${\tilde{P}}\in\mathbb{C}^{g}$. For our purposes, this ${\tilde{P}}$ will
either be an approximation of some point $P$ by means of conic approximation,
or a computed version of a point $P\in C$. We let $\tilde{M}$ be the matrix
whose rows are ${\tilde{P}}^{t}Q_{i}$ with the computed ${\tilde{P}}$. As
usual, we denote the $j$th singular value of a matrix by $\sigma_{j}$.

Recall that the curve $C$ is a local complete intersection. This statement can
be made effective as follows: The row vectors $P^{t}Q_{i}$s span the gradients
of the quadrics in $I_{2}(C)$ at a point $P$. Hence, if we consider the matrix
built of these row vectors then its top $g-2$ singular values are the only
non-trivial ones. Also, the corresponding singular vectors of this matrix are
the gradients of quadrics in $I_{2}(C)$ which span the normal space of the
curve at $P$. We use this observation to set up the proof of Lemma 4.6 as
follows. Considering a point $\tilde{P}$, let
$\tilde{Q}_{1},\ldots,\tilde{Q}_{g-2}$ be the quadrics corresponding to the
top $g-2$ singular vectors of $\tilde{M}$ at the point $\tilde{P}$. We also
set $\tilde{q}_{m}:=\sqrt{\sum_{i=1}^{g-2}|\tilde{Q}_{i}|^{2}}$, where the
norm is the operator norm.

######  Lemma 4.6.

In the setting above, let $\delta_{\tilde{P}}$ be the vector whose coordinates
are $({\tilde{P}}^{t}\tilde{Q}_{i}\tilde{P})_{i}$, then there is a point $P\in
C$ such that $|\tilde{P}-P|<r$ where $r\in\mathbb{R}$ is such that
$r<(|\delta_{\tilde{P}}|-\tilde{q}_{m}r^{2})/2\sigma_{g-2}(\tilde{M})$.

The proof emerges in projective coordinates where $\tilde{P}$ is normalized
but without normalizing $P$. So the distance measured is slightly bigger from
the Fubini-Study metric, which is defined as ${\arccos(|P_{1}\cdot P_{2}|)}$
for points $P_{1},P_{2}$ in
$(\mathbb{C}^{n}\smallsetminus\\{0\\})/\mathbb{R}^{+}$. The latter is an
$\mathbb{S}^{1}$ bundle over $\mathbb{P}^{n-1}$. Instead, we use the value
$1-|P_{1}\cdot P_{2}|$. This is diffeomorphic to the actual Fubini-Study, and
is faster to compute, as we minimize the distance in the $\mathbb{C}^{*}$
bundle over the projective space, and the points $|P|$ are taken to be of unit
norm.

###### Proof.

Consider the map

$\displaystyle B(0,r)^{g-2}$ $\displaystyle\to\mathbb{C}^{g-2}$ $\displaystyle
z$
$\displaystyle\mapsto(\tilde{Q}_{1}(z),\ldots,\tilde{Q}_{g-2}(z))=((\tilde{P}+z)^{t}\tilde{Q}_{i}(\tilde{P}+z))_{i}=\delta_{\tilde{P}}+2\tilde{M}z+(z^{t}\tilde{Q}_{i}z)_{i}.$

Thus, we wish to solve the equation
$0=\delta_{\tilde{P}}+2\tilde{M}z+(z^{t}\tilde{Q}_{i}z)_{i}$. As our map is
open, and each of the coefficients in the degree $2$ tail is bounded by
$r^{2}|\tilde{Q}_{i}|$, the result follows. ∎

Note that in the proof above we merely bound the norm, and thus we do not have
to move to real coordinates, which we would have had to do if we wanted to
explicitly refer to the closest point. We also remark that we make use of this
technique of using effective forms of the implicit and inverse function
theorem later in this section; e.g. in Proposition 4.21 and Lemma 4.22.

######  Definition 4.7.

Considering complex parameter $x$, we denote $P_{x}:=P+xT+x^{2}R,\quad
T_{x}:=T+2xR$, let $M_{x}$ to be the matrix whose rows are $P_{x}^{t}Q_{i}$,
let $M_{T,x}:=\left(\begin{smallmatrix}M_{x}\\\
P_{x}^{\dagger}\end{smallmatrix}\right)$, and
$M_{R,x}:=\left(\begin{smallmatrix}M_{T,x}\\\
T_{x}^{\dagger}\end{smallmatrix}\right)$. Let $v_{x}$ be the closest point on
$C$ to $P_{x}$. Finally, we denote the vector whose entries are $T^{t}Q_{i}R$
by $W_{T}$, and the vector whose entries are $R^{t}Q_{i}R$ by $W_{R}$. We
refer to Figure 1 for intuitive views of $P,P_{x},v_{x}$ and the tangency
point. Finally we set $q_{m}:=\sqrt{\sum_{i=1}^{\dim I_{2}(C)}|Q_{i}|^{2}}$.
Note that we have $\tilde{q}_{m}|M|<q_{m}$ comparing $q_{m}$ with
$\tilde{q}_{m}$ defined in Definition 4.5.

Figure 1. An illustration of the conic approximation argument for a plane
quartic. Note that both $P,v_{x}$ and the tangency point are on the quartic
curve, and that the conic approximation is at $P$ – not at the tangency point.

######  Proposition 4.8.

$||M_{x}-M||_{F}^{2}$ is bounded above by
$q_{m}^{2}(|x|+|R|\cdot|x|^{2})^{2}$.

###### Proof.

First note that the norm of $T$ is 1. The statement follows from the
definition of the Frobenius norm. ∎

The first of the two approximation tasks mentioned before Lemma 4.6 shows that
the difference between the curve and its conic approximation is small. In
essence, it means that we want to prove that for sufficiently small $|x|$, the
magnitude $|P_{x}-v_{x}|$ is $O(|x|^{3})$, whereas $\frac{d}{dx}|P_{x}-v_{x}|$
is $O(|x|^{2})$. Moreover, we have to give actual constants to achieve a
certified result, and not simply be done with “sufficiently small” and $O$
notation. In practice, this entails finding a radius in which the expressions
are dominated by the leading term in $x$ for each of the expressions involved
and computing the coefficient of this term. Sadly, in this process, the
following complicated expressions are unavoidable.

######  Definition 4.9.

We define the following real numbers

$r_{P,x}:=\min\left(\frac{1}{2\sqrt{|R|}},\frac{\sigma_{g-2}(M)}{2q_{m}},\frac{|W_{T}|}{2|W_{R}|},\left(\frac{\sigma_{g-2}(M)^{2}}{\tilde{q}_{m}|W_{T}|}\right)^{1/3}\right),\quad\epsilon_{P,x}:=|x^{3}W_{T}|/2\sigma_{g-2}(M).$

######  Corollary 4.10.

For $|x|<r_{P,x}$, the distance from $P_{x}$ to $v_{x}$ is bounded by
$\epsilon_{P,x}$.

###### Proof.

By the definition of $P,T,R$ we have:

$P_{x}^{t}Q_{i}P_{x}=(P+xT+x^{2}R)^{t}Q_{i}(P+xT+x^{2}R)=2x^{3}T^{t}Q_{i}R+x^{4}R^{t}Q_{i}R=2x^{3}W_{T}+x^{4}W_{R}.$

We now want to find a requirement on $x$ so that
$|P_{x}-v_{x}|<|x^{3}W_{T}|/2\sigma_{g-2}(M_{x})$. Per Lemma 4.6 it suffices
to find a bound on $x$ so that:

$|x^{3}W_{T}/2\sigma_{g-2}(M_{x})|<\left(|2x^{3}W_{T}|-|x^{4}W_{R}|-|\tilde{q}_{m}(x^{3}W_{T}/2\sigma_{g-2}(M_{x}))^{2})|\right)/2\sigma_{g-2}(M_{x}).$

Subtracting the LHS we get

$0<\frac{|x|^{3}}{2\sigma_{g-2}(M_{x})}\left(|W_{T}|-|xW_{R}|-\tilde{q}_{m}|x|^{3}(|W_{T}|/2\sigma_{g-2}(M_{x}))^{2}\right),$

which holds if

$\frac{|W_{T}|}{2}>|x|\cdot|W_{R}|,\
|x|^{3}\frac{\tilde{q}_{m}|W_{T}|^{2}}{(2\sigma_{g-2}(M_{x}))^{2}}\qquad\Leftrightarrow\qquad|x|<\frac{|W_{T}|}{2|W_{R}|},\left(\frac{2\sigma_{g-2}(M_{x})^{2}}{\tilde{q}_{m}|W_{T}|}\right)^{1/3}.$

The result now follows by bounding $\sigma_{g-2}(M_{x})$: By Theorem 4.1,
together with the first two terms in LABEL:dfn:r_{P and LABEL:x} of $r_{P,x}$
and the estimate for $|M_{x}-M|$ in Proposition 4.8, we have

$|\sigma_{g-2}(M_{x})-\sigma_{g-2}(M)|<q_{m}|x|\sqrt{5}/2<\sigma_{g-2}(M)\sqrt{5}/4.$

By the last two terms in the definition of $r_{P,x}$, we are done. ∎

######  Proposition 4.11.

For $|x|<r_{P,x}$, we have
$||M_{T,x}-M_{T}||_{F}^{2}<q_{m}^{2}(|x|+|R|\cdot|x|^{2})^{2}+\epsilon_{P,x}$.

###### Proof.

It follows from Proposition 4.8 and Corollary 4.10. ∎

######  Definition 4.12.

We reformulate Definition 4.5 for the point $v_{x}$ in the role of $P$: let
$M_{v_{x}}$ be the matrix whose rows are $v_{x}^{t}Q_{i}$, and
$M_{T,v_{x}}:=\left(\begin{smallmatrix}M_{v_{x}}\\\
v_{x}^{\dagger}\end{smallmatrix}\right)$, and
$M_{R,v_{x}}:=\left(\begin{smallmatrix}M_{T,{v_{x}}}\\\
T_{v_{x}}^{\dagger}\end{smallmatrix}\right)$. We denote by $T_{v_{x}}$ a
normalized vector that spans the kernel of $M_{v_{x}}$, i.e., $T_{v,x}$ spans
the tangent space to $C$ at $v_{x}$. Next, we denote by $R_{v_{x}}$ the
solution of the system $2M_{R,v_{x}}\cdot X=-M_{T,v_{x}}T_{v_{x}}$, as well as
the two orthogonality requirements $R_{v_{x}}\perp T_{v_{x}},v_{x}$ similar to
Equation 14.

######  Definition 4.13.

We define the following real numbers

$\displaystyle r_{T,x}$
$\displaystyle:=\min\left(\frac{\sigma_{g-1}(M_{T})}{8q_{m}},\
\left(\frac{\sigma_{g-1}(M_{T})\sigma_{g-2}(M)}{2(q_{m}+2)|W_{T}|}\right)^{1/3},\
\frac{\sigma_{g-2}(M)}{q_{m}+1}\right),$ $\displaystyle\epsilon_{M_{x}}$
$\displaystyle:=q_{m}\epsilon_{P,x},\qquad\epsilon_{M_{T,x}}:=\sqrt{1+q_{m}^{2}}\epsilon_{P,x},\qquad\epsilon_{T,x}:=|x|^{2}\cdot\frac{9|W_{T}|}{\sigma_{g-1}(M_{T})}.$

######  Corollary 4.14.

For $|x|<r_{P,x}$, we have $||M_{T,v_{x}}-M_{T,x}||_{F}<\epsilon_{M_{T,x}}$.

###### Proof.

It follows from the definition of the Frobenius norm together with Corollary
4.10. ∎

######  Corollary 4.15.

For $|x|<r_{P,x},r_{T,x}$, we have $|T_{v_{x}}C-T_{x}|<\epsilon_{T,x}$.

###### Proof.

We start by bounding the distance of $T_{x}$ from the kernel of $M_{T,x}$. Per
the definition of $P,T,R$, we have

$M_{T,x}\cdot
T_{x}=\sum_{i=1}^{g-2}(T+2xR)^{t}Q_{i}(P+xT+x^{2}R)=2x^{2}\sum_{i=1}^{g-2}(T+xR)^{t}Q_{i}R=2x^{2}(W_{T}+xW_{R})$

Hence, the distance of $T_{x}$ from the kernel is bounded by the magnitude of
this quantity divided by $\sigma_{g-1}(M_{T,x}$). Second, the Frobenius norm
of the difference between $M_{T,v_{x}}$ and the matrix $M_{T,x}$ is bounded by
$\epsilon_{M_{T,x}}$. Using Proposition 4.11 and LABEL:dfn:r_{P and LABEL:x}
of $r_{P,x}$ we have
$||M_{T}-M_{T,x}||_{F}^{2}<4q_{m}^{2}|x|^{2}+\epsilon_{P,x}^{2}$ as well as
$|W_{T}+xW_{R}|<2|W_{T}|$ when $|x|<r_{P,x}$. Now we have

$\displaystyle|T_{v_{x}}C-T_{x}|<$
$\displaystyle\epsilon_{M_{T,x}}/(\sigma_{1}(M_{T_{,}x})-\epsilon_{M_{T,x}})+2|x|^{2}\cdot|W_{T}+xW_{R}|/\sigma_{g-1}(M_{T,x})$
$\displaystyle<$
$\displaystyle\frac{\sqrt{1+q_{m}^{2}}\epsilon_{P,x}}{\sigma_{1}(M_{T})-\sqrt{4q_{m}^{2}|x|^{2}+\epsilon_{P,x}^{2}}-\sqrt{1+q_{m}^{2}}\epsilon_{P,x}}+\frac{4|x|^{2}\cdot|W_{T}|}{\sigma_{g-1}(M_{T})-\sqrt{4q_{m}^{2}|x|^{2}+\epsilon_{P,x}^{2}}}$
$\displaystyle<$
$\displaystyle\frac{(1+q_{m})\epsilon_{P,x}}{\sigma_{1}(M_{T})-2q_{m}|x|-2\epsilon_{P,x}-q_{m}\epsilon_{P,x}}+\frac{4|x|^{2}\cdot|W_{T}|}{\sigma_{g-1}(M_{T})-2q_{m}|x|-\epsilon_{P,x}}$
$\displaystyle<$
$\displaystyle\frac{(1+q_{m})\epsilon_{P,x}+4|x|^{2}\cdot|W_{T}|}{\sigma_{g-1}(M_{T})-2q_{m}|x|-(q_{m}+2)\epsilon_{P,x}}.$

If we now assume that $|x|$ satisfies the following inequalities

$\sigma_{g-1}(M_{T})/4\ >\ 2q_{m}|x|,\
(q_{m}+2)\epsilon_{P,x}\qquad\text{and}\qquad(q_{m}+1)\epsilon_{P,x}\ <\
|x|^{2}|W_{T}|/2,$

where the former controls the denominator, and the latter limits
$\epsilon_{P,x}$ in the numerator, then we have
$|T_{v_{x}}C-T_{x}|<\frac{(\frac{1}{2}+4)|x|^{2}|W_{T}|}{\sigma_{g-1}(M_{T})(1-\frac{1}{4}-\frac{1}{4})}=\epsilon_{T,x}$
. However, as $\epsilon_{P,x}=|x^{3}W_{T}|/2\sigma_{g-2}(M)$, these
inequalities are exactly the first, second, and third terms in LABEL:dfn:r_{T
and LABEL:x} of $r_{T,x}$ respectively. ∎

######  Definition 4.16.

We define the following real numbers

$r_{R,x}:=\min\left(\frac{\sigma_{g}(M_{R})}{3\cdot
4q_{m}},\sqrt{\frac{\sigma_{g}(M_{R})\sigma_{g-1}(M_{T})}{3\cdot
36|W_{T}|}},\left(\frac{\sigma_{g}(M_{R})\sigma_{g-2}(M)}{3(2+q_{m})|W_{T}|}\right)^{1/3},\frac{\sigma_{g-1}(M_{T})}{9},\sqrt{\frac{\sigma_{g}(M_{R})\sigma_{g-2}(M)}{2\cdot
2(1+q_{m})}},\frac{\sigma_{g}(M_{R})\sigma_{g-1}(M_{T})}{2\cdot 18}\right)$

and

$\epsilon_{M_{R,x}}:=\sqrt{\epsilon_{M_{T,x}}^{2}+\epsilon_{T,x}^{2}}=\sqrt{(1+q_{m}^{2})\epsilon_{P,x}^{2}+\epsilon_{T,x}^{2}},\qquad\epsilon_{R,x}:=\frac{28|x|\cdot|W_{T}|}{\sigma_{g}(M_{R})}.$

######  Corollary 4.17.

For $|x|<r_{P,x},r_{T,x}$, we have
$||M_{R}-M_{R,x}||_{F}^{2}<4q_{m}^{2}|x|^{2}+\epsilon_{T,x}^{2}+\epsilon_{P,x}^{2}$.

###### Proof.

It follows from the definition of the Frobenius norm together with Corollary
4.10 and Corollary 4.15. ∎

######  Corollary 4.18.

For $|x|<r_{P,x},r_{T,x},r_{R,x}$, we have $|R_{v_{x}}-R_{x}|<\epsilon_{R,x}$.

###### Proof.

As in the proof of Corollary 4.15, we compute two quantities. First, we have
to see how far $R$ is from solving $2M_{R,x}\cdot X=-M_{T,x}T_{x}$. The
difference for the $i$’th entry is

$2R^{t}Q_{i}(P+xT+xR^{2})+(T+2xR)^{t}Q_{i}(T+2xR)=6x(xR+T)^{t}Q_{i}R=6x(W_{T}+xW_{R}).$

Next, $R_{v_{x}}$ solves the system $2M_{R,v_{x}}\cdot
X=-M_{T,v_{x}}T_{v_{x}}$, as well as the two orthogonality requirements
$R_{v_{x}}\perp T_{v_{x}},v_{x}$. Hence, by Proposition 4.3 the difference
between $R_{v_{x}}$ and the solution is bounded by

$\frac{\epsilon_{M_{R,x}}}{\sigma_{g}(M_{R,x})(\sigma_{g}(M_{R,x})-\epsilon_{M_{R},x})}+\frac{\epsilon_{T,x}}{\sigma_{g}(M_{R,x})-\epsilon_{M_{R},x}}.$

Hence, for $|x|<r_{T,x}$ the bound resulting from the two expressions is

$|R_{v_{x}}-R_{x}|<\frac{\epsilon_{M_{R,x}/\sigma_{g}(M_{R,x})}}{(\sigma_{g}(M_{R,x})-\epsilon_{M_{R,x}})}+\frac{\epsilon_{T,x}}{\sigma_{g}(M_{R,x})-\epsilon_{M_{R},x}}+\frac{2\cdot
6|x|\cdot|W_{T}|}{\sigma_{g}(M_{R,x})}.$

As in the two preceding cases for $\epsilon_{P,x}$ and $\epsilon_{T,x}$, we
have to bound all the summands here in terms of the leading coefficients. To
this end we will find bounds on $|x|$ so that all the denominators above are
controlled by the inequality:

$\displaystyle\sigma_{g}(M_{R})/2$
$\displaystyle>\sqrt{4q_{m}^{2}|x|^{2}}+2\sqrt{\epsilon_{T,x}^{2}}+2\sqrt{\epsilon_{P,x}^{2}}+\sqrt{q_{m}^{2}\epsilon_{P,x}^{2}}=2q_{m}|x|+2\epsilon_{T,x}+2\epsilon_{P,x}+q_{m}\epsilon_{P,x}$
$\displaystyle=2q_{m}|x|+2|x|^{2}\cdot\frac{9|W_{T}|}{\sigma_{g-1}(M_{T})}+(2+q_{m})|x|^{3}\cdot|W_{T}|/\sigma_{g-2}(M),$

the quantity $\epsilon_{T,x}$ above is controlled by the inequality

$|x|\cdot|W_{T}|>\epsilon_{T,x}=|x|^{2}\cdot\frac{9|W_{T}|}{\sigma_{g-1}(M_{T})},$

and the quantity $\epsilon_{M_{R,x}}/\sigma_{g}(M_{R})$ above is controlled by
the inequality

$|x|\cdot|W_{T}|\cdot\sigma_{g}(M_{R})>\epsilon_{M_{R,x}}=(1+q_{m})\epsilon_{P,x}+\epsilon_{T,x}=(1+q_{m})|x|^{3}\cdot\frac{|W_{T}|}{2\sigma_{g-2}(M)}+|x|^{2}\cdot\frac{9|W_{T}|}{\sigma_{g-1}(M_{T})}.$

The first of these three inequalities holds by the left three terms in
LABEL:dfn:r_{R and LABEL:x} of $r_{R,x}$, the second inequality by the forth
term, and the third inequality by the last two terms. Hence, we have the
following inequalities

$|M_{R}-M_{R,x}|_{F}+\epsilon_{M_{R,x}}<\frac{1}{2\sigma_{g}(M_{R})}\qquad\text{and}\qquad\frac{2\epsilon_{M_{R,x}}}{\sigma_{g}(M_{R})},\epsilon_{T,x}<|x|\cdot|W_{T}|.$

Thus, under the requirement $|x|<r_{R,x}$, we have
$|R_{v_{x}}-R_{x}|<\frac{(1+1+12)|x|\cdot|W_{T}|}{\sigma_{g}(M_{R})(1-3\cdot\frac{1}{2}\cdot\frac{1}{3})}=\epsilon_{R,x}$.
∎

######  Lemma 4.19.

Assume $C_{1},C_{2}$ are curves in affine complex space, where
$q_{1}(x_{1}),q_{2}(x_{2})$ are local parametrizations of them about
$x_{1,0},x_{2,0}$ respectively. Assume also that the function $x_{2}(x_{1})$
minimizes $|q_{1}(x_{1})-q_{2}(x_{2}(x_{1}))|^{2}$ about
$x_{1,0},x_{2,0}=x_{2}(x_{1,0})$. Then $x_{2}(x_{1})$ is holomorphic and

$|q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime}|^{2}<|q_{1}^{\prime\prime}-x_{2}^{\prime\prime}q_{2}^{\prime}-x_{2}^{\prime}q_{2}^{\prime\prime}|\cdot|q_{1}-q_{2}|.$

###### Proof.

We consider the $0$ contour line of
$\frac{d}{dx_{1},dx_{1}^{\dagger}}|q_{1}(x_{1})-q_{2}(x_{2}(x_{1}))|^{2}$ at
the point $(x_{1,0},x_{2,0})$. As all of the functions here are dependent on
just one variable and its conjugate. We will use the prime symbol to denote
the derivation. Since
$|q_{1}-q_{2}|^{2}=q_{1}^{\dagger}q_{1}+q_{2}^{\dagger}q_{2}-q_{2}^{\dagger}q_{1}-q_{1}^{\dagger}q_{2}$,
and since both $\partial\circ\mathrm{conj}$ and $\overline{\partial}$ are
trivial on holomorphic function, we have

$\frac{d}{dx_{1},dx_{1}^{\dagger}}(|q_{1}(x_{1})-q_{2}(x_{2}(x_{1}))|^{2})=((q_{1}-q_{2})^{\dagger}(q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime}),\text{the
conjugate of that}).$

Hence, $x_{2}$ is indeed holomorphic in $x_{1}$, and the Hessian of the
distance $|q_{1}-q_{2}|$ is

$\mathrm{Hess}:=\left(\begin{matrix}(q_{1}-q_{2})^{\dagger}(q_{1}^{\prime\prime}-x_{2}^{\prime\prime}q_{2}^{\prime}-x_{2}^{\prime}q_{2}^{\prime\prime})\
&(q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime})^{\dagger}(q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime})\\\
\text{conjugate of top right}\ &\text{conjugate of top
left}\end{matrix}\right)$

i.e., to solve $x_{2}(x_{1})$ we have to solve the following system of
differential equations:

$\displaystyle 0=$
$\displaystyle(q_{1}-q_{2})^{\dagger}(q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime}),$
$\displaystyle 0=$
$\displaystyle\det\mathrm{Hess}=|(q_{1}-q_{2})^{\dagger}(q_{1}^{\prime\prime}-x_{2}^{\prime\prime}q_{2}^{\prime}-x_{2}^{\prime}q_{2}^{\prime\prime})|^{2}-|q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime}|^{4}.$

Unless the situation is degenerate, the second equation gives us the
following, as desired:

$|(q_{1}-q_{2})^{\dagger}(q_{1}^{\prime\prime}-x_{2}^{\prime\prime}q_{2}^{\prime}-x_{2}^{\prime}q_{2}^{\prime\prime})|=|q_{1}^{\prime}-x_{2}^{\prime}q_{2}^{\prime}|^{2}.$

∎

######  Corollary 4.20.

For $|x|<r_{P,x},r_{T,x},r_{R,x}$, we have

$\frac{d}{dx}|v_{x}-P_{x}|<\sqrt{\epsilon_{R,x}\epsilon_{P,x}}=\sqrt{14}|x|^{2}|W_{T}|/\sqrt{\sigma_{g-2}(M)\sigma_{g}(M_{R})}.$

###### Proof.

It follows from Lemma 4.19. ∎

### From an approximate theta hyperplane to a hyperplane

In what follows, we apply the conic approximation for several points
simultaneously. Recall the settings related to $P,T,R$ from Definition 4.5,
Definition 4.7 and Definition 4.12 comparing with Figure 1. We normalize $R$
by setting $c:=|R|$ and $\tilde{R}:=R/c$. We consider $n:=g-1$ analytic
$0$-neighborhoods $\mathbb{B}_{i}:=B(0,r_{i})\subset\mathbb{C}$. We define
$g-1$ analytic functions

$\displaystyle f_{i}:\mathbb{B}_{i}$ $\displaystyle\to\mathbb{C}^{n+1},\qquad$
$\displaystyle u_{i}:\mathbb{B}_{i}$ $\displaystyle\to\mathbb{C}^{n+1},$
$\displaystyle x$ $\displaystyle\mapsto
P_{i}+xT_{i}+c_{i}x^{2}\tilde{R}_{i}+u_{i}(x)\qquad$ $\displaystyle x$
$\displaystyle\mapsto v_{i,x}-P_{i,x}$

where $P_{i},T_{i},\tilde{R}_{i}\in\mathbb{C}^{n+1},\ P_{i,x}$ and $v_{i,x}$
are as in Definition 4.5, Definition 4.7 and Definition 4.12 respectively, the
constant $c_{i}=|R_{i}|$, and the function $u_{i}$ satisfies
$e_{i}:=|u_{i}(x)|<a_{i}|x^{3}|$,
$e^{\prime}_{i}:=|du_{i}(x)/dx|<b_{i}|x|^{2}$ for $a_{i},b_{i}$ determined by
Corollary 4.10 and Corollary 4.20 for
$|x|<r_{P_{i},x},r_{T_{i},x},r_{R_{i},x}$. We also define:

(15) $\displaystyle F:\prod^{n}_{i=1}\mathbb{B}_{i}$
$\displaystyle\to\mathbb{C}^{n}$ $\displaystyle(x_{1},\ldots,x_{n})$
$\displaystyle\mapsto\left(\nabla_{x_{j}}f_{j}(x_{j})\wedge\bigwedge_{i=1}^{n}f_{i}(x_{i})\right)_{j}=(\det(\nabla_{x_{j}}f_{j}(x_{j}),f_{1}(x_{1}),\ldots,f_{n}(x_{n})))_{j}.$

The wedge product $\wedge_{i}f_{i}$ is an equation of the hyperplane which
passes through small perturbations of the points $P_{i}$. On the other hand,
the $j$th row of the matrix $(\nabla_{x_{j}}f_{j})_{j}$ spans the tangent
space to the curve at the perturbation of $P_{j}$. Now, if we focus on the
$j$th coordinate of the function $F$, it essentially checks whether the
tangent direction lies inside the hyperplane determined by the wedge product
or deviates from it. We can make the following claim in terms of numerical
bounds: If the value of $F$ at the origin is close to $0$, and the
differential of $F$ is sufficiently far from being singular, and $F$ is well
approximated by its differential, then $F$ will attain a value close to $0$
near the origin. In order to give the needed quantitative bounds for the
Jacobian matrix of $F$, denoted ${\mathcal{J}}_{F}$, we set the following
notations:

$\displaystyle\mu_{i}:\mathbb{R}^{+}$ $\displaystyle\to\mathbb{R}^{+},\qquad$
$\displaystyle\eta_{i}:(\mathbb{R}^{+})^{n-1}$
$\displaystyle\to\mathbb{R}^{+}.$ $\displaystyle r$ $\displaystyle\mapsto
1+2|c_{i}|r+b_{i}r^{2}\qquad$
$\displaystyle(r_{1},\ldots\hat{r_{i}}\ldots,r_{n})$
$\displaystyle\mapsto\prod_{j\neq
i}(1+r_{j}+|c_{j}|r_{j}^{2}+a_{j}r_{j}^{3}),$

where the hat symbol $\hat{}$ is to remove the entry in the corresponding
position. We will customize the notation
$\eta_{i}(r_{1},\ldots\hat{r_{i}}\ldots,r_{n})$ by using $\eta_{i}(r_{i})$
instead. For a function $\psi:\mathbb{R}^{N}\to\mathbb{R}^{M}$ or
$\psi:\mathbb{C}^{N}\to\mathbb{C}^{M}$, we have the following notation for the
expansion $\psi$ around $0$:

(16)
$\mathrm{Term}_{1}(\psi):=\sum_{i}x_{i}{d\psi(x)/dx_{i}}\restriction_{x_{i}=0}={{\mathcal{J}}_{\psi}\restriction_{x_{i}=0}}x,\quad\mathrm{Tail}_{2}(\psi):=\psi-\psi(0)-\mathrm{Term}_{1}(\psi).$

######  Proposition 4.21.

Under the setting above, the diagonal $i$th entry of ${\mathcal{J}}_{F}$ is
$2c_{i}\det(\tilde{R}_{i},P_{1},\ldots,P_{n})$, and the non-diagonal entry in
$(i,j)$ is $\det(T_{i},P_{1},\ldots,P_{j-1},T_{j},P_{j+1},\ldots,P_{n}).$
Moreover,

$\displaystyle|\mathrm{Tail}_{2}(F)\restriction_{\prod\mathbb{B}_{i}}|\leq|\mathrm{Tail}_{2}\left(\Sigma_{i}\mu_{i}(r_{i})\eta_{i}(r_{i})\right)|,$

where the norm is the symmetric tensor norm.

###### Proof.

The $k$th entry of the image of $F$ is

$\det(T_{k}+2x_{k}c_{k}\tilde{R}_{k}+du_{k}/dx_{k},P_{1}+x_{1}T_{1}+c_{1}x_{1}^{2}\tilde{R}_{1}+u_{1},\ldots,P_{n}+x_{n}T_{n}+c_{n}x_{n}^{2}\tilde{R}_{n}+u_{n}).$

Expanding the expression above by their $x_{j},u_{j},du_{j}/dx_{j}$ terms , we
get $T_{k}\cdot\wedge_{i}P_{i}$ plus

$2c_{k}x_{k}\det(\tilde{R}_{k},P_{1},\ldots,P_{n})+\sum_{j\neq
k}x_{j}\det(T_{k},P_{1},\ldots P_{j-1},T_{j},P_{j+1}\ldots P_{n})$

plus all the terms in $\mathrm{Tail}_{2}$ which either include
$u_{j},du_{j}/dx_{j}$, or are nonlinear in the $x_{j}$. All the determinants
coming as coefficients to these tail terms have norms bounded by $1$, as being
determinants of matrices whose rows are of norm $1$. The bound on the
symmetric tensor norm of the tail follows since the spectral norm is bounded
by the norm of the sum of the entries, and from the bounds on
$u_{j},du_{j}/dx_{j}$. ∎

We are now prepared to establish criteria for certifying that an approximated
theta hyperplane corresponds to a true theta hyperplane. In the following
lemma, for ease of notation, we will denote by
$H^{\prime}\in{\mathbb{C}^{n+1}}^{*}$ an $L_{2}$ normalization of
$\wedge_{i}P_{i}=\wedge_{i}f_{i}(x_{i})\restriction_{x_{i}=0}$, determined up
to a norm 1 scalar. We write $H^{\prime}$ and not $H$; signifying the fact
that it is merely a putative approximation of a theta hyperplane. In general,
we will put a _prime symbol_ as a superscript for the computed values.

######  Lemma 4.22.

With notation introduced in Equation 16, if

(17)
$\left|\langle\wedge_{i}f_{i}(x_{i})\restriction_{0},H^{\prime}\rangle\cdot{\mathcal{J}}_{F}\restriction_{0}^{-1}\left((\nabla_{x_{j}}{f_{j}}\restriction_{0})\wedge
H^{\prime}\right)_{j}\right|<\min_{i}r_{i}-\mathrm{Tail}_{2}(\Sigma_{i}{\mu_{i}(r_{i})\eta_{i}(r_{i})})/\sigma_{n}({\mathcal{J}}_{F}\restriction_{0})$

then there is a solution for $F=0$ on $\prod_{i}\mathbb{B}_{i}$.

###### Proof.

Expanding $F$ around $x=0$ we get
$F=F(0)+({\mathcal{J}}_{F}\restriction_{0})\cdot x+\mathrm{Tail}_{2}(F)$. That
is, the image of $F$ contains an open ball of radius
$\sigma_{n}({\mathcal{J}}_{F}\restriction_{0})\min_{i}r_{i}-\sup_{\prod\mathbb{B}_{i}}\mathrm{Tail}_{2}F$
around $F(0)$, and the linear approximation and the tail are as in the
formula. The result follows. ∎

### Dealing with inaccuracies of the input points

We face a problem which is very similar to the problem of locally
approximating the curve with a conic; namely, we approximate points on the
curve with approximate computations. The error-bounding computations are
similar in spirit, but technically different: on the one hand they are harder
as one has to account not only for error accumulation but also for errors
introduced by the computations themselves, on other they are easier as one
does not have to deal with bounding polynomials expressions in $|x|$ as in
e.g. Corollary 4.20.

######  Definition 4.23 (Notations for computing the numerical stabilitiy).

Recall that $H^{\prime}$ denotes an approximated theta hyperplane. We
normalize $H^{\prime}$ so that $|H^{\prime}|=1$. Let $P,P^{\prime}$ be
representation in $\mathbb{C}^{g}$ points $\mathbb{P}^{g-1}$ such that

$\mathbb{P}P\in C\cap
Z(H^{\prime}),\qquad|P|=|P^{\prime}|=1,\qquad|P-P^{\prime}|<\epsilon_{0},\qquad\sum_{i}|{P^{\prime}}^{t}Q_{i}P^{\prime}|^{2}<\epsilon_{0}^{2}.$

Let $T,R,M_{T},M_{R}$ be the vectors in $\mathbb{C}^{g}$ and complex matrices,
that are defined in Lemma 4.4. We denote by
$T^{\prime},R^{\prime},M_{T^{\prime}},M_{R^{\prime}}$ the machine computed
versions of them represented with floating-point numbers. More precisely,
$T^{\prime}$ is the computed lowest singular vector of $M_{T}^{\prime}$ of
norm 1 whereas $R^{\prime}$ is computed via Equation 14, but reformulated in
terms of the computed $P^{\prime}$ and $T^{\prime}$. The requirement on
$|P-P^{\prime}|$ follows from the requirement on
$\sum_{i}|{P^{\prime}}^{t}Q_{i}P^{\prime}|^{2}$ for an appropriate choice of
$\epsilon_{0}$ by Lemma 4.6. We refer again to Figure 1 for the relation
between $P^{\prime}$ and $P$.

In the following propositions, we present bounds between values and their
machine computed versions. In Section 7, the explicit bounds that arise from
our key examples will appear. For our computations, we use the linear algebra
library `LAPACK`. Let us start with a remark that gives details about the
methods in `LAPACK` about error computing.

######  Remark 4.24 (Using LAPACK as a certification accuracy oracle).

Denote `E` the machine numerical epsilon. We trust `LAPACK` accuracy bounds:
`SERRBD`$(A)$, `VERRBD`$(A,{\verb|i|})$, `UERRBD`$(A,{\verb|i|})$, which
denote the errors in singular values of the matrix $A$, and in the angles of
singular vectors for the matrix $A$ respectively. The values and the vectors
are computed by the function `zgesvd`. The values of `VERRBD` and `UERRBD` are
computed by the function `ddisna`, whereas `SERRBD` is the top computed
singular value times `E`. We note that if several consecutive singular values
are close then the values given by `ddisna` are meaningless. Finally, we also
trust the least square solution function `zgelss` in the sense that the error
is the reported error in the over-defined case. That being mentioned, the
reported error is the square root of a sum of numbers, so we have to account
for the error in computing this sum as well.

We continue with bounding the singular values of machine-computed matrices.

######  Proposition 4.25.

Let $A$ be an $n\times m$ matrix. If we consider the computed version
$A^{\prime}$ of $A$ we have the following bound

(18)
$|\sigma^{\prime}_{i}(A^{\prime})-\sigma_{i}(A)|\leq\sigma^{\prime}_{1}(A^{\prime})\cdot\,{\verb|E|}+\sqrt{\sum_{j}\epsilon_{\mathrm{row}_{j}}^{2}}\cdot(1+(2n+1){\verb|E|}),$

where $\epsilon_{\mathrm{row}_{j}}$ is the bound between the $L_{2}$-norm of
the difference between the $j$th row of $A^{\prime}$ and the $j$th row of $A$.

###### Proof.

We first write the following difference as follows

$|\sigma^{\prime}_{i}(A^{\prime})-\sigma_{i}(A)|\leq|\sigma^{\prime}_{i}(A^{\prime})-\sigma_{i}(A^{\prime})|+|\sigma_{i}(A^{\prime})-\sigma_{i}(A)|.$

Following the `LAPACK` documentation, the first summand is bounded by
$\sigma^{\prime}_{1}(A^{\prime})\,\cdot$ `E`. The second summand is bounded by
the second summand in Equation 18, which follows from the definition of the
Frobenius form and Theorem 4.1. ∎

######  Definition 4.26.

We denote the bound on the singular values of $A$ above and below by
$\sigma_{i}(A)^{+}$ and $\sigma_{i}(A)^{-}$ respectively.

######  Remark 4.27 (Ignoring the error in estimating the error).

We will compute our error estimates using a computer. The acute reader will
notice that this computation typically involves adding and sometimes
multiplying positive real numbers. Thus this computation also introduces
errors. However, since we are talking here about positive real numbers the
error introduced in the computation will have a size of the error time
$1+O({\verb|E|})$. This extra error is meaningless and will be ignored in the
presentation.

We now bound the difference between the singular values of $M_{T}$ and the
computed matrix $M_{T^{\prime}}$.

######  Proposition 4.28.

We have the following bound for the difference between the computed and the
actual singular values:

(19)
$|\sigma^{\prime}_{i}(M_{T^{\prime}})-\sigma_{i}(M_{T})|<\sigma^{\prime}_{1}(M_{T^{\prime}})\cdot{\verb|E|}+{\epsilon_{0}\sqrt{(1+\sum_{j}(|Q_{j}|+g{\verb|E|})}},$

where $\epsilon_{0}$ is the bound for $|P^{\prime}-P|$ as introduced in
Definition 4.23.

###### Proof.

Per Proposition 4.25, we have

$|\sigma^{\prime}_{i}(M_{T^{\prime}})-\sigma_{i}(M_{T})|<\sigma^{\prime}_{1}(M_{T^{\prime}})\cdot{\verb|E|}+|\sigma_{i}(M_{T^{\prime}})-\sigma_{i}(M_{T})|.$

It follows from Theorem 4.1, the norm
$|\sigma_{i}(M_{T^{\prime}})-\sigma_{i}(M_{T})|$ is bounded by the norm
$|M_{T^{\prime}}-M_{T}|$. We compute the error bound for the latter by
bounding the Frobenius norm $||M_{T^{\prime}}-M_{T}||_{F}$ as follows. The
difference between the last line of $M_{T^{\prime}}$ and the last row of
$M_{T}$ is simply $({P^{\prime}}-P)^{\dagger}$, and the difference between all
the other rows is given by $Q_{i}(P^{\prime}-P)$. The quantity $g$`E` accounts
for the $g$ operations involved in computing each $Q_{i}P^{\prime}$ as the
entries of each $Q_{i}$ are small integers and the norm of the entries of
$P^{\prime}$ are $<1$. ∎

We now estimate $|T^{\prime}-T|$ and obtain a bound for it.

######  Proposition 4.29.

Let $T^{\prime}$ and $T$ be the vectors defined in Definition 4.23. Then

(20)
$|T^{\prime}-T|<\epsilon_{1}:=1+|T^{\prime}|-2\cos(\verb|E|\cdot\sigma_{1}^{\prime}(M_{T^{\prime}})/{\verb|ddisna|}(M_{T^{\prime}},g))+\sigma_{g}(M_{T})^{+}/\sigma_{g-1}(M_{T})^{-}.$

###### Proof.

For the first three summands in Equation 20, per `LAPACK`’s documentation of
the subroutine `ddsina` as discussed in Remark 4.24, the angle between the
true $i$th singular vector and the computed singular vector of a matrix $A$ is
given by $\verb|E|\cdot\sigma_{1}^{\prime}(A)/{\verb|ddisna|}(A,i)$. The last
summand follows from Proposition 4.2. ∎

Set
$\epsilon_{M_{T^{\prime}}}:=\epsilon_{0}\sqrt{(1+\sum_{j}(|Q_{j}|+g{\verb|E|})}$,
which is the quantity bounding $||M_{T^{\prime}}-M_{T}||_{F}$ in Proposition
4.28.

######  Proposition 4.30.

We have the following bounds considering the vectors $R$ and $R^{\prime}$:

  1. (1)

$||M_{R^{\prime}}-M_{R}||_{F}<(1+2{\verb|E|})\sqrt{\epsilon_{M_{T^{\prime}}}^{2}+|T-T^{\prime}|^{2}}$,
which we denote by $\epsilon_{M_{R^{\prime}}}$,

  2. (2)

$|\sigma^{\prime}_{i}(M_{R^{\prime}})-\sigma_{i}(M_{R})|<\sigma^{\prime}_{1}(M_{R^{\prime}})\cdot\verb|E|+\epsilon_{M_{R^{\prime}}}.$

###### Proof.

Item 1 follows from the definition of the matrix $M_{R}$, which has been
formed from the matrix $M$ with the additional two equations as introduced in
Definition 4.5. For Item 2 we proceed as in the proof of Proposition 4.28. We
have

(21)
$|\sigma^{\prime}_{i}(M_{R^{\prime}})-\sigma_{i}(M_{R})|<\sigma^{\prime}_{1}(M_{R^{\prime}})+||M_{R^{\prime}}-M_{R}||_{F}.$

This concludes a proof after combining Equation 21 with Item 1. ∎

Note that the computed version of Equation 14 can be written in terms of the
computed vector $R^{\prime}$s as follows:

(22) $\quad
M_{R^{\prime}}R^{\prime}=-\frac{1}{2}\left(\begin{smallmatrix}{T^{\prime}}^{t}Q_{1}T^{\prime}\\\
\vdots\\\ {T^{\prime}}^{t}Q_{\dim I_{2}(C)}T^{\prime}\\\ 0\\\
0\end{smallmatrix}\right)$

Let $\delta_{R}$ be the bound on the solution accuracy of this system, given
by

$\delta_{R}:=(\text{sum of last $\dim I_{2}(C)+2-g$ ``}{\verb+zgelss+}\text{
solution" entries})\cdot(1+(\dim I_{2}(C)+2-g)\verb|E|),$

per Remark 4.24. We denote by $\xi_{T}$ the upper bound on the top singular
value of the matrix whose $j$th row is given by ${T^{\prime}}^{t}Q_{j}$. Then
we have the following bound

$|({T^{\prime}}^{t}Q_{i}T^{\prime})_{i}-(T^{t}Q_{i}T)_{i}|<\gamma_{R}:=\frac{1}{2}(2\xi_{T}\epsilon_{1}+\epsilon_{1}^{2}q_{m}),$

where $q_{m}$ is the norm defined in Definition 4.7.

We can now make the following error estimation:

######  Proposition 4.31.

With the notation $\sigma_{g}(M_{R})^{-}$ from Definition 4.26, the difference
$|R-R^{\prime}|$ is bounded by

(23)
$\epsilon_{2}:=\frac{1}{2}|({T^{\prime}}^{t}Q_{i}T^{\prime})_{i}|\cdot\left(\frac{\epsilon_{M_{R}}}{\sigma_{g}{M_{R}}^{-}(\sigma_{g}(M_{R})^{-}-\epsilon_{M_{R}})}+(\gamma_{R}+\delta_{R})/(\sigma_{g}(M_{R})^{-}-\epsilon_{M_{R}})\right).$

###### Proof.

This follows from Proposition 4.3. ∎

In what follows, we work out the version of Lemma 4.22 which accounts for the
numerical errors $|P_{ij}-P^{\prime}_{ij}|$, and the resulting errors in the
$T_{ij},R_{ij}$. We do not elaborate here on the accuracy of
$r_{P,x},r_{T,x},r_{R,x}$ or the coefficients of the monomials in $|x|$
appearing in $\epsilon_{P,x}$ and $\epsilon_{R,x}$, which are both used in
bounding the $u_{i}$. For our purposes, on which we will expand in Section 7,
it will suffice to show that all these numbers are $O(1)$, whereas the
numerical errors in $P_{ij},T_{ij},R_{ij}$ are orders of magnitude smaller. So
there is no reason to delve into precise computations to show that one is
greater than the other. This is not the case for Lemma 4.22 for the simple
reason that the expressions appearing there, e.g. the top singular value of
${\mathcal{J}}_{F}$, are complicated, and so the accumulated error analysis is
complicated.

Suppose we are given $g-1$ many points $P^{\prime}_{i}$ that are estimates of
the true tangency points. Let $P_{i},P^{\prime}_{i}$ be representations in
$\mathbb{C}^{g}$ of points in $\mathbb{P}^{g-1}$ such that

$\mathbb{P}P_{i}\in C\cap
Z(H^{\prime}),\qquad|P_{i}|=|P^{\prime}_{i}|=1,\qquad|P_{i}-P^{\prime}_{i}|<\epsilon_{0,i},\qquad\sum_{j}|{P^{\prime}_{i}}^{t}Q_{j}P^{\prime}_{i}|^{2}<\epsilon_{0,i}^{2}.$

Let $T_{i},T^{\prime}_{i},\tilde{R}_{i},\tilde{R}^{\prime}_{i}$ etc. be all
the relevant objects mentioned in Definition 4.23 with the extra subscript $i$
which refers to the $i$th putative tangency points $P_{i}$ and
$P^{\prime}_{i}$. So e.g. $\epsilon_{1,i}$ is the bound for
$|T_{i}^{\prime}-T_{i}|$ considering the corresponding $T_{i}$ and
$T_{i}^{\prime}$ and similarly, $\epsilon_{2,i}$ concerns the bound on
$|R_{i}-R^{\prime}_{i}|$.

We now rewrite the function defined in Equation 15 in its computed version,
using the prime symbol for denoting a computed version as usual.

(24) $\displaystyle F^{\prime}:\prod_{i}\mathbb{B}_{i}$
$\displaystyle\to\mathbb{C}^{n}$ $\displaystyle(x_{1},\ldots,x_{n})$
$\displaystyle\mapsto\left(\nabla_{x_{j}}f^{\prime}_{j}(x_{j})\wedge\bigwedge_{i=1}^{n}f^{\prime}_{i}(x_{i})\right)_{j}=(\det(\nabla_{x_{j}}f^{\prime}_{j}(x_{j}),f^{\prime}_{1}(x_{1}),\ldots,f^{\prime}_{n}(x_{n})))_{j},$

where $f^{\prime}_{i}:\mathbb{B}_{i}\to\mathbb{C}^{n}$ is given by
$f^{\prime}_{i}(x_{i})=P^{\prime}_{i}+x_{i}T^{\prime}_{i}+c^{\prime}x^{2}\tilde{R_{i}^{\prime}}$.

######  Lemma 4.32.

The Frobenius norm of the difference
${\mathcal{J}_{F}}\restriction_{0}-{\mathcal{J}_{F^{\prime}}}\restriction_{0}$
is bounded by

$\displaystyle\epsilon_{{\mathcal{J}}_{F}}:=$
$\displaystyle(1+\verb|E|)^{g^{2}+3g}\Big{(}2\sum_{i}\epsilon_{2,i}(1+\epsilon_{2,i})\prod_{k\neq
i}(1+\epsilon_{0,k})$ $\displaystyle+$ $\displaystyle
2\sum_{i}(|R^{\prime}_{i}|+\epsilon_{2,i})((1+\epsilon_{2,i})\prod_{k\neq
i}(1+\epsilon_{0,k})-1)+\sum_{i,j}\big{(}(1+\epsilon_{1,i})(1+\epsilon_{1,j})\prod_{k\neq
i,j}(1+\epsilon_{0,k})-1\big{)}\Big{)}.$

###### Proof.

The norm estimation is based on the computations in the proof of Proposition
4.21. We denote

$\delta_{P_{i}}:=P^{\prime}_{i}-P_{i},\delta_{T_{i}}:={T^{\prime}}_{i}-T_{i},\delta_{R_{i}}:={R^{\prime}}_{i}-R_{i},\delta_{c_{i}}:=c^{\prime}_{i}-c_{i}.$

We also set $X_{I}$ to be the matrix whose first row is $\delta_{R_{i}}$ if
$i\in I$, and $R_{i}$ otherwise, and whose $k>1$ row is $\delta_{P_{k}}$ if
$k\in I$ and $P_{i}$ otherwise. The differences between the $g$ diagonal
elements in Proposition 4.21 can be written as

$\displaystyle
2(c_{i}+\delta_{c_{i}})\det(R_{i}+\delta_{R_{i}},P_{1}+\delta_{P_{1}},\ldots,P_{n}+\delta_{P_{n}})-2c_{i}\det(R_{i},P_{1},\ldots,P_{n})$
$\displaystyle=2\delta_{c_{i}}\det(R_{i}+\delta_{R_{i}},P_{1}+\delta_{P_{1}},\ldots,P_{n}+\delta_{P_{n}})+2c_{i}\sum_{\begin{subarray}{c}I\subset\\{1,\dots,n+1\\}\\\
\\#I>1\end{subarray}}\det(X_{I}).$

The norms of these summands are bounded by the following quantities
respectively:

$2|\delta_{c_{i}}|(1+|\delta_{c_{i}}|)\prod_{i\neq k}(1+\epsilon_{0,k}),\qquad
2|c_{i}|((1+|\delta_{c_{i}}|)\prod_{i\neq k}(1+\epsilon_{0,k})-1).$

Likewise, we now compute the differences between the non-diagonal elements in
Proposition 4.21, of which there are $g(g-1)$, are given by

$\det(T^{\prime}_{i},P^{\prime}_{1},\ldots,P^{\prime}_{j-1},T^{\prime}_{j},P^{\prime}_{j+1},\ldots,P^{\prime}_{n})-\det(T_{i},P_{1},\ldots,P_{j-1},T_{j},P_{j+1},\ldots,P_{n}).$

A similar argument for the diagonal elements gives us the bound

$(1+\epsilon_{1,i})(1+\epsilon_{1,j})\prod_{k\neq i,j}(1+\epsilon_{0,k})-1.$

To conclude the proof we observe that $|\delta_{c,i}|<\epsilon_{2,i}$, and
that $|c_{i}|=|R_{i}|$, then we count the number of arithmetic operations
involved in the computation. ∎

######  Proposition 4.33.

Suppose that $x$ and $x^{\prime}$ are the solutions of $F(x)=0$ and
$F^{\prime}(x^{\prime})=0$ of the smallest norms. For any $\epsilon_{x}$ which
satisfies

(25) $\displaystyle\epsilon_{x}<$
$\displaystyle\frac{\mathrm{Tail}_{2}(\sum_{i}\mu_{i}((|x^{\prime}_{i}|+\epsilon_{x})_{i})\eta_{i}(|x^{\prime}_{i}|+\epsilon_{x}))}{\sigma_{g-1}({\mathcal{J}}_{F}^{\prime}\restriction_{0})-\epsilon_{{\mathcal{J}}_{F}}}+\frac{\sum_{i}3|c^{\prime}_{i}|\cdot|x^{\prime}_{i}|^{2}+2|c^{\prime}_{i}|^{2}|x_{i}^{\prime}|^{3}}{\sigma_{g-1}({\mathcal{J}}_{F}^{\prime}\restriction_{0})}$
$\displaystyle+$
$\displaystyle\frac{\epsilon_{{\mathcal{J}}_{F}}}{\sigma_{g-1}({\mathcal{J}}_{F^{\prime}}\restriction_{0})(\sigma_{g-1}({\mathcal{J}}_{F^{\prime}}\restriction_{0})-\epsilon_{{\mathcal{J}}_{F}})}+\frac{|((T_{j}\wedge\bigwedge_{i}P_{i})-(T^{\prime}_{j}\wedge\bigwedge_{i}P^{\prime}_{i}))_{j}|}{\sigma_{g-1}({\mathcal{J}}_{F^{\prime}}\restriction_{0})-\epsilon_{{\mathcal{J}}_{F}}}$

and $|x^{\prime}_{i}|+\epsilon_{x}<r_{P_{i},x},r_{T_{i},x},r_{R_{i},x}$, the
magnitude $|x^{\prime}-x|$ is bounded by $\epsilon_{x}$.

###### Proof.

The first two terms are tail components which come Lemma 4.22 for $F$ and
$F^{\prime}$ respectively. The last two terms come from Proposition 4.3 when
comparing the solutions of $F(0)=-{\mathcal{J}}_{F}\restriction_{0}\cdot x$,
and $F^{\prime}(0)=-{\mathcal{J}}_{F^{\prime}}\restriction_{0}\cdot
x^{\prime}$. ∎

######  Remark 4.34.

We have several computational oriented remarks regarding Proposition 4.33:

  1. (1)

Realistically (i.e., this is what happens in practice), the biggest term by
far in Equation 25 is the third. The numerator in that term is far bigger than
the other numerators, and the denominator is far smaller.

  2. (2)

As in the case of $r$ in Lemma 4.6 and its utilization in Corollary 4.10, the
quantity $\epsilon_{x}$ appears in both sides of the inequality in Equation
25. On the RHS it appears as a higher degree term. For the same reason as in
Lemma 4.6, the higher degree term is inconsequential.

  3. (3)

One computes $x^{\prime}$ by solving the linear equation for
${\mathcal{J}}_{F^{\prime}}$. This raises two issues in theory. However, in
practice, we will see in Section 7 that in our cases, the quantities appearing
in the inequalities here differ by several orders of magnitude.

    1. (a)

There is an error in computing ${\mathcal{J}}_{F^{\prime}}$ given
$P^{\prime}_{i},T^{\prime}_{i},R^{\prime}_{i}$. This error is identical to the
error from $\epsilon_{{\mathcal{J}}_{F}}$ introduced in Lemma 4.32, when one
substitutes all the $\epsilon_{j,i}$ for `E`.

    2. (b)

One has to control the inequality analog to Equation 17, with
${\mathcal{J}}_{F^{\prime}},x^{\prime}$ instead of ${\mathcal{J}}_{F},x$, to
ensure its validity.

  4. (4)

Like any gradient descent method, the method at hand can serve as a basis for
refining an existing solution, not merely bounding its accuracy. We believe
this is a worthy effort, but it involves a lot of parameter tuning and
interaction with the gradient descent of Lemma 4.6. Thus, it lies beyond the
scope of the current work. That being said, we did run one iteration with
hand-tailored parameters for the descent. See Section 7.6 Item 3.

######  Theorem 4.35 (Bounding $\sqrt{\sum_{i}|P^{\prime}_{i}-(\text{the
closest point to }P^{\prime}_{i}\text{ in }H\cap C)|^{2}}$).

With $x$ and $x^{\prime}$ as in Proposition 4.33, we have

$\displaystyle\sum_{i}$ $\displaystyle|P^{\prime}_{i}-\text{closest point to
}P^{\prime}_{i}\text{ in }H\cap
C|^{2}<\sum_{i}|f_{i}(x_{i})-f^{\prime}_{i}(x^{\prime}_{i})|^{2}+\sum_{i}|P^{\prime}_{i}-f^{\prime}_{i}(x^{\prime}_{i})|^{2}$
$\displaystyle<$
$\displaystyle\sum_{i}|(P_{i}-P^{\prime}_{i})+(x_{i}-x^{\prime}_{i})T_{i}+x^{\prime}(T_{i}-T^{\prime}_{i})+(x_{i}^{2}-{x^{\prime}_{i}}^{2})R_{i}+{x^{\prime}_{i}}^{2}(R_{i}-R^{\prime}_{i})+u_{i}(x)|^{2}+\sum_{i}|P^{\prime}_{i}-f^{\prime}_{i}(x^{\prime}_{i})|^{2}$
$\displaystyle<$
$\displaystyle\left(\epsilon_{0,i}+\epsilon_{x}+|x^{\prime}|\epsilon_{1,i}+(2\epsilon_{x}(|x^{\prime}|+\epsilon_{x})|x^{\prime}|+\epsilon_{x}^{2})(|R^{\prime}_{i}|+\epsilon_{2,i})+|x^{\prime}|^{2}\epsilon_{2,i}+a_{i}(|x_{i}^{\prime}|+\epsilon_{x})^{3}\right)^{2}+\left(|x^{\prime}|+|R^{\prime}_{i}|\cdot|x^{\prime}_{i}|^{2}\right)^{2}.$

###### Proof.

The first inequality is immediate from the definition, and then we simply
expand. ∎

When computing the quantity above, to bound the estimation error, one has to
multiply it by $(1+$`E`$)^{N}$ where $N$ is the number of machine operations
involved in the computation.

Intuitively, in Lemma 4.22, we perturb points along $C$, which means that if
the $H^{\prime}\wedge T_{P_{i}}C$s are small, then the magnitude
$|H-H^{\prime}|$ should be much smaller than the distance bound from Theorem
4.35. Below we make this intuitive view precise.

######  Theorem 4.36 (Bounding $|H-H^{\prime}|$).

Let $n_{H^{\prime}}$ be the unit normal to $H^{\prime}$. Then under the
assumption of Proposition 4.33, the magnitude $|H-H^{\prime}|$ is bounded by

$(|x^{\prime}|+\epsilon_{x})\Big{(}\sum_{i}|\langle
T^{\prime}_{i},n_{H^{\prime}}\rangle|+\epsilon_{1,i}\Big{)}+\mathrm{Tail}_{2}\Big{(}\prod_{i=1}^{n}(1+(|x_{i}^{\prime}|+\epsilon_{x})+|c_{i}^{\prime}|(|x_{i}^{\prime}|+\epsilon_{x})^{2}+a_{i}^{\prime}(|x_{i}^{\prime}|+\epsilon_{x})^{3})\Big{)}.$

###### Proof.

Since $H^{\prime}=\wedge_{i=1}^{n}P_{i}$, the points $P_{1},\ldots,P_{n+1}$
together $n_{H^{\prime}}$ form a basis to $\mathbb{C}^{n+1}$. Moreover, since
$P_{i}\perp T_{i}$, the coefficient of $P_{i}$ in the representation of
$T_{i}$ in this basis is trivial, that is to say $|T_{i}\wedge\bigwedge_{k\neq
i}P_{i}|<|\langle T_{i},n_{H^{\prime}}\rangle|$. Hence:

$\displaystyle|H-H^{\prime}|=\wedge_{i}(P_{i}+x_{i}T_{+}x_{i}^{2}c_{i}R_{i}+u_{i}(x_{i}))-\wedge_{i}P_{i}|$
$\displaystyle<\sum_{i}|x_{i}||\langle
T_{i},n_{H^{\prime}}\rangle|+\mathrm{Tail}_{2}\Big{(}\prod_{i=1}^{n}(1+|x_{i}|+|c_{i}|\cdot|x_{i}|+a_{i}|x_{i}|)\Big{)}$
$\displaystyle<(|x^{\prime}|+\epsilon_{x})\Big{(}\sum_{i}|\langle
T^{\prime}_{i},n_{H^{\prime}}\rangle|+\epsilon_{1,i}\Big{)}+\mathrm{Tail}_{2}\Big{(}\prod_{i=1}^{n}(1+(|x_{i}^{\prime}|+\epsilon_{x})+|c_{i}^{\prime}|(|x_{i}^{\prime}|+\epsilon_{x})^{2}+a_{i}^{\prime}(|x_{i}^{\prime}|+\epsilon_{x})^{3})\Big{)}.$

∎

We close this section with the following important remark.

######  Proposition 4.37.

Assuming that $\dim\ker M_{T_{i}^{\prime}}<2$ for all the points
$P^{\prime}_{i}$ of a certified hyperplane, the corresponding theta
characteristic $\theta$ is odd i.e., $\dim H^{0}(C,\theta)=1$, where
$M_{T^{\prime}}$ is the computed version of the matrix $M_{T}$ defined in
Definition 4.5.

##  5\. Numerically certifying Steiner sets

In our study, we are provided with putative approximate theta hyperplanes,
which are numerically computed using theta functions, and thus come with
putative level 2 structure as being labelled with half-characteristics, i.e.,
picking an even two torsion point in the Jacobian we identify
$\mathrm{Pic}^{g-1}$ with $\mathrm{Pic}^{0}$. Thus, considering these theta
hyperplanes as points in an odd affine $SP_{2}(2g)$ homogeneous space, we give
them linear homogeneous space coordinates by picking a distinguished $0$ among
the points in the even affine space. The syzygetic quadruples, and so the
Steiner sets, can be built just by means of these coordinates. These are the
objects we require for our main theorem. In the previous section, we have
shown how to certify an approximate theta hyperplane representing an odd theta
characteristic. However, this does not certify the coordinates that encode the
homogeneous level 2 structure. In this section, we aim to certify the level 2
structure data, however, we merely certify the part we need for our purposes,
namely the subsets of a small number of Steiner sets which suffice for the
dimension computations.

We recall from Equation 1 that the Steiner set of a non trivial 2-torsion
point $\alpha\in\mathrm{Jac}(C)[2]$ is

$\Sigma_{\alpha}=\\{\\{\theta,\theta+\alpha\\}:{\dim}H^{0}(C,\theta)={\dim}H^{0}(C,\theta+\alpha)=1\\}.$

In other words, a pair of such pairs of distinct and odd theta characteristics
$\\{\\{\theta_{1},\theta_{2}\\},\\{\theta_{3},\theta_{4}\\}\\}$ are in the
same (resp. not in the same) Steiner set if and only if they are syzygetic
(resp. azygetic); which is if
$\theta_{1}+\theta_{2}+\theta_{3}+\theta_{4}=2K_{C}$ (resp. $\neq 2K_{C}$).
Note that in the azygetic case, the sum of the four theta characteristics is
$2K_{C}$ plus some nonzero 2-torsion point.

We will certify Steiner sets by certifying putative syzygetic and azygetic
quadruples. There are two ways to tackle this problem. The first is a combined
algebraic and analytic approach. Given a quadruple of theta characteristics
that are suspected of being syzygetic, considered as $4(g-1)$ points in the
canonical system $|K_{C}|^{*}$, we seek to find a quadric in
$\mathrm{Sym}^{2}H^{0}(K_{C})$ which is close to zero on all of them. One may
then use this computation to bound the distance of these points from the image
of $2K_{C}\in\mathrm{Pic}^{4(g-1)}$, and show that it is smaller than the
distance of any nonzero 2-torsion point plus $2K_{C}$. For a general curve,
this would require us to “cross the transcendental divide”. However, with the
specific curves we use for our certification, whose Jacobians are isogenous to
a product of elliptic curves, one may project the divisor at hand to a product
of elliptic curves after identifying the full automorphism group of the curve
and analyze the kernel of the isogeny. The alternative approach, which is the
one we use in this work is an algebraic combinatorial approach. It requires
more computation, but it works for general curves, i.e., not only just ones
with many automorphisms. In this approach, the analytical side is slightly
simpler. In fact, we do not aim to directly prove the syzygecy inside a
putative Steiner set but rather prove the azygecy between elements which are
not in the same Steiner set. Thus instead of showing a quantity that is $0$,
we have to show that a sufficient number of quantities are not $0$.

######  Proposition 5.1.

Let $S,\overline{S}$ be two subsets of the set of unordered pairs of odd theta
characteristics of $C$ such that $\\#\overline{S}=2^{2g}-2$ and any pair of
elements in either $\binom{\overline{S}}{2}$ or $\overline{S}\times S$ forms
an azygetic quadruple then $S$ is a constituent of a Steiner set.

###### Proof.

Each element of $\overline{S}$ is in a different Steiner set, and each element
of $S$ is in a Steiner set where none of the elements of $\overline{S}$ are,
but there is only one Steiner set left. ∎

To reduce the number of putative azygetic quadruples we check during the
certification process, we also employ approximations of the quadrics derived
from the non-certified syzygetic quadruples of theta characteristics through
numerical computations. Specifically, we need a numerical version of the
following proposition.

######  Proposition 5.2.

Let $D_{1},D_{2},D_{3},D_{4}$ be effective divisors of degree $2g-2$ such that
$D_{1}+D_{2}\sim D_{3}+D_{4}\sim D_{1}+D_{3}\sim 2K_{C}$ then $D_{2}+D_{4}\sim
2K_{C}$.

###### Proof.

It follows from $D_{2}+D_{4}\sim D_{1}+D_{2}+D_{3}+D_{4}-(D_{1}+D_{3})$. ∎

To work with Proposition 5.2 in numerical terms for the certification
analysis, we reformulate it as follows. Before that we recall that $I_{n}(C)$
is the ideal of the homogeneous forms of degree $n$ containing
$C\subset\mathbb{P}^{g-1}$.

######  Proposition 5.3.

Given the divisors $D_{1},D_{2},D_{3},D_{4}$, we consider the set of quadrics

$\overline{Q}_{ij}:=\\{Q\in\mathrm{Sym}^{2}H^{0}(K)\,|\,Q\restriction_{{D_{i}\cup
D_{j}}}=0\\}$

for $i,j=1,\dots,4$ with $i\neq j$. We assume that $\dim\overline{Q}_{ij}>\dim
I_{2}(C)$ for $(i,j)=(1,2),(3,4),(1,3)$. If we pick representatives
$Q_{ij}\in\overline{Q}_{ij}\smallsetminus I_{2}(C)$ for $(i,j)=(1,2),(3,4)$
then there is a quadric $Q_{13}\in\overline{Q}_{13}$ and a quadric
$Q_{24}\in\mathbb{P}^{g-1}$ so that $Q_{12}Q_{34}-Q_{13}Q_{24}\in I_{4}(C)$.

###### Proof.

We consider the 1-dimensional subspace

$L:=\mathrm{span}(Q_{12},I_{2}(C))/I_{2}(C)\otimes\mathrm{span}(Q_{34},I_{2}(C))/I_{2}(C)\subset(\mathrm{Sym}^{2}H^{0}(K)/I_{2}(C))^{\otimes
2}\cong H^{0}(2K)^{\otimes 2}.$

Since $Z(Q_{12}Q_{34})\cap C\geq D_{1}+D_{3}$, the projection of $L$ to a
point in $|4K|$ sits on the image of

$\mathrm{Sym}^{2}H^{0}(K)/I_{2}(C)\otimes\overline{Q}_{13}/I_{2}(C)$

in $|4K|$. Considering $L$ as a subspace of the preimage of this point in
$H^{0}(2K)^{\otimes 2}$, we may project it to $H^{0}(2K)^{\otimes
2}/(H^{0}(2K)\otimes\overline{Q}_{13}/I_{2}(C))$. A preimage of the projection
is a desired $Q_{24}$. The difference $Q_{12}Q_{34}-Q_{13}Q_{24}$ is in
$I_{4}(C)$ as the two products of quadrics are represented by $L$. ∎

Unlike the proof of Proposition 5.2, this one has a tractable numerical
analogue, achieved by means of computing the SVD of the matrices representing
the vector spaces $\overline{Q}_{ij}$. We now introduce notation to realize
it.

######  Definition 5.4.

Fixing a positive integer $k$, let $D$ be an effective divisor of degree
$N>0$, which is supported on distinct points. Assume that the points in the
support of $D\in(|K_{C}|^{*})^{N}$ are ordered, and their representatives in
$H^{0}(K_{C})^{*}$ are of norm $1$. Let $M_{D}$ be the

$\deg D\times\dim\mathrm{Sym}^{k}H^{0}(K_{C})=\deg D\times\binom{k+g-1}{k}$

complex matrix whose entry at the row $r$, the column $c$ is the evaluation of
the $r$th point of $D$ at the $c$th degree $k$ monomial on $H^{0}(K_{C})$ in
lexicographic order.

Note that the lower $\dim I_{k}(C)$ singular values of $M_{D}$ are trivial. In
the most interesting case for us, when $k=2$, the divisors $D_{1},D_{2}$ are
of degree $2g-2$. We abuse the notation and write $M_{ij}:=M_{D_{i}+D_{j}}$.
On the other hand, the complex unitary norm on $\mathrm{Sym}^{k}H^{0}(K_{C})$
under the monomial basis induces an isomorphism $I_{k}(C)^{\perp}\cong
H^{0}(kK_{C})$. Moreover, this norm induces the Fubini-Study metric on
$|kK_{C}|$. In what follows, we establish a distance-like notion between the
elements of $\mathrm{Sym}^{k}H^{0}(K_{C})$ and the divisors, which measures
the proximity. The objects of interest will be considered normalized to norm 1
with respect to the relevant basis.

######  Definition 5.5.

Given a polynomial $p\in\mathrm{Sym}^{k}H^{0}(K_{C})$ of norm $1$, we define
the _gap_ between the divisor $D$ and the polynomial $p$ as:

(26) ${\mathfrak{g}}(D,p):=|M_{D}\cdot
V_{p}|=\sqrt{\sum_{\begin{subarray}{c}{x\in D}\\\
{|x|=1}\end{subarray}}|p(x)|^{2}},$

where $V_{p}$ is the vector whose entries are the coefficients of $p$.

######  Remark 5.6.

The gap between the divisor $D$ in ${|K|^{*}}^{2k(g-1)}$ and $p$ is bounded
from below by the maximum of the Fubini Study distances between the points of
$Z(p)\cap C$ and $D$, and from above by the sum of these distances.

###### Proof.

This follows from the definition of the gap ${\mathfrak{g}}$. ∎

######  Definition 5.7.

We say that $D_{i},D_{j}$ are $\epsilon$-approximately syzygetic (resp.
$A$-numerically azygetic) if the $(\dim I_{2}(C)+1)$-th lower singular value
of $M_{ij}$ is $<\epsilon$ (resp. $>A$). Let $D_{1},D_{2},D_{3},D_{4}$ be
effective divisors of degree $2g-2$, supported on distinct points. Assume that
the pairs $(D_{1},D_{2})$, $(D_{3},D_{4})$ are $\epsilon$-approximately
syzygetic, and the pair $(D_{1},D_{3})$ is syzygetic.

Recalling the notation introduced in Proposition 5.3, we choose $Q_{13}$ to be
a norm 1 vector in the 1-dimensional vector space
$\mathrm{span}({\overline{Q}_{13}}^{\perp},I_{2}(C))^{\perp}$, where
orthogonality is in the monomial basis for $\mathrm{Sym}^{2}H^{0}(K_{C})$.
Similarly we take $H^{0}(2K_{C})$ to be the orthogonal complement of
$I_{2}(C)$ in $\mathrm{Sym}^{2}H^{0}(K)$, and set
$V:=\mathrm{span}(H^{0}(2K)\otimes[Q_{13}],I_{2}(C))^{\perp}$, where
orthogonality is in the monomial basis on $\mathrm{Sym}^{4}H^{0}(K_{C})$. This
is a $\dim H^{0}(4K_{C})-\dim H^{0}(2K_{C})=4(g-1)$ dimensional subspace of
$\mathrm{Sym}^{4}H^{0}(K)$. Evaluating these quartics on the $4(g-1)$ points
of $D_{1}+D_{3}$ we get a square $4(g-1)\times 4(g-1)$ matrix. We denote its
lowest singular value by $B$.

It is worth noting that the numerical azygecy guarantees the exact azygecy.

######  Proposition 5.8.

With the settings above, $(D_{2},D_{4})$ is

$\frac{\sqrt{g-1}(2+2\sqrt{2})}{\sqrt{2}-\sqrt{4g-4}\max(|M_{12}|,|M_{34}|)\frac{\epsilon}{B}}\max(|M_{12}|,|M_{34}|)\epsilon$

approximately syzygetic.

###### Proof.

For $(i,j)=(1,2),(3,4)$, denote by $Q^{\prime}_{ij}$ the quadric representing
the lower $\dim I_{2}(C)+1$ singular vector of $M_{ij}$. Then
$Q^{\prime}_{ij}$ represents a hyperplane in $|2K|$ whose gap from the points
of $D_{i}+D_{j}$ is at most $\epsilon$ by definition. If we look at the gap
between the sum of the divisors and the quartic
$Q^{\prime}_{12}Q^{\prime}_{34}$ we have that

$\displaystyle{\mathfrak{g}}$
$\displaystyle(D_{1}+D_{2}+D_{3}+D_{4},Q^{\prime}_{12}Q^{\prime}_{34})^{2}=\sum_{x\in
D_{1}+D_{2}+D_{3}+D_{4}}|Q^{\prime}_{12}(x)|^{2}\cdot|Q^{\prime}_{34}(x)|^{2}$
$\displaystyle=$ $\displaystyle\sum_{x\in
D_{1}+D_{2}}|Q^{\prime}_{12}(x)|^{2}\cdot|Q^{\prime}_{34}(x)|^{2}+\sum_{x\in
D_{3}+D_{4}}|Q^{\prime}_{12}(x)|^{2}\cdot|Q^{\prime}_{34}(x)|^{2}$
$\displaystyle<$
$\displaystyle\,\,\epsilon^{2}\deg(D_{1}+D_{2})|M_{34}|^{2}+\epsilon^{2}\deg(D_{3}+D_{4})|M_{12}|^{2}.$

Similarly, we also have

$\displaystyle{\mathfrak{g}}(D_{1}+D_{3},Q^{\prime}_{12}Q^{\prime}_{34})^{2}=\sum_{x\in
D_{1}+D_{3}}|Q^{\prime}_{12}(x)|^{2}\cdot|Q^{\prime}_{34}(x)|^{2}<\,\,\epsilon^{2}\deg(D_{1})|M_{34}|^{2}+\epsilon^{2}\deg(D_{3})|M_{12}|^{2}.$

Recalling $V$ and $B$ from Definition 5.7, the norm of the projection of
$[Q^{\prime}_{12}Q^{\prime}_{34}]$ on $V$ is bounded by
$\epsilon_{R}:={\mathfrak{g}}(D_{1}+D_{3},Q^{\prime}_{12}Q^{\prime}_{34})/B$.
i.e., there exists a quadric $Q$ and a quartic $R$ so that
$Q^{\prime}_{12}Q^{\prime}_{34}-QQ_{13}+R\in I_{4}(C)$. The norm of $R$ is
bounded by $e_{R}$, and
${\mathfrak{g}}(D_{1}+D_{3},R)<{\mathfrak{g}}(D_{1}+D_{3},Q^{\prime}_{12}Q^{\prime}_{34})$.
This gives the following bound for the gap $D_{1}+D_{2}+D_{3}+D_{4}$ and
$QQ_{13}$:

$\displaystyle{\mathfrak{g}}$
$\displaystyle(D_{1}+D_{2}+D_{3}+D_{4},QQ_{13})<{\mathfrak{g}}(D_{1}+D_{3},Q^{\prime}_{12}Q^{\prime}_{34})+{\mathfrak{g}}(D_{1}+D_{2}+D_{3}+D_{4},Q^{\prime}_{12}Q^{\prime}_{34})$
$\displaystyle<$
$\displaystyle(\sqrt{4g-4}+\sqrt{8g-8})\cdot\max(|M_{12}|,|M_{34}|)\epsilon=\sqrt{g-1}(2+2\sqrt{2})\max(|M_{12}|,|M_{34}|).$

We now observe that

$|Q^{\prime}_{12}|=|Q^{\prime}_{34}|=1\quad\Rightarrow\quad
1\geq|Q^{\prime}_{12}Q^{\prime}_{34}|\geq\frac{1}{\sqrt{2}}\quad\Rightarrow\quad\geq|QQ_{13}|\geq\frac{1}{\sqrt{2}}-\epsilon_{R},$

and so we have $1\geq|Q|\geq\frac{1}{\sqrt{2}}-\epsilon_{R}$ since
$|Q_{13}|=1$. The result follows. ∎

######  Theorem 5.9.

Let $S,\overline{S}$ be two subsets of the set of unordered pairs of odd theta
characteristics of $C$ so that $\\#\overline{S}=2^{2g}-1$. Assume that there
are distinguished elements $\overline{s}\in\overline{S}$, and $s\in S$ such
that

  1. (1)

each pair in $\binom{\overline{S}}{2}$ either intersects in a theta
characteristic, or forms an $A$-numerically azygetic quadruple,

  2. (2)

for any $s^{\prime}\in S\smallsetminus\\{s\\}$, the theta characteristics in
$s,s^{\prime}$ are $\epsilon$-approximately syzygetic,

  3. (3)

for any
${\overline{s}}^{\prime}\in\overline{S}\smallsetminus\\{\overline{s}\\}$, the
theta characteristics in $s,{\overline{s}}^{\prime}$ are $A$-numerically
azygetic,

  4. (4)

$A/e$ is bigger than the bound from Proposition 5.8

then $S$ is a constituent of a Steiner set.

###### Proof.

It follows from Proposition 5.1 and Proposition 5.8. ∎

##  6\. Numerically certifying the dimensions of $V_{2,\alpha}$ and their
intersection

In this section, we assume that we are given theta hyperplanes up to bounded
accuracy as points in $|K_{C}|$, and we show how to bound from below
$\mathrm{codim}V_{2,\alpha}$ and
$\cap_{\alpha\in\mathrm{Jac}(C)[2]\smallsetminus\\{0\\}}V_{2,\alpha}$, which
are the fundamental quantities of the present work, introduced in Definition
2.2. Assuming the certification of theta hyperplanes and Steiner sets, these
bounds are the only remaining obstacles in the proof of Theorem 2.7.

We collect some notation for the following propositions. Recall the definition
of a Steiner set $\Sigma_{\alpha}$ for a nonzero 2-torsion point $\alpha$ from
Equation 1. Assume that the linear forms $l_{\theta},l_{\theta+\alpha}$
representing theta hyperplanes in the Steiner set $\Sigma_{\alpha}$ are
normalized and known up to accuracy $\epsilon$ in $L_{2}$-norm in the monomial
basis, and let $q_{\theta,\theta+\alpha}=l_{\theta}\cdot l_{\theta+\alpha}$.
Fixing a Steiner set $\Sigma_{\alpha}$, we denote by $M_{\alpha}$ the matrix
whose rows are the coordinate vectors of the quartic forms
$q^{2}_{\theta,\theta+\alpha}$ in the monomial basis of
$\mathrm{Sym}^{2}\mathrm{Sym}^{2}H_{0}(K)$. Let $\\{u_{\alpha,i}\\}_{i}$ be an
orthonormal basis for the kernel of $M_{\alpha}$. Last but not least, recall
from Equation 6 that $d_{\alpha}$ denotes the expected codimension of
$V_{2,\alpha}$ inside $\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})$. We set

$r:=d_{\alpha}+\dim\mathrm{Sym}^{n}\mathrm{Sym}^{2}H^{0}(K_{C}+\alpha)=\dim
I_{n}(\Delta(H^{0}(K_{C}+\alpha))+\dim\mathrm{Sym}^{n}H^{0}(2K_{C}).$

######  Proposition 6.1.

With the notation above, we make the following statements.

  1. (1)

The quadratic forms $q_{\theta,\theta+\alpha}=l_{\theta}\cdot
l_{\theta+\alpha}$ are accurate up to $2\epsilon+\epsilon^{2}$.

  2. (2)

The matrix $M_{\alpha}$ is accurate up to $\epsilon^{\prime}:=\sqrt{2}\cdot
2^{g-1}(2(2\epsilon+\epsilon^{2})+(2\epsilon+\epsilon^{2})^{2})$ in Frobenius
norm.

  3. (3)

If the $r$th singular value of $M_{\alpha}$ is greater than
$\epsilon^{\prime}$ then the codimension of $V_{2,\alpha}$ is $d_{\alpha}$.

  4. (4)

If the hypothesis of Item 3 holds then the basis $\\{u_{\alpha,i}\\}_{i}$ is
accurate up to
$\epsilon^{\prime\prime}:=\frac{\sigma_{r+1}+\epsilon^{\prime}}{\sigma_{r}-\epsilon^{\prime}}$.

###### Proof.

To prove Item 1 note that the error in the norm of $l_{\theta}\otimes
l_{\theta+\alpha}$ is bounded by $2\epsilon\cdot 1+\epsilon^{2}$. Thus, when
taking the symmetric product without normalizing, the norm is between this
quantity and half of it. Item 2 follows from bounding the operator norm of
$M_{\alpha}$ by its Frobenius norm. By Lemma 2.3, we already have an
inequality in the other direction, so the equality in Item 3 holds. Finally,
Item 4 follows from Proposition 4.3. ∎

######  Proposition 6.2.

With the notations above, assume that $\\{u_{{\alpha_{k}},i}\\}_{i,k}$ is
accurate up to $\epsilon^{\prime\prime}$ for the Steiner set
$\Sigma_{\alpha_{k}}$ for each $1<k<m$. Assume that the
$(\dim\ker\pi^{2}_{K_{C}})$th singular value of the matrix whose rows are
$\\{u_{{\alpha_{k}},i}\\}$ is greater than
$\epsilon^{\prime\prime}\sqrt{d_{\alpha}\cdot m}$ then
$\dim\cap_{k=1}^{m}V_{2,\alpha_{k}}C$ is at most $\dim\ker\pi^{2}_{K_{C}}$.
Therefore these spaces are equal.

###### Proof.

This follows from the definition of singular value. ∎

##  7\. A small matter of programming: Concluding the proof of the main
theorem

As explained in the introduction, the object of this section is fourfold. We
discuss the three sources of numerical errors in our certification process:
first, inherent errors in using numerical linear algebra, which we treat as an
external certification oracle, next inaccurate input data, and finally
numerical errors incurred during the certification. Addressing these sources
of errors concludes the proof of the main theorem – Theorem 2.7. Finally, we
discuss testing the certification code. Before we delve into its details, we
give a bird’s-eye review of the structure of the certification code and its
auxiliary output.

We have a computer program for each genus that we work out. The genus 6 case
is certified by a program running on Wiman’s sextic, whereas the genus 7 case
is certified by a program running on the Fricke-Macbeath curve. The programs
are identical except for the parts which handle the automorphism groups of the
curves and the levels of accuracy. There are three types of functions in both
programs:

  * •

Service functions: e.g. `read_array_cmplx` which reads an array of complex
numbers from a file, or `det_mat` which computes the determinant of a matrix.

  * •

Functions which deal with the group action: e.g. `S4linear_action` which takes
an element in the symmetric group $S_{4}$ and gives its projective linear
action on the canonical system of Wiman’s sextic, and similarly for
`D14linear_action` and the action of the dihedral group $D_{14}$ on the
canonical system of the Fricke-Macbeath curve.

  * •

Functions in the certification code: e.g.
`certify_pt_on_H_cap_C_and_compute_M_T`, which does exactly what its name
suggests, namely certifying that an input point $P^{\prime}$ is close both to
a hyperplane $H^{\prime}$ and to the curve, and computing the matrix
$M_{T^{\prime}}$.

In this section, we will mostly cover the last item. We will only cover the
first two when discussing the testing of our code. As the code computes
numerical estimations, there are a lot of numbers involved:

  * •

The dimensions of spaces involved, such as `DIM_H0_2K`, which is $\dim
H^{0}(K_{C})$.

  * •

Bounds used in hard asserts such as `PT_HYPERPLANE_INCIDENCE_EPSILON` which
all of the input data have to pass.

  * •

Soft bounds used to decide whether some input datum is accurate enough, such
as
`GOOD_MAX_EPSILON_0`. We discard data which does not pass these bounds. These
constants are preceded by the prefix `GOOD_`.

  * •

Bounds such as `SERVICE_DISTINCT_PTS_COORDINATES_MIN_DIFF` which is used to
check whether points are distinct. These later pass the hard assert
`IDENTICAL_PTS_EPSILON` which is several orders of magnitudes smaller.
Similarly, we have bounds like
`SERVICE_FOUR_THETAS_GAP_LOOK_UP`, which indicates the spectral gap we look
for. These constants are preceded by the prefix `SERVICE_`.

  * •

The maximal error in computing the singular values of the matrices such as
$M_{T^{\prime}}$ or $J_{F^{\prime}}$, which is gathered in the class
`VerifiedZSVD` and its member functions.

  * •

The top singular value, the lowest numerically non-trivial singular value and
the one following it, for the matrix arising in Proposition 6.2, and the
maximum of these values for the matrices involved in all incarnations in
Proposition 6.1. We expect the location of the last numerically non-trivial
singular value to agree with the theoretical rank computations.

The code outputs some auxiliary information:

  * •

The group orbits of hyperplanes.

  * •

The points together with their hyperplanes.

  * •

The inverse of the above: the theta hyperplanes that are accurate enough with
their points.

  * •

The orbit of the automorphism group action on the two-torsion points. We use
only certain part of the group: the symmetric group $S_{4}$ for Wiman’s
sextic, and the dihedral group $D_{14}$ for the Fricke-Macbeath curve.

###  7.1. Errors inherent to using SVD

Per `LAPACK` documentation, the error in the singular values is the machine
precision `E` (for IEEE-754 double `E`=$2^{-53}\sim$`1.1e-16`) times the top
computed singular value. When bounding the norm of an operator, the maximal
error is exactly this value as well. There are two components in the
certification process that we care about more than the singular values. The
first one is in Proposition 4.31, where $R$ is a solution to the linear system
in Equation 22. The second one is Proposition 6.2, which uses the left pseudo
kernel of the matrices from Proposition 6.1. See Remark 4.24 for the details
of the accuracy bounds in these computations.

###  7.2. Numerical errors in the input data

Part of our input is a generating set for $I_{2}(C)$, which is comprised of
quadrics $Q_{k}$ with integer coefficients. Being small integers, they do not
come with any inherent errors.

Another part of our input is a set of putative theta hyperplanes
$H^{\prime}_{i}$ and their approximate intersection points with the curve $C$,
given by $P^{\prime}_{ij}$. See Definition 4.23 for the notation. The
following quantities dependant on $H^{\prime}_{i}$ and $P^{\prime}_{ij}$ have
computable bounds:

$H^{\prime}_{i}\cdot
P^{\prime}_{ij},\quad\sqrt{\sum_{k}|{P^{\prime}_{ij}}^{t}Q_{k}P^{\prime}_{ij}|^{2}},\quad\sqrt{\sum_{j}|H^{\prime}_{i}T^{\prime}_{ij}|^{2}}.$

Note that the last one is not part of the input data per se. We compute and
certify the first two bounds in the function
`certify_pt_on_H_cap_C_and_compute_M_T`. These bounds are
`PT_HYPERPLANE_INCIDENCE_EPSILON` and `POINT_ON_QUADRIC_EPSILON` respectively,
which attain the bounds `1e-11`, `1e-10` respectively. Some points are
discarded at an earlier step in the function `read_tangency_pts` after failing
the `GOOD_INPUT_PT_HYPERPLANE_INCIDENCE_EPSILON` bound, which is `1e-11`. The
third quantity is only computed in the function
`verify_conic_approximation_1st_order_invertibility`; in practice, the numbers
we get for this quantity are in the segment $[10^{-5},10^{-9}]$, which are by
far the biggest inaccuracies in the input data. We discard the hyperplanes
that give errors bigger than `GOOD_MAX_SQRT_SUM_i_Hk_DOT_Tki_SQ`, which is set
to be `1e-6`. This error dominates all the other input errors, which in turn
dominate the errors introduced by our SVD computation, which dominate errors
introduced by any other numerical computations (which only introduces errors
of the order of magnitude of `E`).

###  7.3. The numerics of Section 4, and the proof of the main theorem

For most objects considered in this section, we do not have only two but three
incarnations as listed in the following.

  * •

The _Platonic incarnation_ : the theta hyperplanes $H_{i}$ and the tangency
points.

  * •

The points $P_{ij}$, which are true intersection points of $H^{\prime}_{i}$
with the curve.

  * •

The numerically computed objects we actually have at our disposal – given by
$P_{ij}^{\prime}$ and $H^{\prime}_{i}$.

The technical part of Section 4 is roughly divided into two parts, bounding
the differences between the Platonic and the second incarnations, and then the
second and the last incarnations. Although one cannot talk about certification
without considering the combination of the two parts, they are somewhat
different in flavour. Loosely speaking, in the first part we establish claims
whose condition is that
$r_{P_{i},x},r_{T_{i},x},r_{R_{i},x},\epsilon_{P_{i},x}/|x|^{3},\epsilon_{R_{i},x}/|x|,c_{i},\sigma_{g-1}({\mathcal{J}}_{F})$
are all “$O(1)$”, where the $r_{\bullet,x}$ and the
$\epsilon_{\bullet,x}/|x|^{\bullet}$ are minima of short multiplicative (this
includes divisions and taking roots) expressions in
$|R|,|W_{R}|,|W_{T}|,q_{m},\sigma_{g-2}(M),\sigma_{g-1}(M_{T}),\sigma_{g}(M_{R})$.
The second part shows how one computes the error in estimating these
quantities. Again, broadly speaking, the claims in the first part would hold
if the errors computed in the second part were small compared to the
quantities in the first part.

The following theorem gives bounds related to the computed hyperplanes
$H_{i}^{\prime}$ and the points $P_{ij}^{\prime}$. The bounds which appear
below do not bound all the input objects of the corresponding types. Rather,
we only consider the number of hyperplanes and the respective points that
suffice to conclude the proof of the main theorem. In particular, the number
of objects should suffice to get the relevant space dimensions.

######  Theorem 7.1.

The computed hyperplanes $H^{\prime}_{i}$ and points $P^{\prime}_{ij}$ which
we used for the certification satisfy the following bounds.

  * •

For any used $i$, we have $\sqrt{\sum_{j}|P^{\prime}_{ij}-\text{closest point
to }P^{\prime}_{ij}\text{ in }H_{i}\cap C|^{2}}<$`1e-6`.

  * •

For any used $i$, we have $|H^{\prime}_{i}-H_{i}|<$`1e-10`.

###### Proof.

Section 4 gives exact bounds on what the relations between “$O(1)$” and
“small” should be. However, the bounds we get on the numbers involved are so
good, that there is no point in performing the meticulous computation. The
functions involved in the certification of these two claims are as follows.

    certify_pt_on_H_cap_C_and_compute_M_T
    compute_T_and_extend_M_T_to_M_R
    compute_R
    compute_norm_W_T_W_R
    verify_conic_approximation_1st_order_invertibility

Which are all called from the function `verify_and_uniqify_pts`. The numerical
bounds involved are:

    MATRICES_M_M_T_MIN_NON_TRIVIAL_SING_VAL .1
    MATRICES_M_M_T_MAX_NON_TRIVIAL_SING_VAL 50 (for Wiman) 10 (for Fricke-Macbeath)
    MATRIX_M_R_MIN_NON_TRIVIAL_SING_VAL .3 (for Wiman) .7 (for Fricke-Macbeath)
    MATRICES_M_M_T_M_R_MAX_TRIVIAL_SING_VAL 1e-11 (for Wiman) 1e-10 (for Fricke-Macbeath)
    MATRIX_M_Qi_T_MAX_TOP_SING_VAL 10 (for Wiman) 5 (for Fricke-Macbeath)
    MIN_NORM_W_T_W_R .01
    MAX_NORM_W_T_W_R 10
    GOOD_MAX_NORM_R 5
    GOOD_MIN_NORM_R .2
    GOOD_MAX_EPSILON_012 1e-10 (for Wiman) 1e-11 (for Fricke-Macbeath)
    MAX_EPSILON_J_F 1e-9 (for Wiman) 5e-11 (for Fricke-Macbeath)
    GOOD_MAX_EPSILON_0 1e-13
    GOOD_MAX_SQRT_SUM_i_Hk_DOT_Tki_SQ 1e-6
    MAX_X_PRIME_NORM 3e-7
    CONIC_APPROXIMATION_J_F_MAX_SING_VAL 2
    CONIC_APPROXIMATION_J_F_MIN_SING_VAL .05

Whereas the maximal error in computing the singular values of the matrices
involved are as follows for the Fricke-Macbeath curve and Wiman’s sextic
respectively:

    conic approx: 6.90502e-17
    M_Qi_T: 2.71928e-16
    M_R: 4.07084e-16
    M: 4.07084e-16
    M_T: 4.07084e-16

    conic approx: 2.54872e-16
    M: 2.34172e-15
    M_Qi_T: 9.98561e-16
    M_T: 2.34172e-15
    M_R: 2.34172e-15

One easily verifies that the claim in Remark 4.34 Item 1 around Equation 25
holds. Thus, by Proposition 4.33 and Theorem 4.35, the quantity
$\sqrt{\sum_{i}|P^{\prime}_{i}-\text{closest point to }P^{\prime}_{i}\text{ in
}H\cap C|^{2}}$ is bounded by
`MAX_X_PRIME_NORM`+`MAX_EPSILON_J_F`/`CONIC_APPROXIMATION_J_F_MIN_SING_VAL`$\
{}^{2}$ \+ quantities which are several orders of magnitude smaller, which is
smaller than `1e-6`.

Turning our attention to $|H^{\prime}_{i}-H_{i}|$, by Theorem 4.36, the bound
is given by `1e-6` $\cdot(\sum_{i}|\langle
T^{\prime}_{i},n_{H^{\prime}}\rangle|$
+${\small\verb|GOOD_MAX_EPSILON_012|}$)+$\mathrm{Tail}_{2}\prod_{i}\eta_{i}($`1e-6`$)$,
where the bound on each $|\langle T^{\prime}_{i},n_{H^{\prime}}\rangle|$ is
given by
`GOOD_MAX_SQRT_SUM_i_Hk_DOT_Tki_SQ`, which is `1e-6`. Thus the highest terms
in the expression are `1e-6`$(g-1)$`1e-6` from the first term, and
$(g-1)$`GOOD_MAX_NORM_R`$\cdot$`1e-6`$\ {}^{2}$ from the latter. ∎

###  7.4. Errors incurred by certifying the Steiner systems in Section 5

The ingredients involved in the certification of the Steiner sets are:

  * •

The accuracy of the approximated tangency points, for which we have bounds
presented above.

  * •

The operator norms of the matrices $M_{ij}$ defined in Definition 5.4.

  * •

The maximal (resp. minimal) computed $\epsilon$ (resp. $A$) for
$\epsilon$-approximately syzygetic (resp. $A$-numerically azygetic) pairs of
pairs of theta characteristics.

  * •

The smallest singular value $B$ defined in Definition 5.7.

Note that the operator norms $|M_{ij}|$ incur computational errors which we
bound via the bound on the computational errors from the tangency points.
Specifically, if a point is known with accuracy $\epsilon$ then the
corresponding line in the matrix $M_{ij}$ is known with accuracy
$2\epsilon+\epsilon^{2}$. Thus by Theorem 7.1 the matrices $M_{ij}$s are known
with accuracy $\sqrt{\frac{4(g-1)}{g-1}}\cdot(2\cdot 10^{-6}+10^{-12})<$`5e-6`
in terms of the Frobenius norm. Referring to the last item, one can bound the
error in computing $B$ by tracing the steps in Definition 5.7. This is done by
bounding the norm of the product of two quadrics between $1/2$ and exactly the
product of the norms, bounding the Frobenius norms, and applying Proposition
4.2. The reason we do not explicitly trace these steps is that we did similar
computations several times before, and the quotient $\epsilon/B$ we get is in
practice much smaller than anything that errors in computing $B$ might
contribute.

######  Proposition 7.2.

For the subset of hyperplanes which pass our asserts, the labeling of the
Steiner sets obtained from the non-certified level 2 structure is correct.

###### Proof.

The combinatorial and numerical frameworks of Proposition 5.1 and Theorem 5.9
are verified in the functions `four_thetas_constraints`,
`certify_one_steiner`, `fill_all_good_pairs`,
`four_theta_on_quartics` and
`fill_torsion_action_table_and_one_rep_per_steiner`. Recall that by Theorem
5.9 there are four numerical estimations at play here: the $\dim I_{2}(C)$-th
singular value of $\epsilon$-approximately syzygetic sets measured exactly
with $\epsilon$, and $B$ from Definition 5.7, the $A$-numerically azygecy
coefficient, and $|M_{ij}|$. They are connected by comparing $A/\epsilon$ with

$\frac{\sqrt{g-1}(2+2\sqrt{2})}{\sqrt{2}-\sqrt{4g-4}\max(|M_{12}|,|M_{34}|)\frac{\epsilon}{B}}\max(|M_{12}|,|M_{34}|).$

The bounds on these coefficients are:

$\epsilon<$`1e-8` | `EPSILON_SYZYGECY_TRIVIAL_BOUND`
---|---
$|M_{ij}|<$`3` | `FOUR_THETAS_MIN_NON_TRIVIAL_SING_VAL`
$B>$`2e-4` | `GOOD_MIN_B`
$A>$`.2` | `A_AZYGECY`

.

The maximal errors incurred in the singular value computations of $M_{ij}$,
(i.e., in all the constants beside $B$) are `6.90502e-17` for the Fricke-
Macbeath curve, and `2.54872e-16` for Wiman’s sextic. As explained above and
in more detail in Definition 5.7, there are several singular value
computations involved in computing $B$; the biggest error in these
computations is `2e-14` for Wiman’s sextic and `5e-15` for the Fricke-Macbeath
curve. ∎

###  7.5. Errors incurred by certifying the Section 6

Recalling the proof outline given in Section 2, we use the fact that by Lemma
2.8 it suffices to prove the claim for one curve in each genus. In Lemma 2.3
we have established an upper bound on the dimension $\dim V_{2,\alpha}$. Under
the assumption that this upper bound is an equality, we would have
$V_{2,\alpha}\supset\ker\pi^{2}_{K_{C}}$, which shows that
$\cap_{\alpha}V_{2,\alpha}\supset\ker\pi^{2}_{K_{C}}$. Thus this gives us a
lower bound on the dimension of the intersection. We now complete the proof of
the main theorem Theorem 2.7.

###### Proof of Theorem 2.7.

We first give a lower bound on $\dim V_{2,\alpha}$, which turns out to be
identical to the upper bound displayed in Equation 6. Next, we give a lower
bound on
$\dim\mathrm{span}\cup_{\alpha}V_{2,\alpha}^{\perp}=\dim(\cap_{\alpha}V_{2,\alpha})^{\perp}$.
Thus, we get an upper bound on the $\dim(\cap_{\alpha}V_{2,\alpha})$, which
then coincides with $\dim\ker\pi^{2}_{K_{C}}$, which is the lower bound we
have for the dimension of this intersection. By Proposition 7.2 our Steiner
sets are correct, so we can compute the lower bound using the hyperplanes we
certified. The accuracy of each such hyperplane is bounded by `1e-10`. Hence
the accuracy of each line of the matrix representing the Steiner set in
$\mathrm{Sym}^{2}\mathrm{Sym}^{2}H^{0}(K_{C})$, where each line is the
$\mathrm{Sym}^{2}$ of the product of hyperplanes, is bounded by `4e-10`. We
here are ignoring terms of magnitude `1e-20`. The matrices in question have at
most $2^{g-2}(2^{g-1}-1)<2^{2g-2}$ lines, and so the Frobenius norm of the
error in computing these matrices is bounded by $64\cdot 4\cdot
10^{-10}<3\cdot 10^{-8}$. To begin with, this gives us a bound on the error
computed singular values, which are

    one_steiner_max_sing_val = 6.02764,
    one_steiner_expected_min_non_trivial_sing_val = 0.0379421,
    one_steiner_expected_max_trivial_sing_val = 1.03606e-13

for the Fricke-Macbeath curve, where the maximal error incurred in this
computation is `6.69202e-16`;

    one_steiner_max_sing_val = 7.78234,
    one_steiner_expected_min_non_trivial_sing_val = 0.00833475,
    one_steiner_expected_max_trivial_sing_val = 6.46229e-14

for Wiman’s sextic, where the maximal error incurred in this computation is
`8.64014e-16`.

We have now proven the dimension claim for the $V_{2,\alpha}$. Next, by
Proposition 4.2, the error in the bases of the kernels is bounded by
$\frac{|1.04\cdot 10^{-13}-3\cdot 10^{-8}|}{.04-3\cdot 10^{-7}}<2\cdot
10^{-7}$ for the Fricke-Macbeath curve, and $\frac{|6.5\cdot 10^{-14}-3\cdot
10^{-8}|}{.01-3\cdot 10^{-8}}<1.1\cdot 10^{-6}$ for Wiman’s sextic.

For the Fricke-Macbeath curve, we intersect only $96$ many $V_{2,\alpha}$,
whereas for Wiman’s sextic, we intersect only $56$ of them. In both cases,
this means the Frobenius norm of the error is bounded by $\sqrt{100}\cdot
1.1\cdot 10^{-6}<1.1\cdot 10^{-5}$. The singular values we compute for the
corresponding matrices are

    sing[0] = 7.26053 sing[170]  = 1.17418 svd.s_[171] = 1.86973e-12

for the Fricke-Macbeath curve, where the error incurred in this computation is
`8.06081e-16`; and

    sing[0] = 9.21738 sing[119]  = 1.07752 svd.s_[120] = 4.10144e-12

for Wiman’s sextic, where the error incurred in this computation is
`1.02333e-15`. ∎

###  7.6. Testing

Although testing cannot prove that a program is correct, it is an important
part of any programming endeavour, as it may provide a probabilistic assurance
for the correctness of the code. First, it is important to note that the
entire point of Section 4 and Section 5 is to certify the computations which
were done in a completely different context. This full separation between the
way the input that was generated and our certification code is an important
test in itself.

Next, there are two types of numerical tests we perform. The first of them is
verifying known facts, like the essential dimensions appearing in the matrices
$M^{\prime},M^{\prime}_{T},M^{\prime}_{R}$ being $g-2,g-1,g$, or the ones in
Section 6 being smaller than the dimension bounds given in Section 2. An
inequality of these numbers in the reverse direction would have just meant
that our theorem is wrong. Another verification is via the fact that the
point-hyperplane (approximate) incidence relation is equivariant under the
symmetry group action. The second type is an incarnation of the rule of thumb
“numbers behave randomly unless they have a reason not to”. The incarnations
of this rule in our code are:

  1. (1)

The spectral gaps on the essential dimensions of the spans of the rows or
columns of $M^{\prime},M^{\prime}_{T},M^{\prime}_{R}$, or the matrices
checking azygecy, or the matrices certifying the computations in Section 5 are
big. That is to say: if we denote this dimension by $d_{\text{ess}}$ for each
of the cases then in these cases $\sigma_{d_{\text{ess}}+1}$ is bounded from
above by `1e-10`, `1e-12`, `1e-11`,`1e-12` respectively, whereas
$\sigma_{d_{\text{ess}}}$ is bounded from below by `.3`,`.01`, `0.00833475`,
`1.07752` respectively, and $\sigma_{1}$ is bounded from above by `50`, `10`,
`7.78234`, `9.21738` respectively. This is in contrast to the case for the
spectral gap for checking syzgecy, where $d_{\text{ess}}=\dim I_{2}(C)-1$ and
$\sigma_{d_{\text{ess}}+1}$ is bounded from above merely by `1e-8`, which is
far bigger a bound than the corresponding bound for azygecy. The fact that
both this bound and the bound `GOOD_MAX_SQRT_SUM_i_Hk_DOT_Tki_SQ` on
$\sqrt{\sum_{i}|H^{\prime}_{k}\wedge T^{\prime}_{ki}|^{2}}$ are relatively big
reflects the fact that the distance of $\sum_{j}P^{\prime}_{ij}$ from being
half canonical is big, due to the way the $P^{\prime}_{ij}$s are computed –
see the discussion above Theorem 4.36.

  2. (2)

Points and hyperplanes are in practice far off from one another:

        SERVICE_MIN_FUBINI_STUDY_HYPERPLANES_DISTANCE 1e-5
    SERVICE_DISTINCT_PTS_COORDINATES_MIN_DIFF 1e-3 (for Wiman) 5e-4 (for Fricke-Macbeath)

We use the `SERVICE_` prefix here since `IDENTICAL_PTS_EPSILON` and
`FUBINI_STUDY_HYPERPLANES_EPSILON` are reserved for points appearing multiple
times due to the `Julia` code or the group action.

  3. (3)

One step of the iterative gradient descent (via inverse function theorem)
algorithm alluded to in Remark 4.34 Item 4 is performed, and it improves the
solution by at least a factor of `GRAD_DESCENT_MIN_IMPROVEMENT` which is
$0.85$.

## References

  * [AC21] Daniele Agostini and Lynn Chua.  Computing Theta Functions with Julia.  Journal of Software for Algebra and Geometry, 11:41–51, 2021.
  * [ACGH85] Enrico Arbarello, Maurizio Cornalba, Phillip Griffiths, and Joe Harris.  Geometry of algebraic curves: volume I.  Springer-Verlag, New York, 1985.
  * [Aro64] Siegfried Heinrich Aronhold.  Über den gegenseitigen Zusamemmenhang der 28 Doppeltangenten einer allgemeiner Curve 4ten Grades.  Monatberichter der Akademie der Wissenschaften zu Berlin, pages 499–523, 1864.
  * [BT17] Paul Breiding and Sascha Timme.  HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia.  In International Congress on Mathematical Software, 2017.
  * [CEFS13] Alessandro Chiodo, David Eisenbud, Gavril Farkas, and Frank-Olaf Schreyer.  Syzygies of torsion bundles and the geometry of the level $l$ modular variety over $\overline{\mathcal{M}}_{g}$.  Inventiones mathematicae, 194(1):73–118, 2013.
  * [ÇFM23] Türkü Özlüm Çelik, Samantha Fairchild, and Yelena Mandelshtam.  Crossing the Transcendental Divide: From Translation Surfaces to Algebraic Curves.  Experimental Mathematics, 0(0):1–19, 2023.
  * [ÇKRS19] Türkü Özlüm Çelik, Avinash Kulkarni, Yue Ren, and Mahsa Sayyary Namin.  Tritangents and their space sextics.  Journal of Algebra, 538:290–311, 2019.
  * [Cob29] Arthur B. Coble.  Algebraic Geometry and Theta Functions.  American Mathematical Society, Colloquium Publications Series, 1929.
  * [CS03] Lucia Caporaso and Edoardo Sernesi.  Characterizing curves by their odd theta-characteristics.  Journal für die reine und angewandte Mathematik, 2003(562):101–135, 2003.
  * [DHB+04] Bernard Deconinck, Matthias Heil, Alexander Bobenko, Mark van Hoeij, and Marcus Schmies.  Computing Riemann theta functions.  Mathematics of Computation, 73:1417–1442, 2004.
  * [Dol03] Igor Dolgachev.  Lectures on Invariant Theory.  London Mathematical Society Lecture Note Series. Cambridge University Press, 2003.
  * [Dol12] Igor V. Dolgachev.  Classical Algebraic Geometry: A Modern View.  Cambridge University Press, 2012.
  * [Dv01] Bernard Deconinck and Mark van Hoeij.  Computing Riemann matrices of algebraic curves.  Physica D: Nonlinear Phenomena, 152-153:28–46, 2001.  Advances in Nonlinear Mathematics and Science: A Special Issue to Honor Vladimir Zakharov.
  * [FJK19] Jörg Frauendiener, Carine Jaber, and Christian Klein.  Efficient computation of multidimensional theta functions.  Journal of Geometry and Physics, 141:147–158, 2019.
  * [FK11] Jörg Frauendiener and Christian Klein.  Algebraic Curves and Riemann Surfaces in Matlab.  In Alexander I. Bobenko and Christian Klein, editors, Computational Approach to Riemann Surfaces, pages 125–162. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.
  * [Fri99] Robert Fricke.  Ueber eine einfache Gruppe von 504 Operationen.  Mathematische Annalen, 52(2–3):321–339, 1899.
  * [GB21] Sohrab Ganjian and Nils Bruin.  RiemannTheta: A SageMath package for evaluating Riemann theta functions, 2021.
  * [Guà02] Jordi Guàrdia.  Jacobian Nullwerte and algebraic equations.  Journal of Algebra, 253:112–132, 2002.
  * [Hid18] Ruben A. Hidalgo.  About the Fricke–Macbeath curve.  European Journal of Mathematics, 4(1):313–325, 2018.
  * [Leh15] David Lehavi.  Effective Reconstruction of Generic Genus 4 Curves from Their Theta Hyperplanes.  International Mathematics Research Notices, 2015(19):9472–9485, 2015.
  * [Leh22] David Lehavi.  Effective Reconstruction of Generic Genus 5 Curves from their Theta Hyperplanes.  Experimental Mathematics, 0(0):1–15, 2022.
  * [LS96] H. Lange and E Sernesi.  Quadrics containing a Prym-canonical curve.  Journal of Algebraic Geometry, 5(2):387–399, 1996.
  * [Mac65] Alexander Murray Macbeath.  On a Curve of Genus 7.  Proceedings of the London Mathematical Society, 15(1):527–542, 1965\.
  * [MN19] Pascal Molin and Christian Neurohr.  Computing period matrices and the Abel-Jacobi map of superelliptic curves.  Mathematics of Computation, 88:847–888, 2019.
  * [Muk88] Shigeru Mukai.  Curves, K3 Surfaces and Fano $3$-folds of Genus $\leq$ 10.  In Algebraic Geometry and Commutative Algebra, pages 357–377. Academic Press, 1988.
  * [Mum07] David Mumford.  Tata Lectures on Theta I (reprint of 1983 edition).  Birkhäuser, Boston, 2007.
  * [NBZ19] Jeroen Sijsling Nils Bruin and Alexandre Zotine.  Numerical computation of endomorphism rings.  In ANTS XIII: Proceedings of the Thirteenth Algorithmic Number Theory Symposium, The Open Book Series 2-1, pages 155–171, 2019.
  * [Sch] Frank-Olaf Schreyer.  KoszulDivisorOnPic14M8 – Construction of random curves of degree 14 and genus 8.
  * [vW06] Paul B. van Wamelen.  Computing with the analytic Jacobian of a genus 2 curve.  In Wieb Bosma and John Cannon, editors, Discovering Mathematics with Magma: Reducing the Abstract to the Concrete, pages 117–135. Springer Berlin Heidelberg, Berlin, Heidelberg, 2006.
  * [Wim96] Anders Wiman.  Zur Theorie der endlichen Gruppen von birationalen Transformationen in der Ebene.  Mathematische Annalen, 48(1):195–240, 1896.