Computing Riemann matrices of algebraic curves

Bernard Deconinck and Mark van Hoeij

November 30, 1999

Abstract

A black-box program for the explicit calculation of Riemann matrices of arbitrary compact connected Riemann surfaces is presented. All such Riemann surfaces are represented as plane algebraic curves. These algebraic curves are allowed to have arbitrary singularities. The method of calculation of the Riemann matrix is essentially its definition: we numerically integrate the holomorphic differentials of the Riemann surface over the cycles of a canonical basis of the homology of the Riemann surface. Both the holomorphic differentials and the canonical basis of the homology of the Riemann surface are obtained exactly through symbolic calculations. This program is included in MapleV.6, as part of the algcurves package.

1 Motivation

Integrable partial differential equations such as the Korteweg-deVries (KdV) equation and the Nonlinear Schrödinger (NLS) equation have been widely used in the last thirty years for the description of various physical phenomena. Especially the soliton solutions of these equations have claimed a well-deserved niche in both theory and experiment. It is however fair to say that their popularity has not been shared with their periodic and quasiperiodic counterparts.

There are many issues one could point at to explain this discrepancy, but the main one must be that the theory of the periodic and quasiperiodic solutions of these equations invariably is connected with algebraic geometry and the theory of Riemann surfaces [2, 8]. If one manages to obtain explicit formulas for the solutions, they usually involve Riemann theta functions, parametrized by some Riemann surface. This lack of explicitness should be compared with the usually much simpler explicit formulas one obtains for one- and twosoliton solutions, in terms of exponential and rational functions.

A Riemann theta function with g phases is given by

ϑ(ϕ_{1}, …, ϕ_{g}|B) = ∑_{m∈ℤ^{g}} exp (2πi (\frac{1}{2}m · B · m + m · ϕ)) ,  (1)

with ϕ = (ϕ_{1}, …, ϕ_{g}). Here B is a g × g Riemann matrix, i.e., a symmetric matrix whose imaginary part is positive definite. Not all such theta functions appear in the solutions of integrable partial differential equations. Only those for which the Riemann matrix originates from a Riemann surface (as described in the next section; this is a much smaller set) result in solutions. Essentially this means that the solutions given in terms of theta functions are implicitly parametrized in terms of Riemann surfaces. As stated, the Riemann theta function is very explicit. The non-explicitness arises because of the complicated nature of the map from the Riemann surface to the Riemann matrix, which is described in the next section.

This paper presents a program to turn this map into an effective computation. The paper is not intended for the specialist in Riemann surfaces. Rather, it is aimed at the interested applied mathematician or physicist. For that reason, terminology is kept to a minimum and hopefully all concepts are made sufficiently clear. The details of the program, including the technical mathematical issues one encounters, are presented elsewhere [7].

2 Terminology and problem formulation

In this section, the required ingredients from the theory of Riemann surfaces are introduced. More details can be found in the standard references [10, 22, 23, 24, 25] and in the review paper [8]. An excellent place to read up on Riemann surfaces and how they relate to plane algebraic curves is the book by Brieskorn and Knörrer [5].

Consider a plane algebraic curve, defined over the complex numbers ℂ, i.e., consider the subset of ℂ^{2} consisting of all points (x, y) satisfying a polynomial relation in two variables x and y with complex coefficients:

F (x, y) = a_{n}(x)y^{n} + a_{n−1}(x)y^{n−1} + … + a_{1}(x)y + a_{0}(x) = 0.  (2)

Here a_{j}(x), j = 0, …, n are polynomials in x. Write a_{j}(x) = ∑_{i} a_{ij}x^{i}, where the coefficients a_{ij} are complex numbers. Assuming a_{n}(x) ≢ 0, n is the degree of F (x, y) considered as a polynomial in y. We only consider irreducible algebraic curves, so F (x, y) cannot be written as the product of two nonconstant polynomials with complex coefficients.

Let d denote the degree of F (x, y) as a polynomial in x and y, i.e. d is the largest i + j for which the coefficient a_{ij} of x^{i}y^{j} in F (x, y) is non-zero. The behavior at infinity of the algebraic curve is examined by homogenizing F (x, y) = 0 by introducing an auxiliary variable z: then F (x, y, z) = z^{d}F (x/z, y/z) = 0 is a homogeneous polynomial equation of degree d. Finite points (x/z, y/z) ∈ ℂ^{2} on Γ correspond to triples (x : y : z) with z ≠ 0. Since for these points (x : y : z) = (x/z : y/z : 1), we may take z = 1, so finite points can be denoted by (x, y) instead of (x/z, y/z). Points at infinity correspond to triples (x : y : z), with z = 0, hence, at a point at infinity, at least one of the two coordinate functions x or y (which are the functions x/z and y/z in homogeneous coordinates) is infinite. Because F (x, y, 0) is a homogeneous polynomial of degree d, there are at most d points at infinity. In what follows, Γ is used to denote the set of points on the algebraic curve given by F (x, y) = 0, consisting of both finite points (x, y), as well as points (x, y, 0) at infinity.

The algebraic curve can have singular points. Finite singular points on the algebraic curve satisfy F (x, y) = ∂_{x}F (x, y) = ∂_{y}F (x, y) = 0. Points at infinity can also be singular. Singular points at infinity satisfy ∂_{x}F (x, y, z) = ∂_{y}F (x, y, z) = ∂_{z}F (x, y, z) = 0 (then also F (x, y, z) = 0, by Euler’s theorem for homogeneous functions). Desingularizing Γ results in a Riemann surface, i.e., a one-dimensional complex-analytic manifold (so it is two-dimensional over the real numbers; it is a surface). There are various ways of desingularizing algebraic curves. The method used here is explained in Section 5. Each nonsingular point on Γ corresponds to one point on the Riemann surface, whereas a singular point on the algebraic curve can correspond to more points on the Riemann surface. Therefore, as long as singular points are avoided, as for instance in Section 3, it is convenient to identify Γ with the Riemann surface. All Riemann surfaces obtained this way are connected (because F (x, y) is irreducible) and compact (because the points at infinity are included). Conversely, it is known [2, 25] that every compact connected Riemann surface can be obtained this way. In what follows, all Riemann surfaces considered are understood to be connected and compact.

An algebraic curve (2) defines an n-sheeted algebraic covering y(x) of the extended complex x-plane. For all but a finite number of values of x in the extended complex plane ℂ ⋃{∞} there are n values for y(x) in ℂ ⋃{∞}. A value for x corresponds to a singularity or a branch point if and only if there are fewer than n values for y(x). A branch point of this n-sheeted covering is defined as an x-value x = b where y(x) does not return to its original value when one analytically continues y(x) along a small circle around x = b. Labeling the sheets of the algebraic covering y(x) with numbers 1, 2, …, n defines a permutation σ_{b}, acting on the sheet labels, giving how they interchange when y(x) is analytically continued counterclockwise around the branch point x = b. The collection of such permutations around all branch points determines the monodromy group of the algebraic curve (2). More details are found in [9] or [26]. In Section 3, the ingredients for a program for the calculation^{1} of the above permutations and hence of a representation of the monodromy group are presented.

A Riemann surface is a one-dimensional complex-analytic manifold. Hence it is a twodimensional manifold over the real numbers and its only topological invariant is its genus g, since the surface is connected. The genus of the Riemann surface represented by the algebraic covering is easily found from the Hurwitz formula [10, 25]: let e_{P} denote the ramification index of y(x) at a point P on the Riemann surface, which is defined as the number of sheets of the algebraic covering y(x) that meet at P. The total branching number W of y(x) is then defined to be the sum of e_{P} −1, taken over all branch points P on the Riemann surface. Then the Hurwitz formula gives

g = \frac{W}{2} − n + 1.  (3)

1Throughout this paper, “calculation” is used when exact results are obtained, whereas “computation”

is used for numerical results.

Define a cycle to be a closed, oriented, smooth or piecewise smooth curve. For the purposes of integrating holomorphic functions in the complex plane, there are no nontrivial cycles, because Cauchy’s theorem allows the deformation of such cycles to points, reducing the integral around any cycle of a holomorphic function to zero. On the other hand, a Riemann surface of genus g is topologically equivalent to a sphere with g handles; it has g holes. On such a surface there are many cycles which cannot be deformed to points. A cycle going around a handle or a cycle encircling a hole cannot be deformed to a point. If two cycles can be continuously deformed into one another, they are homotopic. The concept of homology is weaker than the concept of homotopy. A cycle is homologous to zero if it bounds a piece of the surface Γ. The sum of two cycles is defined to be the cycle obtained by traversing one cycle after the other following the prescribed orientation (hence a cycle can have multiple components). The negative of a cycle is the cycle traversed in the opposite direction. Finally, two cycles are homologous to each other if their difference is homologous to zero.

