Skip to content

Commit

Permalink
Allow all=1 in several functions and tests; fix old bug in ctx_z_shif…
Browse files Browse the repository at this point in the history 
…t_a0
  • Loading branch information
@j-kieffer
j-kieffer committed Aug 30, 2024
1 parent e08ca3d commit 6a892af
Show file tree
Hide file tree
Showing 15 changed files with 466 additions and 169 deletions.
15 changes: 9 additions & 6 deletions src/acb_theta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,15 +313,18 @@ void acb_theta_sum_jet_all(acb_ptr th, const acb_theta_ctx_z_struct * vec, slong
/* Quasilinear algorithms */

int acb_theta_ql_nb_steps(slong * pattern, const arb_mat_t cho, slong prec);
int acb_theta_ql_setup(acb_ptr rts, acb_ptr t, slong * guard, slong * easy_steps,
int acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong ** pts,
slong * nb, arf_t err, slong * fullprec, acb_srcptr z, const acb_mat_t tau,
arb_srcptr distances, slong s, ulong a, slong prec);
void acb_theta_ql_recombine(acb_ptr th, acb_srcptr th0, acb_srcptr cofactors,
const slong * pts, slong nb, const arf_t err, slong fullprec,
slong s, ulong a, int all, slong g, slong prec);
int acb_theta_ql_setup(acb_ptr rts, acb_ptr rts0, acb_ptr t, slong * guard, slong * easy_steps,
acb_srcptr zs, slong nb, const acb_mat_t tau, arb_srcptr distances,
slong nb_steps, int all, slong prec);
int acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
arf_t err, slong * fullprec, acb_srcptr z, const acb_mat_t tau,
arb_srcptr distances, slong s, ulong a, slong prec);
void acb_theta_ql_steps(acb_ptr th, acb_ptr th_init, acb_srcptr rts,
slong nb, slong nb_steps, arb_srcptr distances, const slong * easy_steps,
slong g, slong prec);
acb_srcptr rts_all, slong nb, slong nb_steps, arb_srcptr distances,
const slong * easy_steps, int all, slong g, slong prec);
int acb_theta_ql_exact(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
const slong * pattern, int all, int shifted_prec, slong prec);

Expand Down
5 changes: 2 additions & 3 deletions src/acb_theta/ctx_z_shift_a0.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@ acb_theta_aj_is_zero(ulong a, slong j, slong g)
return !((a >> (g - 1 - j)) & 1);
}

/* We ignore uinv here. */
void
acb_theta_ctx_z_shift_a0(acb_theta_ctx_z_t res, const acb_theta_ctx_z_t ctx,
const acb_theta_ctx_tau_t ctx_tau, ulong a, slong prec)
Expand Down Expand Up @@ -62,10 +61,10 @@ acb_theta_ctx_z_shift_a0(acb_theta_ctx_z_t res, const acb_theta_ctx_z_t ctx,
}
acb_mul(acb_theta_ctx_c(res), c, acb_theta_ctx_exp_a_tau_a_div_4(ctx_tau, a), prec);

/* Compute c, u, r, v */
/* Compute c, r, v; u and uinv must be left invariant */
_arb_vec_set(acb_theta_ctx_r(res), acb_theta_ctx_r(res), g);
acb_abs(abs, acb_theta_ctx_c(res), prec);
arb_mul(acb_theta_ctx_u(res), acb_theta_ctx_u(ctx), abs, prec);
arb_set(acb_theta_ctx_u(res), acb_theta_ctx_u(ctx));
acb_mul(acb_theta_ctx_c(res), acb_theta_ctx_c(res), acb_theta_ctx_c(ctx), prec);

acb_theta_char_get_arb(v_shift, a, g);
Expand Down
6 changes: 3 additions & 3 deletions src/acb_theta/ql_all_new.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
#include "acb_theta.h"

