Implementation of elliptic functions/integrals

[bluescarni]

Not really an issue, more like a hopefully polite request :)

I am an enthusiastic user of the elliptic functions/integrals module in mpmath. Any plans for implementing them in arb as well?

[fredrik-johansson]

Yes, these are on the horizon.

Jacobi theta functions and the Dedekind eta function are the easiest ones to implement (the theta series are fast, and error bounds are trivial). One might also want to implement transformations to reduce to the smallest/fastest possible argument in every case.

Elliptic integrals are a bit more complicated. The complete integrals can be computed using arithmetic-geometric mean iteration. I think most of the necessary error analysis exists in the literature. For incomplete elliptic integrals, I'm not so sure. I hope Carlson's symmetric iteration extends to balls/intervals in an obvious way, but don't remember the details.

I think working out the interface is a serious concern. Do you use the variable q, tau, k, m, etc? Avoiding singularities and branch cuts is very important. With balls, you can't really compute a function that includes a q^(1/24) singularity and then expect to cancel it out later, like mpmath often does. Generally, it makes sense to implement the functions that are most computationally convenient, not necessarily the ones that are mathematically most standard.

I would also like to have many of these functions for formal power series, but that can obviously be postponed.

[bluescarni]

@fredrik-johansson Cheers!

Personally, I am working a lot with the Weierstrassian formalism for the elliptic functions, so what I did was essentially leveraging the Jacobi elliptic functions and Legendre integrals to implement numerically \wp, \sigma, \zeta, etc. (see A+S chapter 18).

You can find one instance here:

https://github.com/bluescarni/stark_weierstrass/blob/master/weierstrass_elliptic.py

The interface of those implementations is kinda different from mpmath, as they were built to perform calculations for a specific paper. I tinkered with the idea of submitting them to mpmath but never gotten around doing it. I have no complaints against the interface in mpmath, but if you want to take a look a possible alternative is the one in the file above (kind of object-orienty).

There are better ways of computing these of course, see for instance the general schema here for the Weierstrassian elliptic and related:

http://imajna.oxfordjournals.org/content/10/1/119.abstract

And this for the Weierstrass inverse (basically, an elliptic integral):

http://ir.library.oregonstate.edu/xmlui/handle/1957/17420

@darioizzo and I have done some work towards the implementation of the first paper in C/C++ (but I don't think it's in a usable state yet).

[fredrik-johansson]

I've started working on a module for modular forms, and I will be adding the Jacobi theta functions soon.

I looked at Carlson's paper (http://arxiv.org/abs/math/9409227v1) again, and fortunately, the incomplete elliptic integrals are easy. One can basically just apply the symmetric iteration until the arguments overlap and then take their union.

In fact the specific algorithms given by Carlson (which finish with a power series expansion for faster termination) should work too. I just hope his error bounds are correct. The way it is written makes me a bit worried that "the error decreases by a factor 4" actually means something like "the error decreases by a factor 4 - O(1/n)". And there is perhaps some subtlety when extending the algorithms to balls. So I would like to see a formal proof.

[fredrik-johansson]

By the way, I've decided that theta functions definitely will be functions of tau i.e. theta_3(z; tau) = sum_n exp(pi i (n^2 tau + 2 n z)), etc.

[bluescarni]

@fredrik-johansson Great to hear that! I am assuming that the implementation is going to allow complex arguments, like mpmath's?

Any thought about the implementation of the Weierstrassian functions? They can be implemented in terms of the theta functions, but I don't know if there are any pitfalls when moving over to ball/interval arithmetic.

[fredrik-johansson]

Yeah, everything will be for complex arguments.

I might not implement derivatives of theta functions just yet, though. If you need those for the Weierstrassian functions, you will have to settle for finite differences, for the time being...

[bluescarni]

@fredrik-johansson Ok that seems reasonable, thanks for working on this!

Are there specific arb functions dealing with finite differences or otherwise how would one provide a guarantee on the error with respect to the exact derivative?

[fredrik-johansson]

