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 approximation. (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 considered[14]. Suchthetafunctionswithcharacteristicsareuptoanexponential factor Riemann theta functions (1) evaluated at a shifted argument. Thus the computation of the Riemann theta function (1) also allowsthe computation 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 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

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

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