Modulo some borderline cases (and unforseen bugs) abelfunctions is getting close to being able to compute period matrices of algebraic curves. The final ingredient is an algorithm for computing a basis of holomorphic differentials on the curve. In this article, we review the theory behind the algorithm and discuss how SymPy is used in the implementation of the algorithm.
Given a Riemann surface defined by the locus of a complex plane algebraic curve there is a space of differentials which are holomorphic on . This space, denoted , has dimension where is the genus of the Riemann surface.
Using a theorem of Mñuk , a basis for this space of holomorphic differentials can be computed using the integral basis of the algebraic function field
The holomorphic differerentials are all of the form where are the so-called adjoint polynomials of and is the (projective) degree of . Therefore, finding a basis of differentials reduces to finding conditions on the coefficients of the adjoint polynomials. Mñuk proves that the the space of adjoint polynomials of the Riemann surface, , is determined completely by the integral basis of the algebraic function field.
That is, the adjoint polynomials are those such that for all . Since these requirements impose linear conditions on the coefficients of .
abelfunctions uses SymPy for all of its symbolic computaitons. However, personal experience shows that it can be difficult keeping halfway decent performance in the symbolic portions of the code. For example, some manual manipulation trickery is used in both puiseux.py:_new_polynomial() and integral_basis.py:(several functions) as opposed to using built-in SymPy functions. On a personal note, these manual implementations were chosen because I don't understand SymPy well enough and need to take more time to explore its intricacies.
The SymPy Wiki's FAQ Page admits that its current implementation is approximately 10x slower than most other computer algebra systems. This is mostly due to the fact that it is implemented purely in Python. Ondřej Čertík has been working on an experimental SymPy core written in C++ that will hopefully boost performance!
Thankfully, we can determine the aforementioned linear conditions on the coefficients in a reasonable amount of time without having to resort to "manual" approaches.
Let's look at the example . Here, so the adjoint polynomials of are of the form
We need to keep track of these undetermined coefficients so we construct like so:
P(x,y) = 2 3 c_0_0 + c_0_1*y + c_0_2*y + c_0_3*y + c_1_0*x + c_1_1*x*y + 2 3 2 2 2 2 c_1_2*x*y + c_1_3*x*y + c_2_0*x + c_2_1*x *y + c_2_2*x *y 3 3 + c_3_0*x + c_3_1*x *y
The integral basis can be computed using abelfunctions.
beta = 2 y y [1, -, --] x 3 x
Now, at this point, we could do the "manual approach" and compute product of with , do a sympy.expand(), and then use some expression string parsing functions to figure out which terms contain negative exponents. Then, we would have to parse the expression string again to find the coefficients and set the appropriate expressions equal to zero. This has already been attempted and my head hurt after a while with all of the spaghetti code and bugs that cropped up.
Instead, let's take a stab at using SymPy's polynomial division algorithms. What's nice here, as opposed to the Puiseux series and integral basis contexts, is that we're only working with symbolic expressions that have integral exponents. (Rational exponents have been difficult to deal with so far.) We cast the appropriate expressions into polynomials in since we are only dividing by such creatures.
Remainder: 2 3 4 5 2 / 2 c_0_0*y + c_0_1*y + c_0_2*y + c_0_3*y + x *\c_2_0*y + 3 4\ / 2 3 4 c_2_1*y + c_2_2*y / + x*\c_1_0*y + c_1_1*y + c_1_2*y 5\ + c_1_3*y / Resulting Equations: [c_2_2, c_2_1, c_2_0, c_1_3, c_1_2, c_1_1, c_1_0, c_0_3, c _0_2, c_0_1, c_0_0]
Note that none of these terms are divisible by . Therefore, the set of linear equations we obtain from this particular integral basis element is . We repeat this process for each integral basis element and solve the resulting system of linear equations on to obtain our adjoint polynomials; which, in this case, is the polynomial
In the event where the denominator of contains several factors of the form we might have to repeat the above process for every such factor. (Note that since we are working over the denominator splits completely.) This is because we have to impose regularity at all possible singular points, not just those at as in the above example. Whether or not we can simply divide by and easily arrive at the desired conditions on is something to experiment with. One approach would be to shift and then divide by . I'm not sure, though, if this is at all necessary nor am I sure what is the best way in SymPy to factor polynomials over . One approach wold be to find the roots over using Poly.all_roots().
My hope is to have a decenly-perfoming implementation of the Mñuk-based algorithm by the end of the week. Once this is done it's time to compute period matrices!
 Mñuk, M. Computing adjoint curves. J. Symb. Comput. 23, 229–240 (1997)