From 297c04cb7e85488ff63d762cbfd7bd49cd28b9c4 Mon Sep 17 00:00:00 2001 From: Jean Kieffer Date: Mon, 10 Feb 2025 13:07:55 +0100 Subject: [PATCH] Allow sqr parameter in acb_theta_jet; call acb_theta_jet in acb_theta_all --- src/acb_theta.h | 2 +- src/acb_theta/all.c | 122 +------------------------------ src/acb_theta/jet.c | 28 +++++-- src/acb_theta/test/main.c | 2 - src/acb_theta/test/t-all.c | 68 ----------------- src/acb_theta/test/t-g2_chi3_6.c | 2 +- src/acb_theta/test/t-jet.c | 7 +- 7 files changed, 32 insertions(+), 199 deletions(-) delete mode 100644 src/acb_theta/test/t-all.c diff --git a/src/acb_theta.h b/src/acb_theta.h index f3b677b1f7..8a3dc50ac3 100644 --- a/src/acb_theta.h +++ b/src/acb_theta.h @@ -206,7 +206,7 @@ int acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec); void acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, - const acb_mat_t tau, slong ord, int all, slong prec); + const acb_mat_t tau, slong ord, int all, int sqr, slong prec); void acb_theta_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec); void acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, diff --git a/src/acb_theta/all.c b/src/acb_theta/all.c index 54a5f38e5a..d0a5e0aebf 100644 --- a/src/acb_theta/all.c +++ b/src/acb_theta/all.c @@ -9,131 +9,11 @@ (at your option) any later version. See . */ -#include "acb.h" -#include "fmpz_mat.h" -#include "acb_mat.h" #include "acb_theta.h" void acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, int sqr, slong prec) { - slong g = acb_mat_nrows(tau); - slong n2 = 1 << (2 * g); - fmpz_mat_t mat; - acb_mat_t new_tau, N, ct; - acb_ptr new_zs, exps, cs, aux, units; - arb_ptr rs; - acb_t s; - ulong * ch; - slong * e; - slong kappa; - slong j, ab; - int res; - - if (nb <= 0) - { - return; - } - - fmpz_mat_init(mat, 2 * g, 2 * g); - acb_mat_init(new_tau, g, g); - acb_mat_init(N, g, g); - acb_mat_init(ct, g, g); - new_zs = _acb_vec_init(nb * g); - exps = _acb_vec_init(nb); - cs = _acb_vec_init(nb); - aux = _acb_vec_init(n2 * nb); - units = _acb_vec_init(8); - rs = _arb_vec_init(nb * g); - ch = flint_malloc(n2 * sizeof(ulong)); - e = flint_malloc(n2 * sizeof(slong)); - acb_init(s); - - /* Reduce tau then z */ - res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, ct, exps, zs, nb, tau, prec); - if (res) - { - res = acb_theta_reduce_z(new_zs, rs, cs, new_zs, nb, new_tau, prec); - } - - if (res) - { - /* Setup */ - _acb_vec_unit_roots(units, 8, 8, prec); - acb_theta_char_table(ch, e, mat, -1); - if (sqr) - { - kappa = acb_siegel_kappa2(mat); - acb_mat_det(s, ct, prec); - } - else - { - kappa = acb_siegel_kappa(s, mat, new_tau, prec); - } - - acb_theta_jet_notransform(aux, new_zs, nb, new_tau, 0, 0, 1, sqr, prec); - - /* Account for reduce_z */ - for (j = 0; j < nb; j++) - { - if (sqr) - { - acb_sqr(&cs[j], &cs[j], prec); - } - _acb_vec_scalar_mul(aux + j * n2, aux + j * n2, n2, &cs[j], prec); - } - - /* Account for reduce_tau */ - for (j = 0; j < nb; j++) - { - if (sqr) - { - acb_sqr(&exps[j], &exps[j], prec); - } - for (ab = 0; ab < n2; ab++) - { - acb_mul(&th[j * n2 + ab], &aux[j * n2 + ch[ab]], &exps[j], prec); - acb_mul(&th[j * n2 + ab], &th[j * n2 + ab], s, prec); - acb_mul(&th[j * n2 + ab], &th[j * n2 + ab], - &units[((sqr ? 2 : 1) * (kappa + e[ab])) % 8], prec); - } - } - } - else - { - /* Use local_bound to avoid returning NaN */ - arb_t c, rho; - arb_init(c); - arb_init(rho); - - for (j = 0; j < nb; j++) - { - acb_theta_local_bound(c, rho, zs + j * g, tau, 0); - for (ab = 0; ab < n2; ab++) - { - arb_zero_pm_one(acb_realref(&th[j * n2 + ab])); - arb_zero_pm_one(acb_imagref(&th[j * n2 + ab])); - } - _acb_vec_scalar_mul_arb(th + j * n2, th + j * n2, n2, c, prec); - } - - arb_clear(c); - arb_clear(rho); - } - - - fmpz_mat_clear(mat); - acb_mat_clear(new_tau); - acb_mat_clear(ct); - acb_mat_clear(N); - _acb_vec_clear(new_zs, nb * g); - _acb_vec_clear(exps, nb); - _acb_vec_clear(cs, nb); - _acb_vec_clear(aux, n2 * nb); - _acb_vec_clear(units, 8); - _arb_vec_clear(rs, nb * g); - acb_clear(s); - flint_free(e); - flint_free(ch); + acb_theta_jet(th, zs, nb, tau, 0, 1, sqr, prec); } diff --git a/src/acb_theta/jet.c b/src/acb_theta/jet.c index c047cca19c..ba6610e928 100644 --- a/src/acb_theta/jet.c +++ b/src/acb_theta/jet.c @@ -16,7 +16,7 @@ void acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, - slong ord, int all, slong prec) + slong ord, int all, int sqr, slong prec) { slong g = acb_mat_nrows(tau); slong n = 1 << g; @@ -60,19 +60,33 @@ acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, { res = acb_theta_reduce_z(new_zs, rs, cs, new_zs, nb, new_tau, prec); } + sqr = sqr && (ord == 0); if (res) { /* Setup */ - _acb_vec_unit_roots(units, 8, 8, prec); - kappa = acb_siegel_kappa(s, mat, new_tau, prec); + if (sqr) + { + _acb_vec_unit_roots(units, 4, 4, prec); + kappa = acb_siegel_kappa2(mat); + acb_mat_det(s, ct, prec); + } + else + { + _acb_vec_unit_roots(units, 8, 8, prec); + kappa = acb_siegel_kappa(s, mat, new_tau, prec); + } acb_theta_char_table(ch, e, mat, (all ? -1 : 0)); - acb_theta_jet_notransform(aux, new_zs, nb, new_tau, ord, *ch, all, 0, prec); + acb_theta_jet_notransform(aux, new_zs, nb, new_tau, ord, *ch, all, sqr, prec); /* Account for reduce_z */ for (j = 0; j < nb; j++) { + if (sqr) + { + acb_sqr(&cs[j], &cs[j], prec); + } _acb_vec_scalar_mul(aux + j * nbth * nbjet, aux + j * nbth * nbjet, nbth * nbjet, &cs[j], prec); _arb_vec_neg(r, rs + j * g, g); @@ -90,10 +104,14 @@ acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, for (j = 0; j < nb; j++) { acb_theta_jet_exp_qf(jet, zs + j * g, N, ord, prec); + if (sqr) + { + acb_sqr(&jet[0], &jet[0], prec); + } for (ab = 0; ab < nbth; ab++) { - acb_mul(t, s, &units[(kappa + e[ab]) % 8], prec); + acb_mul(t, s, &units[(kappa + e[ab]) % (sqr ? 4 : 8)], prec); _acb_vec_scalar_mul(th + j * nbth * nbjet + ab * nbjet, aux + j * nbth * nbjet + (all ? ch[ab] : 0) * nbjet, nbjet, t, prec); diff --git a/src/acb_theta/test/main.c b/src/acb_theta/test/main.c index d8ee598320..8498319bae 100644 --- a/src/acb_theta/test/main.c +++ b/src/acb_theta/test/main.c @@ -15,7 +15,6 @@ #include "t-agm_mul.c" #include "t-agm_mul_tight.c" #include "t-agm_sqrt.c" -#include "t-all.c" #include "t-char_dot.c" #include "t-char_is_even.c" #include "t-char_is_goepel.c" @@ -79,7 +78,6 @@ test_struct tests[] = TEST_FUNCTION(acb_theta_agm_mul), TEST_FUNCTION(acb_theta_agm_mul_tight), TEST_FUNCTION(acb_theta_agm_sqrt), - TEST_FUNCTION(acb_theta_all), TEST_FUNCTION(acb_theta_char_dot), TEST_FUNCTION(acb_theta_char_is_even), TEST_FUNCTION(acb_theta_char_is_goepel), diff --git a/src/acb_theta/test/t-all.c b/src/acb_theta/test/t-all.c deleted file mode 100644 index 089377382d..0000000000 --- a/src/acb_theta/test/t-all.c +++ /dev/null @@ -1,68 +0,0 @@ -/* - Copyright (C) 2024 Jean Kieffer - - This file is part of FLINT. - - FLINT 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 3 of the License, or - (at your option) any later version. See . -*/ - -#include "test_helpers.h" -#include "acb_mat.h" -#include "acb_theta.h" - -TEST_FUNCTION_START(acb_theta_all, state) -{ - slong iter; - - /* Test: agrees with jet_notransform */ - for (iter = 0; iter < 20 * flint_test_multiplier(); iter++) - { - slong g = 1 + n_randint(state, 2); - slong n2 = 1 << (2 * g); - slong nb = n_randint(state, 3); - slong mprec = 50 + n_randint(state, 400); - slong prec = mprec + 50; - slong bits = n_randint(state, 4); - int sqr = n_randint(state, 2); - acb_mat_t tau; - acb_ptr z; - acb_ptr th, test; - - acb_mat_init(tau, g, g); - z = _acb_vec_init(nb * g); - th = _acb_vec_init(nb * n2); - test = _acb_vec_init(nb * n2); - - /* Sample tau not too far from reduced domain */ - acb_siegel_randtest_reduced(tau, state, prec, bits); - acb_mat_scalar_mul_2exp_si(tau, tau, -1); - acb_siegel_randtest_vec_reduced(z, state, nb, tau, 0, prec); - - /* Call theta_all at precision mprec, test against jet_notransform */ - acb_theta_all(th, z, nb, tau, sqr, mprec); - acb_theta_jet_notransform(test, z, nb, tau, 0, 0, 1, sqr, prec); - - if (!_acb_vec_overlaps(th, test, nb * n2)) - { - flint_printf("FAIL\n"); - flint_printf("g = %wd, mprec = %wd, prec = %wd, nb = %wd, sqr = %wd, tau, z:\n", - g, mprec, prec, nb, sqr); - acb_mat_printd(tau, 5); - _acb_vec_printd(z, nb * g, 5); - flint_printf("th, test:\n"); - _acb_vec_printd(th, nb * n2, 5); - _acb_vec_printd(test, nb * n2, 5); - flint_abort(); - } - - acb_mat_clear(tau); - _acb_vec_clear(z, nb * g); - _acb_vec_clear(th, nb * n2); - _acb_vec_clear(test, nb * n2); - } - - TEST_FUNCTION_END(state); -} diff --git a/src/acb_theta/test/t-g2_chi3_6.c b/src/acb_theta/test/t-g2_chi3_6.c index 5d4903ed64..b1ecf8a051 100644 --- a/src/acb_theta/test/t-g2_chi3_6.c +++ b/src/acb_theta/test/t-g2_chi3_6.c @@ -25,7 +25,7 @@ acb_theta_g2_chi8_6(acb_poly_t res, const acb_mat_t tau, slong prec) z = _acb_vec_init(2); acb_init(c); - acb_theta_jet(dth, z, 1, tau, 1, 1, prec); + acb_theta_jet(dth, z, 1, tau, 1, 1, 0, prec); acb_theta_g2_chi3_6(res, dth, prec); for (k = 0; k < 16; k++) { diff --git a/src/acb_theta/test/t-jet.c b/src/acb_theta/test/t-jet.c index 9ed05e344a..5dd8cf7643 100644 --- a/src/acb_theta/test/t-jet.c +++ b/src/acb_theta/test/t-jet.c @@ -30,6 +30,7 @@ TEST_FUNCTION_START(acb_theta_jet, state) slong nbjet = acb_theta_jet_nb(ord, g); int all = iter % 2; slong nbth = (all ? n * n : 1); + int sqr = n_randint(state, 2); acb_mat_t tau; acb_ptr z; acb_ptr th, test; @@ -46,8 +47,12 @@ TEST_FUNCTION_START(acb_theta_jet, state) _acb_vec_scalar_mul_2exp_si(z, z, nb * g, 1); /* Call jet at precision mprec, test against jet_notransform */ - acb_theta_jet(th, z, nb, tau, ord, all, mprec); + acb_theta_jet(th, z, nb, tau, ord, all, sqr, mprec); acb_theta_jet_notransform(test, z, nb, tau, ord, 0, all, 0, prec); + if (ord == 0 && sqr) + { + _acb_vec_sqr(test, test, nb * nbth, prec); + } if (!_acb_vec_overlaps(th, test, nb * nbth * nbjet)) {