Many of the cycles around holes and handles are deformable into each other, but one easily sees that there are 2g cycles which can not be deformed into each other and which are not homologous. Fig. 1 illustrates this for the case of a g = 2 surface.

Figure 1: A genus 2 surface with a canonical basis of cycles.

The intersection index of two cycles counts the number of intersections, taking the orientation of the cycles into account. If the cycles a and b do not intersect, their intersection index a ∘ b is zero; if the cycles a and b intersect once, a ∘ b = 1 if t_{a} × t_{b} points out of the surface, where t_{a} (t_{b}) is the tangent vector to a (b) at the intersection point. If t_{a} × t_{b} points into the surface, a ∘ b = −1. For more details, see [19].

On a Riemann surface Γ of genus g, it is possible to choose the 2g nonhomologous cycles such that their intersection indices are as follows:

a_{i} ∘ a_{j} = 0, b_{i} ∘ b_{j} = 0, a_{i} ∘ b_{j} = δ_{ij}, i, j = 0, …, g,  (4)

where δ_{ij} is the Krönecker delta: it is one if i = j and zero otherwise. A basis for the homology of the Riemann surface Γ with these intersection indices is called a canonical basis of cycles. A canonical basis for the homology of a genus 2 surface is drawn in Fig. 1. Notice that a canonical basis for the homology is not uniquely determined. In 1984, C. L. Tretkoff and M. D. Tretkoff [26] published an algorithm for the explicit construction of a canonical basis of cycles of the homology of a Riemann surface. The input of the algorithm is a monodromy representation of the Riemann surface. The implementation of this algorithm is discussed in Section 4.

A differential on a Riemann surface is a one-form on the Riemann surface. It is globally defined and can be integrated over paths on the Riemann surface. A holomorphic differential is a differential with no poles. Consider the case of the extended complex x-plane, i.e., the Riemann sphere. Let f(x) be a non-zero meromorphic function on the Riemann sphere. Since meromorphic functions on the Riemann sphere are rational functions, f(x) = p_{n}(x)/q_{m}(x), with p_{n}(x) and q_{m}(x) polynomials in x of degrees n and m respectively. Then ω = f(x) dx is by definition a meromorphic differential on the Riemann sphere, i.e., a differential whose only singularities are poles. This differential has poles at the roots of q_{m}(x). Hence, in order for ω to be holomorphic, q_{m}(x) is constant. Without loss of generality, let q_{m}(x) = 1. Furthermore, using τ = 1/x as a local parameter at x = ∞, we can write ω = −p_{n}(1/τ) dτ/τ^{2} at infinity. Hence ω has a pole at infinity unless p_{n}(1/τ) has at least a double root at τ = 0, but this is impossible. Hence, on the Riemann sphere no non-zero holomorphic differentials exist. However, for genus greater than zero, the situation is different; non-zero holomorphic differentials do exist, and furthermore cycles that are not homologous to zero exist as well. In fact, a cycle is homologous to zero if and only if the integral of every holomorphic differential along that cycle is zero. The set of all holomorphic differentials is a ℂ-vector space of dimension g, where g is the genus. A basis for the cohomology of the Riemann surface is a set of g linearly independent holomorphic differentials. Its calculation is discussed in Section 5.

The values of the integrals of the holomorphic differentials along the cycles of the homology are closely related to the geometry and analytic structure of the Riemann surface under consideration [22]. The periods of the holomorphic differentials are defined as the value of their integrals around a canonical basis of the homology. Given a canonical basis of both the homology {a_{i}, b_{i}, i = 1, …, g} and the cohomology {ω_{i}, i = 1, …, g}, a period matrix Ω of the Riemann surface Γ is

Ω = (A B) ,  (5)

which is a g × 2g matrix, consisting of two g × g blocks:

A = (A_{ij})_{i,j=1}^{g} , A_{ij} = ∮_{a_{j}} ω_{i},  (6)

B = (B_{ij})_{i,j=1}^{g} , B_{ij} = ∮_{b_{j}} ω_{i}.  (7)

As it is defined, the period matrix Ω depends on the choices of the bases of both the homology and the cohomology. Choosing different bases results in an equivalent period matrix for the same Riemann surface. Up to this equivalence, the period matrix completely determines the Riemann surface (Torelli’s theorem, see [12]). Often one chooses a specific basis of the cohomology such that

∮_{a_{j}} ω̂_{i} = δ_{ij}.  (8)

With this basis of the cohomology , A ≡ I, the g × g identity matrix. Then B is called a Riemann matrix. One easily sees that from any period matrix a Riemann matrix B̂ for the Riemann surface is determined by

B̂ = A^{−1}B.  (9)

The Riemann matrix by definition does not depend on the choice of a basis for the cohomology. However, it still depends on the choice of the canonical basis of the homology. The Riemann matrix is symmetric and the eigenvalues of its imaginary part are positive definite. These properties are not obvious from the definition. They are consequences of the Riemann conditions on the periods of differentials on Riemann surfaces. See [8] or [25] for proofs. The computation of a period matrix is discussed in Section 6. From there, (9) easily gives a Riemann matrix for the Riemann surface Γ. This Riemann matrix now defines a Riemann theta function, using (1).

Notice that all concepts above are introduced using existence statements. The way these concepts were introduced here mirrors the standard texts, such as [10, 25] or [8].

The mathematics in this paper is not new. Methods for the homology, differentials, their integrals, and hence the period matrix are known in the literature. However, to the best of our knowledge a complete implementation of all these ingredients had not been accomplished. The black-box program implemented by the authors for the computation of the Riemann matrix of a Riemann surface is written in Maple and it is included in MapleV.6, as part of the algcurves package, as are the other programs that are presented.

Some remarks are in order:

• Other methods for the computation of Riemann matrices are found in the literature, see [21, 11] and Chapter 5 of [2] and the references therein. Although the final output of these methods is identical to ours, the point of view is very different. The input of our program is a plane irreducible algebraic curve with complex coefficients. This is a representation of a Riemann surface which arises naturally in the theory of integrable partial differential equations. See, for instance, Section 4.6 of [2] or [6, 14]. The homology and the cohomology of the Riemann surface are then inferred from this representation.

• Allowing singular plane algebraic curves is important for the theory of integrable partial differential equations: in [6], a representation of the Riemann surface for solutions of the Kadomtsev-Petviashvili (KP) equation is obtained which is generically singular. Not allowing these singularities [20], results in the omission of many solutions which can be expressed in terms of Riemann theta functions.

• Riemann surfaces corresponding to an algebraic curve of the form y^{2} = P (x), with P (x) a polynomial in x, are called hyperelliptic. The class of hyperelliptic Riemann surfaces is an important subclass of Riemann surfaces because it is a large class for which bases for both the homology and the cohomology can be written down explicitly. It should be stressed that the program discussed here is in no way restricted to hyperelliptic Riemann surfaces. Arbitrary irreducible plane algebraic curves (2) are allowed. The genus of the Riemann surface Γ and the covering degree n are not restricted. Due to the complexity of the algorithm, practical (machine-dependent) limits on n or g do, of course, exist.

In the theory of the KdV and the NLS equations only hyperelliptic Riemann surfaces arise. Nevertheless, it is important from the point of view of integrable differential equations to consider more general Riemann surfaces, since they arise in the theory of other integrable equations, such as the Boussinesq equation, where cubic surfaces appear. In the theory of periodic and quasiperiodic solutions of the Kadomtsev-Petviashvili (KP) equations, arbitrary Riemann surfaces appear [13].

• For now, the coefficients of (2) are allowed to be algebraic numbers. These coefficients must be given exactly, floating point coefficients are at the moment (this may change in later implementations) not allowed. The reason for not allowing floating point numbers is that the geometry of the Riemann surface is highly dependent on the accuracy of the coefficients in (2): if an algebraic curve has singularities, than almost surely, the nature of this singularity will be affected by inaccuracies in the coefficients of the curve. This may affect the genus, homology and cohomology of the Riemann surface. For the KdV equation, it is known that periodic solutions are stable with respect to small changes in the Riemann surface [17]. Hence, although the geometry of the Riemann surface may be severely affected by small changes in the coefficients, this ultimately has only a small effect on the solution of the KdV equation. A similar statement is known to hold for the KP2 equation [15].

