From 151193e52dec57b69ce157911afd52eeb1fd40bd Mon Sep 17 00:00:00 2001 From: Daniel Schultz Date: Sat, 22 May 2021 18:50:30 +0200 Subject: [PATCH 1/2] add acb_elliptic_p_prime --- acb_elliptic.h | 2 + acb_elliptic/p_prime.c | 58 ++++++++++++++++++++++++++ acb_elliptic/test/t-p_p_prime.c | 73 +++++++++++++++++++++++++++++++++ doc/source/acb_elliptic.rst | 4 ++ 4 files changed, 137 insertions(+) create mode 100644 acb_elliptic/p_prime.c create mode 100644 acb_elliptic/test/t-p_p_prime.c diff --git a/acb_elliptic.h b/acb_elliptic.h index c0ebcea6..b19ef371 100644 --- a/acb_elliptic.h +++ b/acb_elliptic.h @@ -50,6 +50,8 @@ void acb_elliptic_pi_inc(acb_t res, const acb_t n, const acb_t phi, const acb_t void acb_elliptic_p(acb_t r, const acb_t z, const acb_t tau, slong prec); +void acb_elliptic_p_prime(acb_t r, const acb_t z, const acb_t tau, slong prec); + void acb_elliptic_p_jet(acb_ptr r, const acb_t z, const acb_t tau, slong len, slong prec); void _acb_elliptic_p_series(acb_ptr res, acb_srcptr z, slong zlen, const acb_t tau, slong len, slong prec); diff --git a/acb_elliptic/p_prime.c b/acb_elliptic/p_prime.c new file mode 100644 index 00000000..920b6343 --- /dev/null +++ b/acb_elliptic/p_prime.c @@ -0,0 +1,58 @@ +/* + Copyright (C) 2021 Daniel Schultz + + 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 . +*/ + +#include "acb_elliptic.h" +#include "acb_modular.h" + +void +acb_elliptic_p_prime(acb_t r, const acb_t z, const acb_t tau, slong prec) +{ + acb_struct tz[4]; + acb_t t1, t2, t3; + int i, real; + + real = acb_is_real(z) && arb_is_int_2exp_si(acb_realref(tau), -1) && + arb_is_positive(acb_imagref(tau)); + + acb_init(t1); + acb_init(t2); + acb_init(t3); + + for (i = 0; i < 4; i++) + acb_init(tz + i); + + acb_modular_theta(tz + 0, tz + 1, tz + 2, tz + 3, z, tau, prec); + + /* (-2*pi*eta^2/theta1)^3*theta2*theta3*theta4 */ + acb_const_pi(t2, prec); + acb_mul_2exp_si(t2, t2, 1); + acb_neg(t2, t2); + acb_modular_eta(t3, tau, prec); + acb_mul(t1, t3, t3, prec); + acb_mul(t3, t1, t2, prec); + acb_div(t1, t3, tz + 0, prec); + acb_mul(t2, t1, t1, prec); + acb_mul(t3, t1, t2, prec); + acb_mul(t1, tz + 1, tz + 2, prec); + acb_mul(t2, t1, tz + 3, prec); + acb_mul(r, t3, t2, prec); + + if (real) + arb_zero(acb_imagref(r)); + + acb_clear(t1); + acb_clear(t2); + acb_clear(t3); + + for (i = 0; i < 4; i++) + acb_clear(tz + i); +} + diff --git a/acb_elliptic/test/t-p_p_prime.c b/acb_elliptic/test/t-p_p_prime.c new file mode 100644 index 00000000..69ccbdb0 --- /dev/null +++ b/acb_elliptic/test/t-p_p_prime.c @@ -0,0 +1,73 @@ +/* + Copyright (C) 2021 Daniel Schultz + + 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 . +*/ + +#include "acb_elliptic.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("p_p_prime...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) + { + acb_struct pj[2]; + acb_t tau, z, p, pp; + slong prec; + + acb_init(tau); + acb_init(z); + acb_init(p); + acb_init(pp); + acb_init(pj + 0); + acb_init(pj + 1); + + prec = 2 + n_randint(state, 400); + + acb_randtest(z, 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_p(p, z, tau, prec); + acb_elliptic_p_prime(pp, z, tau, prec); + acb_elliptic_p_jet(pj, z, tau, 2, prec); + + if (!acb_overlaps(p, pj + 0) || !acb_overlaps(pp, pj + 1)) + { + flint_printf("FAIL (overlap)\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("p = "); acb_printd(p, 30); flint_printf("\n\n"); + flint_printf("pp = "); acb_printd(pp, 30); flint_printf("\n\n"); + flint_printf("pj0 = "); acb_printd(pj + 0, 30); flint_printf("\n\n"); + flint_printf("pj1 = "); acb_printd(pj + 1, 30); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(tau); + acb_clear(z); + acb_clear(p); + acb_clear(pp); + acb_clear(pj + 0); + acb_clear(pj + 1); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/acb_elliptic.rst b/doc/source/acb_elliptic.rst index f474d441..7ad0a0c5 100644 --- a/doc/source/acb_elliptic.rst +++ b/doc/source/acb_elliptic.rst @@ -299,6 +299,10 @@ The main reference is chapter 23 in [NIST2012]_. \frac{\theta_4^2(z,\tau)}{\theta_1^2(z,\tau)} - \frac{\pi^2}{3} \left[ \theta_2^4(0,\tau) + \theta_3^4(0,\tau)\right]. +.. function:: void acb_elliptic_p_prime(acb_t res, const acb_t z, const acb_t tau, slong prec) + + Computes the derivative `\wp'(z, \tau)` of Weierstrass's elliptic function `\wp(z, \tau)`. + .. function:: void acb_elliptic_p_jet(acb_ptr res, const acb_t z, const acb_t tau, slong len, slong prec) Computes the formal power series `\wp(z + x, \tau) \in \mathbb{C}[[x]]`, From 43e74b1f57db6ca4a6bf72eae683c3469e8094c2 Mon Sep 17 00:00:00 2001 From: Daniel Schultz Date: Sat, 22 May 2021 19:01:30 +0200 Subject: [PATCH 2/2] check relation --- acb_elliptic/test/t-p_p_prime.c | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/acb_elliptic/test/t-p_p_prime.c b/acb_elliptic/test/t-p_p_prime.c index 69ccbdb0..dc5ad294 100644 --- a/acb_elliptic/test/t-p_p_prime.c +++ b/acb_elliptic/test/t-p_p_prime.c @@ -24,7 +24,7 @@ int main() for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) { acb_struct pj[2]; - acb_t tau, z, p, pp; + acb_t tau, z, p, pp, g2, g3, t; slong prec; acb_init(tau); @@ -33,6 +33,9 @@ int main() acb_init(pp); acb_init(pj + 0); acb_init(pj + 1); + acb_init(g2); + acb_init(g3); + acb_init(t); prec = 2 + n_randint(state, 400); @@ -57,12 +60,34 @@ int main() flint_abort(); } + acb_elliptic_invariants(g2, g3, tau, prec); + acb_pow_ui(pj + 0, pp, 2, prec); + + acb_mul(t, p, g2, prec); + acb_add(t, t, g3, prec); + acb_pow_ui(pj + 1, p, 3, prec); + acb_mul_ui(pj + 1, pj + 1, 4, prec); + acb_sub(pj + 1, pj + 1, t, prec); + + if (!acb_overlaps(pj + 0, pj + 1)) + { + flint_printf("FAIL (check pp^2 = 4p^3-g2*p-g3)\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("p = "); acb_printd(p, 30); flint_printf("\n\n"); + flint_printf("pp = "); acb_printd(pp, 30); flint_printf("\n\n"); + flint_abort(); + } + acb_clear(tau); acb_clear(z); acb_clear(p); acb_clear(pp); acb_clear(pj + 0); acb_clear(pj + 1); + acb_clear(g2); + acb_clear(g3); + acb_clear(t); } flint_randclear(state);