Skip to content

Commit

Permalink
Allow sqr parameter in acb_theta_jet; call acb_theta_jet in acb_theta…
Browse files Browse the repository at this point in the history 
…_all
  • Loading branch information
@j-kieffer
j-kieffer committed Feb 10, 2025
1 parent 5d6c6ba commit 297c04c
Show file tree
Hide file tree
Showing 7 changed files with 32 additions and 199 deletions.
2 changes: 1 addition & 1 deletion src/acb_theta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
122 changes: 1 addition & 121 deletions src/acb_theta/all.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,131 +9,11 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#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);
}
28 changes: 23 additions & 5 deletions src/acb_theta/jet.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down
2 changes: 0 additions & 2 deletions src/acb_theta/test/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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),
Expand Down
68 changes: 0 additions & 68 deletions src/acb_theta/test/t-all.c
View file

This file was deleted.

2 changes: 1 addition & 1 deletion src/acb_theta/test/t-g2_chi3_6.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++)
{
Expand Down
7 changes: 6 additions & 1 deletion src/acb_theta/test/t-jet.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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))
{
Expand Down

0 comments on commit 297c04c

Please sign in to comment.