Users can consider floating point coefficients, but these need to be converted to a different form (rational, for instance), before the program will accept the input.

• The program treats some classes of Riemann surfaces differently: shortcuts are used for the case of algebraic curves with real coefficients, for hyperelliptic Riemann surfaces given as a quadratic polynomial in y, etc. Also, throughout the program, many checks for the correctness of the answer are executed. These details will be given in [7].

3 Calculation of the monodromy group of a plane algebraic curve

The monodromy group of an algebraic covering y(x) requires several ingredients. One starts by selecting a base point x = a in the complex x-plane. This base point is not allowed to be a discriminant point of the algebraic covering. In other words, for x = a, n distinct y-values exist. These n y-values are now assigned an order, (y_{1}, y_{2}, …, y_{n}). This ordering of the n y-values in effect labels the sheets of the algebraic covering y(x). For each branch point b one chooses a path γ_{b} in the complex x-plane which starts and ends at x = a and encircles only branch point x = b, counterclockwise. The n-tuple (y_{1}, y_{2}, …, y_{n}) is then analytically continued around this path γ_{b}. When one returns to x = a, a new n-tuple is found, which has the same entries as (y_{1}, y_{2}, …, y_{n}), but ordered differently: (y_{σ_{b}(1)}, y_{σ_{b}(2)}, …, y_{σ_{b}(n)}), where σ_{b} is a permutation acting on the set of labels {1, 2, …, n}. If x = b is indeed a branch point, then by definition σ_{b} is not the identity.

We now discuss the steps of this procedure more systematically.

1. The problem points of the analytic continuation: To avoid difficulties during the analytical continuation of (y_{1}, y_{2}, …, y_{n}), the paths γ_{b} should avoid passing through certain points, which will be referred to as problem points {b_{1}, b_{2}, …, b_{m}}. The roots of the polynomial a_{n}(x) are elements of this set. As before, a_{n}(x) is the coefficient of y^{n} in (2). The reason for including these roots in the set of problem points (i.e., excluding these points from the paths) is that all y-values are finite if and only if x is not a root of a_{n}(x).

The remaining problem points are the discriminant points of the algebraic covering y(x), which are the x-values for which y(x) has fewer than n elements [27]. These points are called discriminant points because they are roots of the discriminant of F (x, y) with respect to y, which is a polynomial in x. All branch points of an algebraic covering are discriminant points, but not all discriminant points are branch points. Discriminant points which are not branch points are singular points. Branch points can also be singular points.

In principle, only the knowledge of the branch points is required for the purposes of computing the permutations. However, to avoid numerical problems during the analytical continuation of (y_{1}, y_{2}, …, y_{n}), the other problem points, such as singularities, are avoided as well. The monodromy program treats them on the same footing as the branch points, although they are removed from the output of the program, since their permutations σ_{b} are all equal to the identity.

2. The choice of the base point: In order to compute the monodromy group of the algebraic covering y(x), y(x) is analyticly continued along paths around the problem points of the analytic continuation. For numerical accuracy, it is advantageous to keep some distance from the problem points. To this end, with every problem point b_{i}, a radius r(b_{i}) is associated: (ρ denotes the distance)

r(b_{i}) = \frac{2}{5}ρ(b_{i}, \{b_{1}, b_{2}, …, b_{m}\} − \{b_{i}\}),  (10)

or, r(b_{i}) is two fifths of the distance of b_{i} to the next-nearest problem point. The number 2/5 is somewhat arbitrary; other numbers between 0 and 1/2 might be used. Let C(b_{i}, r(b_{i})) denote the circle with center b_{i} and radius r(b_{i}). Then the circles C(b_{i}, r(b_{i})) do not intersect each other.

Now a base point a is chosen, such that the real part of a is smaller than the real parts of any of the b_{i}. By this choice, the arguments of b_{i} − a are between −π/2 and π/2.

3. The labeling of the sheets: At the base point x = a there are n distinct finite y-values. These are determined numerically as the solutions of F (a, y) = 0. Let these n y-values be assigned an order (y_{1}, y_{2}, …, y_{n}), which is denoted as y(a). Assigning such an order to these y-values labels the sheets of the algebraic covering y(x): sheet one is the sheet containing y_{1}, sheet 2 is the sheet containing y_{2}, and so on.

4. The ordering of the problem points of the analytic continuation: An ordering needs to be imposed on the problem points. We order these points according to their argument with respect to the base point: if arg(b_{i} −a) < arg(b_{j} −a), then b_{i} precedes b_{j} in the ordering, where arg(·) denotes the argument function. If arg(b_{i}−a) = arg(b_{j} −a), then b_{i} precedes b_{j} in the ordering if |b_{i} − a| < |b_{j} − a|. This ordering results in an mtuple of problem points: (b_{1}, b_{2}, …, b_{m}). For convenience, the same notation as above is used for the ordered problem points as for the elements of the non-ordered set.

5. The choice of the paths: Now paths are chosen for the analytic continuation. The paths chosen are composed of line segments and semi-circles. The simplest path L(b_{i}) around b_{i} consists of one line segment from a to b_{i}−r(b_{i}). This is followed by C(b_{i}, r(b_{i})), starting at b_{i} − r(b_{i}). Successively, a line segment is followed from b_{i} − r(b_{i}), back to a. However, in many cases, this path will intersect one of the circles C(b_{j}, r(b_{j})), j ≠ i. This indicates that the path comes close to the problem point b_{j}. To avoid accuracy issues during the analytic continuation (see below), this should be avoided. This is remedied as indicated in Fig. 2: the path takes a detour along a semi-circle around b_{j}. Whether this semi-circle goes above or below b_{j} depends on the relative positions of a, b_{i} and b_{j}. The semi-circle is chosen such that the new path is deformable to L(b_{i}), without crossing any problem points of the analytic continuation.

(a)  (b)  (c)

Figure 2: Choosing the path from x = a to x = b_{i}. The path around b_{i} is indicated in a thick black line. (a) The simplest path intersects C(b_{s}, r(b_{s})). (b) This is remedied by a new path which is a deformation of the previous one. The new path intersects C(b_{j}, r(b_{j})). (c) This is remedied by another path, which is a deformation of both previous paths

This process is interated, until a path is obtained, which stays at least r(b_{j}) away from b_{j}, for j = 1, 2, …, m. In itself, the iteration of this process is not sufficient to ensure that the chosen path is deformable to the straight-line path from the base point a to b_{i} − r(b_{i}). To ensure that a correct path is chosen, the program explicitly checks for the presence of problem points between the chosen path and the straight-line path. If such points are present, the path is modified to go around them, after which the check procedure is re-iterated.

6. The analytic continuation: Let two non-problem points x = x_{1} and x = x_{2} be given. Corresponding to x_{1} is an ordered n-tuple y(x_{1}). When a path is followed in the complex x-plane from x_{1} to x_{2}, the entries of y(x_{1}) follow paths on the Riemann surface to the roots of F (x_{2}, y) = 0, which gives rise to an n-tuple y(x_{2}), whose ordering is induced by the ordering of y(x_{1}). If the path between them deviates little from a straight-line segment and does not pass through or near any problem points, then

y(x_{2}) = y(x_{1}) + y′(x_{1})(x_{2} − x_{1}) + ϑ(|x_{2} − x_{1}|^{2}),  (11)

and the last term is small when x_{2} and x_{1} are suffiently close (to make this precise one needs to bound the second derivative of y(x) to find a bound for ϑ(|x_{2} − x_{1}|^{2})). Here y′(x_{1}) is the n-tuple of derivatives to the algebraic curve at y(x_{1}). By implicit differentiation,

y′(x_{1}) = − (\frac{F_{x}(x_{1}, y_{1}(x_{1}))}{F_{y}(x_{1}, y_{1}(x_{1}))}, \frac{F_{x}(x_{1}, y_{2}(x_{1}))}{F_{y}(x_{1}, y_{2}(x_{1}))}, …, \frac{F_{x}(x_{1}, y_{n}(x_{1}))}{F_{y}(x_{1}, y_{n}(x_{1}))}) ,  (12)

