COMPUTING RIEMANN THETA FUNCTIONS 1. Motivation An Abelian function is a 2g-fold ...

26 Pages·2004·3.8 MB·English

MATHEMATICSOFCOMPUTATION
Volume73,Number247,Pages1417{1442
S0025-5718(03)01609-0
ArticleelectronicallypublishedonDecember19,2003
COMPUTING RIEMANN THETA FUNCTIONS
BERNARDDECONINCK,MATTHIASHEIL,ALEXANDERBOBENKO,
MARKVANHOEIJ,ANDMARCUSSCHMIES
Abstract. The Riemann theta function is a complex-valued function of g
complexvariables. Itappearsintheconstructionofmany(quasi-)periodicso-
lutionsofvariousequationsofmathematicalphysics. Inthispaper,algorithms
for its computation are given. First, a formula is derived allowing the point-
wiseapproximationofRiemannthetafunctions,witharbitrary,user-speci(cid:12)ed
precision. Thisformulaisusedtoconstructauniformapproximationformula,
againwitharbitraryprecision.
1. Motivation
An Abelian function is a 2g-fold periodic, meromorphic function of g complex
variables. Thus, Abelian functions are generalizations of elliptic functions to more
thanonevariable. Justasellipticfunctionsareassociatedwithellipticsurfaces,i.e.,
Riemannsurfacesofgenus1,Abelianfunctions areassociatedto Riemannsurfaces
of higher genus. As in the elliptic case, any Abelian function can be expressed as
a ratio of homogeneous polynomials of an auxiliary function, the Riemann theta
function [11]. Many di(cid:11)erential equations of mathematical physics have solutions
thatarewrittenintermsofAbelianfunctions, andthus intermsofRiemanntheta
functions (see [3, 7] and references therein). Thus, to compute these solutions,
one needs to compute either Abelian functions or theta functions. Whittaker and
Watson[20, x21.1]state,\When itis desiredto obtainde(cid:12)nite numericalresultsin
problems involving Elliptic functions, the calculations are most simply performed
with the aid of certain auxiliary functions known as Theta-functions." Whittaker
and Watson are referring to Jacobi’s theta functions, but the same can be said for
Abelian functions, using the Riemann theta function. This paper addresses the
problem of computing values of the Riemann theta function and its derivatives.
2. Definition
The Riemann theta function is de(cid:12)ned by
X
(1) (cid:18)(zj(cid:10))= e2(cid:25)i(12n(cid:1)(cid:10)(cid:1)n+n(cid:1)z);
n2Zg
ReceivedbytheeditorJune7,2002.
2000 Mathematics Subject Classi(cid:12)cation. Primary14K25,30E10,33F05, 65D20.
Key words and phrases. Riemanntheta function, pointwise approximation, uniform approxi-
mation.
(cid:13)c2003 American Mathematical Society
1417
1418 B.DECONINCK,M.HEIL,A.BOBENKO,M.VANHOEIJ,ANDM.SCHMIES
where z 2 Cg, (cid:10) 2 Cg(cid:2)g, such that (cid:10) is symmetric ((cid:10)T = (cid:10)) and the imaginary
part of (cid:10), Im((cid:10)), is strictly positive de(cid:12)nite. Such an (cid:10) is called a Riemann
matrix. Also,
Xg
n(cid:1)z = n z ;
i i
i=1
the scalar product of the integer vector n with z; and
Xg
n(cid:1)(cid:10)(cid:1)n= (cid:10) n n :
ij i j
i;j=1
The positive de(cid:12)niteness of Im((cid:10)) guaranteesthe convergenceof (1), for all values
of z. Then the series (1) converges absolutely in both z and (cid:10) and uniformly
on compact sets. Thus, it de(cid:12)nes a holomorphic function of both z and (cid:10). The
function de(cid:12)ned by (1) is also known as a multidimensional theta function, or as
a theta function of several variables. There are as many di(cid:11)erent conventions for
writing the Riemann theta function as there are names for it. These di(cid:11)erent
conventions di(cid:11)er from (1) by at worst a complex scaling transformation on the
arguments. An extensive overviewofthe wealthofproperties of(cid:18)(zj(cid:10)) is found in
[14, 15, 16].
Remarks. (cid:15) The Riemann theta function was devised by Riemann as a gen-
eralization of Jacobi’s theta functions of one variable [12] for solving the
Jacobi inversion problem on general compact connected Riemann surfaces
[17]. For these purposes Riemann considered only theta functions associ-
ated with Riemann surfaces. Riemann theta functions in their full gen-
erality, as de(cid:12)ned in (1), were considered (cid:12)rst by Wirtinger [21], whose
convention for the arguments is adopted here.
(cid:15) Often, so-called Riemann theta functions with characteristics are consid-
ered[14]. Suchthetafunctionswithcharacteristicsareuptoanexponential
factor Riemann theta functions (1) evaluated at a shifted argument. Thus
the computation of the Riemann theta function (1) also allowsthe compu-
tation of theta functions with arbitrary characteristics.
(cid:15) In many applications, the Riemann theta function (1) originates from a
speci(cid:12)c Riemann surface, i.e., the Riemann matrix (cid:10) is the normalized
periodmatrix of the Riemann surface. Obtaining this periodmatrix is a
nontrivial problem, which is now also reduced to a black-box program.
This issue was addressed in [6].
(cid:15) Someofthe algorithmsgiveninthispaperhavebeenusedalreadybysome
of the authors (for instance, tori with constant mean curvature, and Will-
more tori with umbilic lines were constructed in [10] using the results of
[2, 4]; multiphase solutions of the Kadomtsev-Petviashvili equation were
constructedusing Schottky uniformizationin [3,5]) but they werenot eas-
ily available for use by others. Now, these algorithms are implemented as
black-box programs in Maple and Java. The maple implementation is in-
cluded in the Maple distribution as of Maple 8. It is also available from
http://www.math.fsu.edu/~hoeij/RiemannTheta/. TheJavaimplemen-
tationisavailablefromwww-sfb288.math.tu-berlin.de/~jem. Theseim-
plementations are discussed in the appendices.
COMPUTING RIEMANN THETA FUNCTIONS 1419
3. Rewriting the Riemann theta function
The Fourier series representation (1) is the starting point for the computation
of Riemann theta functions. Separating both z and (cid:10) in their real and imaginary
parts, z =x+iy, (cid:10)=X +iY, and using Y =YT, we obtain
X
(cid:18)(zj(cid:10)) = e2(cid:25)i(12n(cid:1)(cid:10)(cid:1)n+n(cid:1)z)
n2Zg
X
= e2(cid:25)i(12n(cid:1)X(cid:1)n+n(cid:1)x+i(12n(cid:1)Y(cid:1)n+n(cid:1)y))
n2Zg
X
= e2(cid:25)i(12n(cid:1)X(cid:1)n+n(cid:1)x)e(cid:0)2(cid:25)(12n(cid:1)Y(cid:1)n+n(cid:1)y)
n2Zg
X (cid:16) (cid:17)
= e2(cid:25)i(12n(cid:1)X(cid:1)n+n(cid:1)x)e(cid:0)2(cid:25) 12(n+Y(cid:0)1y)(cid:1)Y(cid:1)(n+Y(cid:0)1y)(cid:0)12y(cid:1)Y(cid:0)1(cid:1)y
n2Zg
X (cid:16) (cid:17) (cid:16) (cid:17)
(2) = e(cid:25)y(cid:1)Y(cid:0)1(cid:1)y e2(cid:25)i(12n(cid:1)X(cid:1)n+n(cid:1)x)e(cid:0)(cid:25) n+Y(cid:0)1y (cid:1)Y(cid:1) n+Y(cid:0)1y
n2Zg
At this point, all exponential growth has been isolated: the factor multiplying the
sum grows double-exponentially as the components of z leave the real line, due to
thefactthatY(cid:0)1isstrictlypositivede(cid:12)nite. ThisfollowssinceY isstrictlypositive
de(cid:12)nite, which also guarantees the existence of Y(cid:0)1. All terms remaining in the
sum are either oscillating (the (cid:12)rst factor of every term) or a damped exponential
(the second factor of every term).
Since for most applications, the exponential growth term cancels out, it will
be disregarded in almost all that follows: any statements of approximation,
pointwise or uniform, of the Riemann theta function pertain to the in-
(cid:12)nite sum in (2), without considering the exponential growth.
Let [V] be the vector with integer component closest to V, and let [[V]] =
V (cid:0)[V]. Continuing the rewriting of the theta function,
X
(cid:18)(zj(cid:10))=e(cid:25)y(cid:1)Y(cid:0)1(cid:1)y e2(cid:25)i(12n(cid:1)X(cid:1)n+n(cid:1)x)
n2Zg
(cid:16) h i hh ii(cid:17) (cid:16) h i hh ii(cid:17)
(cid:2)e(cid:0)(cid:25) n+ Y(cid:0)1y + Y(cid:0)1y (cid:1)Y(cid:1) n+ Y(cid:0)1y + Y(cid:0)1y
X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
=e(cid:25)y(cid:1)Y(cid:0)1(cid:1)y e2(cid:25)i 12 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x
n2Zg
(cid:16) hh ii(cid:17) (cid:16) hh ii(cid:17)
(3) (cid:2)e(cid:0)(cid:25) n+ Y(cid:0)1y (cid:1)Y(cid:1) n+ Y(cid:0)1y :
Thelaststepisachievedbyshiftingthesummationindexn. Fromthisformulation,
it is clear that the size of each one of the terms in the in(cid:12)nite sum is controlled by
its second exponential factor.
Similarly,formulaeareworkedoutforthederivativesofthetafunctions. Denote
theN-thorderdirectionalderivativeoftheg-variableRiemannthetafunctionalong
the g-dimensional vectors k(1), :::, k(N) as
Xg @N(cid:18)(zj(cid:10))
(4) D(k(1);:::;k(N))(cid:18)(zj(cid:10))= k(1):::k(N) :
i1;:::;iN=1 i1 iN @zi1:::@ziN
1420 B.DECONINCK,M.HEIL,A.BOBENKO,M.VANHOEIJ,ANDM.SCHMIES
Then
D(k(1);:::;k(N))(cid:18)(zj(cid:10))
X
=(2(cid:25)i)N (k(1)(cid:1)n):::(k(N)(cid:1)n)e2(cid:25)i(12n(cid:1)(cid:10)(cid:1)n+n(cid:1)z)
n2Zg X(cid:16) (cid:0) (cid:2) (cid:3)(cid:1)(cid:17) (cid:16) (cid:0) (cid:2) (cid:3)(cid:1)(cid:17)
=(2(cid:25)i)Ne(cid:25)y(cid:1)Y(cid:0)1(cid:1)y k(N)(cid:1) n(cid:0) Y(cid:0)1y ::: k(N)(cid:1) n(cid:0) Y(cid:0)1y
n2Zg
(cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
(cid:2)e2(cid:25)i 12 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x
(cid:16) hh ii(cid:17) (cid:16) hh ii(cid:17)
(5) (cid:2)e(cid:0)(cid:25) n+ Y(cid:0)1y (cid:1)Y(cid:1) n+ Y(cid:0)1y ;
using similar steps as before. The exponential growth is still factored out, but the
remainingin(cid:12)nite sumnowcontainstermsthatgrowalgebraicallyasthe argument
oftheRiemannthetafunctionleavestherealline. Theorderofthisgrowthisequal
to the order of the derivative. However, the size of individual terms within the
in(cid:12)nite sum is againdeterminedby the lastfactorwhichis exponentially decaying.
Remark. Riemann theta functions of genus greater than one have been computed
before, for instance in [8], where genus three was considered. There four distinct
representations of these Riemann theta functions were used, based on whether
the di(cid:11)erent phases behaved as trigonometric or hyperbolic functions (two limits
of elliptic functions), depending on the parameters in the Riemann matrix. These
distinctrepresentationswererequiredtohaveamanageablenumberofcontributing
terms to the series. Such a variety of representations is not needed here: all limit
cases are incorporated in the form (3). This is obtained by the separation of the
real and imaginary parts of both the argument z and the Riemann matrix (cid:10).
4. Pointwise approximation
The formulae (3) and (5) are the basis for the approximationof theta functions
andtheir derivatives. All approximationsarebasedon determining whichterms of
the in(cid:12)nite sum are dominant.
For the Riemann theta function,
(cid:18)(zj(cid:10))e(cid:0)(cid:25)y(cid:1)Y(cid:0)1(cid:1)y
X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
2(cid:25)i 1 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x
= e 2
n2Zg
(cid:16) hh ii(cid:17) (cid:16) hh ii(cid:17)
(cid:2)e(cid:0)(cid:25) n+ Y(cid:0)1y (cid:1)Y(cid:1) n+ Y(cid:0)1y
X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
2(cid:25)i 1 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x
= lim e 2
R!1
(cid:16)SR hh ii(cid:17) (cid:16) hh ii(cid:17)
(6) (cid:2)e(cid:0)(cid:25) n+ Y(cid:0)1y (cid:1)Y(cid:1) n+ Y(cid:0)1y ;
where
(cid:8) (cid:0) (cid:2)(cid:2) (cid:3)(cid:3)(cid:1) (cid:0) (cid:2)(cid:2) (cid:3)(cid:3)(cid:1) (cid:9)
(7) S = n2Zgj(cid:25) n+ Y(cid:0)1y (cid:1)Y (cid:1) n+ Y(cid:0)1y <R2 :
R
Below, we show this limit exists, hence proving that the Riemann theta function
is well de(cid:12)ned. This proof is di(cid:11)erent from that found in many references (see
[14], for example) in that it also allows the estimation of the error made in the
COMPUTING RIEMANN THETA FUNCTIONS 1421
approximation of the truncation of the series by only considering a (cid:12)nite value of
R.
SinceY isstrictlypositivede(cid:12)nite,ithasaCholeskydecopmpos(cid:0)ition,(cid:2)Y(cid:2) =TT(cid:3)(cid:3)T(cid:1).
Let (cid:3) be the lattice of all vectors v(n) of the form v(n)= (cid:25)T n+ Y(cid:0)1y ,
for n2Zg. Then S is rewritten as
R n (cid:12) o
(cid:12)
(8) S = v(n)2(cid:3)(cid:12)jjv(n)jj<R :
R
Thus,the approximationrequires (cid:12)nding alllattice points in (cid:3) that are alsoinside
the g-dimensional sphere jjv(n)jj=R. Let
(9)
X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
(cid:18)R( j(cid:10))=e(cid:25)y(cid:1)Y(cid:0)1(cid:1)y e2(cid:25)i 12 n(cid:0)Y(cid:0)1y (cid:1)X(cid:1) n(cid:0)Y(cid:0)1y + n(cid:0)Y(cid:0)1y (cid:1)x e(cid:0)jjv(n)jj2;
z
SR
then
(10) (cid:15)(R)=j(cid:18)(zj(cid:10))(cid:0)(cid:18) (zj(cid:10))je(cid:0)(cid:25)y(cid:1)Y(cid:0)1(cid:1)y
(cid:12) R (cid:12)
(cid:12) (cid:12)
(cid:12)X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17) (cid:12)
=(cid:12)(cid:12) e2(cid:25)i 12 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x e(cid:0)jjv(n)jj2(cid:12)(cid:12)
(cid:12) (cid:12)
(cid:3)nSR(cid:12) (cid:12)
X (cid:12) (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17) (cid:12)
(cid:20) (cid:12)(cid:12)e2(cid:25)i 12 n(cid:0) Y(cid:0)1y (cid:1)X(cid:1) n(cid:0) Y(cid:0)1y + n(cid:0) Y(cid:0)1y (cid:1)x e(cid:0)jjv(n)jj2(cid:12)(cid:12)
(cid:3)XnSR
(11) = e(cid:0)jjv(n)jj2:
(cid:3)nSR
Thus (cid:15)(R) is the absolute error of the oscillatory part of the Riemann theta
function, due to the truncation to a (cid:12)nite number of terms. In practice, this
oscillatory part is of order 1, so (cid:15)(R) is also a measure of the relative error.
In orderto estimate this last sum, the following theorems and lemmas are used.
Theorem 1 (Mean-value theoremfor subharmonic functions). Let (cid:10) be a domain
in Rg. Let u be a twice continuously di(cid:11)erentiable function on (cid:10), continuous on
the boundary of @(cid:10) of (cid:10), which is subharmonic on (cid:10). Then for any ball Bg(y) in
R
(cid:10) of radius R and center y
Z
1
(12) u(y)(cid:20) u(x)dx:
RgVolB1g(0) Bg(y)
R
A proof of this theorem is found in [9].
Lemma 1. For every integer p(cid:21)0 and any dimension g >0 the function u(y)=
e(cid:0)jjyjj2jjyjjp with y 2Rg is subharmonic on RgnBg(0), with
q R
p
1
R= g+2p+ g2+8p:
2
Proof. Sinceu(y)dependsonlyonr =jjyjj,thiscalculationisdoneing-dimensional
spherical coordinates.
@2u g(cid:0)1@u
(cid:1)u(y) = +
@r2 (cid:18) r @r (cid:19)
p2+(g(cid:0)2)p
= 2u(y) 2r2(cid:0)(g+2p)+ :
2r2
1422 B.DECONINCK,M.HEIL,A.BOBENKO,M.VANHOEIJ,ANDM.SCHMIES
q
p
This last factor is positive if r > 1 g+2p+ g2+8p, from which the result
2
follows. (cid:3)
Lemma 2. Let (cid:3) be a g-dimensional a(cid:14)ne lattice in Rg, i.e., (cid:3) = fXn+xj
n2Zgg, with X 2Gl(n;R), x2Rg, and let p2Z, positive. Then
(cid:18) (cid:19) (cid:18) (cid:19)
(13) X jjyjjpe(cid:0)jjyjj2 (cid:20) g 2 g(cid:0) g+p;(R(cid:0)(cid:26)=2)2 ;
2 (cid:26) 2
y2(cid:3);jjyjj(cid:21)R
R
where (cid:0)(z;d) = 1tz(cid:0)1e(cid:0)zdz, the incomplete Gamma function [1],
d q
p
1
R> g+2p+ g2+8p+(cid:26)=2;
2
and
(cid:26)=minfjjx(cid:0)yjj jx;y 2(cid:3);x6=yg:
Proof. By the previous lemmas,
Z
X X
jjyjjpe(cid:0)jjyjj2 (cid:20) 1 jjxjjpe(cid:0)jjxjj2dx
y2(cid:3);jjyjj(cid:21)R y2(cid:3);jjyjj(cid:21)R(cid:18)((cid:26)(cid:19)=2)gVolB1g(0) ZB(cid:26)g=2(y)
(14) = 1 2 g X jjxjjpe(cid:0)jjxjj2dx
VolB1g(0)(cid:18)(cid:26)(cid:19) yZ2(cid:3);jjyjj(cid:21)R B(cid:26)g=2(y)
(cid:20) 1 2 g jjxjjpe(cid:0)jjxjj2dx
VolB1g(0) (cid:26)(cid:18) (cid:19) RZgnBRg(cid:0)(cid:26)=2(0)
VolSg(cid:0)1(0) 2 g 1
= 1 e(cid:0)r2rg(cid:0)1+pdr
VolBg(0) (cid:26)
(cid:18) (cid:19)1 Z R(cid:0)(cid:26)=2
g 2 g 1
= e(cid:0)tt(g+p)=2(cid:0)1dt
2 (cid:26)
(cid:18) (cid:19) ((cid:18)R(cid:0)(cid:26)=2)2 (cid:19)
g 2 g g+p
= (cid:0) ;(R(cid:0)(cid:26)=2)2 :
2 (cid:26) 2
Here VolSg(cid:0)1(0)=VolVg(0)=g=r was used. (cid:3)
r r
Remark. The inequality following (14) is quite crude. Speci(cid:12)cally, it adds all the
space in between the balls around the lattice points. Taking this into account,
another estimate is as follows:
X
jjyjjpe(cid:0)jjyjj2
y2(cid:3);jjyjj(cid:21)R
(cid:18) (cid:19) Z
(cid:20) 1 2 g X jjxjjpe(cid:0)jjxjj2dx
VolB1g(0)(cid:18)(cid:26)(cid:19) y2(cid:3);jjyjj(cid:21)ZR B(cid:26)g=2(y)
(cid:25) 1 2 gc((cid:3);(cid:26)=2) jjxjjpe(cid:0)jjxjj2dx:
VolB1g(0) (cid:26) RgnBRg(cid:0)(cid:26)=2(0)
Here c((cid:3);(cid:26)=2) is the \(cid:12)ll factor" of the lattice (cid:3) with balls of radius (cid:26)=2: it is
the fraction of the lattice volume that is (cid:12)lled if a ball of radius (cid:26)=2 is placed at
every lattice point. This estimate is justi(cid:12)ed, since the function being integrated
only depends on the radial variable. Thus, within thin radial shells the function is
COMPUTING RIEMANN THETA FUNCTIONS 1423
Figure 1. Fill factor c((cid:3);(cid:26)=2) for the case of a square lattice.
almost constant, independent of the proximity to a lattice point. The crudeness of
theinequalityfollowing(14)isillustratedinFigure1,wherethe(cid:12)llfactorc((cid:3);(cid:26)=2)
isshownforthescenarioofminimalj(cid:3)j=(cid:26)g,i.e.,asquarelattice. Thenc((cid:3);(cid:26)=2)=
(cid:25)g=2=(2g(cid:0)1g(cid:0)(g=2)),whichdecaysfastas g !1. Figure 1showsthatfora square
lattice,notusing the (cid:12)ll factorresultsin the errorbeing overestimatedbyonly one
digitforg =6. Usingthe(cid:12)llfactorforgenericlattices,where(cid:3)>(cid:26)g,usuallyleads
to a more signi(cid:12)cant improvement.
2(cid:25)g=2 ((cid:26)=2)g
Thusc((cid:3);(cid:26)=2)=VolBg (0)=j(cid:3)j= ,wherej(cid:3)jisthedeterminant
(cid:26)=2 g(cid:0)(g=2) j(cid:3)j
ofthe lattice. It is the volumeofa cellof(cid:3) spannedby asetofgeneratingvectors.
This gives
X
jjyjjpe(cid:0)jjyjj2
y2(cid:3);jjyjj(cid:21)R
(cid:18) (cid:19) Z
(cid:25)< VolS1g(cid:0)1(0) 2 g 2(cid:25)g=2 ((cid:26)=2)g 1 e(cid:0)r2rg(cid:0)1+pdr
VolBg(0) (cid:26) g(cid:0)(g=2) j(cid:3)j
1 Z R(cid:0)(cid:26)=2
2(cid:25)g=2 1 1 1
=g e(cid:0)tt(g+p)=2(cid:0)1dt
g(cid:0)(g=2)j(cid:3)j2
(cid:18) (R(cid:0)(cid:26)=2)2 (cid:19)
(cid:25)g=2 g+p
(15) = (cid:0) ;(R(cid:0)(cid:26)=2)2 :
j(cid:3)j(cid:0)(g=2) 2
Using (13) with p=0, (11) becomes
(cid:18) (cid:19)
(16) (cid:15)(R)(cid:20) X e(cid:0)jjv(n)jj2 (cid:20) g 2 g(cid:0)(cid:16)g;(R(cid:0)(cid:26)=2)2(cid:17);
2 (cid:26) 2
(cid:3)nSR
wherep(cid:26) is the lengthofthe shortestlattice vectorin(cid:3): (cid:26)=minfjjxjj;x2(cid:3)g, and
R>( 2g+(cid:26))=2(tosatisfythesubharmonicityconditionofLemma1). Asclaimed,
thisprovesthattheRiemannthetafunctioniswellde(cid:12)ned,sincelimR!1(cid:0)(z;R)=
0 [1, eq. 6.5.3].
1424 B.DECONINCK,M.HEIL,A.BOBENKO,M.VANHOEIJ,ANDM.SCHMIES
Thus, in order to approximate the oscillatory part of a Riemann theta function
of g complex variables with a predetermined error (cid:15), one solves the equation
(cid:18) (cid:19)
(cid:16) (cid:17)
g 2 g g
(17) (cid:15)= (cid:0) ;(R(cid:0)(cid:26)=2)2 ;
2 (cid:26) 2
p p
for R, with R > ( 2g+(cid:26))=2, real. If no solution exists, then R = ( 2g +(cid:26))=2
su(cid:14)ces. These results combine to give
Theorem 2 (Pointwise approximation). The Riemann theta function is approxi-
mated by
(18)
X (cid:16) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:16) h i(cid:17) (cid:17)
(cid:18)( j(cid:10))(cid:25)e(cid:25) (cid:1) (cid:0)1(cid:1) e2(cid:25)i 12 (cid:0) (cid:0)1 (cid:1) (cid:1) (cid:0) (cid:0)1 + (cid:0) (cid:0)1 (cid:1) e(cid:0)jj ( )jj2;
yY y n Y y X n Y y n Y y x v n
z
SR
(cid:8) (cid:12) (cid:9)
with absolute error (cid:15) on the sum. Here S = v(n)2(cid:3)(cid:12) jjv(n)jj<R , (cid:3) =
(cid:8)p (cid:9) R
(cid:25)T(n+[[Y(cid:0)1y]]) jn2Zg . The shortest distance between any two points of
p
(cid:3) is denoted by (cid:26). Then the radius R is dete(cid:0)rmined as the g(cid:1)reater of ( 2g+(cid:26))=2
and the real positive solution of (cid:15)=g2g(cid:0)1(cid:0) g=2;(R(cid:0)(cid:26)=2)2 =(cid:26)g.
This theorem gives a pointwise approximation to the Riemann theta function:
for every z at which the Riemann theta function is approximated, S is di(cid:11)erent,
R
although R is unchanged as long as g and (cid:15) remain the same. The error used in
Theorem 2 is referred to as the 100% Error (100%E).
pAnother good estimate for R is obtained from (15). Then R is the greater of
( 2g+(cid:26))=2 and the real positive solution of
(cid:0) (cid:1)
(cid:15)=(cid:25)g=2(cid:0) g=2;(R(cid:0)(cid:26)=2)2 =(j(cid:3)j(cid:0)(g=2)):
This error is referred to as the Fill Factor Error (FFE).
Figure 2 compares the size R of the ellipsoid, required to obtain a prescribed
error(cid:15), as a result of using the 100%E or the FFE, for the cases g =2 and g =16.
In these (cid:12)gures, the worst-case scenario j(cid:3)j = (cid:26)g was assumed. As expected, the
di(cid:11)erence between the 100%E and the FFE is small for small genus but becomes
signi(cid:12)cant for larger genus.
Two examples are shown in Table 1, for various values of (cid:15). Both of these use
the Riemann matrix (cid:10) with (cid:10) = i, j = 1;:::;g, and (cid:10) = (cid:0)1=2, j 6= k. The
g jj jk
(cid:12)rstexample(leftsideofthe table)correspondstog =2,thesecond(rightside)to
g =6. For both examples, the Riemann theta function is evaluated at z =0, thus
in this case the oscillatory partof the Riemann theta function equals the Riemann
theta function. The value of the Riemann theta function is given in the (cid:12)rst line.
The (cid:12)rstcolumnofanexample gives(cid:15),used forthe pointwiseapproximation,with
either the 100%E or the FFE. In the table N(100%E) and N(FFE) denote the
number of terms that are used to compute the oscillatory part of the Riemann
theta function, i.e., the number of elements of S , using the 100%E or the FFE,
R
respectively. AE denotes Actual Error,i.e., the di(cid:11)erence betweenthe actualvalue
of the oscillatory part of the Riemann theta function (computed using 30 digits
of accuracy and (cid:15) = 10(cid:0)30) and the computed value, using the two di(cid:11)erent error
formulae.
COMPUTING RIEMANN THETA FUNCTIONS 1425
5
4
R
3
2
0 2 4 6 8 10
(a) g =2
8
R 6
4
0 2 4 6 8 10
(b) g =16
Figure 2. R,assumingthe100%Error(solidline)andFillFactor
Error(dashedline)asfunctions of(cid:0)log ((cid:15)), whichis the number
10
of accurate digits.
Using (3), we see that the (cid:12)rst example, with g =2, computes
X1
cos(n1n2(cid:25))e(cid:0)(cid:25)(n21+n22):
n1;n2=(cid:0)1
Since for this example g = 2, there should be little di(cid:11)erence between the compu-
tation using the 100%E or the FFE. This is con(cid:12)rmed: although the value of R is
slightly di(cid:11)erent for both computations, the number of elements of S is identical,
R
resulting in both computations having the same accuracy.
For the second example, with g = 6, there is a di(cid:11)erence between the two
computations, although it is not as outspoken as for g =16 in Figure 2. It is clear
that the Actual Error using the FFE computation is closer to (cid:15) than using the
100%E.Eveninthiscase,the computationdoesseveralordersofmagnitudebetter
thanprescribed. This,ofcourse,is due to the inequalitiesthatwereusedto obtain
these error estimates.
1426 B.DECONINCK,M.HEIL,A.BOBENKO,M.VANHOEIJ,ANDM.SCHMIES
Table 1. Two examples of using the pointwise approximation to
compute the oscillatory part of a Riemann theta function.
g=2 g=6
(cid:18)(0;0j(cid:10)2)e(cid:0)(cid:25)y(cid:1)Y(cid:0)21y =1:1654010572 (cid:18)(0;0;0;0;0;0j(cid:10)6)e(cid:0)(cid:25)y(cid:1)Y(cid:0)61y =1:3945305616
(cid:15) N(100%E) AE(100%FE) (cid:15) N(100%E) AE(100%FE) N(FFE) AE(FFE)
=N(FFE) =AE(FFE)
1E-1 5 7.5E-3 1E-1 485 4.3E-5 233 9.2E-4
1E-2 9 1.5E-5 1E-2 797 3.8E-6 485 4.3E-5
1E-3 13 1.2E-6 1E-3 1341 2.6E-7 797 3.8E-6
1E-4 21 5.1E-11 1E-4 2301 1.3E-8 1341 2.6E-7
1E-5 21 5.1E-11 1E-5 3321 4.2E-10 2301 1.3E-8
1E-6 21 5.1E-11 1E-6 4197 3.8E-11 3321 4.2E-10
1E-7 21 5.1E-11 1E-7 5757 2.3E-12 4197 3.8E-11
1E-8 25 1.9E-12 1E-8 8157 9.2E-14 5757 2.3E-12
1E-9 29 1.8E-13 1E-9 10237 3.5E-15 8157 9.2E-14
1E-10 37 1.5E-17 1E-10 12277 2.7E-16 10237 3.5E-15
Remarks. (cid:15) ThemethodforapproximatingtheoscillatorypartofaRiemann
theta function requires the determination of the elements of S . This is
R
the set of all elements of (cid:3) that lie inside the g-dimensional sphere de-
termined by jjv(n)jj = R. It is easier to consider the equivalent problem
of determining the points with integer coordinates that lie inside the g-
dimensional ellipsoid determined by (cid:25)(n(cid:0)c)(cid:1)Y (cid:1)(n(cid:0)c) = R2, where c
is the center of the ellipsoid. This is done recursively, by remarking that
every (g (cid:0) 1)-dimensional plane section of a g-dimensional ellipsoid is a
(g(cid:0)1)-dimensional ellipsoid. A (cid:12)nite number (due to the discrete nature
of the lattice) of such (g(cid:0)1)-dimensional sections are taken. All of these
are(g(cid:0)1)-dimensionalellipsoids,onwhichthissectionprocessisrepeated,
until one arrives at a one-dimensional ellipsoid, i.e., a line piece. At this
level it is easy to determine which integer points are in this piece. This
method takes full advantageof the triangularformofthe Choleskydecom-
position of Y =TTT. This is done as follows: all integer points satisfying
p
(19) kT(n(cid:0)c)k<R= (cid:25)
arepsought for. Since T is upper triangular, this implies jTgg(ng (cid:0)cg)j <
R= (cid:25), or
R R
c (cid:0) p <n <c + p :
g g g
(cid:25)T (cid:25)T
gg gg
This gives a set of allowed values of n . For each one of these, write
g
n=(n^;n )T, c=(c^;c )T. With
g g
(cid:18) (cid:19)
T^ ^t
T = ;
0 T
gg
(19) becomes kT^(n^ (cid:0)c^)+^t(n (cid:0)c )k2+T2 (n (cid:0)c )2 < R2=(cid:25), which is
g g gg g g
rewritten as
q
kT^(n^ (cid:0)c^+T^(cid:0)1^t(n (cid:0)c ))k< R2=(cid:25)(cid:0)T2 (n (cid:0)c )2:
g g gg g g
This is a similar problem to (19) but in dimension g(cid:0)1 instead of g and
with a di(cid:11)erent R, c and T. The procedure is now repeated.