implement inverse Weierstrass elliptic function and lattice invariants

acb_elliptic.h

@@ -58,6 +58,12 @@ void acb_elliptic_zeta(acb_t res, const acb_t z, const acb_t tau, slong prec);
 
 void acb_elliptic_sigma(acb_t res, const acb_t z, const acb_t tau, slong prec);
 
+void acb_elliptic_roots(acb_t e1, acb_t e2, acb_t e3, const acb_t tau, slong prec);
+
+void acb_elliptic_invariants(acb_t g2, acb_t g3, const acb_t tau, slong prec);
+
+void acb_elliptic_inv_p(acb_t res, const acb_t z, const acb_t tau, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_elliptic/inv_p.c

@@ -0,0 +1,35 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+void
+acb_elliptic_inv_p(acb_t res, const acb_t z, const acb_t tau, slong prec)
+{
+    acb_t e1, e2, e3;
+
+    acb_init(e1);
+    acb_init(e2);
+    acb_init(e3);
+
+    acb_elliptic_roots(e1, e2, e3, tau, prec);
+
+    acb_sub(e1, z, e1, prec);
+    acb_sub(e2, z, e2, prec);
+    acb_sub(e3, z, e3, prec);
+
+    acb_elliptic_rf(res, e1, e2, e3, 0, prec);
+
+    acb_clear(e1);
+    acb_clear(e2);
+    acb_clear(e3);
+}
+

acb_elliptic/invariants.c

@@ -0,0 +1,31 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+#include "acb_modular.h"
+
+void
+acb_elliptic_invariants(acb_t g2, acb_t g3, const acb_t tau, slong prec)
+{
+    acb_struct t[2];
+
+    acb_init(t);
+    acb_init(t + 1);
+
+    acb_modular_eisenstein(t, tau, 2, prec);
+
+    acb_mul_ui(g2, t, 60, prec);
+    acb_mul_ui(g3, t + 1, 140, prec);
+
+    acb_clear(t);
+    acb_clear(t + 1);
+}
+

acb_elliptic/roots.c

@@ -0,0 +1,50 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+#include "acb_modular.h"
+
+void
+acb_elliptic_roots(acb_t e1, acb_t e2, acb_t e3, const acb_t tau, slong prec)
+{
+    acb_t t1, t2, t3, t4;
+
+    acb_init(t1);
+    acb_init(t2);
+    acb_init(t3);
+    acb_init(t4);
+
+    acb_modular_theta(t1, t2, t3, t4, t1, tau, prec);
+
+    acb_pow_ui(t2, t2, 4, prec);
+    acb_pow_ui(t4, t4, 4, prec);
+
+    acb_sub(e2, t2, t4, prec);
+    acb_mul_2exp_si(t3, t4, 1);
+    acb_add(e1, t2, t3, prec);
+    acb_mul_2exp_si(t3, t2, 1);
+    acb_add(e3, t3, t4, prec);
+
+    acb_const_pi(t3, prec);
+    acb_mul(t3, t3, t3, prec);
+    acb_div_ui(t3, t3, 3, prec);
+
+    acb_mul(e1, e1, t3, prec);
+    acb_mul(e2, e2, t3, prec);
+    acb_mul(e3, e3, t3, prec);
+    acb_neg(e3, e3);
+
+    acb_clear(t1);
+    acb_clear(t2);
+    acb_clear(t3);
+    acb_clear(t4);
+}
+

acb_elliptic/test/t-inv_p.c

@@ -0,0 +1,66 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("inv_p....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t tau, z, w, pw;
+        slong prec;
+
+        acb_init(tau);
+        acb_init(z);
+        acb_init(w);
+        acb_init(pw);
+
+        prec = 2 + n_randint(state, 400);
+
+        acb_randtest(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
+        acb_randtest(w, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
+        acb_randtest(tau, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
+        if (arf_sgn(arb_midref(acb_imagref(tau))) < 0)
+            acb_neg(tau, tau);
+
+        acb_elliptic_inv_p(w, z, tau, prec);
+        acb_elliptic_p(pw, w, tau, prec);
+
+        if (!acb_contains(pw, z))
+        {
+            flint_printf("FAIL (containment)\n");
+            flint_printf("tau = "); acb_printd(tau, 30); flint_printf("\n\n");
+            flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
+            flint_printf("w = "); acb_printd(w, 30); flint_printf("\n\n");
+            flint_printf("pw = "); acb_printd(pw, 30); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_clear(tau);
+        acb_clear(z);
+        acb_clear(w);
+        acb_clear(pw);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_elliptic/test/t-invariants.c

@@ -0,0 +1,82 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("invariants....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t tau, e1, e2, e3, g2, g3, h2, h3;
+        slong prec;
+
+        acb_init(tau);
+        acb_init(e1);
+        acb_init(e2);
+        acb_init(e3);
+        acb_init(g2);
+        acb_init(g3);
+        acb_init(h2);
+        acb_init(h3);
+
+        prec = 2 + n_randint(state, 400);
+
+        acb_randtest(tau, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
+        if (arf_sgn(arb_midref(acb_imagref(tau))) < 0)
+            acb_neg(tau, tau);
+
+        acb_elliptic_roots(e1, e2, e3, tau, prec);
+        acb_elliptic_invariants(g2, g3, tau, prec);
+
+        acb_mul(h2, e1, e1, prec);
+        acb_addmul(h2, e2, e2, prec);
+        acb_addmul(h2, e3, e3, prec);
+        acb_mul_2exp_si(h2, h2, 1);
+
+        acb_mul(h3, e1, e2, prec);
+        acb_mul(h3, h3, e3, prec);
+        acb_mul_2exp_si(h3, h3, 2);
+
+        if (!acb_overlaps(g2, h2) || !acb_overlaps(g3, h3))
+        {
+            flint_printf("FAIL (overlap)\n");
+            flint_printf("tau = "); acb_printd(tau, 30); flint_printf("\n\n");
+            flint_printf("g2 = "); acb_printd(g2, 30); flint_printf("\n\n");
+            flint_printf("g3 = "); acb_printd(g3, 30); flint_printf("\n\n");
+            flint_printf("h2 = "); acb_printd(h2, 30); flint_printf("\n\n");
+            flint_printf("h3 = "); acb_printd(h3, 30); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_clear(tau);
+        acb_clear(e1);
+        acb_clear(e2);
+        acb_clear(e3);
+        acb_clear(g2);
+        acb_clear(g3);
+        acb_clear(h2);
+        acb_clear(h3);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_elliptic.rst

@@ -3,6 +3,11 @@
 **acb_elliptic.h** -- elliptic integrals and functions of complex variables
 ===============================================================================
 
+This module supports computation of elliptic (doubly periodic) functions,
+and their inverses, elliptic integrals.
+See :ref:`acb_modular.h <acb-modular>` for the closely related modular forms
+and Jacobi theta functions.
+
 Warning: incomplete elliptic integrals have very complicated
 branch structure when extended to complex variables.
 For some functions in this module, branch cuts may be
@@ -252,6 +257,15 @@ in [Car1995]_ and chapter 19 in [NIST2012]_.
 Weierstrass elliptic functions
 -------------------------------------------------------------------------------
 
+Elliptic functions may be defined on a general lattice
+`\Lambda = \{m 2\omega_1 + n 2\omega_2\ : m, n \in \mathbb{Z}\}`
+with half-periods `\omega_1, \omega_2`.
+We simplify by setting
+`2\omega_1 = 1, 2\omega_2 = \tau` with `\operatorname{im}(\tau) > 0`.
+To evaluate the functions on a general lattice, it is enough to make a
+linear change of variables.
+The main reference is chapter 23 in [NIST2012]_.
+
 .. function:: void acb_elliptic_p(acb_t res, const acb_t z, const acb_t tau, slong prec)
 
     Computes Weierstrass's elliptic function
@@ -284,6 +298,31 @@ Weierstrass elliptic functions
     Sets *res* to the Weierstrass elliptic function of the power series *z*,
     with periods 1 and *tau*, truncated to length *len*.
 
+
+.. function:: void acb_elliptic_invariants(acb_t g2, acb_t g3, const acb_t tau, slong prec)
+
+    Computes the lattice invariants `g_2, g_3`. The Weierstrass elliptic
+    function satisfies the differential equation
+    `[\wp'(z, \tau)]^2 = 4 [\wp(z,\tau)]^3 - g_2 \wp(z,\tau) - g_3`.
+    Up to constant factors, the lattice invariants are the first two
+    Eisenstein series (see :func:`acb_modular_eisenstein`).
+
+.. function:: void acb_elliptic_roots(acb_t e1, acb_t e2, acb_t e3, const acb_t tau, slong prec)
+
+    Computes the lattice roots `e_1, e_2, e_3`, which are the roots of
+    the polynomial `4z^3 - g_2 z - g_3`.
+
+.. function:: void acb_elliptic_inv_p(acb_t res, const acb_t z, const acb_t tau, slong prec)
+
+    Computes the inverse of the Weierstrass elliptic function, which
+    satisfies `\wp(\wp^{-1}(z, \tau), \tau) = z`. This function is given
+    by the elliptic integral
+
+    .. math ::
+
+        \wp^{-1}(z, \tau) = \frac{1}{2} \int_z^{\infty} \frac{dt}{\sqrt{(t-e_1)(t-e_2)(t-e_3)}}
+            = R_F(z-e_1,z-e_2,z-e_3).
+
 .. function:: void acb_elliptic_zeta(acb_t res, const acb_t z, const acb_t tau, slong prec)
 
     Computes the Weierstrass zeta function