where a subindex x or y denotes partial differentiation and y_{i}(x_{1}), i = 1, …, n denotes the i-th component of y(x_{1}). Hence, under the above conditions, the first two terms of (11) give a good approximation to y(x_{2}). Having the unordered entries of y(x_{2}) and comparing them with the ordered approximation y(x_{1}) + y′(x_{1})(x_{2} − x_{1}) determines the ordering of these entries, resulting in the ordered n-tuple y(x_{2}). Clearly, in order to avoid matching up the entries of y(x_{2}) with the wrong entries of the ordered approximation, the accepted size of |x_{2} − x_{1}| depends on the separation of the entries of y(x_{2}). This process is illustrated in Fig. 3 with n = 3.

Figure 3: The analytic continuation of y(x) from x_{1} to x_{2}, using a linear approximation, with n = 3. Since the path from x_{1} to x_{2} avoids problem points, the paths on the Riemann surface from y(x_{1}) to y(x_{2}) never intersect (even though they appear to intersect in this picture).

If |x_{2} − x_{1}| is not small, or if the path connecting them deviates significantly from a straight-line segment, then an analytic continuation from y(x_{1}) to y(x_{2}) is obtained by iterating the above process along small segments of the path, such that the necessary conditions above are satisfied. Note that y(x_{2}) is dependent on the path chosen from x_{1} to x_{2}. For brevity of notation, this dependence is not made explicit.

7. The monodromy group: Consider a closed path starting from x = a and returning there after encircling one branch point x = b. After analytic continuation of y(a) along this path, the entries of y(a) are recovered, but they are shuffled by the permutation σ_{b}.

The collection of all σ_{b} generates the monodromy group, which is represented here as a subgroup of S_{n}, the group of permutations of {1, 2, …, n}. Note that this representation depends on the choice of the labeling of the y-values at x = a, so it is only unique up to conjugation. More details are found in [9].

This representation of the monodromy group is discrete in character. Because of this, it is obtained exactly, despite the numerical nature of the analytic continuation.

8. Infinity: The point x = ∞ might also be a branch point. The corresponding permutation σ_{∞} does not need to be computed by analytic continuation, since it can be determined from the other σ_{b} using the following relation.

σ_{∞} ∘ σ_{b_{m}} ∘ σ_{b_{m−1}} ∘ … ∘ σ_{b_{2}} ∘ σ_{b_{1}} = 1.  (13)

This states that analytical continuation along a closed path in the extended complex x-plane that encircles all branch points will act as the identity permutation. Such a path is deformable to a point and analytic continuation along this path does not permute the entries of y(a).

The following example computes a representation of the monodromy group for the algebraic covering y(x), corresponding to F (x, y) = y^{3} − x^{7} + 2x^{3}y = 0: the command monodromy(f,x,y) gives a list with three entries. The first entry is the choice of the basepoint x = a. The second entry is y(a), which is a list of n elements. The third entry is a list of the branch points b_{i} with their permutations σ_{b_{i}}, given in disjoint cycle^{2} notation.

>with(algcurves):         # load the algcurves package
>f:=y^3-x^7+2*x^3*y:      # define the algebraic curve
>genus(f,x,y);            # calculate the genus of the algebraic curve
                        2
># calculate the monodromy representation for y(x)
>m:=monodromy(f,x,y,showpaths):
>m[1];                    # the base point x=a
                  -1.44838920232
>m[2];                    # the sheets y(a)
          [-3.20203812254, 1.60101906127 - 1.26997391750 I,
              1.60101906127 + 1.26997391750 I]
>m[3];                    # the branch points with their permutations
[[-.319697769990 - .983928563571 I, [[1, 3]]],
[.836979627962 - .608101294789 I, [[2, 3]]],
[-1.03456371594, [[2, 3]]],
[0, [[1, 3]]],
[.836979627962 + .608101294789 I, [[2, 3]]],
[-.319697769990 + .983928563571 I, [[1, 2]]],
[infinity, [[1, 3, 2]]]]

2We apologize for any confusion that may arise because of the dual use of the term “cycle”. Unfortunately, this term is standard in combinatorics as well as in the theory of Riemann surfaces.

Thus, starting from the ordered sheets 1, 2 and 3 above the base point x = −1.44838920232 (i.e., from the ordered y-values −3.20203812254, 1.60101906127 −1.26997391750i, 1.60101906127 +1.26997391750i), and encircling x = −.319697769990 − .983928563571i, one finds that sheet one has become sheet 3 and sheet 3 has become sheet 1.

The optional argument showpaths produces Fig. 4. This shows the paths followed in the complex x-plane for the analytic continuation of y(x).

Figure 4: The complex x-plane and the paths followed in it for the analytic continuation of y(x), with y^{3} − x^{7} + 2x^{3}y = 0. The base point x = a is at the left of the first circle from the left.

4 Calculation of the homology of a Riemann surface

Choosing a base point x = a results in a choice of paths around the branch points, avoiding all other problem points. The representation of the monodromy calculated above then allows the explicit construction of a canonical basis for the homology for the Riemann surface Γ. The algorithm to obtain such a basis is due to Tretkoff and Tretkoff [26]. Since our program implements this algorithm essentially without changes, the algorithm is not discussed here in much detail. More details are available in [26] or [7]. A short geometrical sketch of the algorithm is as follows:

1. On the Riemann surface Γ, indicate all points with x = a and y contained in y(a). This gives n points A_{i} = (a, y_{i}(a)), i = 1, 2, …, n on the Riemann surface, which label the sheets of Γ.

2. Let x = b denote a branch point in the complex x-plane. Denote one of the disjoint cycles of σ_{b} by τ. Then the sheets labeled by A_{k}, with k ∈ τ meet at a branch point B on the Riemann surface, with x-value b. Similarly, each one of the disjoint cycles of any σ_{b} corresponds one-to-one with a branch point on the Riemann surface, with x-value x = b. In particular, the number of branch points on the Riemann surface with x = b is the number of disjoint cycles in σ_{b}, not including fixed points. The total number of points on the Riemann surface with x = b is the number of disjoint cycles in σ_{b}, including fixed points. In what follows, we denote the disjoint cycle in σ_{b} corresponding to B by σ_{B}.

Now, On the Riemann surface, indicate all branch points, including branch points at infinity. Let t denote the total number of branch points and B_{i}, i = 1, 2, …, t the branch points on the Riemann surface.

3. Now join every branch point B_{i} to each point A_{j}, for which j ∈ σ_{B_{i}}, i = 1, …, t, using paths which only meet at the points B_{i} and A_{j}. Thus, every branch point B_{i} is connected to all points A_{j} which can be reached by paths emanating from B_{i} without passing through other branch points B_{j}, j ≠ i.

This creates a non-directed graph on the Riemann surface, with n + t vertices.

4. Reduce this graph to a spanning tree, by removing a number of edges, say r edges. Denote these edges by e_{i}, i = 1, …, r.

5. This spanning tree contains no closed paths, by definition. Adding to it the removed edge e_{1} gives rise to a unique closed path on the Riemann surface. Fix an orientation on this closed path, thus defining a cycle c_{1}.

6. Similarly, every other removed edge e_{k} gives rise to a closed path on the Riemann surface. If this path has any edges in common with the cycles c_{1}, c_{2}, …, c_{k−1}, then an orientation is induced from these cycles on the cycle c_{k}. Otherwise, an orientation is chosen. This way, a collection of r cycles c_{1}, …, c_{r} is obtained on the Riemann surface.

Tretkoff and Tretkoff [26] show that the cycles constructed above are all nontrivial, i.e., they can not be contracted to a point. Furthermore,

r = 2g + n − 1.  (14)

Since n > 1 (otherwise y(x) is a single-valued function), r > 2g. The above construction hence results in more cycles than are required for a basis of the homology: r − 2g = n − 1 of these cycles are dependent in the homology of the Riemann surface.

In [26], Tretkoff and Tretkoff present an algorithmic way to cut the above-mentioned graph on the Riemann surface across the edges e_{1}, …, e_{r}. This results in a planar graph which contains cut copies of the cycles c_{1}, …, c_{r}. Tretkoff and Tretkoff [26] then show how this planar graph is used to find the intersection numbers K_{ij} = c_{i} ∘ c_{j}, i, j = 1, …, r, resulting in an r × r intersection matrix K = (K_{ij})_{i,j=1}^{g} = (c_{i} ∘ c_{j})_{i,j=1}^{g}. Because only 2g of the cycles c_{1}, …, c_{r} are independent, the rank of this matrix is 2g. Furthermore, an r × r matrix α with integer entries and determinant ±1 exists such that