There are no functions for this currently. I have not thought this through, but I think you could derive a generic error bound for finite differences in terms of the Cauchy integral formula. Which means you could use acb_calc_cauchy_bound. Of course this will be horribly slow compared to proper implementations of the derivatives.

[argriffing]

Speaking of differentiation and arb, what about 'automatic' or 'algorithmic' differentiation? If it is implemented it could be better than finite differences but worse than custom coded derivatives of functions.

[fredrik-johansson]

Automatic differentiation is just power series arithmetic, which Arb already supports for most of its functions. Theta function derivatives will be done the same way.

[fredrik-johansson]

I've added some code for evaluating theta functions. It's not fully optimized right now except for computing theta constants (z = 0), and it assumes |q| << 1; I will implement argument reduction for large |q| later (this is a bit tricky).

[bluescarni]

@fredrik-johansson Cheers!

Just FYI, have you seen these papers

http://www.sciencedirect.com/science/article/pii/S0377042711001270 http://www.researchgate.net/publication/226331661_Fast_computation_of_complete_elliptic_integrals_and_Jacobian_elliptic_functions

and related? They claim to be faster than Carlson's method, but I do not know if they offer mathematically rigorous guarantees.

[fredrik-johansson]

@bluescarni They may be faster, but the error analysis is probably harder.

I've implemented the Jacobi theta functions complete with the modular transformations.

This should be enough to compute elliptic functions, at least if one is willing to make them functions of the half-period ratio tau instead of insisting on m or k.

I guess for an elliptic function f(z,tau), you'd usually want to evaluate it for several different values of z with tau fixed. This could be made a lot faster with some precomputation. But I'm not sure what the interface for that should look like.

[bluescarni]

@fredrik-johansson Regarding the interface, I think that you are right that one often wants to evaluate the functions with fixed tau and varying z. With "precomputation", do you mean to evaluate "offline" part of the series expansion in order to reduce the runtime cost, or are you referring to transformations on tau to be performed before actually doing the computation?

[fredrik-johansson]

Both: you could both precompute the transformation information for a given tau and precompute the series coefficients.

[fredrik-johansson]

I've implemented Weierstrass's elliptic function \wp.

It's definitely far from optimal right now, even for a single evaluation. I don't really want to try to micro-optimize this code yet though.

[bluescarni]

@fredrik-johansson Cheers, that sounds awesome! I imagine the implementation of the other weierstrassian functions (zeta, sigma and P') would follow the same pattern?

Seems like I need to start writing the acb C++ wrapper :)

[fredrik-johansson]

I think you really need derivatives of theta functions for zeta, sigma and P'.

I should really do the derivatives...

[bluescarni]

Ahh yes right, forgot about that part.

[fredrik-johansson]

I've added computation of theta function derivatives and derivatives of \wp, or rather their power series expansions (as it turns out, computing the first few theta function derivatives simultaneously is barely more expensive than computing the zeroth derivative alone).

Right now the tau transformation isn't supported, but that shouldn't be too hard. I should have time to do it next week.

I still need to do lots of cleanup, rename some functions, and add test code for the derivatives.

One obvious optimization I haven't done yet is argument reduction with respect to z (subtracting the largest possible multiple of tau). This will probably have to wait.

[fredrik-johansson]

This technically means that arb can evaluate all elliptic functions, since all elliptic functions are rational functions in \wp and \wp'.

I should perhaps add a function for 1 / \wp also, avoiding division by zero at the origin.

[bluescarni]

@fredrik-johansson Thanks, this looks awesome! Looking forward to test them out.

[fredrik-johansson]

I've implemented incomplete elliptic integrals. Also, the various elliptic integrals and functions are now in their own module: http://arblib.org/acb_elliptic.html

[fredrik-johansson]

Added the Weierstrass zeta and sigma functions for good measure.

[fredrik-johansson]

And finally, I implemented the inverse Weierstrass function. This issue seems safe to close now (though we could still add more notational versions, e.g. the various Jacobi elliptic functions).