\documentclass{amsart}

\input{preamble}

\begin{document}

\title{A polynomial-time algorithm to compute theta functions of Schottky groups}

\author{Marc Masdeu}
	\address{Departament de Matem\`atiques\\Universitat Aut\`onoma de
		Barcelona\\Edifici C, 08193 Bellaterra, Barcelona, Catalonia}
    \address{Centre de Recerca Matemàtica\\Bellaterra, Barcelona, Catalonia}

\email{marc.masdeu@uab.cat}


\subjclass[2020]{14Q99, 14H42 (14K25, 11Y99, 30F10)}

\date{}


\begin{abstract}
  Abelian integrals, holomorphic differentials, period matrices and the Schottky--Klein prime function
  of a Riemann surface uniformized by a classical Schottky group are given by Poincar\'e theta series, which
  are sums over the elements of a free group. The algorithms that evaluate these series sum a truncated, possibly
  adaptively pruned, portion of the Cayley tree of the group, and the number of terms needed grows exponentially with
  the number of correct digits. We observe that the sum over all group elements starting with a fixed generator is a
  holomorphic function on the complement of a disc, and that these functions satisfy an exact linear functional
  equation coming from the self-similarity of the Cayley tree, given by a transfer operator built from the generators.
  Representing them by truncated Taylor expansions turns the evaluation of the series into a finite-dimensional
  linear problem, which we solve iteratively. The resulting algorithm computes the integrals $\Omega_n$, the period
  matrix and the prime function to $P$ bits of precision; for a fixed Schottky datum the number of arithmetic operations is
  cubic in $P$ and in the genus. We prove that it converges geometrically under a standard Schottky-type hypothesis on the
  group, which is closely related to the absolute convergence of the series. Along the way we give a complete proof of the cross-ratio formula for the period matrix, including a term
  that is missing from one version of it in the literature. We compare with the existing algorithms, and illustrate the method by
  reproducing a published value of the prime function and computing it to $30$ digits, and by computing period
  matrices of genus up to six to about $25$ digits.
\end{abstract}

\maketitle

\section*{Disclosure of AI use}
{\sloppy
This work was prepared with substantial assistance from AI language models, used through the Claude Code
command-line agent: Claude Opus~5.5, Claude Sonnet~5.5
and Claude Haiku~5.5, all developed by Anthropic.
They were used to draft and edit the text of this article, to survey the literature and to suggest references,
to propose and check proofs and calculations, to write and run the Sage code and the numerical experiments
reported here (including the benchmarks and the examples of Section~\ref{sec:applications}), and, in separate sessions
acting as referees, to criticise earlier versions of the manuscript. No AI system is listed as an author.
The author conceived the project and the algorithm's origin in the $p$-adic setting, directed all the work, and takes full
responsibility for the content of the article, including the correctness of its statements and the accuracy
of its references. The software is available at \url{https://github.com/mmasdeu/archimedean_theta}.\par}

\section{Introduction}

Throughout we follow the notations from~\cite{Schmies2011}, which we now recall.
Let $G$ be a (classical) Schottky group, with loxodromic generators
\[
\sigma_1, \ldots, \sigma_N \in  \PGL_2(\CC).
\]
Denote by $A_n$ and $B_n$ the repelling and attracting fixed points of $\sigma_n$, respectively, and by
$\mu_n = \sigma_n'(B_n)$, with $|\mu_n|<1$, its loxodromic factor; thus
\begin{equation}\label{eq:multiplier}
\frac{\sigma_n(z) - B_n}{\sigma_n(z) - A_n} = \mu_n\,\frac{z - B_n}{z - A_n}.
\end{equation}
For convenience, we will also write $\sigma_{-i} = \sigma_i^{-1}$.

Since $G$ is freely generated by the $\sigma_i$, each $\sigma \in G$ has a unique representation of the form
\[
\sigma = \sigma_{i_1} \cdots \sigma_{i_k},\quad i_j \in \{\pm 1, \ldots, \pm N\},
\]
where we require $i_j \neq -i_{j+1}$ for all $j$. The list of integers $[i_1,\ldots, i_k]$ will be called the
\emph{Tietze word} associated to $\sigma$.

Denote by $H_n=\langle\sigma_n\rangle$ the subgroup generated by $\sigma_n$. Every coset in $G/H_n$ has a unique
representative whose Tietze word does not end in $\pm n$, that is, a representative with Tietze word of the form
\[
[i_1,\ldots, i_k],\quad i_j \in \{\pm 1, \ldots, \pm N\},\quad i_k \neq \pm n.
\]

Consider now the series
\[
\omega_n(z) = \sum_{\sigma \in G/H_n}\left( \frac{1}{z - \sigma(B_n)} -  \frac{1}{z - \sigma(A_n)}\right)dz,
\]
which is a holomorphic differential provided that the sum converges absolutely. The series $\omega_n(z)$ can be integrated
to give
\[
\Omega_n(z) = \int_\infty^z \omega_n = \sum_{\sigma \in G/H_n}\log\frac{z - \sigma(B_n)}{z - \sigma(A_n)}.
\]
One may also obtain the period matrix $(B_{mn})_{1\leq m,n\leq N}$, where $B_{mn}=\int_z^{\sigma_m(z)}\omega_n$ (this
is independent of $z$ and of the path, up to $2\pi i\ZZ$) is given by
\begin{equation}\label{eq:periodmatrix}
B_{mn} = \delta_{mn} \log \mu_n + \sum_{\substack{\sigma \in H_m\backslash G / H_n\\ \sigma \notin H_n \text{ if } m=n}}
\log \left(B_m:A_m:\sigma(B_n):\sigma(A_n)\right),
\end{equation}
with $\delta_{mn}$ being Kronecker's delta function and $(a:b:c:d)$ the cross-ratio
\[
(a:b:c:d) = \frac{(a-c)(b-d)}{(a-d)(b-c)}
\]
of the points $a,b,c,d$. The sum runs over the double cosets, represented by the reduced words which neither begin with
$\sigma_{\pm m}$ nor end with $\sigma_{\pm n}$, \emph{together with the trivial double coset when $m\neq n$}, which
contributes $\log(B_m:A_m:B_n:A_n)$. This is the formula of~\cite[Theorem~30]{BobenkoIntro2011}
(see also~\cite[(5.14)]{CrowdyMarshallMcurves2011}); the term of the trivial double coset for $m\neq n$ is easy
to overlook, and is missing in~\cite[(4.4)]{Schmies2011} (see Remark~\ref{rem:identity}). We give a complete proof in
Proposition~\ref{prop:periods}.

These series go back to Schottky, Burnside and Baker~\cite{Burnside1891,Baker1897} (see
also~\cite{Hejhal1972,BobenkoIntro2011}), and algorithms to evaluate them have been proposed in the literature; we
review them in Section~\ref{sec:existing-algorithms}. For example,
\cite{Schmies2011} studies the convergence of the series $\omega_n$ and gives a priori error bounds, \cite{Lyamaev2022}
sharpens the estimates and the summation strategy, and~\cite{YamamotoMiyoshi2026} does the same for the product
defining the Schottky--Klein prime function. All the algorithms which evaluate the series rely on the summation of a
certain finite subset of their terms, and the number of terms which are needed grows exponentially with the target
precision (see Section~\ref{sec:existing-algorithms}). For the prime function of a planar circular domain there are
also methods which avoid the series altogether, by collocation~\cite{CrowdyMarshall2007} or by solving boundary value
problems~\cite{CKGN2016}; they are specific to that setting (Schottky doubles), and for the first of them no convergence proof is
available~\cite{CrowdyMarshallMcurves2011}.

The goal of this note is to describe a polynomial-time algorithm to compute the above quantities,
by means of an iterative procedure. In fact, we will give an algorithm to compute $\Omega_n(z)$ as a
function of $z$. By computing its derivative we will recover $\omega_n(z)$. We will also give an algorithm
to compute the period matrix, and the Schottky--Klein prime function.

\subsection*{Origin, main ideas and results}
The algorithm is the archimedean analogue of the algorithm of Masdeu and Xarles~\cite{masdeu-xarles} for $p$-adic theta
functions of Mumford curves, which is where this work originated (see Section~\ref{sec:applications}).

The Cayley tree of $G$ is self-similar: the subtree of descendants of a vertex only depends on the last letter of the
corresponding word. Consequently, the sum of the terms of the series over all words \emph{beginning} with a given
generator $\sigma_i$ is one holomorphic function $F^{(i)}_n$ on the complement of a closed disc, and these $2N$ functions
satisfy the linear system
\[
F^{(i)}_n = E^{(i)}_n + \sum_{j\neq -i} T_{ij} F^{(j)}_n,
\]
where $E_n^{(i)}$ is an explicit logarithm of a ratio of linear functions, and $T_{ij}$ is a ``transfer operator'':
the composition with the M\"obius transformation $\sigma_i^{-1}$ followed by a normalization
(Theorem~\ref{thm:functional-equation}). Expanding each $F^{(i)}_n$ in a Taylor series in a coordinate that
is natural to the disc, the operators $T_{ij}$ become explicit matrices, which can be computed in $O(p^2)$ operations
via a three-term recurrence (Lemma~\ref{lem:recurrence}). Our main results are:
\begin{itemize}
  \item an algorithm (Algorithm~\ref{alg:main}) which computes $\Omega_n$, for all $n$, to $P$ bits of precision
  with $O(N^3 P^3)$ arithmetic operations (Theorem~\ref{thm:convergence}),
  whereas the term-by-term summation needs $2^{\Theta(P)}$ terms, together with a proof that it converges
  geometrically under hypothesis (S) of Section~\ref{sec:preliminaries}, which is closely related to the absolute
  convergence of the series (Remark~\ref{rem:S-equiv}); the constants depend on the geometry of the Schottky datum;
  \item a corresponding algorithm for the period matrix (Section~\ref{sec:periods}), together with a careful proof of the
  classical cross-ratio formula~\eqref{eq:periodmatrix}, and for the Schottky--Klein
  prime function (Section~\ref{sec:prime}), the latter being the basic tool for the explicit solution of
  boundary value problems in multiply connected planar domains~\cite{CKGN2016,Crowdy2020};
  \item an implementation in Sage, and numerical experiments (Section~\ref{sec:experiments}), including
  a comparison with the truncation by word length and with a pruned summation of Bogatyrev type.
\end{itemize}

The rest of the article is organized as follows. In Section~\ref{sec:preliminaries} we recall the necessary facts
on Schottky groups, state the convergence hypothesis, and prove the recursive structure of the series.
Section~\ref{sec:existing-algorithms} reviews the state of the art. Our algorithm and its analysis occupy
Section~\ref{sec:algorithm}. Section~\ref{sec:periods} and Section~\ref{sec:prime} treat the period matrix and the
prime function. Section~\ref{sec:experiments} reports on numerical experiments, and Section~\ref{sec:applications}
discusses applications and open directions.


\section{Preliminaries}
\label{sec:preliminaries}

\subsection{Schottky data}
Let $\Lambda$ be the closure of the set of limit points for $G$, and set
$\Omega = \Ch \setminus \Lambda$, on which $G$ acts discontinuously. The quotient
$\Omega/G$ is a Riemann surface of genus $N$.

We work with \emph{classical} Schottky groups, given by the following data.
\begin{definition}\label{def:schottky-data}
  A \emph{Schottky datum} consists of open discs $\cD_i\subset\CC$, for $i\in I=\{\pm 1,\ldots,\pm N\}$, with centers
  $a_i$ and radii $r_i$, and loxodromic $\sigma_1,\ldots,\sigma_N\in\PGL_2(\CC)$ such that
  \begin{enumerate}
    \item The closed discs $\overline{\cD_i}$ are mutually disjoint,
    \item For all $i\in I$, we have $\sigma_i(\Ch \setminus \cD_{-i}) = \overline{\cD_i}$.
  \end{enumerate}
  The set
  \[
  \cF = \Ch \setminus \bigcup_{i\in I} \overline{\cD_i}
  \]
  is a fundamental domain for $G=\langle\sigma_1,\ldots,\sigma_N\rangle$.
\end{definition}

We will assume that all the data defining this fundamental domain is available to us, and (after a change of
coordinates) that $\infty\in\cF$. The attracting fixed point $B_n$ of $\sigma_n$ lies in $\cD_n$, and the repelling
one $A_n$ lies in $\cD_{-n}$. The group $G$ is not determined by the discs alone: the transformation $\sigma_i$ is
only determined up to composition with a M\"obius transformation preserving $\overline{\cD_i}$, and which one is chosen changes the group.

\begin{remark}
  Schmies~\cite{Schmies2011} works with isometric circles (the \emph{iso-classical} case), for which the discs
  $\cD_{\pm i}$ are the interiors of the isometric circles of $\sigma_{\mp i}$. Nothing in this article requires this.
\end{remark}

For $i\in I$ we use the coordinate
\begin{equation}\label{eq:coordinate}
t_i(z) = \frac{r_i}{z-a_i},
\end{equation}
which identifies $\Ch\setminus\cD_i$ with the closed unit disc $\{|t|\leq 1\}$, with $t_i(\infty)=0$.

\subsection{Convergence}
For $\sigma\neq 1$ with Tietze word $[i_1,\ldots,i_k]$, let
\[
O(\sigma) = \sigma\bigl(\Ch\setminus\cD_{-i_k}\bigr) = \sigma_{i_1}\cdots\sigma_{i_{k-1}}\bigl(\overline{\cD_{i_k}}\bigr)
\]
(a closed disc: the ``outer disc of level $k$'' in the terminology of~\cite{Lyamaev2022}). Note that
$O(\sigma)\subset \overline{\cD_{i_1}}$. We will assume the following hypothesis, which is what makes the series of
interest converge absolutely.

\begin{enumerate}
  \item[\textbf{(S)}] There are constants $C_{\mathrm S}>0$ and $0<\theta<1$ such that for all $k\geq 1$,
  \[
  \sum_{\sigma\in\cW_k}\diam O(\sigma)\ \leq\ C_{\mathrm S}\theta^k,
  \]
  where $\cW_k=\{\sigma\in G: |\sigma|=k\}$.
\end{enumerate}

Sufficient conditions for (S) go back to Burnside~\cite{Burnside1891} (sufficiently small discs), and were made
quantitative by Schmies~\cite{Schmies2011} and, for real hyperelliptic curves, by Bogatyrev~\cite{Bogatyrev1999,Bogatyrev2012}. For instance, in the setting of
Theorem~1 in~\cite{Lyamaev2022} the sum of the diameters of the inner circles of the domain $\sigma\cF$ is at most
$\frac{\gamma-1}{\gamma+1}$ times the diameter of its outer circle, where $\gamma=\max_j\gamma_j$; these inner circles are the $O(\sigma\sigma_l)$, so this gives (S) with $\theta=\frac{\gamma-1}{\gamma+1}$.
Our results are stated and proved under (S). We now recall why (S) is essentially the absolute convergence of the
series $\omega_n$; the facts involved are standard, we only sketch the argument, and it is not used elsewhere.

\begin{remark}[sketch]\label{rem:S-equiv}
  Let $\delta$ be the Hausdorff dimension of $\Lambda$, which is also the critical exponent of the Poincar\'e series
  $\sum_{\sigma\in G}|\sigma'(z)|^s$~\cite{Sullivan1979}. By the Koebe distortion theorem, $\diam O(\sigma)$ is comparable, with
  constants that only depend on the datum and on $z$, to $|\sigma'(z)|$ for any fixed $z\in\cF$, so that
  $\sum_{\sigma\in\cW_k}\diam O(\sigma)^s$ grows like $e^{k\Pr(s)+o(k)}$, where $\Pr$ is the pressure function of the
  associated expanding map. By Bowen's formula~\cite{Bowen1979}, $\Pr$ is strictly decreasing and vanishes at $s=\delta$. Hence (S)
  holds if and only if $\delta<1$. On the other hand, for a classical Schottky group the Poincar\'e series diverges
  at $s=\delta$~\cite{Sullivan1979}, so the series $\omega_n$, whose terms are of size $\asymp|\sigma'(z)|$, converge
  absolutely if and only if $\delta<1$ (compare~\cite[\S1.8.2]{BobenkoIntro2011}). Sufficient geometric conditions are Schottky's
  ``circle decomposability'' (which holds for instance for the Schottky uniformizations of real Riemann surfaces of decomposing type
  by Fuchsian groups of the second kind, see~\cite[Theorems 31, 32]{BobenkoIntro2011}) and the criteria
  of~\cite{Burnside1891,Schmies2011,Bogatyrev2012} quoted above; there are classical Schottky groups with $\delta>1$, for which
  the series diverge (see the references in~\cite[\S1.8.2]{BobenkoIntro2011}). For a given datum, $\delta$ can be computed
  (rigorously) with the methods of~\cite{McMullen1998,PollicottVytnova2022}.
\end{remark}

\subsection{A recursive definition of \texorpdfstring{$G$}{G}}
\label{sec:recurs-defin-g}
For each $\sigma \in G$, denote by $|\sigma|$ its length, namely the number of terms in its Tietze word representation.
Recall that $\cW_k$ denotes the set of elements of length $k$.
Given $\sigma\neq 1$, define $h(\sigma) = \sigma_j$ if $|\sigma_{-j}\sigma| < |\sigma|$ (the \emph{head} of $\sigma$, its first letter). As an extra piece of notation,
write
\[
\cW_k^{(i)} = \{ \sigma \in \cW_k ~|~ h(\sigma)=\sigma_i\}
\]
for the subset of elements in $G$ of length $k$ that start with the generator $\sigma_i$.

\begin{lemma}\label{lem:tree}
  For all $k\geq 1$,
  \[
  \cW_{k+1}^{(i)} = \coprod_{j\neq -i} \sigma_i \cW_k^{(j)}.
  \]
\end{lemma}
\begin{proof}
  Every $\sigma\in\cW_{k+1}^{(i)}$ is of the form $\sigma_i\tilde\sigma$ for a unique $\tilde\sigma\in \cW_k$, and since
  the Tietze word of $\sigma$ is reduced, $h(\tilde\sigma)\neq\sigma_{-i}$. Hence $\tilde\sigma\in\cW_k^{(j)}$ for some
  $j\neq -i$. Conversely, if $j\ne -i$ and $\tilde\sigma\in \cW_k^{(j)}$, then $\sigma_i\tilde\sigma$ is reduced and
  belongs to $\cW_{k+1}^{(i)}$.
\end{proof}

For a fixed $n\in\{1,\ldots,N\}$, let $\cR_n\subset G$ be the set of elements whose Tietze word does not end in
$\pm n$ (together with the identity). It is a set of representatives of $G/H_n$, and we write
$\cR_n^{(i)}=\cR_n\cap\cW^{(i)}$ for those that start with $\sigma_i$, with $\cW^{(i)}=\bigcup_k \cW_k^{(i)}$.
The analogue of Lemma~\ref{lem:tree} is
\begin{equation}\label{eq:cosetrec}
\cR_n^{(i)} = \{\sigma_i : |i|\neq n\}\ \sqcup \coprod_{j\neq -i}\sigma_i\cR_n^{(j)}.
\end{equation}
Indeed, $\sigma_i\tilde\sigma$ with $\tilde\sigma\neq 1$ ends in $\pm n$ if and only if $\tilde\sigma$ does, while the
word $\sigma_i$ itself ends in $\pm n$ if and only if $|i|=n$.

\subsection{A recursive definition of \texorpdfstring{$\Omega_n$}{Omega\_n}}
\label{sec:recurs-defin-omegan}

Note that $\Omega_n(z)$ can be written as
\[
\Omega_n(z)= \log \prod_{\sigma \in G / H_n} \frac{z - \sigma(B_n)}{z - \sigma(A_n)}.
\]
For $\sigma\in\cR_n$, the points $\sigma(A_n)$ and $\sigma(B_n)$ lie in $O(\sigma)$ (because $A_n,B_n\in\cD_{\pm n}$ and
$\sigma$ does not end in $\pm n$). We deduce from (S) that
\[
\sum_{\sigma\in\cR_n}\bigl|\sigma(B_n)-\sigma(A_n)\bigr|\ \leq\ \sum_{k\geq 0}\sum_{\sigma\in\cW_k}\diam O(\sigma)<\infty
\]
(the term $k=0$ being the single point $\sigma=1$), so the series defining $\Omega_n$ and $\omega_n$ converge absolutely and
uniformly on compact subsets of $\Omega$; note that the points $\sigma(A_n),\sigma(B_n)$ are fixed points of loxodromic
elements of $G$, hence lie in $\Lambda$.

For $i\in I$, let
\begin{equation}\label{eq:Fi}
F_n^{(i)}(z) = \sum_{\sigma \in \cR_n^{(i)}} \log\frac{z - \sigma(B_n)}{z - \sigma(A_n)},
\end{equation}
where each summand is the branch of the logarithm which vanishes at $z=\infty$ (it is a function without zeros or poles on
the simply connected set $\Ch\setminus\cD_i$ which contains $\infty$, because both $\sigma(B_n)$ and $\sigma(A_n)$ lie in
$\cD_i$). Thus
\begin{equation}\label{eq:Omega-decomp}
\Omega_n(z) = \log\frac{z-B_n}{z-A_n} + \sum_{i\in I}F_n^{(i)}(z),\qquad z\in\cF,
\end{equation}
up to the choice of the branch of the first logarithm. Each $F_n^{(i)}$ is holomorphic on $\Ch\setminus\cD_i$ and vanishes at
$\infty$. Moreover, the singularities of~\eqref{eq:Fi} lie in the compact set
$Z_i=\sigma_i\bigl(\bigcup_{j\neq-i}\overline{\cD_j}\bigr)$, which is contained in $\cD_i$, because
$\overline{\cD_j}\cap\overline{\cD_{-i}}=\emptyset$ for $j\ne -i$. Hence
\begin{equation}\label{eq:Ri}
R_i = \min_{z\in Z_i}\,|t_i(z)|\ >\ 1,
\end{equation}
and $F_n^{(i)}$, expressed in the coordinate $t_i$, is given by a power series $\sum_{m\geq 1}c_m t_i^m$
converging on the disc $|t_i|<R_i$, which strictly contains the closed unit disc. Let $R=\min_i R_i>1$.

We now describe how the functions $F_n^{(i)}$ are related to each other. Let $\cH_i$ denote the space of functions
holomorphic on a neighborhood of $\Ch\setminus\cD_i$ and vanishing at $\infty$. For $j\neq -i$ we have
$\sigma_i^{-1}(\Ch\setminus\cD_i)=\overline{\cD_{-i}}\subset\Ch\setminus\overline{\cD_j}$, so that we may define
\[
T_{ij}:\cH_j\to\cH_i,\qquad (T_{ij}f)(z) = f\bigl(\sigma_i^{-1}z\bigr)-f\bigl(\sigma_i^{-1}\infty\bigr).
\]

\begin{theorem}[Functional equation]\label{thm:functional-equation}
  For every $i\in I$,
  \[
  F^{(i)}_n \;=\; E^{(i)}_n \;+\; \sum_{j\neq -i}T_{ij}F_n^{(j)},
  \]
  where $E_n^{(i)}(z)=\log\frac{z-\sigma_i(B_n)}{z-\sigma_i(A_n)}$ if $|i|\neq n$, and $E_n^{(i)}=0$ if $|i|=n$.
\end{theorem}
\begin{proof}
  By~\eqref{eq:cosetrec}, $F_n^{(i)}$ is the sum of $E_n^{(i)}$ and of the sum over $j\neq -i$ and $\tilde\sigma\in\cR_n^{(j)}$ of
  $\log\frac{z-\sigma_i\tilde\sigma(B_n)}{z-\sigma_i\tilde\sigma(A_n)}$. If $\sigma_i=\smtx{a}{b}{c}{d}$, then for
  all $z,w$ one has $\sigma_i(z)-\sigma_i(w)=\frac{(ad-bc)(z-w)}{(cz+d)(cw+d)}$. Applying this with $z$ replaced by $u=\sigma_i^{-1}(z)$ gives
  \[
  \frac{z-\sigma_i(w_1)}{z-\sigma_i(w_2)} = \frac{u-w_1}{u-w_2}\cdot\frac{cw_2+d}{cw_1+d}.
  \]
  Taking $w_1=\tilde\sigma(B_n)$ and $w_2=\tilde\sigma(A_n)$ we find that the $\tilde\sigma$-summand equals
  the corresponding summand of $F_n^{(j)}$ evaluated at $\sigma_i^{-1}z$, plus a constant. Since both sides vanish at $z=\infty$, the
  constant is minus the value of the first term at $z=\infty$, which is exactly the normalization in $T_{ij}$.
  Summing over $\tilde\sigma$ (the sums converge absolutely) gives the result.
\end{proof}

Write $\cH=\bigoplus_{i\in I}\cH_i$, $F_n=(F_n^{(i)})_i$, $E_n=(E_n^{(i)})_i$ and let $\cT:\cH\to\cH$ be the operator
with components $T_{ij}$ ($j\neq -i$) and $0$ otherwise. Then $F_n=E_n+\cT F_n$. Unwinding using Lemma~\ref{lem:tree}, we
see that
\[
(\cT^kE_n)^{(i)}(z)=\sum_{\sigma\in\cR_n^{(i)},\ |\sigma|=k+1}\log\frac{z-\sigma(B_n)}{z-\sigma(A_n)},
\]
that is, $\cT^kE_n$ is precisely the contribution of the level $k+1$ of the Cayley tree, and
\begin{equation}\label{eq:neumann}
F_n=\sum_{k\geq 0}\cT^kE_n=(1-\cT)^{-1}E_n,
\end{equation}
where the series converges and $1-\cT$ is invertible on a suitable space of functions (see the proof of Theorem~\ref{thm:convergence}).
The point of this formula is that each level is not computed term by term: a level is a vector of $2N$ functions, and
passing to the next level is one application of the fixed operator $\cT$.


\section{Existing algorithms}
\label{sec:existing-algorithms}

In this section we review the existing algorithms to compute theta functions of Schottky groups. All of them evaluate the
series in the form
\[
\sum_{\sigma\in\cR_n} \bigl(\text{explicit term attached to }\sigma\bigr)
\]
over a finite subtree $\cS_0\subset\cR_n$ of the Cayley tree, and differ in the way $\cS_0$ is chosen and in the error
estimate for the discarded part. Note that if we normalize so that the evaluation point $z$ is in $\cF$, the term attached
to $\sigma$ is of size $\asymp |\sigma(B_n)-\sigma(A_n)|\leq\diam O(\sigma)$.

\subsubsection*{Truncation by word length}
The simplest method~\cite{FairchildRiosOrtiz2024} (see also~\cite{Schmies2011}) keeps all words of length at most $L$, a set
of $1+\sum_{k=1}^L 2N(2N-1)^{k-1}$ elements. By hypothesis (S) the neglected part is $O(\theta^L)$, and
explicit values of $\theta$ in terms of cross-ratios of the circles are given by Bogatyrev's lemma
(Lemma~6.1 in~\cite{FairchildRiosOrtiz2024}, based on~\cite[Lemma 6.1]{Bogatyrev2012}). The cost is exponential in the
precision: to achieve error $2^{-P}$ one needs $L\approx P\log 2/\log(1/\theta)$ and hence about $(2N-1)^{L}$ terms.

\subsubsection*{Bogatyrev's and Burnside's estimates, and adaptive summation}
Bogatyrev~\cite{Bogatyrev1999} observed that the rates at which the terms decay along different branches of the tree differ
widely, and proposed to decide separately for each branch whether to stop: one descends from a vertex $\sigma$ as long as an
estimate $K(\sigma)|\sigma(B_n)-\sigma(A_n)|$ for the total contribution of the descendants of $\sigma$ is larger than a tolerance $\nu$
(depth-first traversal, a posteriori error estimate). Schmies~\cite{Schmies2005,Schmies2011} uses the sharper
Burnside-type estimates~\cite{Burnside1891}, valid under a criterion on the Schottky data, and orders the vertices by a
weight, always expanding the vertex of largest weight (a priority queue), which yields an a priori error estimate. Both are
described in detail in~\cite[\S2]{Lyamaev2022}. Lyamaev~\cite{Lyamaev2022,Lyamaev2022b} proves a new
estimate for the sum over a subtree in terms of the term at its root, valid in situations where Burnside's estimate
is not applicable, and uses it to reduce the number of terms needed ``by several times in the case of slow
convergence and by tens of percent in the usual situations'' (for real hyperelliptic curves, in the setting of the
\emph{Schottky model} used by Bogatyrev for extremal polynomials).

Yamamoto and Miyoshi~\cite{YamamotoMiyoshi2026} apply the same philosophy to the infinite product defining the
Schottky--Klein prime function: each group element is assigned a cross-ratio ``potential'' bounding its contribution,
uniform bounds on how the potential changes when a generator is prepended allow an efficient enumeration of the
elements above a threshold, and the relative error decays exponentially in the threshold at a rate governed by the
Hausdorff dimension of the limit set. Burnside's method is also the one used to compute period matrices in the numerical
Schottky uniformizations of~\cite{Hidalgo2001,HidalgoFigueroa2005}.

These methods are very effective at moderate accuracy, for which they were designed.
However, their cost grows exponentially with the number of digits, with an exponent which we now determine. Let
$\delta\in(0,1)$ be as in Remark~\ref{rem:S-equiv} and let $m_\Lambda$ be the Patterson--Sullivan measure of $G$, a probability
measure on $\Lambda$ which is conformal of dimension $\delta$~\cite{Sullivan1979}, so that $m_\Lambda(O(\sigma))\asymp\diam O(\sigma)^\delta$ by bounded
distortion. A pruned summation stops at a \emph{cut} $\cS$ of the tree, that is, a set of vertices whose discs $O(\sigma)$ cover $\Lambda$ without overlaps;
the part of the series below a vertex $\sigma$ is $\asymp\diam O(\sigma)$, so the discarded part, bounded by the sum of the absolute values of its terms,
is $\asymp\sum_{\sigma\in\cS}\diam O(\sigma)$, while $\sum_{\sigma\in\cS}\diam O(\sigma)^\delta\asymp m_\Lambda(\Lambda)=1$. By H\"older's inequality,
$\sum_{\cS}\diam O(\sigma)^\delta\leq|\cS|^{1-\delta}\bigl(\sum_{\cS}\diam O(\sigma)\bigr)^\delta$, so that a cut whose discarded part is at most $\varepsilon$
has at least a constant times $\varepsilon^{-\delta/(1-\delta)}$ elements, and this is attained, up to constants, by cutting at the level set $\diam O(\sigma)\approx\varepsilon^{1/(1-\delta)}$.
Hence any method which sums the terms individually and bounds the discarded part by the sum of the absolute values of its terms
needs about $2^{\delta P/(1-\delta)}$ terms to reach the accuracy $2^{-P}$: pruning improves the constant in the exponent, not the exponential growth,
and the exponent $\delta/(1-\delta)$ tends to infinity as $\delta\to1$, that is, as the data approach the boundary of the range of convergence (compare the
dimension-dependent rate in~\cite{YamamotoMiyoshi2026}). In our experiments (Section~\ref{sec:experiments}) the number of terms visited by a pruned
depth-first summation grows like $\nu^{-0.38}$ and its error decays like $\nu^{0.62}$ in the tolerance $\nu$, consistent with $\delta\approx0.38$ for that
datum: the number of terms is multiplied by about $10^{\delta/(1-\delta)}\approx4$ for each additional digit.

\subsubsection*{Methods which avoid the series}
For the Schottky--Klein prime function of a circular domain in the plane, Crowdy and Marshall~\cite{CrowdyMarshall2007}
proposed an algorithm which does not use any sum or product over the group: the prime function (more precisely its
square $X=\omega^2$) and some auxiliary functions are expanded in Fourier--Laurent series about the centres of the
circles, and the coefficients are found by imposing the defining functional relations at collocation points on the
circles and solving the resulting overdetermined linear system by least squares (building on a code of Trefethen).
In their example, $15$ Laurent coefficients give $15$ correct digits, while the product truncated at level $9$ gives
about $12$; no convergence proof is given (as noted in~\cite{CrowdyMarshallMcurves2011}).
Their cost is polynomial in the number of circles. Write $N$ for the number of pairs of circles and $N_1$ for the truncation order.
Step~1 of~\cite[\S4]{CrowdyMarshall2007} has $4NN_1+2N$ real unknowns and at least $2N(2N_1+1)$ real equations, the same
matrix serving for the $N$ right-hand sides, and Step~2 has $2NN_1$ complex unknowns; a dense least-squares solution therefore
costs $O(N^3N_1^3)$ operations, and each evaluation of the prime function afterwards costs $O(NN_1)$. (Adding a circle adds only $O(N_1)$ unknowns,
as remarked in~\cite{CrowdyMarshall2007}.) This is the same exponent as ours, with $N_1$ playing the r\^ole of the number $p$ of Taylor coefficients.
What is not known is how large $N_1$ has to be for $P$ bits: the Laurent coefficients should decay geometrically, so that
$N_1=O(P)$ for a fixed domain, as in the observed $N_1=15$ for $15$ digits, but neither this nor the conditioning of the least-squares problems
is proved, and the constants depend on the geometry. We implemented the method of~\cite{CrowdyMarshall2007} to measure it (Section~\ref{sec:collocation}): the observed cost is cubic in $N_1$ and in $N$, and about $0.7$ digit is gained per coefficient. We have not analysed the cost of the
boundary value formulation of~\cite{CKGN2016}, whose implementation can use fast multipole acceleration. This was superseded by the
boundary value formulation of Crowdy, Kropf, Green and Nasser~\cite{CKGN2016} (see also~\cite{Crowdy2020}), with a
public \textsc{Matlab} implementation (\texttt{SKPrime}). These methods have the important advantage of not needing
the convergence of the series, but they are specific to Schottky doubles of planar domains. Our algorithm is also based on
expansions at the circles, but the linear system relating the coefficients is not obtained by collocation: it is
derived \emph{exactly} from the group action (Theorem~\ref{thm:functional-equation}), applies to any classical
Schottky group satisfying (S), and the truncation order and the number of iterations are the only sources of error,
which we analyze in Theorem~\ref{thm:convergence}.

\subsubsection*{Functional equations}
The idea of replacing a Poincar\'e series by a system of functional equations for functions analytic in discs is not new.
In the \emph{method of functional equations} of Mityushev and Rogosin~\cite{MityushevRogosin2000}, boundary value problems in
circular multiply connected domains are reduced to systems of functional equations, in which the unknown functions are analytic in
the discs and are related through compositions with the inversions in the boundary circles; these systems are solved by successive
approximations, which amounts to summing Poincar\'e series level by level, and the method has been used to study the convergence of
Poincar\'e series~\cite{Mityushev1998} and to compute them efficiently, also for nearly touching discs~\cite{MityushevRylko2013}.
Our equation $F_n=E_n+\cT F_n$ is of the same nature. What we add, as far as we know, is its formulation for the Abelian integrals of
an arbitrary classical Schottky group (not necessarily generated by inversions in the circles of a planar domain), an explicit
discretization of the transfer operator by Taylor expansions computed with the recurrence of Lemma~\ref{lem:recurrence}, and a
complexity analysis with a proof of geometric convergence of the discretized iteration (Theorem~\ref{thm:convergence}). We have not
compared our implementation with theirs.

\subsubsection*{Other approaches}
Starting from an algebraic model of the curve one may compute periods by numerical integration, which is the
standard route to Riemann matrices and the setting of the efficient evaluation of Riemann theta
functions~\cite{DHBHS2004}; in our setting the algebraic model is exactly what is unknown. See
also~\cite{FairchildRiosOrtiz2024} for the computation of an algebraic curve from a Schottky group, where the
truncation of the Poincar\'e theta series limits the attainable precision, and which was one of the motivations of the present work.

\medskip
In Table~\ref{tab:comparison} we summarize the comparison.
\begin{table}[t]
\centering
\small
\renewcommand{\arraystretch}{1.25}
\begin{tabularx}{\textwidth}{@{}>{\raggedright\arraybackslash}p{4.6cm}>{\raggedright\arraybackslash}X>{\raggedright\arraybackslash}p{3.3cm}@{}}
\toprule
Method & Cost for $P$ bits & Error control \\
\midrule
Truncation by word length~\cite{FairchildRiosOrtiz2024} & $(2N-1)^{\Theta(P)}$ terms & a priori (Bogatyrev) \\
Pruned DFS / priority queue~\cite{Bogatyrev1999,Schmies2011,Lyamaev2022,YamamotoMiyoshi2026} & about $2^{\delta P/(1-\delta)}$ terms & a priori or a posteriori \\
Collocation~\cite{CrowdyMarshall2007} (prime function, planar domains) & $O(N^3N_1^3)$ operations (dense least squares), $N_1$ coefficients per circle, with $N_1=O(P)$ expected but not proved (Section~\ref{sec:collocation}) & a posteriori \\
Boundary value problems~\cite{CKGN2016} (prime function, planar domains) & not analysed here; fast in practice & a posteriori \\
\midrule
\textbf{This article} & $O(N^3P^3)$ operations, constants depending on the datum & geometric rate, non-explicit constants (Thm.~\ref{thm:convergence}) \\
\bottomrule
\end{tabularx}
\caption{Comparison of the cost of the algorithms to compute $\Omega_n(z)$ with $2^{-P}$ accuracy.}
\label{tab:comparison}
\end{table}


\section{The algorithm}
\label{sec:algorithm}

\subsection{Transfer matrices}
For $j\neq -i$, let $\varphi_{ij}=t_j\circ\sigma_i^{-1}\circ t_i^{-1}$, a M\"obius transformation. It maps the closed unit
disc into the open unit disc (because $\sigma_i^{-1}t_i^{-1}$ maps it onto $\overline{\cD_{-i}}$, which is disjoint
from $\overline{\cD_j}$). In the coordinates $t_i$ and $t_j$, the operator $T_{ij}$ reads
\[
T_{ij}\Bigl(\sum_{k\geq 1}f_kt^k\Bigr)=\sum_{k\geq 1}f_k\bigl(\varphi_{ij}(u)^k-\varphi_{ij}(0)^k\bigr).
\]
Let $M^{(ij)}\in\CC^{p\times p}$ be the matrix whose entry $M^{(ij)}_{km}$ is the coefficient of $u^m$ in $\varphi_{ij}(u)^k$,
for $0\leq k,m<p$. Then $T_{ij}$ acts, on the coefficient vector $f=(f_1,\ldots,f_{p-1})$ of a truncated series, by
$f\mapsto \bigl(\sum_{k=1}^{p-1}f_kM^{(ij)}_{km}\bigr)_{m=1,\ldots,p-1}$.

\begin{lemma}\label{lem:recurrence}
  Let $\varphi(u)=\frac{au+b}{cu+d}$ with $bd\neq0$, and for $k\geq0$ write $\varphi(u)^k=\sum_{m\geq0}M_{km}u^m$.
  Then $M_{k0}=(b/d)^k$ and, for $m\geq 0$,
  \[
  (m+1)\,bd\,M_{k,m+1}=\bigl(k(ad-bc)-(ad+bc)m\bigr)M_{km}-ac\,(m-1)\,M_{k,m-1}.
  \]
  In particular, all the entries of $M$ can be computed with $O(p^2)$ arithmetic operations.
\end{lemma}
\begin{proof}
  Put $f=\varphi^k$. Then $f'/f=k\bigl(\tfrac{a}{au+b}-\tfrac{c}{cu+d}\bigr)=\tfrac{k(ad-bc)}{(au+b)(cu+d)}$, i.e.
  $(acu^2+(ad+bc)u+bd)f'=k(ad-bc)f$. Comparing the coefficients of $u^m$ gives the recurrence.
\end{proof}

For $\varphi=\varphi_{ij}$ the hypothesis $bd\neq0$ holds automatically: $b/d=\varphi(0)=t_j(\sigma_i^{-1}\infty)$ is finite and nonzero, because
$\sigma_i^{-1}\infty\in\cD_{-i}$ is a finite point different from the centre $a_j$. (We checked that the recurrence agrees with the
straightforward computation of powers of power series to the digits displayed below.)

Finally, the seed $E_n^{(i)}$ is expanded in the coordinate $t_i$ as follows: for a point $w\in\cD_i$, we have
$z-w=\frac{r_i}{t}\bigl(1+\frac{a_i-w}{r_i}t\bigr)$ with $t=t_i(z)$, so
\begin{equation}\label{eq:seed}
\log\frac{z-w}{z-w'}=\log\Bigl(1+\tfrac{a_i-w}{r_i}t\Bigr)-\log\Bigl(1+\tfrac{a_i-w'}{r_i}t\Bigr)
=\sum_{m\geq1}\frac{(-1)^{m+1}}{m\,r_i^m}\bigl((a_i-w)^m-(a_i-w')^m\bigr)t^m.
\end{equation}

\subsection{The iteration}

\begin{algorithm}[h]
\caption{Computation of $\Omega_n$}\label{alg:main}
\KwIn{A Schottky datum (Definition~\ref{def:schottky-data}), a truncation order $p$, a number of iterations $K$.}
\KwOut{Coefficient vectors $f^{(i)}\in\CC^{p-1}$, $i\in I$, approximating $F_n^{(i)}$ in the coordinate $t_i$.}
\textbf{Precompute:} for every pair $(i,j)$ with $j\neq-i$, the matrix $M^{(ij)}$ (Lemma~\ref{lem:recurrence}) and its
truncation $\widehat M^{(ij)}=\bigl(M^{(ij)}_{km}\bigr)_{1\leq k,m<p}$\;
\textbf{Seeds:} for $|i|\neq n$ compute $e^{(i)}$ from~\eqref{eq:seed} with $w=\sigma_i(B_n)$, $w'=\sigma_i(A_n)$, and $e^{(i)}=0$ for $|i|=n$\;
$f^{(i)}\gets e^{(i)}$ for all $i$\;
\For{$s=1,\ldots,K$}{
  \For{$i\in I$}{
    $g^{(i)}\gets e^{(i)}+\sum_{j\neq-i}f^{(j)}\,\widehat M^{(ij)}$\;
  }
  $f\gets g$\;
}
\Return{$f$}
\end{algorithm}

The transfer matrices do not depend on $n$; only the seeds do. Hence for all the functions
$\Omega_1,\ldots,\Omega_N$ the precomputation is shared. Given the output, $\Omega_n(z)$ is evaluated for $z\in\overline\cF$ as
\begin{equation}\label{eq:eval}
\Omega_n(z)\approx\log\frac{z-B_n}{z-A_n}+\sum_{i\in I}\sum_{m=1}^{p-1}f^{(i)}_m\,t_i(z)^m,
\end{equation}
and $\omega_n(z)=\Omega_n'(z)$ by differentiating termwise ($\frac{d}{dz}t_i=-t_i^2/r_i$). For $z\notin\overline\cF$ one first moves
$z$ into $\overline\cF$ using the group action, via $\Omega_n(\sigma_m z)=\Omega_n(z)+B_{mn}$
(see Section~\ref{sec:periods}); this is cheap whenever the word length of the element needed is small.

\begin{remark}
  In an implementation, we do not need to exponentiate: the operators $T_{ij}$ are linear, and act directly on the
  coefficients of the logarithms $F_n^{(j)}$. (Our first implementation, see Section~\ref{sec:experiments}, performed the
  iteration on the power series of $\exp F_n^{(j)}$, which requires an exponential and a logarithm of a power series per step.)
\end{remark}

\begin{remark}\label{rem:direct}
  Since the algorithm just computes a Neumann series~\eqref{eq:neumann} for a matrix $\widehat M$ of size $2N(p-1)$, one may
  instead solve the linear system $(1-\widehat M)f=e$ directly, for all $n$ at once (with $N$ right-hand sides). This
  has the same asymptotic cost $O(N^3p^3)$ and avoids having to choose $K$. Both variants are implemented: the
  cost comparison below uses the direct solver, and the experiments on the genus the iteration.
\end{remark}

\subsection{Convergence and complexity}

\begin{theorem}\label{thm:convergence}
  Assume (S), and let $\lambda\in(\theta,1)$. There are constants $C>0$, $p_0\geq1$ and $0<\kappa<1$, where $\kappa$
  only depends on the Schottky datum, such that for all $p\ge p_0$ and all $K\geq0$ the output of
  Algorithm~\ref{alg:main} satisfies
  \[
  \sup_{z\in\overline\cF}\Bigl|\Omega_n(z)-\Omega_n^{(p,K)}(z)\Bigr|\ \leq\ C\bigl(\lambda^K+\kappa^{p}\bigr)
  \]
  (modulo $2\pi i$), where $\Omega_n^{(p,K)}$ is the right hand side of~\eqref{eq:eval}.
  In particular, if $p,K\geq c\,P$ for a constant $c$ which only depends on the datum, the error is $O(2^{-P})$, and the
  number of arithmetic operations is $O(N^2p^2K)$ for each $n$, that is, $O(N^3P^3)$ for all $n$ together.
\end{theorem}

Here the operations are counted in exact arithmetic, and the constants depend on the datum (see Remark~\ref{rem:transfer}). We have
not analysed the propagation of rounding errors, and we make no claim on the bit complexity. In the experiments of
Section~\ref{sec:experiments} we work with $40$ to $64$ bits more than the target accuracy, and the results do not change when the precision is
increased. For the matrices of Lemma~\ref{lem:recurrence}, we compared the recurrence with the computation by repeated multiplication of
power series: for the data of Section~\ref{sec:experiments} and of Example~\ref{ex:fro}, with $p=120$ at $200$ bits, the relative difference
is below $2\cdot10^{-62}$, that is, below the working precision. In all our examples $ac=0$ in Lemma~\ref{lem:recurrence}, so that the recurrence has only two
terms; we have not tested the general three-term case.

\begin{proof}
  \emph{Function spaces.} For $\rho\geq1$ and $i\in I$ let $U_i^\rho=\{z\in\Ch:|t_i(z)|<\rho\}$, the complement of
  the closed disc of centre $a_i$ and radius $r_i/\rho$; it is an open disc in $\Ch$ containing $\infty$, and
  $U_i^1=\Ch\setminus\overline{\cD_i}$. Let $Y_\rho$ be the Banach space of tuples $f=(f^{(i)})_{i\in I}$ of bounded holomorphic
  functions $f^{(i)}$ on $U_i^\rho$ vanishing at $\infty$, with the norm $\|f\|_\rho=\max_i\sup_{U_i^\rho}|f^{(i)}|$. We identify
  $f^{(i)}$ with its Taylor series in $t_i$, and denote by $\Pi_p$ the truncation of each component to degree $<p$.

  \emph{Choice of radii.} For $j\neq-i$ we have $\sigma_i^{-1}(\Ch\setminus\cD_i)=\overline{\cD_{-i}}\subset U_j^1$. Moreover,
  the compact set $Q_i=\sigma_i\bigl(\{\infty\}\cup\bigcup_{j\neq -i}\overline{\cD_j}\bigr)$ contains $Z_i$ and is
  contained in $\cD_i$, so $R'=\min_i\min_{Q_i}|t_i|>1$ and $R'\leq R$. By compactness we may choose $1<\rho<\rho'<R'$ such that
  \begin{equation}\label{eq:radii}
    W_{ij}:=\sigma_i^{-1}\bigl(\overline{U_i^{\rho'}}\bigr)\subset U_j^1\quad\text{for all } i\in I,\ j\neq-i.
  \end{equation}
  We set $\kappa=\rho/\rho'$, and write $r=\min_i r_i$.

  \emph{Step 1: $\cT$ is bounded from $Y_\rho$ to $Y_{\rho'}$.} If $z\in U_i^{\rho'}$, then $\sigma_i^{-1}(z)$ and $\sigma_i^{-1}(\infty)$
  lie in $W_{ij}\subset U_j^1\subset U_j^\rho$, so $(T_{ij}f)(z)$ is defined, holomorphic in $z\in U_i^{\rho'}$, vanishes at $\infty$ and is bounded by
  $2\|f\|_\rho$. Hence $\|\cT f\|_{\rho'}\leq2(2N-1)\|f\|_\rho$. By Cauchy's estimates, the Taylor coefficients of $g\in Y_{\rho'}$ satisfy
  $|g_m|\leq\|g\|_{\rho'}\rho'^{-m}$, so that $\|(1-\Pi_p)g\|_\rho\leq\|g\|_{\rho'}\kappa^p/(1-\kappa)$. Consequently
  \begin{equation}\label{eq:trunc}
    \|(1-\Pi_p)\cT\|_{\rho\to\rho}\leq C_1\kappa^p,\qquad C_1=\frac{2(2N-1)}{1-\kappa}.
  \end{equation}

  \emph{Step 2: $\|\cT^k\|_{\rho\to\rho}\leq C_2\theta^k$.} By induction on $k$,
  \[
  (\cT^kf)^{(i)}(z)=\sum f^{(j_k)}\bigl(\sigma^{-1}z\bigr)-f^{(j_k)}\bigl(\sigma^{-1}\infty\bigr),
  \qquad \sigma=\sigma_i\sigma_{j_1}\cdots\sigma_{j_{k-1}},
  \]
  where the sum runs over the sequences $(j_1,\ldots,j_k)$ with $j_1\neq-i$ and $j_{l+1}\neq-j_l$; so $\sigma$ runs over the
  words of length $k$ with head $\sigma_i$, each one appearing $2N-1$ times. Applying~\eqref{eq:radii} letter by letter, and using
  $U^1_j\subset U^{\rho'}_j$, we see that $V_\sigma=\sigma^{-1}(U_i^\rho)$ is contained in
  $W_{j_{k-1}j_k}$ (with $j_0=i$), where $|t_{j_k}|<1$. There $\bigl|\tfrac{d}{dz}f^{(j_k)}\bigr|\leq\|f\|_\rho/((\rho-1)r)$ by Cauchy's
  estimate, since $\bigl|\tfrac{dt}{dz}\bigr|=|t|^2/r_{j_k}$. As $V_\sigma$ is a bounded round disc containing $\sigma^{-1}z$ and $\sigma^{-1}\infty$,
  each summand is bounded by $\diam(V_\sigma)\,\|f\|_\rho/((\rho-1)r)$.

  It remains to compare $\diam V_\sigma$ with $\diam O(\sigma^{-1})$, where $O(\sigma^{-1})=\sigma^{-1}(\Ch\setminus\cD_i)$ because the last
  letter of $\sigma^{-1}$ is $\sigma_{-i}$. The pole of $\sigma^{-1}$ is $q=\sigma(\infty)$, and $q\in Q_i$ (if $\sigma=\sigma_i\sigma''$, then $\sigma''(\infty)$ is
  $\infty$ or lies in $O(\sigma'')\subset\overline{\cD_{j_1}}$), so $x=|q-a_i|\leq r_i/R'$. For a M\"obius transformation
  $\smtx{a'}{b'}{c'}{d'}$ with pole $q$ inside the disc $|z-a_i|<s$, the image of the complement of that disc is a disc of diameter
  $2s|a'd'-b'c'|/\bigl(|c'|^2(s^2-x^2)\bigr)$. Applying this with $s=r_i/\rho$ and $s=r_i$,
  \[
  \frac{\diam V_\sigma}{\diam O(\sigma^{-1})}=\frac{r_i^2-x^2}{\rho\,(r_i^2/\rho^2-x^2)}\leq\frac{1}{\rho\,(\rho^{-2}-R'^{-2})}=:C_3.
  \]
  Since $\sigma\mapsto\sigma^{-1}$ is a bijection of $\cW_k$, hypothesis (S) gives
  $\|\cT^kf\|_\rho\leq C_2\theta^k\|f\|_\rho$ with $C_2=(2N-1)C_3C_{\mathrm S}/((\rho-1)r)$.

  \emph{Step 3: an adapted norm.} Fix $\lambda_0\in(\theta,\lambda)$. By Step~2, $c_*=\sum_{k\geq0}\lambda_0^{-k}\|\cT^k\|_{\rho\to\rho}<\infty$, and
  $\|f\|_*=\sum_{k\geq0}\lambda_0^{-k}\|\cT^kf\|_\rho$ is a norm on $Y_\rho$ with $\|f\|_\rho\leq\|f\|_*\leq c_*\|f\|_\rho$ and
  $\|\cT f\|_*\leq\lambda_0\|f\|_*$. Let $\mathcal{A}_p=\Pi_p\cT=\cT-(1-\Pi_p)\cT$. By~\eqref{eq:trunc},
  $\|\mathcal{A}_p\|_*\leq\lambda_0+c_*C_1\kappa^p\leq\lambda$ for $p\geq p_0$.

  \emph{Step 4: the algorithm.} By the discussion following Lemma~\ref{lem:recurrence}, the vectors computed by
  Algorithm~\ref{alg:main} are the coefficients of $f_0=\Pi_pE_n$ and $f_{s+1}=\Pi_pE_n+\mathcal{A}_pf_s$: indeed, for a polynomial $f$
  the truncation $\widehat M^{(ij)}$ computes exactly the coefficients of degrees $1,\ldots,p-1$ of $T_{ij}f$. Since $E_n,F_n\in Y_{\rho'}$ (their
  singularities lie in the sets $Z_i$, and $\rho'<R$), the fixed point $F_n^{(p)}=(1-\mathcal{A}_p)^{-1}\Pi_pE_n$ exists, and
  $\|f_K-F_n^{(p)}\|_*\leq\lambda^K\|f_0-F_n^{(p)}\|_*$, where $\|f_0-F_n^{(p)}\|_*\leq\frac{2c_*}{1-\lambda}\|\Pi_pE_n\|_\rho\leq\frac{2c_*\|E_n\|_{\rho'}}{(1-\lambda)(1-\kappa)}$ is bounded independently of $p$ (by $\|(1-\mathcal{A}_p)^{-1}\|_*\leq(1-\lambda)^{-1}$, and by Cauchy's estimates, as in the proof of~\eqref{eq:trunc}). Finally, $F_n=E_n+\cT F_n$
  gives $(1-\mathcal{A}_p)F_n=\Pi_pE_n+(1-\Pi_p)F_n$, whence
  \[
  \|F_n-F_n^{(p)}\|_*=\bigl\|(1-\mathcal{A}_p)^{-1}(1-\Pi_p)F_n\bigr\|_*\leq\frac{c_*}{1-\lambda}\,\frac{\|F_n\|_{\rho'}\kappa^p}{1-\kappa}.
  \]
  For $z\in\overline\cF$ and every $i$ we have $z\in\overline{U_i^1}\subset U_i^\rho$, so
  $|\Omega_n(z)-\Omega_n^{(p,K)}(z)|\leq 2N\|F_n-f_K\|_\rho\leq 2N\|F_n-f_K\|_*$, and the bound follows.

  \emph{Cost.} The matrices $\widehat M^{(ij)}$ are computed with $O(p^2)$ operations each (Lemma~\ref{lem:recurrence}),
  and there are $2N(2N-1)$ of them; each iteration performs the same number of products of a matrix of size $p-1$ by a vector,
  i.e., $O(N^2p^2)$ operations; and one needs $K=O(P)$ iterations for each of the $N$ values of $n$. The seeds
  and the final evaluation need $O(Np)$ operations each.
\end{proof}

\begin{remark}\label{rem:transfer}
  The constant $\lambda$ can be taken arbitrarily close to $\theta$, while $\kappa$ is determined by how much the discs can be
  enlarged in~\eqref{eq:radii}, that is, by the gaps between the discs; $C$ and $p_0$ also depend on $N$ (through the factor $2N-1$ in $C_1$
  and $C_2$). To reach $2^{-P}$ one needs $p\approx P\log2/\log(1/\kappa)$ and $K\approx P\log2/\log(1/\lambda)$, so the cost is
  $O\bigl(N^3P^3/(\log(1/\kappa)^2\log(1/\lambda))\bigr)$ operations: polynomial for data whose rates stay bounded away from $1$, but not uniformly in
  families in which the discs approach each other or $\delta\to1$ (compare Figure~\ref{fig:geometry}). The operator $\cT$ is a transfer operator, in the sense of
  thermodynamic formalism, for the expanding map defined by the generators on the discs, acting on spaces of holomorphic functions;
  such operators are compact (indeed nuclear) and have exponentially decaying eigenvalues, see~\cite{BandtlowJenkinson2008}
  for explicit bounds. Similar operators, acting on $|\sigma'|^s$ instead of on our logarithms, are the basis of the algorithms
  for the Hausdorff dimension of limit sets~\cite{McMullen1998,PollicottVytnova2022}, and Algorithm~\ref{alg:main} is a finite
  section (Galerkin) approximation of the equation $(1-\cT)F_n=E_n$. The techniques of~\cite{BandtlowJenkinson2008} could be used to
  make the constants of Theorem~\ref{thm:convergence} explicit. In practice we choose $p$ and $K$ by the
  rule ``increase $p$ until two consecutive values of $p$ give the same result to the desired digits'', with $K=p$,
  which is reliable in our experiments (Section~\ref{sec:experiments}).
\end{remark}


\section{The period matrix}
\label{sec:periods}

The functions $\Omega_n$ are not single valued on $\Omega$, but their differentials are $G$-invariant. Therefore, for $\sigma\in G$ the
difference $\Omega_n(\sigma z)-\Omega_n(z)$ is a constant, which, for $\sigma=\sigma_m$ is the \emph{$B$-period}
\[
B_{mn}=\int_z^{\sigma_m(z)}\omega_n=\Omega_n(\sigma_mz)-\Omega_n(z),
\]
defined modulo $2\pi i$. The $A$-periods are normalized: $\oint_{\partial\cD_m}\omega_n=2\pi i\,\delta_{mn}$ for the positively oriented circle,
since the series converges uniformly on $\partial\cD_m$ and the term attached to $\sigma$ contributes $2\pi i$ times the number of the points
$\sigma(B_n),\sigma(A_n)$ in $\cD_m$, counted with signs $+1$ and $-1$: for $\sigma\neq1$ both points lie in $O(\sigma)$, so they are
both inside or both outside $\cD_m$, while for $\sigma=1$ we have $B_n\in\cD_n$ and $A_n\in\cD_{-n}$.

\begin{algorithm}[h]
\caption{Period matrix}\label{alg:periods}
\KwIn{A Schottky datum, $p$, $K$.}
\KwOut{The matrix $(B_{mn})$, up to $2\pi i\ZZ$.}
Run Algorithm~\ref{alg:main} for $n=1,\ldots,N$, sharing the precomputation\;
\For{$m=1,\ldots,N$}{
  Choose $z\in\partial\cD_{-m}$, and let $w=\sigma_m(z)\in\partial\cD_m$ (any $z$ will do)\;
  \For{$n=1,\ldots,N$}{$B_{mn}\gets\Omega_n(w)-\Omega_n(z)$, using~\eqref{eq:eval}\;}
}
\end{algorithm}

Since $z$ and $w$ lie on the circles where $|t_i|=1<R$, the evaluation by~\eqref{eq:eval} converges geometrically, with the same rate as in
Theorem~\ref{thm:convergence}. The total cost is $O(N^3p^2K)+O(N^3p)$, i.e., $O(N^3P^3)$.

\begin{proposition}\label{prop:periods}
  For all $m,n$, the period $B_{mn}=\Omega_n(\sigma_mz)-\Omega_n(z)$ is given, modulo $2\pi i$, by~\eqref{eq:periodmatrix}.
\end{proposition}
\begin{proof}
  Fix $z\in\cF$, and write $\sigma_m=\smtx{a}{b}{c}{d}$; we have $c\neq0$, since $\infty\in\cF$ is not a fixed point of $\sigma_m$. The
  difference $\Omega_n(\sigma_mz)-\Omega_n(z)=\sum_{\sigma\in G/H_n}\bigl(g_\sigma(\sigma_mz)-g_\sigma(z)\bigr)$, with
  $g_\sigma(x)=\log\frac{x-\sigma(B_n)}{x-\sigma(A_n)}$, is a difference of two absolutely convergent series, so we may
  group its terms along the orbits of $H_m$ acting on $G/H_n$ by left multiplication.

  If $m=n$, the coset $H_n$ is fixed, and its contribution is $\log\mu_n$ by~\eqref{eq:multiplier}. Every other orbit is of the form
  $\{\sigma_m^k\gamma H_n\}_{k\in\ZZ}$ with all its elements distinct: if $\sigma_m^k\gamma H_n=\gamma H_n$ with $k\neq0$, then
  $\gamma^{-1}\sigma_m^k\gamma\in H_n$, which in the free group $G$ forces $m=n$ and $\gamma\in H_n$. These orbits correspond bijectively to the
  double cosets $H_m\gamma H_n$ appearing in~\eqref{eq:periodmatrix}. Fix such a $\gamma$, and put $b_k=\sigma_m^k\gamma(B_n)$, $a_k=\sigma_m^k\gamma(A_n)$ and
  $g_k=g_{\sigma_m^k\gamma}$. The identity $\sigma_m(x)-\sigma_m(y)=\frac{(ad-bc)(x-y)}{(cx+d)(cy+d)}$ gives
  $g_k(\sigma_mz)=g_{k-1}(z)+\ell_{k-1}$ with $\ell_k=\log\frac{ca_k+d}{cb_k+d}$, so that
  \[
  \sum_{k=-L}^{L}\bigl(g_k(\sigma_mz)-g_k(z)\bigr)=g_{-L-1}(z)-g_L(z)+\sum_{k=-L-1}^{L-1}\ell_k .
  \]
  As $k\to+\infty$ (resp.\ $-\infty$) both $a_k$ and $b_k$ tend to $B_m$ (resp.\ $A_m$), so $g_{\pm L}(z)\to0$. To compute the
  $\ell_k$, note that $cx+d=-(ad-bc)/(c\,\sigma_m(x)-a)$, so that $e^{\ell_k}=\frac{b_{k+1}-q}{a_{k+1}-q}$ with $q=\sigma_m(\infty)=a/c$. Let
  $u(x)=\frac{x-B_m}{x-A_m}$, so that $x-y=\frac{(B_m-A_m)(u(x)-u(y))}{(1-u(x))(1-u(y))}$ and, by~\eqref{eq:multiplier}, $u(\sigma_mx)=\mu u(x)$ with $\mu=\mu_m$.
  Then $u(b_k)=\mu^k\beta$, $u(a_k)=\mu^k\alpha$ and $u(q)=\mu$, where $\beta=u(\gamma B_n)$ and $\alpha=u(\gamma A_n)$ are finite and nonzero
  (two elements of a Schottky group with a common fixed point generate a cyclic group, and this is excluded by the choice of $\gamma$). Therefore
  \[
  e^{\ell_k}=\frac{(\mu^{k+1}\beta-\mu)(1-\mu^{k+1}\alpha)}{(\mu^{k+1}\alpha-\mu)(1-\mu^{k+1}\beta)}
  =\frac{(1-\mu^{k}\beta)(1-\mu^{k+1}\alpha)}{(1-\mu^{k}\alpha)(1-\mu^{k+1}\beta)},
  \]
  and the product over $-L-1\leq k\leq L-1$ telescopes to
  $\frac{(1-\mu^{-L-1}\beta)(1-\mu^{L}\alpha)}{(1-\mu^{-L-1}\alpha)(1-\mu^{L}\beta)}$, which tends to $\beta/\alpha$ since $|\mu|<1$. Finally
  \[
  \frac{\beta}{\alpha}=\frac{(\gamma B_n-B_m)(\gamma A_n-A_m)}{(\gamma B_n-A_m)(\gamma A_n-B_m)}=(B_m:A_m:\gamma(B_n):\gamma(A_n)).
  \]
  Summing over the orbits gives~\eqref{eq:periodmatrix} modulo $2\pi i$. (All the identities above hold modulo $2\pi i$; for all but
  finitely many terms the arguments of the logarithms are close to $1$ and the principal branches satisfy them exactly, so that only
  finitely many multiples of $2\pi i$ occur.)
\end{proof}

\begin{remark}\label{rem:identity}
  Formula~\eqref{eq:periodmatrix} is~\cite[Theorem~30]{BobenkoIntro2011} (with a different convention for the cross-ratio,
  which gives the same value), and agrees with~\cite[(5.14)]{CrowdyMarshallMcurves2011}. In~\cite[(4.4)]{Schmies2011}, on the other hand, the trivial double coset is excluded
  for all $m,n$. This is correct for $m=n$ (where the trivial double coset contributes $\log\mu_n$ instead)
  but not for $m\neq n$: there, by the proof above, the trivial double coset $H_mH_n$ contributes $\log(B_m:A_m:B_n:A_n)$. We
  also checked this numerically in the genus-$2$ example of Section~\ref{sec:experiments}: the sum over all
  double cosets with a representative of length at most $10$, with the trivial one included for $m\neq n$, agrees with the
  value of $B_{mn}$ computed with Algorithm~\ref{alg:periods} to $5\cdot10^{-10}$, whereas omitting it
  gives an error of $0.56$. The same omission explains the discrepancy with the Riemann matrix reported in~\cite[\S5.4]{FairchildRiosOrtiz2024},
  see Example~\ref{ex:fro}.
\end{remark}

Since $(B_{mn})$ is symmetric (Riemann's bilinear relations) with negative definite real part, the matrix
$\tau=\frac{1}{2\pi i}(B_{mn})$ is a Riemann matrix of the Riemann surface $\Omega/G$. The symmetry
$B_{mn}=B_{nm}$ is a convenient check of the output, since it is not used anywhere in the algorithm and is only satisfied up
to the truncation error: in the experiments below it holds to the accuracy we claim.

\section{The Schottky--Klein prime function}
\label{sec:prime}

Let $A,B\in\cF$ be two points, and consider
\[
\Theta_{A,B}(z)=\prod_{\sigma\in G}\frac{z-\sigma(B)}{z-\sigma(A)},\qquad
\Theta'_{A,B}(z)=\prod_{\sigma\in G\setminus\{1\}}\frac{z-\sigma(B)}{z-\sigma(A)},
\]
which converge absolutely, for $z\in\cF$, under (S). The argument of Section~\ref{sec:preliminaries} applies verbatim,
with $\cR_n$ replaced by $G$ and no exclusion in the seeds: the sums over the words beginning with $\sigma_i$ satisfy the system
of Theorem~\ref{thm:functional-equation} with $E^{(i)}=\log\frac{z-\sigma_iB}{z-\sigma_iA}$ for all $i\in I$.
Algorithm~\ref{alg:main} thus computes $\Theta_{A,B}$ with the same cost, and with an error bound of the same form, whose constants now
also depend on $A$ and $B$ (through the distance from the points $\sigma_i(A),\sigma_i(B)$ to $\partial\cD_i$).

Let $D$ be a circular domain in the plane (the unit disc with $N$ disjoint discs removed), and $G$ the Schottky group of its Schottky double,
generated by the compositions of the reflections in the unit circle and in the inner circles~\cite{CrowdyMarshall2007,CKGN2016}. The Schottky--Klein prime
function is given, when the product converges, by the formula of Baker~\cite[Ch.~12]{Baker1897} (see~\cite[(15), (16)]{CrowdyMarshall2007})
\[
\omega(z,w)=(z-w)\prod_{\sigma\in G''}\frac{(z-\sigma w)(w-\sigma z)}{(z-\sigma z)(w-\sigma w)},
\]
where $G''$ is a set of representatives of $(G\setminus\{1\})/\{\sigma\sim\sigma^{-1}\}$. Since the factors for $\sigma$ and
$\sigma^{-1}$ coincide,
\begin{equation}\label{eq:primesq}
\omega(z,w)^2=(z-w)^2\,\frac{\Theta'_{w,z}(w)}{\Theta'_{w,z}(z)}=(z-w)^2\prod_{\sigma\neq1}\frac{(z-\sigma w)(w-\sigma z)}{(z-\sigma z)(w-\sigma w)}.
\end{equation}
This is the function $X(z,w)=\omega(z,w)^2$ computed in~\cite{CrowdyMarshall2007}, which has the advantage of not depending on
the choice of a square root. It is the quantity computed by the function \texttt{SKprime} of our code, from two evaluations
of $\Theta'_{w,z}$.

\section{Numerical experiments}
\label{sec:experiments}

\subsection{Implementation and code availability}
The code, in Sage~\cite{sagemath}, and scripts reproducing every number of this section and of Examples~\ref{ex:fro} and~\ref{ex:x023}, are available
at \url{https://github.com/mmasdeu/archimedean_theta}. The file \texttt{schottky\_theta.py} is a self-contained implementation of Algorithms~\ref{alg:main}
and~\ref{alg:periods}, together with two baselines: summation of all words of length at most $L$ and the pruned depth-first summation
of Bogatyrev type (stop descending below $\sigma$ when the diameter of the outer disc $O(\sigma)$ is smaller than $\nu$; unlike the estimates in~\cite{Schmies2011,Lyamaev2022}
this has no proven error bound, so the tolerance was lowered until the true error was below the target). The prime function
(Section~\ref{sec:prime}) is computed with a first implementation, \texttt{theta\_arch.py}, which iterates on exponentials of power series and
uses the \texttt{darmonpoints} package; on the examples where both apply, the two implementations agree.
All timings are for a single core; those of Tables~\ref{tab:cost} and~\ref{tab:genus} and Figures~\ref{fig:cost} and~\ref{fig:genus} were measured on an
otherwise idle server, and the others on a desktop computer.

\subsection{Correctness checks}
We use the genus $2$ datum with discs of radius $1$ centred at $\pm2\pm2i$, with $\cD_{\pm1}$ centred at $-2\mp2i$ and
$\cD_{\pm2}$ at $2\mp2i$ (the example \texttt{examples/genus2\_basic.sage} of the repository), and generators determined by $\sigma_i=\bigl(\begin{smallmatrix}1&-a_i\\0&r_i\end{smallmatrix}\bigr)^{-1}
\bigl(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\bigr)\bigl(\begin{smallmatrix}1&-a_{-i}\\0&r_{-i}\end{smallmatrix}\bigr)$.
\begin{itemize}
  \item At $z=0.3+0.1i$, the difference between $\Omega_1(z)$ computed by Algorithm~\ref{alg:main} and
  the sum over all words of length at most $6$ ($729$ terms) is $3.1\cdot 10^{-6}$, and at length $8$ ($6561$ terms) it is $8.7\cdot10^{-8}$, which is the size of the truncation error of the naive sum.
  \item The period matrix computed at $p=70$ and $200$ bits of precision is symmetric to $2.7\cdot10^{-43}$, and the value of $B_{11}$ is independent
  of the chosen point $z$ to $3\cdot10^{-42}$. At $150$ bits the symmetry defect is $6.6\cdot10^{-21}$ for $p=30$ and $2.0\cdot10^{-37}$ for $p=60$. The
  entries are
  \begin{align*}
  B_{11}=B_{22}&=-2.57529955131838968508990935760\ldots,\\
  B_{12}=B_{21}&=-0.566042371936581272131437130\ldots
  \end{align*}
  \item The sums over double cosets~\eqref{eq:periodmatrix} up to length $L=5,7,9,10$ approach the computed matrix with
  errors $4\cdot10^{-6}$, $1\cdot10^{-7}$, $3\cdot10^{-9}$, $5\cdot10^{-10}$, in all four entries (this includes the trivial double coset, see Remark~\ref{rem:identity}).
  \item Theorem~\ref{thm:convergence} bounds the error on all of $\overline\cF$, including the boundary circles, where convergence is slowest. On
  $256$ points of the four circles, the maximal error of $\Omega_1$ is $1.4\cdot10^{-13}$ for $p=20$ and $8.3\cdot10^{-37}$ for $p=60$ (about
  $0.58$ digits per coefficient), against $4.4\cdot10^{-22}$ for $p=20$ at the interior point $0.3+0.1i$ used below (about one digit per coefficient).
\end{itemize}

\subsection{Cost versus precision}
We compare the new algorithm with the pruned depth-first summation on the genus $2$ datum above. The reference value of $\Omega_1(0.3+0.1i)$
is computed by the new algorithm with $p=150$ at $500$ bits (it agrees with the value for $p=130$ to $10^{-129}$). For each number $d$ of digits we report the smallest $p$ for which the new
algorithm, with the direct solver of Remark~\ref{rem:direct} in ball arithmetic of $3.33d+40$ bits, has error below $10^{-d}$, and, for the pruned summation, the largest
tolerance $\nu$ (in powers of $10$) for which the true error is below $10^{-d}$ (in $3.33d+40$-bit floating point arithmetic).
\begin{table}[t]
\centering
\small
\begin{tabular}{@{}r rrr c rrr@{}}
\toprule
 & \multicolumn{3}{c}{new algorithm (all $n$)} & & \multicolumn{3}{c}{pruned summation (one $n$)} \\
\cmidrule(lr){2-4}\cmidrule(lr){6-8}
$d$ & $p$ & time (s) & error & & $\nu$ & terms & time (s) \\
\midrule
5  & 5  & 0.005 & $3.7\cdot10^{-6}$  & & $10^{-5}$  & $567$ & 0.036 \\
10 & 10 & 0.018 & $2.3\cdot10^{-11}$ & & $10^{-13}$ & $560\,769$ & 35.6 \\
11 & 11 & 0.025 & $1.4\cdot10^{-12}$ & & $10^{-15}$ & $3\,178\,725$ & 218 \\
12 & 12 & 0.029 & $1.7\cdot10^{-14}$ & & $10^{-16}$ & $7\,712\,583$ & 440 \\
15 & 15 & 0.050 & $1.1\cdot10^{-16}$ & & $10^{-21}$ & $582\,414\,297$ & $32\,914$ \\
20 & 19 & 0.098 & $7.7\cdot10^{-21}$ & & \multicolumn{3}{c}{\emph{not attempted}} \\
40 & 39 & 1.27  & $4.7\cdot10^{-41}$ & & \multicolumn{3}{c}{\emph{not attempted}} \\
80 & 79 & 12.5  & $7.2\cdot10^{-81}$ & & \multicolumn{3}{c}{\emph{not attempted}} \\
\bottomrule
\end{tabular}
\caption{Cost versus precision for $\Omega_n$ in genus $2$. The new algorithm computes both $\Omega_1$ and $\Omega_2$ (and, from the same data, the period
matrix), including the precomputation (minimum of three runs); the pruned summation computes one function. Timings on an otherwise idle $32$-core machine, one core per run.}
\label{tab:cost}
\end{table}

\begin{figure}[t]
\centering
\begin{tikzpicture}
\begin{axis}[width=0.82\textwidth, height=6.2cm, ymode=log, xlabel={target accuracy $d$ (digits)}, ylabel={time (s)},
  xmin=0, xmax=82, ymin=3e-3, ymax=1e13, legend pos=north west, legend cell align=left,
  ytick={1e-2,1,1e2,1e4,1e6,1e8,1e10,1e12}]
\addplot[blue, mark=*, thick] coordinates {(5,0.005) (10,0.018) (11,0.025) (12,0.029) (15,0.050) (20,0.098) (40,1.27) (80,12.5)};
\addlegendentry{new algorithm (all $n$)}
\addplot[red, mark=square*, thick] coordinates {(5,0.036) (10,35.6) (11,218) (12,440) (15,32914)};
\addlegendentry{pruned summation (one $n$), measured}
\addplot[red, dashed, thick] coordinates {(15,32914) (20,32914*4^5)};
\addlegendentry{extrapolation: $\times4$ per digit}
\end{axis}
\end{tikzpicture}
\caption{Time to compute $\Omega_n$ in genus $2$ as a function of the number of correct digits (data of Table~\ref{tab:cost}). Note the logarithmic vertical axis: the new algorithm is polynomial in $d$, the pruned summation exponential. The dashed line is the extrapolation from $d=15$.}
\label{fig:cost}
\end{figure}
For the new algorithm one needs about one Taylor coefficient per digit for this datum, and the time grows
polynomially. For the pruned summation the number of terms is multiplied by about $4=10^{\delta/(1-\delta)}$, with $\delta\approx0.38$, for each additional digit (Section~\ref{sec:existing-algorithms})
($567$ terms for $5$ digits, $5.6\cdot10^5$ for $10$, $7.7\cdot10^6$ for $12$ and $5.8\cdot10^8$ for $15$, which took $9.1$ hours), so that by extrapolation about $6\cdot10^{11}$ terms, that is, about a year of single-core computation, would be needed for $20$ digits, where the new algorithm needs $0.1$ seconds. Moreover, the tolerance $\nu$ had to be taken
considerably smaller than the target accuracy (by a factor of about $10^{3}$ at $d=10$ and $10^{6}$ at $d=15$; for $d=15$, the tolerances $10^{-19}$, $10^{-20}$ and $10^{-21}$ gave errors $1.0\cdot10^{-14}$, $2.5\cdot10^{-15}$ and $5.7\cdot10^{-16}$), because the pruning criterion has no
error guarantee. The estimates of~\cite{Schmies2011,Lyamaev2022} are designed to provide one; we did not implement them, and we expect the
number of terms they need to obey the same frontier bound of Section~\ref{sec:existing-algorithms}.

This comparison has limitations, which we state explicitly. Both programs are written in Sage, but the pruned summation is a recursive
interpreted loop in multiprecision arithmetic (about $60$ microseconds per term), whereas the new algorithm spends its time in dense
linear algebra; a compiled implementation of the baseline in double or double-double arithmetic would be faster by two or three orders of
magnitude for $d\leq15$, which moves the crossover point to larger $d$ but does not change the exponential growth. For this reason we
report the number of terms, which does not depend on the implementation. Both methods were tuned with the help of a reference value, and the
pruned summation computes one function $\Omega_n$, although a single traversal could serve all $n$. The interior test point is a favourable
case for the new algorithm; on the boundary circles the error is larger (Section~\ref{sec:experiments}, Correctness checks). Finally, for the
prime function the natural competitor is the \textsc{Matlab} package \texttt{SKPrime}~\cite{CKGN2016}, which we could not run (it does not
work under Octave).
In the regime of about $5$ digits where the term-by-term methods are normally used, both approaches
are fast ($0.04$ seconds against $0.005$); the difference becomes decisive from about $10$ digits on.

\subsection{Dependence on the genus}
We use the discs of the family of Fairchild and R\'ios Ortiz~\cite{FairchildRiosOrtiz2024}: $2g$ discs of the same radius
$\tan\frac{\eta_2-\eta_1}{2}$ with centres at $\sec\frac{\eta_2-\eta_1}2\,e^{i((\eta_2+\eta_1)/2+2\pi j/g)}$ and the complex conjugates of these
(the function \texttt{point\_construct} of \texttt{theta\_arch.py}), with $\eta_1=\pi/(6g)$ and $\eta_2=\pi/(2g)$, but paired by the generators
$\sigma_i(z)=a_i+r_ir_{-i}/(z-a_{-i})$ of Section~\ref{sec:experiments}, not by those of~\cite{FairchildRiosOrtiz2024}; we call these the
\emph{ring family}. (The family of~\cite{FairchildRiosOrtiz2024} itself is treated in Example~\ref{ex:fro}.) In Table~\ref{tab:genus} we give the time needed to compute all $N=g$ functions $\Omega_n$ (precomputation
and $K=p$ iterations) at $160$ bits, for $p=20$ and $p=40$, and the accuracy obtained. As a measure of the accuracy
we use the symmetry defect $\max_{m,n}|B_{mn}-B_{nm}|$ of the period matrix, which the algorithm does not impose, and
the distance to the result for a larger $p$: for $g=2$ and $g=3$ the matrices for $p=40$ and $p=64$ differ by $2.8\cdot10^{-25}$ and $6.3\cdot10^{-26}$, respectively, so
that $p=40$ gives about $25$ digits, in agreement with the symmetry defect. For $g\geq4$ the symmetry defect is our only measure of
accuracy; it is a weaker test, since it is blind to errors which preserve the symmetry of the data. The convergence is geometric, at a rate of about
$0.6$ digits per Taylor coefficient; this was checked for $g=2$ and $g=3$ with the direct solver against $p=90$ (errors of the period matrix $3\cdot10^{-14}$, $3\cdot10^{-25}$ and $8\cdot10^{-37}$ for $p=20,40,60$ when $g=2$, and $2\cdot10^{-14}$, $6\cdot10^{-26}$ and $1\cdot10^{-37}$ when $g=3$).
\begin{table}[t]
\centering
\small
\begin{tabular}{@{}c c rr c rr@{}}
\toprule
 & min.\ gap & \multicolumn{2}{c}{$p=20$} & & \multicolumn{2}{c}{$p=40$} \\
\cmidrule(lr){3-4}\cmidrule(lr){6-7}
$g$ & between discs & time (s) & symm.\ defect & & time (s) & symm.\ defect \\
\midrule
2 & 0.499 & 0.4  & $3.9\cdot10^{-14}$ & & 1.9   & $7.7\cdot10^{-26}$ \\
3 & 0.342 & 0.8  & $4.7\cdot10^{-14}$ & & 6.7   & $5.9\cdot10^{-26}$ \\
4 & 0.259 & 1.9  & $2.9\cdot10^{-14}$ & & 17.7  & $6.9\cdot10^{-26}$ \\
5 & 0.208 & 3.8  & $2.2\cdot10^{-14}$ & & 32.3  & $6.2\cdot10^{-26}$ \\
6 & 0.174 & 6.7  & $1.6\cdot10^{-14}$ & & 57.2  & $3.3\cdot10^{-26}$ \\
\bottomrule
\end{tabular}
\caption{Period matrices for the ring family, $160$ bits. Times include the precomputation of the transfer matrices and the iterations for all $n$,
in unoptimized Sage code with dense matrix--vector products, on an idle machine (one core).}
\label{tab:genus}
\end{table}

\begin{figure}[t]
\centering
\begin{tikzpicture}
\begin{axis}[width=0.7\textwidth, height=6cm, xmode=log, ymode=log, xlabel={genus $g$}, ylabel={time (s)},
  xtick={2,3,4,5,6}, xticklabels={2,3,4,5,6}, legend pos=south east, legend cell align=left, xmin=1.8, xmax=6.6]
\addplot[blue, mark=*, thick] coordinates {(2,1.9) (3,6.7) (4,17.7) (5,32.3) (6,57.2)};
\addlegendentry{$p=40$}
\addplot[teal, mark=triangle*, thick] coordinates {(2,0.4) (3,0.8) (4,1.9) (5,3.8) (6,6.7)};
\addlegendentry{$p=20$}
\addplot[black, dashed, domain=2:6] {1.9*(x/2)^3};
\addlegendentry{$\propto g^{3}$}
\end{axis}
\end{tikzpicture}
\caption{Time to compute the period matrix of the ring family as a function of the genus (data of Table~\ref{tab:genus}).}
\label{fig:genus}
\end{figure}
Between $g=2$ and $g=6$ the time grows like $g^{3.1}$ for $p=40$ (and $g^{2.6}$ for $p=20$), in line with the $O(N^3)$ of Theorem~\ref{thm:convergence}, while the
accuracy for fixed $p$ does not degrade as $g$ grows and the discs get closer; a term-by-term summation of the same series to this accuracy is out of reach.

\subsection{Dependence on the geometry}
The rate of convergence in $p$ depends on how close the discs are to each other (through $R$ in~\eqref{eq:Ri}). We keep the genus $2$ datum
above but change the common radius of the four discs to $r$, so that the smallest distance between two discs is $4-2r$. The error of $\Omega_1(0.3+0.1i)$ for the truncation
order $p$ is shown in Figure~\ref{fig:geometry}.
\begin{figure}[t]
\centering
\begin{tikzpicture}
\begin{axis}[width=0.82\textwidth, height=6.4cm, ymode=log, xlabel={truncation order $p$}, ylabel={error of $\Omega_1(0.3+0.1i)$},
  xmin=0, xmax=64, ymin=1e-45, ymax=1e1, legend pos=south west, legend cell align=left, ytick={1e-40,1e-30,1e-20,1e-10,1}]
\addplot[blue, mark=*, thick] coordinates {(10,2.3e-11) (20,4.4e-22) (40,2.9e-42)};
\addlegendentry{$r=1.0$ (gap $2.0$)}
\addplot[teal, mark=square*, thick] coordinates {(10,1.3e-6) (20,7.5e-13) (40,5.9e-23) (60,2.6e-33)};
\addlegendentry{$r=1.6$ (gap $0.8$)}
\addplot[red, mark=triangle*, thick] coordinates {(10,9.4e-5) (20,8.1e-10) (40,5.0e-17) (60,2.6e-24)};
\addlegendentry{$r=1.8$ (gap $0.4$)}
\end{axis}
\end{tikzpicture}
\caption{Error of $\Omega_1(0.3+0.1i)$ as a function of the truncation order $p$ for the genus $2$ datum with four discs of common radius $r$ (smallest gap $4-2r$); the lines show geometric convergence, with slopes of about $1$, $0.5$ and $0.4$ digits per coefficient. All three curves use the direct solver of Remark~\ref{rem:direct} at each $p$; the reference is $p=120$ at $200$ bits. For $r=1$ the values agree with Table~\ref{tab:cost}.}
\label{fig:geometry}
\end{figure}
The convergence remains geometric, but the number of digits per coefficient drops from about $1$ to about $0.5$ and then about $0.4$ as the gap shrinks,
so that the cost for a fixed number of digits increases accordingly (as $p^3$ in our implementation). Thus for data which is close to the boundary of the range in which
the series converge, the break-even point with the adaptive summation moves, though the latter also becomes more expensive there:
we did not determine where it lies.

\subsection{The Schottky--Klein prime function of a triply connected domain}
Crowdy and Marshall~\cite[\S5]{CrowdyMarshall2007} compute $X=\omega^2$ for the triply connected domain obtained by removing from the unit disc the discs bounded by the circles
$|z-\frac12|=\frac1{10}$, $|z-\frac i2|=\frac1{10}$, at the points $z=\frac{-1-i}{2}$, $w=1$, and report the value
$2.39754812221763+1.76164377547306\,i$, obtained with their collocation method (with $15$ or more Laurent coefficients), the product truncated at level $9$ giving
$2.39754812221900+1.76164377547293\,i$. The quantity~\eqref{eq:primesq} for the Schottky datum of the double, computed with $K=p=40$ and $400$ bits of precision, is
\[
2.39754812221762939504895532890+1.76164377547305989716897377676\,i,
\]
which agrees with the published value up to $6.1\cdot10^{-16}$, the rounding of the published value; the values for $p=10$, $20$ and $40$
are within $7\cdot10^{-12}$, $8\cdot10^{-22}$ and $10^{-30}$ of the value for $p=60$, respectively (distances in $\CC$). The word-length
truncation of the product over words of length at most $7$ ($4372$ elements) is only accurate to $10^{-10}$.

\subsection{Comparison with collocation}\label{sec:collocation}
To compare the cost with the method of Crowdy and Marshall~\cite{CrowdyMarshall2007} we implemented it in Sage (\texttt{collocation.py}), in multiprecision
arithmetic: the two least-squares problems of~\cite[\S4]{CrowdyMarshall2007} are solved through the normal equations in ball arithmetic with $200$ bits (twice the
accuracy sought), with $N_2=2N_1+1$ collocation points per circle. With $N_1=10$, $15$, $20$ and $30$ it computes the value of the previous example
with errors $7\cdot10^{-8}$, $1.4\cdot10^{-10}$, $2\cdot10^{-14}$ and $5\cdot10^{-20}$, so it converges to the same value (in~\cite{CrowdyMarshall2007} the value is stable for $N_1=15$ at the digits shown, and our rate is somewhat slower).
We then computed $X=\omega^2$ at $z=\frac{1+i}{10}$, $w=\frac{2-i}{10}$ for the unit disc with $N$ discs of radius $\frac1{10}$ removed, centred at $\frac12e^{2\pi ik/N}$,
$N=2,3,4,6$, comparing with our first implementation (the same $p$ for all discs, $400$ bits) with errors measured against the value for $p=45$ (itself accurate to at least $10^{-46}$).
Both are single-core Sage programs run on the same server, not optimised, so only the trends are meaningful (Table~\ref{tab:colloc}).
\begin{table}[t]
\centering
\small
\begin{tabular}{@{}rrrrrrrrr@{}}
\toprule
 & \multicolumn{4}{c}{collocation~\cite{CrowdyMarshall2007}} & \multicolumn{4}{c}{this article} \\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}
$N$ & \multicolumn{2}{c}{$N_1=20$} & \multicolumn{2}{c}{$N_1=40$} & \multicolumn{2}{c}{$p=20$} & \multicolumn{2}{c}{$p=40$} \\
 & time & error & time & error & time & error & time & error \\
\midrule
$2$ & $5.5$ & $2\cdot10^{-15}$ & $45$ & $2\cdot10^{-27}$ & $0.6$ & $2\cdot10^{-28}$ & $4.1$ & $10^{-51}$ \\
$3$ & $22$ & $1\cdot10^{-15}$ & $168$ & $2\cdot10^{-27}$ & $1.7$ & $2\cdot10^{-28}$ & $11$ & $9\cdot10^{-52}$ \\
$4$ & $46$ & $1\cdot10^{-15}$ & $364$ & $1\cdot10^{-27}$ & $2.7$ & $2\cdot10^{-28}$ & $19$ & $5\cdot10^{-52}$ \\
$6$ & $171$ & $6\cdot10^{-16}$ & $1267$ & $2\cdot10^{-27}$ & $6.8$ & $3\cdot10^{-27}$ & $47$ & $3\cdot10^{-46}$ \\
\bottomrule
\end{tabular}
\caption{Collocation versus the algorithm of this article for the prime function of the unit disc with $N$ discs removed: time in seconds and error.}
\label{tab:colloc}
\end{table}
The cost of collocation is cubic in $N_1$ ($45$ against $5.5$ seconds for $N=2$ when $N_1$ goes from $20$ to $40$) and in $N$ (the time for $N_1=40$ is multiplied by $28$ when $N$ goes from $2$ to $6$, an exponent
of $3.0$), as predicted, whereas ours grows with exponent about $2.2$ in $N$ for $p=40$ in this implementation (which iterates, rather than solving the full linear system).
In the number of correct digits collocation gains about $0.68$ per coefficient, against about $1.3$ for the iteration, which is a factor of about $8$ in the cubic cost for a given accuracy;
for $10^{-27}$ our implementation is $76$, $99$, $136$ and $186$ times faster for $N=2,3,4,6$. The convergence of collocation is geometric, as expected, but its rate and
the conditioning of the least-squares problems are not controlled by any theorem; conversely it does not require the series to converge. The cubic exponent in $N$ is the worst case for both
(a direct solution of our system also costs $O(N^3p^3)$), so for the same linear algebra the difference is one of constants, which depend on the implementation
(we used the normal equations, which require working at twice the precision).

\section{Applications and outlook}
\label{sec:applications}

\subsubsection*{High-precision period matrices and algebraic recognition}
A motivation for this work is the passage from a Schottky group to an algebraic model of the curve
\cite{FairchildRiosOrtiz2024}, where the Riemann matrix $\tau=\frac1{2\pi i}(B_{mn})$ is approximated and the
coefficients of the equations are then recognized as algebraic numbers; numerical Schottky uniformization of real curves
along these lines goes back to~\cite{Hidalgo2001,HidalgoFigueroa2005,HidalgoSeppala2011,Seppala2004}, where the period matrices
are computed with Burnside's method. Recognition requires many more digits than the
moderate accuracy for which the term-by-term methods were designed (algorithms working from an algebraic model, such as~\cite{MolinNeurohr2019},
routinely provide thousands of digits), and the algorithm of this article makes this
regime accessible: for the ring family of Section~\ref{sec:experiments}, the entries of $\tau$ for $g=2,\dots,6$
are obtained to about $25$ digits in under a minute each, and more digits are obtained at the price of a larger $p$.
The algorithm gives $\omega_n$ and its derivatives as explicit expansions, which is what is needed to compute
a canonical embedding. The following example addresses a question left open in~\cite{FairchildRiosOrtiz2024}.

\begin{example}\label{ex:fro}
  For $g\geq2$ and $0<\eta<\theta<\pi/g$, Fairchild and R\'ios Ortiz~\cite[Definition~1.2]{FairchildRiosOrtiz2024} consider the Riemann
  surface $S_{\eta,\theta}$ uniformized by the classical Schottky group $G_{\eta,\theta}$ which pairs the circle $C_1$ orthogonal to the unit circle
  through $e^{i\eta}$ and $e^{i\theta}$ with its complex conjugate, and the rotations of these two circles by $2\pi j/g$, by the M\"obius
  transformations with positive real trace. It is an $M$-curve with an action of the dihedral group of order $2g$~\cite[Theorem~1.3]{FairchildRiosOrtiz2024}.
  They ask~\cite[\S6.3]{FairchildRiosOrtiz2024} whether, for $\eta,\theta$ rational, the algebraic curve $S_{\eta,\theta}$ is defined over $\QQ$. We
  interpret this as rational multiples of $\pi$, as in all their examples (for rational values in radians the generators are not algebraic,
  and the heuristic below does not apply). We keep their notation $\theta$, which should not be confused with the constant of (S).

  For $g=2$ they give an explicit model~\cite[Theorem~1.4]{FairchildRiosOrtiz2024}:
  \[
  S_{\eta,\theta}\colon\ y^2=\prod_{j=1}^6\bigl(x-\varphi(x_j)\bigr),\qquad \varphi=\omega_2/\omega_1,\qquad
  \{x_1,\ldots,x_6\}=\{1,-1,e^{i\theta},e^{i\eta},e^{i(\pi-\theta)},e^{i(\pi-\eta)}\},
  \]
  where the $x_j$ are the Weierstrass points. All the $x_j$ lie in $\overline\cF$, so the values $\varphi(x_j)$ are obtained from
  Algorithm~\ref{alg:main} by evaluating the derivative of~\eqref{eq:eval}. We take their main example $(\eta,\theta)=(\pi/12,\pi/4)$, with
  their generators. With $p=200$ and $p=240$ (at $463$ and $543$ bits, using the iteration with $K=p$, in $12$ and $21$ minutes on one core) we find
  values which agree to $109$ digits. They are real up to $10^{-133}$, and come in reciprocal pairs,
  $\varphi(1)\varphi(-1)=\varphi(e^{i\theta})\varphi(e^{i(\pi-\theta)})=\varphi(e^{i\eta})\varphi(e^{i(\pi-\eta)})=1$, as predicted by~\cite[Theorem~1.4]{FairchildRiosOrtiz2024}:
  \begin{center}
  \small
  \begin{tabular}{@{}lr@{\qquad}lr@{\qquad}lr@{}}
  \toprule
  $x$ & $\varphi(x)$ & $x$ & $\varphi(x)$ & $x$ & $\varphi(x)$ \\
  \midrule
  $1$  & $-0.0871498649464\ldots$ & $e^{i\theta}$ & $0.2448219613034\ldots$ & $e^{i\eta}$ & $-0.0848342452149\ldots$ \\
  $-1$ & $-11.474487087440\ldots$ & $e^{i(\pi-\theta)}$ & $4.0846008857870\ldots$ & $e^{i(\pi-\eta)}$ & $-11.787692546408\ldots$ \\
  \bottomrule
  \end{tabular}
  \end{center}
  The absolute Igusa--Clebsch invariants of this sextic are
  \begin{align*}
  i_1&=I_2^5/I_{10}=1.99201885562106565926941266785215031530687294214131\ldots\cdot10^{18},\\
  i_2&=I_2^3I_4/I_{10}=2.75425460826438289552834739852060109502274709119433\ldots\cdot10^{15},\\
  i_3&=I_2^2I_6/I_{10}=8.83956432659306652162737151388465431238315953340935\ldots\cdot10^{14}.
  \end{align*}
  If $S_{\pi/12,\pi/4}$ were defined over $\QQ$, they would be rational numbers. We searched for integer relations with $100$ digits (\texttt{algdep} of PARI/GP~\cite{PARI}
  on $i_k/10^{m_k}$, $m_k=\lfloor\log_{10}i_k\rfloor$): for every degree $d\leq8$ and each $k$, the integer polynomial
  of degree $d$ found has height about $10^{100/(d+1)}$ (for instance $10^{48.8}$ for $d=1$ and $10^{10.9}$ for $d=8$, for $i_1$), which is the size
  expected for a random real number, and the residual is of the size of the working precision. Hence none of the $i_k$ is a rational
  number of height less than about $10^{48}$, nor an algebraic number of degree at most $8$ whose minimal polynomial has height less
  than about $10^{10}$. This is numerical evidence that the answer to the question of~\cite{FairchildRiosOrtiz2024} is negative
  for $(\eta,\theta)=(\pi/12,\pi/4)$; it does not exclude rational invariants of very large height, for which there is no a priori bound. Since the parameters are rational multiples of $\pi$, we also tested whether $i_k/\pi^n$ is
  algebraic of degree at most $4$ for $|n|\leq30$, and whether $i_k$ lies in $\QQ[\pi]$ (degree at most $8$) or in $\QQ(\pi)$ (numerator and
  denominator of degree at most $4$) with small coefficients. In all cases the relations found by PARI are of the same size as for random
  real numbers of the same precision (forty controls, treated identically; of the $36$ values, one is $0.1$ below the smallest control height in its test, within the spread of the controls), so no such relation exists with small height either.
  With term-by-term summation, $109$ digits are far out of reach (Section~\ref{sec:experiments}).

  A negative answer is also what one expects. For these parameters, the generators of $G_{\eta,\theta}$ have algebraic entries. In genus $1$ the
  corresponding statement is a theorem: a Schottky group of genus $1$ is generated by $z\mapsto qz$, and if $q$ is algebraic, with $0<|q|<1$,
  then $j(q)$ is transcendental~\cite{BDGP1996} (the Mahler--Manin conjecture, which also holds $p$-adically). We are not aware of a
  higher-genus analogue, and~\cite{FairchildRiosOrtiz2024} state the question without conjecturing an answer.

  As a cross-check, for $g=3$ and the same parameters~\cite[\S5.4]{FairchildRiosOrtiz2024} report, using words of length at most $10$, the
  matrix with entries
  \[
  \tau_{nn}\approx0.39980895499614727\,i,\qquad \tau_{mn}\approx-0.004497125148\,i\quad(m\neq n).
  \]
  We obtain (with $p=50$ at $250$ bits, in $30$
  seconds, the result agreeing with $p=30$ to $7\cdot10^{-20}$ and being symmetric to $6\cdot10^{-30}$)
  \begin{align*}
  \tau_{nn}&=0.3998089542813468336080453\ldots\,i,\\
  \tau_{mn}&=-0.0527503085055988846587939\ldots\,i\qquad(m\neq n).
  \end{align*}
  The diagonal agrees with theirs to $7\cdot10^{-10}$, the accuracy of their truncation, but the off-diagonal entries do not. The difference
  is the contribution $\frac{1}{2\pi i}\log(B_m:A_m:B_n:A_n)=-0.04825318274977\ldots\,i$ of the trivial double coset: omitting it we get
  $-0.004497125755826\ldots\,i$, which agrees with their value to $6\cdot10^{-10}$, and the sums over double cosets with representatives of length at most $8$ reproduce
  their value when the trivial double coset is omitted, and ours when it is included. Thus their numerical matrix appears to omit that term (their Theorem~5.2
  includes it), as in Remark~\ref{rem:identity}; we refer to version~1 of~\cite{FairchildRiosOrtiz2024}, and we are not aware of a published version. Since the difference is not an integer matrix, it cannot be explained by a change of the $b$-cycles.
\end{example}

\subsubsection*{Schottky uniformization of modular and Shimura curves}
Every compact Riemann surface has a Schottky uniformization (Koebe), and every curve defined over a number field, in particular
every modular or Shimura curve, has one by a \emph{classical} Schottky group: according to the preprint~\cite{Hidalgo2017classical}, this holds for every Belyi curve, that is, by Belyi's theorem, for every curve defined over a number field. In the $p$-adic world
such uniformizations are arithmetic: Shimura curves at primes dividing the discriminant are Mumford curves uniformized by explicit Schottky
groups coming from definite quaternion algebras~\cite{Cerednik1976,Drinfeld1976} (see~\cite{amoros-milione} for their fundamental domains), and $X_0(p)$
is a Mumford curve over $\QQ_{p^2}$, since its stable reduction is totally degenerate after this extension~\cite{DeligneRapoport1973} (see~\cite{ManinDrinfeld1973} for $p$-adic Schottky groups and their periods); this is the setting of~\cite{masdeu-xarles,guitart-masdeu-xarles}. Over $\CC$
there is no such description, and the Schottky group has to be found by solving a transcendental problem. For a real curve with real
points one may look for it among the Schottky groups with an invariant circle, for which the series converge (Remark~\ref{rem:S-equiv}); Shimura curves
of discriminant $D>1$ have no real points~\cite{Shimura1975}, so this is only available for modular curves. The fast and precise
evaluation of the curve attached to a Schottky datum makes it possible to solve this problem by Newton's method, as the following example shows.

\begin{example}\label{ex:x023}
  The modular curve $X_0(23)$ has genus $2$ and the model
  \[
  X_0(23)\colon\ y^2=(x^3-x+1)(x^3-8x^2+3x-7),
  \]
  of discriminant $2^{12}\cdot23^6$; we checked that its numbers of points over $\mathbb{F}_p$, $3\leq p<60$, $p\neq 23$, agree with the
  coefficients of the newform of level $23$, whose coefficient field is $\QQ(\sqrt5)$. Its Igusa--Clebsch invariants are
  $(I_2,I_4,I_6,I_{10})=(-1144,\,669124,\,-17115576,\,606355001344)$. The sextic has two real roots, so $X_0(23)(\RR)$ has one real
  oval. We use the real Schottky model of Bogatyrev for such curves~\cite{Bogatyrev2012} (see~\cite[\S1]{Lyamaev2022}): two discs
  $\cD_j$ of centre $c_j>0$ and radius $r_j$ on the positive real half-line, $\cD_{-j}=-\cD_j$, and the generators
  \[
  \sigma_j(u)=c_j+\frac{r_j^2}{u+c_j}\qquad(j=1,2),
  \]
  which map $\Ch\setminus\cD_{-j}$ onto $\overline{\cD_j}$. These groups preserve the real line, so the series converge by Schottky's criterion
  (see~\cite[Theorem~31]{BobenkoIntro2011}). Then $u\mapsto-u$ induces the hyperelliptic involution, and the Weierstrass points are
  $0$, $\infty$ and $c_j\pm ir_j$. As in Example~\ref{ex:fro}, the values of $\varphi=\omega_2/\omega_1$ at these six points give a model
  $y^2=\prod_k(x-\varphi_k)$ of the curve, and hence its absolute invariants $(I_4/I_2^2,\ I_2I_4/I_6,\ I_{10}/I_2^5)$ as functions of the
  three real parameters $(r_1,c_2,r_2)$, where we normalize $c_1=1$. We solved the three equations expressing that these invariants equal those of
  $X_0(23)$. A coarse search and a Nelder--Mead minimization at $p=14$ give a starting point, which we refined by Newton's method
  (with the Jacobian computed once, by finite differences, at $p=60$), successively with $p=60$, $120$ and $200$; the whole computation takes
  about $50$ minutes on one core. The discs nearly touch (the gap between $\cD_1$ and $\cD_2$ is about $0.042$), and the convergence is about $0.24$
  digits per Taylor coefficient. We obtain
  \begin{align*}
  r_1&=0.1666951580327760609280827065351772813460\ldots,\\
  c_2&=1.2882213331143925842580331840891030391867\ldots,\\
  r_2&=0.0791274636493090013132263963050174632062\ldots,
  \end{align*}
  with an error of at most about $10^{-40}$ (the final residual is $4\cdot10^{-44}$ and the inverse Jacobian has norm about $150$; the solution at $p=120$ agrees
  to $4\cdot10^{-30}$). The fixed points and multipliers of the generators are
  \begin{align*}
  B_1=-A_1&=1.01379843939097291614505198297\ldots, & \mu_1&=-0.00685194661047901954277103106\ldots,\\
  B_2=-A_2&=1.29064920043929657297740608307\ldots, & \mu_2&=-0.000941445990915403693931876068\ldots
  \end{align*}
  (the multipliers are negative because the $\sigma_j$ reverse the orientation of the real line).

  As an independent check, we computed the Riemann matrix of this Schottky group with Algorithm~\ref{alg:periods}, which was not used in the
  fit:
  \begin{align*}
  \tau_{11}&=-\tfrac12+0.81287551976127870033\ldots\,i, & \tau_{12}=\tau_{21}&=0.68070897251165259288\ldots\,i,\\
  \tau_{22}&=-\tfrac12+1.2292513977736790783\ldots\,i
  \end{align*}
  (with $p=120$; the entries are defined modulo $1$, and we chose them to make $\tau$ symmetric). An integer relation search
  (\texttt{integer\_matrix\_relations} in Sage, at $90$ bits) between $\tau$ and the Riemann matrix of the model above, computed by numerical
  integration with Sage's \texttt{RiemannSurface} (see~\cite{BSZ2018}), finds a $\ZZ$-module of relations of rank $2$, which contains a symplectic matrix of
  determinant $1$. So the two matrices are equivalent under $\mathrm{Sp}_4(\ZZ)$, and the rank $2$ reflects the real multiplication of $J_0(23)$ by
  $\ZZ[\frac{1+\sqrt5}{2}]$.

  Here the equation of $X_0(23)$ was an input, and fitting the parameters of a Schottky group to a given curve is a classical approach to
  numerical uniformization~\cite{Hidalgo2001,HidalgoFigueroa2005,Seppala2004}; what is new is the precision. Obtaining the Schottky group of a modular curve from the modular side alone is possible in
  principle: the covering map $\mathfrak H\to\Omega$ is a quotient of solutions of a differential equation $y''+\frac12 qy=0$, where $q$ is a
  $\Gamma_0(N)$-invariant quadratic differential whose $3g-3$ coefficients are fixed by requiring trivial monodromy along the chosen $a$-cycles.
  With such a computation, Algorithm~\ref{alg:main} would then produce the canonical model of the curve without using any equation. We leave
  this for future work.
\end{example}

\subsubsection*{Conformal mapping and potential theory in multiply connected domains}
The Schottky--Klein prime function is the basic ingredient of the explicit formulas for conformal maps from multiply connected
circular domains, Green's functions and harmonic measures, and of the solutions of boundary value problems arising in fluid
mechanics and elsewhere~\cite{CrowdyMarshall2007,CKGN2016,Crowdy2020}. As Section~\ref{sec:prime} shows, the algorithm applies to
this setting whenever the product converges (for instance for all Schottky groups of Fuchsian type,
see~\cite{CrowdyMarshallMcurves2011}), giving the prime function to $30$ digits in about $4$ seconds with the first implementation (Section~\ref{sec:experiments}).
It is also natural to use it to uniformize real hyperelliptic $M$-curves as in~\cite{CrowdyMarshallMcurves2011}.

\subsubsection*{Other applications}
Schottky models of real hyperelliptic curves are used for Chebyshev polynomials for several intervals and multiband
filters~\cite{Bogatyrev1999,Bogatyrev2012,BGL2017}, the capacity of condensers~\cite{Bezrodnykh2019}, the filtration under a stepped
dam~\cite{BogatyrevGrigorev2020} and planar magnetic nanoelements~\cite{Bogatyrev2017}; Schottky uniformization also appears in finite-gap
integration~\cite{Bobenko1988,BBEIM1994,LiuYouYuan2026}, in dimer models~\cite{BBS2024,BB2024}, and in the construction of higher-genus
polylogarithms~\cite{BBILZ2024,BergerBroedel2026}, whose integration kernels are Poincar\'e series. These computations often involve large genus or
slow convergence, where term-by-term summation is most expensive, but where the rates $\kappa,\lambda$ of Theorem~\ref{thm:convergence} also
deteriorate (Figure~\ref{fig:geometry}); the break-even point should be determined on specific examples. We expect the method to extend to other
Poincar\'e series (of higher weight, or twisted by characters), with $\cT$ replaced by a weighted composition operator.

\subsubsection*{The $p$-adic origin of the algorithm}
This work originated in the $p$-adic setting, where the analogue of $\Omega_n$ is the theta function of a Mumford curve, an infinite
product over the Schottky group~\cite{GvdP,FvdP}. Masdeu and Xarles~\cite{masdeu-xarles} gave an iterative algorithm, polynomial in the
precision, which computes these theta functions as locally analytic functions without enumerating group elements (see also~\cite{MR,negrini-masters}
for earlier algorithms, \cite{guitart-masdeu-xarles,darmon-vonk} for applications to $p$-adic singular moduli, and~\cite{KMMvdP2026} for $p$-adic
Abelian integrals and heights). The algorithm of this article is its archimedean counterpart: the structure which makes both work is the same,
namely that the subtree below a vertex depends only on its last letter, and that the sum over it is an analytic function on the complement of a
disc. What is specific to the archimedean case is that convergence rests on hypothesis~(S), and that the functions are logarithms of products,
which require the normalization at $\infty$.

\subsubsection*{Limitations and further work}
We assume (S), that is, the absolute convergence of the series, as all the algorithms which sum them; the methods
of~\cite{CrowdyMarshall2007,CKGN2016} do not need it, but only apply to the prime function of planar domains. It would be interesting
to know whether the equation $(1-\cT)F=E$, which makes sense without (S), can be used beyond the range of convergence. We have not compared with the methods of Mityushev and coauthors~\cite{MityushevRogosin2000,MityushevRylko2013} or with \texttt{SKPrime}~\cite{CKGN2016}.
The constants in Theorem~\ref{thm:convergence} are not explicit, and we have not analysed rounding errors: it would be desirable to have explicit and rigorous a priori bounds for $p$ and $K$, as
are available in the algorithms of Schmies and Lyamaev, using the methods of~\cite{BandtlowJenkinson2008} and interval arithmetic for the matrix
recurrence (we noted that a naive use of ball arithmetic produces enclosures which are too wide for the larger examples).
The algorithm treats the discs as given: finding Schottky data for a given curve (``numerical uniformization'', see e.g.~\cite{Seppala2004}) is
a separate problem, which is not addressed here.

\printbibliography

\end{document}