αKα^{T} = J = ( \begin{matrix}0_{g} & I_{g} & 0_{g,n−1} \\ −I_{g} & 0_{g} & 0_{g,n−1} \\ 0_{n−1,g} & 0_{n−1,g} & 0_{n−1,n−1} \end{matrix} ) ,  (15)

with 0_{g} the g × g zero matrix, I_{g} the g × g identity matrix and 0_{p,q} the p × q zero matrix. Define the cycles

a_{i} = ∑_{j=1}^{r} α_{ij}c_{j}, b_{i} = ∑_{j=1}^{r} α_{i+g,j}c_{j}, i = 1, …, g.  (16)

It is a straightforward calculation to check that these cycles satisfy (4). Hence the cycles a_{1}, …, a_{g}, b_{1}, …, b_{g} define a canonical basis of cycles for the homology of the Riemann surface. The non-uniqueness of such a basis is then a restatement of the non-uniqueness of the matrix α. This matrix is the transformation matrix from the over-complete basis of cycles c_{1}, …, c_{r} to the canonical basis. Its first 2g rows prescribe the linear combination of the cycles c_{1}, …, c_{r} which results in the canonical basis. Its last n − 1 rows confirm the dependence of the cycles c_{1}, …, c_{r}:

∑_{j=1}^{r} α_{ij}c_{j} = 0, i = 2g + 1, …, 2g + n − 1.  (17)

The following example computes a canonical basis for the homology of the Riemann surface corresponding to y^{3} − x^{7} + 2x^{3}y = 0. The command homology(f,x,y) results in a table. This table has the following entries:

• basepoint: the base point for the analytic continuation of y(x).

• sheets: the ordered n-tuple y(a).

• cycles: the cycles c_{1}, c_{2}, …, c_{r}. The cycle c_{k} is given as a list: the first element specifies the starting sheet y_{i}(a), by giving i. The second element is a branch point x = b in the complex x-plane, together with the disjoint cycle of σ_{b}, which contains i. The third element is a sheet y_{j}(a), given by j. This part of the cycle c_{k} is read as: “From sheet i go to sheet j, by encircling x = b”. It is possible that x = b needs to be encircled more than once, in order to get from sheet i to sheet j. Having arrived at sheet j, this process now repeats. The list is cyclical, meaning that after encircling the last branch point, one arrives again at the initial sheet, so as to obtain a cycle on the Riemann surface.

• linearcombination: the first 2g rows of the matrix α.

• canonicalcycles: the result of combining linearcombination and cycles. Each of

the cycles a_{1}, …, a_{g}, b_{1}, …, b_{g} is given as a list of lists. Adding the lists in the list gives a basis element of the canonical basis of cycles. Usually, there is only one list necessary for each canonical-basis element. Since the canonical-basis elements are obtained from the information in both cycles and linearcombination, their representation is usually more complicated. Also, instead of specifying the disjoint cycle of the permutation of the branch point, the number of times one needs to encircle the branch point in the complex x-plane counterclockwise is given. If this number is negative, the branch point needs to be encircled clockwise as many times as the absolute value of the number.

• genus: this entry gives the genus of the Riemann surface, by halving the dimension of the canonical basis. This topological calculation is completely independent of the one using Puiseux expansions, used by the genus program [29].

>with(algcurves):          # load the algcurves package
>f:=y^3-x^7+2*x^3*y:       # define the algebraic curve
># calculate the homology of the Riemann surface corresponding to f.
>h:=homology(f,x,y):
>h[basepoint];             # the base point x=a
                   -1.44838920232
>h[sheets];                # the sheet labels y(a)
   [-3.20203812254, 1.60101906127 - 1.26997391750 I,
       1.60101906127 + 1.26997391750 I]
>eval(h[cycles]);          # the cycles c1,...,cr
table([
6=[1, %2, 3, [.836979627962 + .608101294789 I, [2, 3]], 2, %1]
1=[1, %2, 3, [0, [1, 3]]]
2=[1, %2, 3, [infinity, [1, 3, 2]]]
3=[1, %1, 2, [infinity, [1, 3, 2]]]
4=[1, %2, 3, [.836979627962 - .608101294789 I, [2, 3]], 2, %1]
5=[1, %2, 3, [-1.03456371594, [2, 3]], 2, %1]
   ])
%1:=[-.319697769990 + .983928563571 I, [1, 2]]

  %2:=[-.319697769990 - .983928563571 I, [1, 3]]
  >eval(h[linearcombination]);    # the first 2g rows of the matrix alpha
                      [1    0    0     0    0    0]
                      [                           ]
                      [0    0    1     0    0    0]
                      [                           ]
                      [0    0    0    -1    0    1]
                      [                           ]
                      [1    0    0    -1    0    0]
  >eval(h[canonicalcycles]);       # the canonical-basis cycles
  table([
   b[2]=[[1, [-.319697769990 + .983928563571 I, 1], 2,
         [.836979627962 - .608101294789 I, 1], 3, [0, -1]]]
   b[1]=[[2, [.836979627962 - .608101294789 I, 1], 3,
         [.836979627962 + .608101294789 I, -1]]]
   a[1]=[[1, [-.319697769990 - .983928563571 I, 1], 3, [0, -1]]]
   a[2]=[[1, [-.319697769990 + .983928563571 I, 1], 2, [infinity, 1]]]
        ])

Thus, the cycle c_{1} is as follows: start on sheet 1; encircle branch point x = −.319697769990− .983928563571i to arrive at sheet 3; encircle branch point x = 0 to arrive at sheet 1. Using the linear combination, the cycle b_{1} is given by b_{1} = −c_{4}+c_{6}, which is rewritten as: start from sheet 2; encircle branch point x = .836979627962 − .608101294789i once, counterclockwise, ending up at sheet 3; encircle branch point x = .836979627962 + .608101294789i one time, counterclockwise and end up back at sheet 2.

The program demonstrated above is not the first program to implement Tretkoff and Tretkoff’s algorithm [26]. Such a program in Turbo Pascal was also announced in [3]. More recently, this program was rewritten in C++ [4]. These programs start from a representation of the monodromy group of a Riemann surface and construct from it a canonical basis for the homology. To the best of our knowledge, the Maple program homology, presented here, is the only program that calculates a canonical basis for the homology of a Riemann surface, starting from the equation of a plane algebraic curve.

5 Calculation of the cohomology of a Riemann surface

The cohomology of a Riemann surface is specified by a basis {ω_{1}, …, ω_{g}} of holomorphic differentials on the surface. These ω_{k} are of the form (see [5] or [18])

ω_{k} = \frac{P_{k}(x, y)}{∂_{y}F (x, y)} dx.  (18)

Here P_{k}(x, y) = ∑_{i+j≤d−3} c_{kij}x^{i}y^{j} is a polynomial in x and y of degree at most d − 3, where as before d is the degree of F (x, y) as a polynomial in x and y. Clearly there are at most (d − 1)(d − 2)/2 linearly independent polynomials P_{k}(x, y) of this form. If Γ is nonsingular, then all polynomials P (x, y) of degree ≤ d − 3 give rise to a holomorphic differential ω = P (x, y)/∂_{y}F (x, y) dx. In other words, the genus of a nonsingular plane algebraic curve of degree d is exactly (d − 1)(d − 2)/2 [5].

The denominator ∂_{y}F (x, y) vanishes at the branch points of y(x) as well as at the singular points, whereas the differential dx vanishes only at the branch points of y(x). Therefore, in order for the differential ω_{k} not to have poles at the singular points, the numerator P_{k}(x, y) must vanish at the singular points. Noether [18] showed that at a singular point P on Γ of multiplicity m_{P} the adjoint polynomial P_{k}(x, y) vanishes at P with multiplicity at least m_{P} − 1. Imposing regularity of the differentials (18) at a point P imposes a number of independent linear conditions on the coefficients c_{kij}. The number of such conditions is called the delta invariant δ_{P}. Once the delta invariants of all points (finite and infinite) are known, the genus is given by (d − 1)(d − 2)/2 minus the sum of all the delta invariants. A point P is a singular point if and only if δ_{P} > 0, and this holds if and only if m_{P} > 1. Note that the Maple command singularities finds {m_{P}, δ_{P} } for all singular points P on the algebraic curve.

