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))
{