static void
acb_theta_ql_dupl(acb_ptr th2, acb_srcptr th0, acb_srcptr th, slong g, slong prec)
acb_theta_ql_dupl_all(acb_ptr th2, acb_srcptr th0, acb_srcptr th, slong g, slong prec)
{
slong n = 1 << g;
acb_ptr v;
Expand Down Expand Up @@ -137,7 +137,7 @@ acb_theta_ql_all_sum(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
acb_theta_sum_a0_tilde(new_th, vec, nb, ctx_tau, distances, prec);
for (j = 0; j < nb; j++)
{
acb_theta_ql_dupl(th + j * n * n, new_th,
acb_theta_ql_dupl_all(th + j * n * n, new_th,
new_th + j * n, g, prec);
}

Expand Down Expand Up @@ -204,7 +204,7 @@ acb_theta_ql_all_mid_err(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t ta
{
for (j = 0; j < nb; j++)
{
acb_theta_ql_dupl(th + j * n * n, new_th,
acb_theta_ql_dupl_all(th + j * n * n, new_th,
new_th + (j + add_zero) * n, g, prec);
}
}
Expand Down
101 changes: 66 additions & 35 deletions src/acb_theta/ql_exact.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,27 @@ acb_theta_ql_exact_sum(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
{
for (j = 0; j < nb; j++)
{
acb_theta_sum_a0_tilde(th + j * n, &vec[j], 1, ctx_tau, distances + j * n, prec);
if (all)
{
acb_theta_sum_all_tilde(th + j * n * n, &vec[j], 1, ctx_tau, distances + j * n, prec);
}
else
{
acb_theta_sum_a0_tilde(th + j * n, &vec[j], 1, ctx_tau, distances + j * n, prec);
}
}
}
else
{
/* distances are all set to zero; use one ellipsoid */
acb_theta_sum_a0_tilde(th, vec, nb, ctx_tau, distances, prec);
if (all)
{
acb_theta_sum_all_tilde(th, vec, nb, ctx_tau, distances, prec);
}
else
{
acb_theta_sum_a0_tilde(th, vec, nb, ctx_tau, distances, prec);
}
}

acb_theta_ctx_tau_clear(ctx_tau);
Expand All @@ -61,28 +75,32 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
slong nba = 1 << (g - s);
slong n0 = 1 << s;
slong n = 1 << g;
slong nbth = (all ? n * n : n);
acb_ptr * new_zs;
acb_ptr * cofactors;
slong * nb_pts;
slong ** pts;
slong * fullprec;
arf_struct * err;
acb_ptr th0, z0s;
slong nb0;
slong nb0, nbth0;
acb_mat_t tau0;
ulong a, a0;
slong j, k, index;
acb_t x;
slong j, a, index;
int res = 1;

new_zs = flint_malloc(nb * nba * sizeof(acb_ptr));
cofactors = flint_malloc(nb * nba * sizeof(acb_ptr));
nb_pts = flint_malloc(nb * nba * sizeof(slong));
err = flint_malloc(nb * nba * sizeof(arf_struct));
pts = flint_malloc(nb * nba * sizeof(slong *));
for (j = 0; j < nb * nba; j++)
{
arf_init(&err[j]);
}
fullprec = flint_malloc(nb * nba * sizeof(arf_struct));
acb_mat_window_init(tau0, tau, 0, 0, s, s);
acb_init(x);

nb0 = 1; /* always start with zero vector */
for (j = 0; j < nb; j++)
Expand All @@ -92,7 +110,7 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
if (res)
{
res = acb_theta_ql_lower_dim(&new_zs[j * nba + a], &cofactors[j * nba + a],
&nb_pts[j * nba + a], &err[j * nba + a], &fullprec[j * nba + a],
&pts[j * nba + a], &nb_pts[j * nba + a], &err[j * nba + a], &fullprec[j * nba + a],
zs + j * g, tau, distances + j * n, s, a, prec);
nb0 += nb_pts[j * nba + a];
}
Expand All @@ -103,6 +121,7 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
cofactors[j * nba + a] = _acb_vec_init(0);
nb_pts[j * nba + a] = 0;
fullprec[j * nba + a] = 0;
pts[j * nba + a] = flint_malloc(0);
}
}
}
Expand All @@ -112,7 +131,8 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
nb0 = 0;
}
z0s = _acb_vec_init(nb0 * s);
th0 = _acb_vec_init(nb0 * n0);
nbth0 = (all ? n0 * n0 : n0);
th0 = _acb_vec_init(nb0 * nbth0);

if (res)
{
Expand All @@ -128,39 +148,40 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
}

/* Call acb_theta_ql_exact in dimension s */
flint_printf("(ql_exact_lower_dim) calling ql_exact on %wd vectors in dimension %wd\n", nb0, s);
res = acb_theta_ql_exact(th0, z0s, nb0, tau0, pattern, 0, shifted_prec, prec);
/* flint_printf("(ql_exact_lower_dim) calling ql_exact on %wd vectors in dimension %wd\n", nb0, s); */
res = acb_theta_ql_exact(th0, z0s, nb0, tau0, pattern, all, shifted_prec, prec);
}

if (res)
{
/* Rescale using cofactors and sum into th */
/* Recombine into th */
index = 1; /* do not use result from zero vector */
_acb_vec_zero(th, nb * n);
_acb_vec_zero(th, nb * nbth);
for (j = 0; j < nb; j++)
{
for (a = 0; a < nba; a++)
{
for (k = 0; k < nb_pts[j * nba + a]; k++)
{
_acb_vec_scalar_mul(th0 + index * n0, th0 + index * n0,
n0, &cofactors[j * nba + a][k], prec);
for (a0 = 0; a0 < n0; a0++)
{
acb_add(&th[j * n + (a0 << (g - s)) + a],
&th[j * n + (a0 << (g - s)) + a],
&th0[index * n0 + a0], fullprec[j * nba + a]);
}
index += 1;
}
for (a0 = 0; a0 < n0; a0++)
{
acb_add_error_arf(&th[j * n + (a0 << (g - s)) + a],
&err[j * nba + a]);
}
acb_theta_ql_recombine(th + j * nbth, th0 + index * nbth0, cofactors[j * nba + a],
pts[j * nba + a], nb_pts[j * nba + a], &err[j * nba + a], fullprec[j * nba + a],
s, a, all, g, prec);
index += nb_pts[j * nba + a];
}
}
}
else
{
flint_printf("WARNING: ql_lower_dim failed, falling back to summation\n");
flint_printf("g = %wd, nb = %wd, tau:\n", g, nb);
acb_mat_printd(tau, 5);
flint_printf("zs:\n");
for (j = 0; j < nb; j++)
{
_acb_vec_printd(zs + j * g, g, 5);
}

acb_theta_ql_exact_sum(th, zs, nb, tau, distances, all, shifted_prec, prec);
res = 1;
}

/* Clear */
for (j = 0; j < nb; j++)
Expand All @@ -177,12 +198,15 @@ acb_theta_ql_exact_lower_dim(acb_ptr th, acb_srcptr zs, slong nb,
for (j = 0; j < nb * nba; j++)
{
arf_clear(&err[j]);
flint_free(pts[j]);
}
flint_free(pts);
flint_free(err);
flint_free(fullprec);
acb_mat_window_clear(tau0);
_acb_vec_clear(z0s, nb0 * s);
_acb_vec_clear(th0, nb0 * n0);
_acb_vec_clear(th0, nb0 * nbth0);
acb_clear(x);
return res;
}

Expand All @@ -194,7 +218,7 @@ acb_theta_ql_exact_steps(acb_ptr th, acb_srcptr zs, slong nb,
slong g = acb_mat_nrows(tau);
slong n = 1 << g;
slong nb_steps = pattern[g - 1];
acb_ptr rts, th_init, t;
acb_ptr rts, th_init, t, rts_all;
slong * easy_steps;
acb_mat_t new_tau;
slong guard, hp;
Expand All @@ -208,11 +232,15 @@ acb_theta_ql_exact_steps(acb_ptr th, acb_srcptr zs, slong nb,
th_init = _acb_vec_init(3 * n * nb);
easy_steps = flint_malloc(nb * sizeof(slong));
acb_mat_init(new_tau, g, g);
if (all)
{
rts_all = _acb_vec_init(nb * n * n);
}

res = acb_theta_ql_setup(rts, t, &guard, easy_steps, zs, nb, tau, distances,
res = acb_theta_ql_setup(rts, rts_all, t, &guard, easy_steps, zs, nb, tau, distances,
nb_steps, all, prec);

flint_printf("(ql_exact_steps) result of setup: %wd, easy_steps[0] = %wd\n", res, easy_steps[0]);
/* flint_printf("(ql_exact_steps) result of setup: %wd, easy_steps[0] = %wd\n", res, easy_steps[0]); */
hp = prec + nb_steps * guard;
acb_mat_scalar_mul_2exp_si(new_tau, tau, nb_steps);

Expand Down Expand Up @@ -374,7 +402,8 @@ acb_theta_ql_exact_steps(acb_ptr th, acb_srcptr zs, slong nb,

if (res) /* th_init is set; now perform steps */
{
acb_theta_ql_steps(th, th_init, rts, nb, nb_steps, distances, easy_steps, g, hp);
acb_theta_ql_steps(th, th_init, rts, rts_all, nb, nb_steps, distances,
easy_steps, all, g, hp);
}
else /* setup did not succeed: fall back to summation */
{
Expand All @@ -396,6 +425,10 @@ acb_theta_ql_exact_steps(acb_ptr th, acb_srcptr zs, slong nb,
_acb_vec_clear(th_init, 3 * n * nb);
flint_free(easy_steps);
acb_mat_clear(new_tau);
if (all)
{
_acb_vec_clear(rts_all, nb * n * n);
}
return res;
}

Expand Down Expand Up @@ -430,8 +463,6 @@ int acb_theta_ql_exact(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
}
}

flint_printf("(ql_exact) using nb_steps = %wd, split = %wd\n", nb_steps, split);

if (nb_steps == 0 && split == 0)
{
acb_theta_ql_exact_sum(th, zs, nb, tau, distances, all, shifted_prec, prec);
Expand Down
10 changes: 4 additions & 6 deletions src/acb_theta/ql_lower_dim.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,8 @@ acb_theta_ql_eld_points(slong ** pts, slong * nb_pts, arb_ptr v,
}

int
acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
arf_t err, slong * fullprec, acb_srcptr z, const acb_mat_t tau,
acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong ** pts,
slong * nb, arf_t err, slong * fullprec, acb_srcptr z, const acb_mat_t tau,
arb_srcptr distances, slong s, ulong a, slong prec)
{
slong g = acb_mat_nrows(tau);
Expand All @@ -79,7 +79,6 @@ acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
arb_ptr y, v, w, new_y, new_w;
acb_ptr u, x;
acb_t f;
slong * pts;
slong j, k;
int res;

Expand Down Expand Up @@ -107,7 +106,7 @@ acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
_acb_vec_get_imag(y, z, g);
arb_mat_vector_mul_col(w, Yinv, y, prec);

res = acb_theta_ql_eld_points(&pts, nb, v, fullprec,
res = acb_theta_ql_eld_points(pts, nb, v, fullprec,
err, distances, a, w, C, C1, prec);
*new_zs = _acb_vec_init((*nb) * s);
*cofactors = _acb_vec_init(*nb);
Expand All @@ -118,7 +117,7 @@ acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
acb_theta_char_get_acb(u, a, g - s);
for (j = 0; j < g - s; j++)
{
acb_add_si(&u[j], &u[j], pts[k * (g - s) + j], prec);
acb_add_si(&u[j], &u[j], (*pts)[k * (g - s) + j], prec);
}

/* Get new_z and cofactor at 0 */
Expand Down Expand Up @@ -152,6 +151,5 @@ acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong * nb,
_arb_vec_clear(new_y, s);
_arb_vec_clear(new_w, s);
acb_clear(f);
flint_free(pts);
return res;
}
Loading

0 comments on commit 6a892af

Please sign in to comment.