For every singular point, there are m_{P} (m_{P} − 1)/2 linear conditions which are easily computed. These arise from the fact that P_{k}(x, y) should vanish at P with multiplicity m_{P} − 1. If δ_{P} = m_{p}(m_{P} − 1)/2 this is a sufficient number of linear equations. Otherwise δ_{P} > m_{p}(m_{P} − 1)/2, and more linear equations are required. Singular points P with δ_{P} > m_{P} (m_{P} − 1)/2 will be called special singularities.

In what follows, an outline is given on how to obtain the linear conditions for the case when the curve has finite special singularities. At the end of this section, we discuss what needs to be done if Γ has special singularities at infinity. As before, the Maple command discussed is demonstrated. The example on which the command is illustrated is also discussed in more detail.

1. Puiseux expansions and desingularization: The polynomials P_{k}(x, y) are called adjoint polynomials of degree d − 3. Several methods are known for computing the conditions on their coefficients.

Since these conditions follow from the behavior of the plane algebraic curve at its singular points, it is no surprise that these conditions are computed using a desingularizing of Γ. There are several ways to desingularize Γ and obtain from it a Riemann surface. A well-known method relies on the use of quadratic transformations to “blow-up” or resolve the singularities [1, 5]. Another method uses Puiseux expansions: near every x-value x = x_{0}, one computes c_{i}, r_{i}, η_{i}, i = 1, …, M such that

F (x = x_{0} + c_{i}t^{r_{i}}, y = η_{i}(t)) = 0,  (19)

here η_{i}(t) is meromorphic in the parameter t: η_{i}(t) = ∑_{j=−N}^{∞} d_{ij}t^{j}, for some integer N. The positive integer r_{i} is called the ramification index of the Puiseux expansion x = x_{0} + c_{i}t^{r_{i}}, y = η_{i}(t). Such a Puiseux expansion corresponds to a point on the desingularization of Γ, hence to a point on the Riemann surface. This point is a branch point if and only if r_{i} > 1, and r_{i} is the branching number. The complex number c_{i} can in principle always be chosen to be one, but for computations other choices are often more convenient. Each Puiseux expansion gives a local description of the r_{i} sheets that meet at one point on the Riemann surface with x-value x_{0}.

Having obtained Puiseux expansions of the algebraic curve at all the singular points, a desingularisation of Γ is obtained as the union of the set Γ with all its singular points removed, with the set of Puiseux expansions at the singular points. Imposing a complex-analytic structure on this set makes it a Riemann surface.

So, there are two ways to compute the adjoint polynomials P_{k}(x, y): using quadratic transformations or (probably easier to implement) using Puiseux expansions. A third way to determine the P_{k}(x, y) is described below. This third way is based on computing an integral basis. We chose to implement this method because an implementation of the integral basis algorithm (which uses Puiseux expansions) was already available.

2. An integral basis: Consider the coordinate functions x and y on the Riemann surface. These two functions are algebraically dependent, by F (x, y) = 0. Denote by A the part of the Riemann surface where both x and y are finite. Also, let O_{A} be the set of all meromorphic functions on the Riemann surface^{3} that have no poles in A. For example, O_{A} contains ℂ[x, y], the set of all polynomials of x and y: since in A both x and y are finite, any polynomial of x and y results in a finite value as well. If Γ has no finite singularities, then every meromorphic function on the Riemann surface without poles in A can be represented as a polynomial in x and y, hence O_{A} = ℂ[x, y] if Γ has no finite singularities. In general, O_{A} is the integral closure of ℂ[x, y] in the meromorphic functions on A: it is the set of all meromorphic functions f on A which satisfy a monic polynomial equation f^{m} + c_{m−1}(x, y)f^{m−1} + … + c_{1}(x, y)f + c_{0}(x, y) = 0, for a certain positive integer m and coefficients c_{i}(x, y), i = 0, 1, …, m−1 which are in ℂ[x, y]. Note that m = 1 implies ℂ[x, y] ⊂ O_{A}, so all polynomials in x and y are in O_{A}. An integral basis {β_{1}, …, β_{n}} of O_{A} can be computed such that every element of O_{A} can be written as a linear combination of β_{1}, …, β_{n} with coefficients which are polynomial in x. Here n is the degree of F (x, y) as a polynomial in y. An efficient method to calculate an integral basis of O_{A}, using Puiseux expansions, is given in [28], and is already available in Maple V.5. It turns out, using the results of [16], that the knowledge of an integral basis reduces the problem of finding the adjoint polynomials to linear algebra.

3All meromorphic functions on the extended complex x-plane, i.e. the Riemann sphere with coordinate function x are rational functions in x, i.e., quotients of polynomials in x. More general, on a Riemann surface with coordinate functions x and y, all meromorphic functions are quotients of polynomials in x and y.

3. Calculation of the holomorphic differentials using a theorem of Mñuk: Denote the set of all adjoint polynomials by Adj(Γ). The elements of Adj(Γ) are polynomials P (x, y) ∈ ℂ[x, y] for which the differential ω = P (x, y)/∂_{y}F (x, y) dx has no poles in A. Then, for any element f ∈ O_{A}, fω also has no poles in A. In fact one shows that in this case fP (x, y) is again a polynomial, and hence is in Adj(Γ), see [16]. Denoting by O_{A} · Adj(Γ) the set of products of elements of O_{A} with elements of Adj(Γ), this statement is written as O_{A} · Adj(Γ) ⊂ Adj(Γ) ⊂ ℂ[x, y], since all elements of Adj(Γ) are by definition polynomials in x and y. Theorem 3.3 of [16] shows that this condition determines the adjoint polynomials completely:

Adj(Γ) = {P (x, y) | O_{A} · P (x, y) ⊂ ℂ[x, y]} .  (20)

From this the linear conditions on the coefficients c_{ij} of P (x, y) = ∑_{i+j≤d−3} c_{ij}x^{i}y^{j} arising from the finite singularities are easily found: having found an integral basis {β_{1}, …, β_{n}}, the above equation is equivalent to demanding that all products β_{j}P (x, y), j = 1, …, n are polynomials in x and y. This is done as follows: using the equation F (x, y) = 0, powers of y^{n} and higher are eliminated from the quantities β_{j}P (x, y). Then these quantities are all reduced to the form G_{j}(x, y)/H_{j}(x), with G_{j}(x, y) a polynomial in x and y, and H_{j}(x) a polynomial in x. This is rewritten as G_{j}(x, y)/H_{j}(x) = Q_{j}(x, y) + R_{j}(x, y)/H_{j}(x), with the degree of R_{j}(x, y) as a polynomial in x less than the degree of H_{j}(x). Condition (20) then states that all coefficients of R_{j}(x, y) as a polynomial in x and y are zero. These coefficients are linear combinations of the c_{ij}, which are equated to zero. After obtaining similar conditions from the singular points at infinity (see below), the total set of linear equations for the coefficients c_{ij} is solved. The solution set of these equations is g-dimensional, because there are g linearly independent holomorphic differentials. By computing a set of g independent solutions and substituting these in P (x, y), a set of g linearly independent adjoint polynomials P_{k}(x, y) is found, and hence by equation (18) a basis ω_{1}, …, ω_{g} of the cohomology is found.

If Γ has special singular points at infinity, then similar reasoning applies, but only after transforming (18) such that it is expressed using the coordinate functions x, z or y, z. Here (x : y : z) are the homogeneous coordinates introduced in Section 2. Recall that for finite points on Γ, z ≠ 0, so that finite points can be denoted by (x, y), with z = 1. Similarly, for infinite points z = 0, but at least one of x or y is non-zero. If at a point at infinity x ≠ 0, then (x : y : z) = (1 : y/x : z/x). In this case, we equate x = 1 and y and z are good local coordinate functions near this point at infinity. Otherwise, if x = 0 but y ≠ 0, then x and z are good local coordinates. In the first case, the differential is transformed to the new coordinate functions using x → 1/z, y → y/z. In the second case, the transformation is x → x/z, y → 1/z. This transformation is now applied to the equation for the plane algebraic curve (2) and the equation for the adjoint polynomial P (x, y) = ∑_{i+j≤d−3} c_{ij}x^{i}y^{j}.

This results in respectively an equation for the algebraic curve in the new coordinates and a polynomial P̃ in the new coordinates, namely the numerator of P (x, y) under the transformation. The coefficients of this new polynomial P̃ are linear combinations of the coefficients c_{ij}. Finding an integral basis for the algebraic curve in the new coordinate functions and applying Mñuk’s result (20) gives linear conditions on the coefficients c_{ij}, in addition to the ones obtained using the coordinate functions x and y.

If for some singular points at infinity x = 0, while for others y = 0, then this process may have to be repeated a total of three times, using all three sets of coordinate functions (x, y), (x, z) and (y, z).

Example: Using F (x, y) = y^{3} − x^{7} + 2x^{3}y = 0, we illustrate the concepts of this section and construct the holomorphic differentials. The interested reader may find it helpful to check our calculations. From (18), all holomorphic differentials are of the form ω = P (x, y)/(3y^{2} + 2x^{3})dx, with P (x, y) = ∑_{i+j≤4} c_{ij}x^{i}y^{j} a polynomial in x and y of degree at most 4. This gives rise to 15 undetermined coefficients c_{ij}. Expressed in homogeneous coordinates (x : y : z), the singular points are P_{1} = (0 : 0 : 1) and P_{2} = (0 : 1 : 0). The second singular point P_{2} is infinite and the conditions it imposes on the coefficients of P (x, y) are derived later.

The multiplicity of P_{1} = (0 : 0 : 1) is m_{P_{1}} = 3, its delta invariant is δ_{P_{1}} = 4. Since δ_{P_{1}} = 4 > 3 = m_{P_{1}}(m_{P_{1}} − 1)/2, the integral basis method is used. The integral basis is found (using the integral_basis command in Maple) to be {1, y/x, y^{2}/x^{3}}. Hence all elements of O_{A} are of the form f = f_{1}(x) + f_{2}(x)y/x + f_{3}(x)y^{2}/x^{3}, where f_{j}(x), j = 1, 2, 3 are polynomials in x. We verify explicitly that y/x and y^{2}/x^{3} satisfy monic polynomial equations with coefficients that are polynomials in x and y. Indeed:

(\frac{y}{x})^{3} + 2x (\frac{y}{x}) − x^{4} = 0, (\frac{y^{2}}{x^{3}} )^{2} + 2 (\frac{y^{2}}{x^{3}} ) − xy = 0.  (21)

A different way to check that y/x and y^{2}/x^{3} have no poles at the two points on the Riemann surface corresponding to P_{1} is to compute the two Puiseux expansions at P_{1}:

\{ \begin{matrix}x = t, \\ y = \frac{1}{2}t^{4} − \frac{1}{16}t^{9} + …, \end{matrix} \{ \begin{matrix}x = −2t^{2}, \\ y = 4t^{3} − 4t^{8} + …, \end{matrix}  (22)

Each of these can be substituted for x and y in y/x and y^{2}/x^{2}. The results have no pole at t = 0, hence y/x and y^{2}/x^{2} do not have a pole at either one of the two points on the Riemann surface corresponding to P_{1}. This integral basis then gives rise to the conditions that P (x, y), yP (x, y)/x and y^{2}P (x, y)/x^{3} are polynomials in x and y. Clearly only the last two of these result in any conditions on the coefficients. Demanding that yP (x, y)/x is a polynomial in x and y, gives c_{00} = 0 = c_{01}. Demanding that y^{2}P (x, y)/x^{3} is a polynomial in x and y gives c_{10} = 0 = c_{20}. As expected, the singular point P_{1} results in δ_{P_{1}} = 4 conditions on the coefficients c_{ij}.

We now turn to the singular point at infinity P_{2} = (0 : 1 : 0). Since y_{P_{2}} ≠ 0, (x, z) are good coordinate functions near this point. After homogenizing y^{3} − x^{7} + 2x^{3}y = 0 and equating y = 1, we find z^{4} − x^{7} + 2x^{3}z^{3} = 0 and this algebraic curve now has a singular point at (x, z) = (0, 0), as expected. The transformed adjoined polynomial is P̃(x, z) = ∑_{i+j≤4} c_{ij}x^{i}z^{4−(i+j)}. Again, the integral basis method is used, since δ_{P_{2}} = 9 > 6 = m_{P_{2}}(m_{P_{2}} − 1)/2. An integral basis is {1, z/x, z^{2}/x^{3}, z^{3}/x^{5}}. Imposing that zP̃(x, z)/x, z^{2}P̃(x, z)/x^{3} and z^{3}P̃(x, z)/x^{5} are polynomials in x and z demands that all c_{ij} = 0, except c_{30} and c_{11}, which are undetermined. Hence the most general adjoint polynomial is

P (x, y) = c_{11}xy + c_{30}x^{3}.  (23)

Thus a basis of holomorphic differentials for the Riemann surface specified by y^{3}−x^{7}+2x^{3}y = 0 is

ω_{1} = \frac{xy}{3y^{2} + 2x^{3}} dx, ω_{2} = \frac{x^{3}}{3y^{2} + 2x^{3}} dx,  (24)

which confirms again that the genus of the Riemann surface considered in this example is g = 2. The calculation of the holomorphic differentials of the Riemann surface specified by y^{3} − x^{7} + 2x^{3}y = 0 using Maple V.6 is given below.

  >with(algcurves):          # load the algcurves package
  >f:=y^3-x^7+2*x^3*y:       # define the algebraic curve
  ># calculate the holomorphic differentials
  >differentials(f,x,y);
                              3
                             x  dx       x y dx
                         [-----------, -----------]
                             2      3     2      3
                          3 y  + 2 x   3 y  + 2 x

6 Computation of a Riemann matrix of a Riemann surface

Having obtained a canonical basis for the homology and the holomorphic differentials of the Riemann surface, a period matrix is found by evaluation of the integrals (6) and (7). Once a period matrix is found, a Riemann matrix for the Riemann surface is found from (9).

From (16),

∮_{a_{j}} ω_{i} = ∑_{k=1}^{r} α_{jk} ∮_{c_{k}} ω_{i}, ∮_{b_{j}} ω_{i} = ∑_{k=1}^{r} α_{g+j,k} ∮_{c_{k}} ω_{i},  (25)

and the computation of a period matrix is reduced to the computation of the integrals ∮_{c_{k}} ω_{i}, k = 1, …, r for every holomorphic differential ω_{i}. But by construction, every one of the cycles c_{k} consists of line segments and semi-circles in the complex x-plane lifted to the Riemann surface. Each one of these line segments or semi-circles is parametrized by x = γ(t), with 0 ≤ t ≤ 1. The lifting of x = γ(t), denoted by y = γ̃(y_{0}, t), is obtained by specifying a starting value of y, y_{0} (essentially the sheet number), and by analytically continuing this value y_{0} along x = γ(t). Hence γ̃(y_{0}, 0) = y_{0} and F (γ(t), γ̃(y_{0}, t)) = 0. This reduces the problem to evaluating integrals of the type

∫_{0}^{1} \frac{P_{i}(γ(t), γ̃(y_{0}, t))}{∂_{y}F (γ(t), γ̃(y_{0}, t))}γ′(t)dt  (26)

These integrals are evaluated numerically using Maple’s numerical integration routine. This has the advantage that the user can specify the number of significant digits to be used in computations. If Digits is the number of significant digits the user specified, the command periodmatrix(f,x,y) attempts to return the periodmatrix of the Riemann surface specified by the plane algebraic curve f = F (x, y) with at least Digits−3 significant digits. The numerical evaluation of these integrals is slow, since for every evaluation of the integrand analytic continuation is required. Of course, requiring more significant digits takes more computer time.

The following example computes a period matrix and a Riemann matrix for the Riemann surface specified by F (x, y) = y^{3} − x^{7} + 2x^{3}y. The accuracy of the output is estimated by the absolute values of the anti-symmetric part of the Riemann matrix (if the error was zero then the Riemann matrix would be symmetric). If the absolute value of an entry of the anti-symmetric part is greater than 10^{Digits−3}, a warning message is printed.

  >with(algcurves):          # load the algcurves package
  >f:=y^3-x^7+2*x^3*y:       # define the algebraic curve
  ># calculate a period matrix for the Riemann surface corresponding to f
  >pm:=periodmatrix(f,x,y):
  >evalf(pm,5)      # use only 5 significant digits, for display purposes
 [-.71618 + .98573 I    -1.1588 - .37652 I    -3.7499    -2.5911 - .37652 I]
 [                                                                         ]
 [-1.8496 - .60096 I    1.1431 + 1.5733 I     -1.4129    -2.5560 + 1.5733 I]
  ># calculate a Riemann matrix for the Riemann matrix corresponding to f
  >rm:=periodmatrix(f,x,y,Riemann);
        [1.690983020 + .9510565349 I    1.500000005 + .3632712707 I]
        [                                                          ]
        [1.500000010 + .3632712735 I    1.309017002 + .9510565221 I]
  >with(linalg):      # load the linalg package
  >evalm(rm-transpose(rm));  # compute the anti-symmetric part of rm
     [                              -8         -8  ]

     [         0              -.5 10   - .28 10   I]
     [                                             ]
     [     -8         -8                           ]
     [.5 10   + .28 10   I              0          ]
  ># increase the significant digits used in computations to 20.
  >Digits:=20;
  ># calculate a Riemann matrix using 20 significant digits.
  >rm:=periodmatrix(f,x,y,Riemann):
  >rmr:=map(Re,rm);          # compute the real part of rm
        [1.6909830056250525744    1.4999999999999999994]
        [                                              ]
        [1.4999999999999999995    1.3090169943749474240]
  >rmi:=map(Im,rm);      # compute the imaginary part of rm
        [.95105651629515357133    .36327126400268044257]
        [                                              ]
        [.36327126400268044298    .95105651629515357255]
  >eigenvals(rmi);           # the eigenvalues of the imaginary part of rm
              .58778525229247312916, 1.3143277802978340147
  ># the imaginary part of the Riemann matrix is indeed positive definite.

Remark: If the roles of the coordinate functions x and y are switched, x is regarded as an algebraic function of y, x = x(y). This changes the monodromy and hence the period matrix computation is very different. On the other hand, the period matrix and the Riemann matrix are merely transformed into equivalent matrices. This equivalence is a consequence of the freedom in choosing a canonical basis of the homology. If a new canonical basis of cycles (a′_{1}, …, a′_{g}, b′_{1}, …, b′_{g}) is chosen as an integer linear combination of the elements of the existing basis by

a′_{i} = ∑_{j=1}^{g} d_{ij}a_{j} + ∑_{j=1}^{g} c_{ij}b_{j}, b′_{i} = ∑_{j=1}^{g} b_{ij}a_{j} + ∑_{j=1}^{g} a_{ij}b_{j},  (27)

then the period matrix Ω transforms according to

Ω → Ω′ = (Ad^{T} + Bc^{T} Ab^{T} + Ba^{T} ) ,  (28)

with a = (a_{ij})_{i,j=1}^{g}, b = (b_{ij})_{i,j=1}^{g}, c = (c_{ij})_{i,j=1}^{g} and d = (d_{ij})_{i,j=1}^{g}. Using the symmetry of the Riemann matrix, we find that it transforms according to

B̂ → B̂′ = (aB̂ + b)(cB̂ + d)^{−1}.  (29)

Because a,b, c and d transform a canonical basis of the homology into a canonical basis of the homology, they must satisfy

( \begin{matrix}a & b \\ c & d \end{matrix} ) ( \begin{matrix}O_{g} & I_{g} \\ −I_{g} & O_{g} \end{matrix} ) ( \begin{matrix}a^{T} & c^{T} \\ b^{T} & d^{T} \end{matrix} ) = ( \begin{matrix}O_{g} & I_{g} \\ −I_{g} & O_{g} \end{matrix} ) ,  (30)

which restates that the intersection indices of the elements of the new basis of the homology are like those of the previous basis.

Acknowledgements

The authors wish to thank B. A. Dubrovin, Y. Kopeliovich, R. MacRae, M. Seppälä and M. D. Tretkoff for useful discussions. This work was initiated while both authors where at the Mathematical Sciences Research Institute during the Fall of 1998 for the Special Program on Symbolic Computation in Geometry and Analysis. It was continued while the first author was at MSRI, the University of Colorado at Boulder, and the University of Washington, and while the second author was at Florida State University. All institutions mentioned are gratefully acknowledged for their hospitality and support. The research presented here was supported through grants DMS-9701755, DMS-9731097 and DMS-9805983.

References

[1] S. S. Abhyankar. Algebraic geometry for scientists and engineers, volume 35 of Mathematical Surveys and Monographs. American Mathematical Society, Providence, RI, 1990.

[2] E. D. Belokolos, A. I. Bobenko, V. Z. Enol’skii, A. R. Its, and V. B. Matveev. Algebrogeometric approach to nonlinear integrable problems. Springer Series in Nonlinear Dynamics. Springer-Verlag, Berlin, 1994.

[3] K. Berry and M. D. Tretkoff. The periodmatrix of macbeath’s curve of genus seven. Contemporary mathematics, 136:31–40, 1992.

[4] K. Berry and M. D. Tretkoff. Private communication, 1999.

[5] E. Brieskorn and H. Knörrer. Plane algebraic curves. Birkhäuser Verlag, Basel, 1986. Translated from the German by John Stillwell.

[6] B. Deconinck and H. Segur. The KP equation with quasiperiodic initial conditions. Physica D, 123:123–152, 1998.

[7] B. Deconinck and M. van Hoeij. An algorithm for computing the riemann matrix of an algebraic curve. in preparation.

[8] B. A. Dubrovin. Theta functions and nonlinear equations. Russian Math. Surveys, 36(2):11–80, 1981.

[9] B. A. Dubrovin, A. T. Fomenko, and S. P. Novikov. Modern geometry—methods and applications. Part II, volume 104 of Graduate Texts in Mathematics. Springer-Verlag, New York, 1985.

[10] H. M. Farkas and I. Kra. Riemann surfaces. Springer-Verlag, New York, second edition, 1992.

[11] P. Gianni, M. Seppälä, R. Silhol, and B. Trager. Riemann surfaces, plane algebraic curves and their period matrices. J. Symbolic Comput., 26(6):789–803, 1998.

[12] P. Griffiths and J. Harris. Principles of algebraic geometry. John Wiley & Sons Inc., New York, 1994. Reprint of the 1978 original.

[13] I. M. Krichever. Integration of nonlinear equations by the methods of algebraic geometry. Funct. Anal. Appl., 11:12–26, 1977.

[14] I. M. Krichever. Nonlinear equations and elliptic curves. J. Soviet Math., 28(1):51–90, 1985.

[15] I. M. Krichever. Spectral theory of two-dimensional periodic operators and its applications. Uspekhi Mat. Nauk, 44(2(266)):121–184, 1989.

[16] M. Mñuk. Computing adjoint curves. J. Symb. Comput., 23:229–240, 1997.

[17] H. P. McKean. Stability for the Korteweg-de Vries equation. Comm. Pure Appl. Math., 30(3):347–353, 1977.

[18] M. Noether. Rationale ausführungen der operationen in der theorie der algebraischen funktionen. Math. Ann., 23:311–358, 1883.

[19] S. P. Novikov, S. V. Manakov, L. P. Pitaevskii, and V. E. Zakharov. Theory of Solitons, volume 1 of Contemporary Soviet Mathematics. Consultants Bureau, New York and London, 1984.

[20] E. Previato. Generalized Weierstrass ℘-functions and KP flows in affine space. Comment. Math. Helv., 62(2):292–310, 1987.

[21] M. Seppälä. Computation of period matrices of real algebraic curves. Discrete Comput. Geom., 11(1):65–81, 1994.

[22] C. L. Siegel. Topics in complex function theory. Vol. I. John Wiley & Sons, Inc., New York, 1988.

[23] C. L. Siegel. Topics in complex function theory. Vol. II. John Wiley & Sons, Inc., New York, 1988.

[24] C. L. Siegel. Topics in complex function theory. Vol. III. John Wiley & Sons, Inc., New York, 1989.

[25] G. Springer. Introduction to Riemann surfaces. Addison-Wesley Publishing Company, Inc., Reading, Mass., 1957.

[26] C.L. Tretkoff and M. D. Tretkoff. Combinatorial group theory, Riemann surfaces and differential equations. Contemporary Mathematics, 33:467–517, 1984.

[27] B. L. van der Waerden. Algebra, Vol. I. Springer-Verlag, New York, 1991.

[28] M. van Hoeij. An algorithm for computing an integral basis in an algebraic function field. J. Symb. Comput., 18:353–363, 1994.

[29] M. van Hoeij. An algorithm for computing the Weierstrass normal form. In ISSAC ’95 Proceedings, pages 90–95, 1995.

Bernard Deconinck Department of Applied Mathematics University of Washington Box 352420 Seattle, WA 98195 deconinc@amath.washington.edu

Mark van Hoeij Department of Mathematics Florida State University 208 James J. Love Building Tallahassee, FL 32306 hoeij@lars.math.fsu.edu