Skip to content

Commit

Permalink
Expose jet_exp_qf; use reduce_tau in jet_00
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-kieffer
j-kieffer committed Sep 4, 2024
1 parent 8c7917a commit 3bb95c5
Show file tree
Hide file tree
Showing 7 changed files with 126 additions and 211 deletions.
3 changes: 2 additions & 1 deletion src/acb_theta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ void acb_theta_jet_mul(acb_ptr res, acb_srcptr v1, acb_srcptr v2, slong ord,
void acb_theta_jet_compose(acb_ptr res, acb_srcptr v, const acb_mat_t N,
slong ord, slong prec);
void acb_theta_jet_exp_pi_i(acb_ptr res, arb_srcptr a, slong ord, slong g, slong prec);
void acb_theta_jet_exp_qf(acb_ptr res, acb_srcptr z, const acb_mat_t N, slong ord, slong prec);

void acb_theta_jet_error_bounds(arb_ptr err, acb_srcptr z, const acb_mat_t tau,
acb_srcptr dth, slong ord, slong prec);
Expand Down Expand Up @@ -280,7 +281,7 @@ void acb_theta_jet_ql_finite_diff(acb_ptr dth, const arf_t eps, const arf_t err,
/* Main functions */

int acb_theta_reduce_tau(acb_ptr new_zs, acb_mat_t new_tau, fmpz_mat_t mat, acb_mat_t N,
acb_ptr exps, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec);
acb_mat_t ct, acb_ptr exps, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec);
int acb_theta_reduce_z(acb_ptr res, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
slong nb, const acb_mat_t tau, slong prec);

Expand Down
6 changes: 4 additions & 2 deletions src/acb_theta/00.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ acb_theta_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong pre
{
slong g = acb_mat_nrows(tau);
fmpz_mat_t mat;
acb_mat_t new_tau, N;
acb_mat_t new_tau, N, ct;
acb_ptr new_zs, exps;
acb_ptr aux, units;
acb_t s;
Expand All @@ -34,13 +34,14 @@ acb_theta_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong pre
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);
aux = _acb_vec_init(nb);
units = _acb_vec_init(8);
acb_init(s);

res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, exps, zs, nb, tau, prec);
res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, ct, exps, zs, nb, tau, prec);

if (res)
{
Expand Down Expand Up @@ -68,6 +69,7 @@ acb_theta_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong pre
fmpz_mat_clear(mat);
acb_mat_clear(new_tau);
acb_mat_clear(N);
acb_mat_clear(ct);
_acb_vec_clear(new_zs, nb * g);
_acb_vec_clear(exps, nb);
_acb_vec_clear(aux, nb);
Expand Down
11 changes: 5 additions & 6 deletions src/acb_theta/all.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
slong g = acb_mat_nrows(tau);
slong n2 = 1 << (2 * g);
fmpz_mat_t mat;
acb_mat_t new_tau, N, c;
acb_mat_t new_tau, N, ct;
acb_ptr new_zs, exps, aux, units;
acb_t s;
ulong * image_ab;
Expand All @@ -37,7 +37,7 @@ acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
fmpz_mat_init(mat, 2 * g, 2 * g);
acb_mat_init(new_tau, g, g);
acb_mat_init(N, g, g);
acb_mat_init(c, g, g);
acb_mat_init(ct, g, g);
new_zs = _acb_vec_init(nb * g);
exps = _acb_vec_init(nb);
aux = _acb_vec_init(n2 * nb);
Expand All @@ -46,7 +46,7 @@ acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
e = flint_malloc(n2 * sizeof(slong));
acb_init(s);

res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, exps, zs, nb, tau, prec);
res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, ct, exps, zs, nb, tau, prec);

if (res)
{
Expand All @@ -56,8 +56,7 @@ acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
if (sqr)
{
kappa = acb_theta_transform_kappa2(mat);
acb_siegel_cocycle(c, mat, new_tau, prec);
acb_mat_det(s, c, prec);
acb_mat_det(s, ct, prec);
}
else
{
Expand Down Expand Up @@ -94,7 +93,7 @@ acb_theta_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,

fmpz_mat_clear(mat);
acb_mat_clear(new_tau);
acb_mat_clear(c);
acb_mat_clear(ct);
acb_mat_clear(N);
_acb_vec_clear(new_zs, nb * g);
_acb_vec_clear(exps, nb);
Expand Down
121 changes: 11 additions & 110 deletions src/acb_theta/jet_00.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,101 +13,19 @@
#include "acb_mat.h"
#include "acb_theta.h"

/* Compute jet of exp (z^T N z) */
static void
acb_theta_jet_exp_qf(acb_ptr res, acb_srcptr z, const acb_mat_t N, slong ord, slong prec)
{
slong g = acb_mat_nrows(N);
slong nb = acb_theta_jet_nb(ord, g);
acb_mat_t tp;
acb_poly_t pol;
acb_ptr aux;
acb_ptr y;
acb_t c;
slong * tup;
slong j, k, l, i;

acb_mat_init(tp, g, g);
acb_poly_init(pol);
aux = _acb_vec_init(nb);
y = _acb_vec_init(g);
acb_init(c);
tup = flint_malloc(g * sizeof(slong));

/* exp((z+h)^T N (z+h)) = exp(z^T N z) exp(z^T (N+N^T) h) exp(h^T N h) */
_acb_vec_zero(res, nb);
acb_mat_vector_mul_col(y, N, z, prec);
acb_dot(&res[0], NULL, 0, z, 1, y, 1, g, prec);
acb_exp(&res[0], &res[0], prec);

acb_mat_transpose(tp, N);
acb_mat_add(tp, tp, N, prec);
acb_mat_vector_mul_row(y, z, tp, prec);
for (j = 0; j < g; j++)
{
_acb_vec_zero(aux, nb);
acb_poly_zero(pol);
acb_poly_set_coeff_acb(pol, 1, &y[j]);
acb_poly_exp_series(pol, pol, ord + 1, prec);
for (l = 0; l <= ord; l++)
{
for (i = 0; i < g; i++)
{
tup[i] = 0;
}
tup[j] = l;
acb_poly_get_coeff_acb(&aux[acb_theta_jet_index(tup, g)], pol, l);
}
acb_theta_jet_mul(res, res, aux, ord, g, prec);
}

for (j = 0; j < g; j++)
{
for (k = j; k < g; k++)
{
_acb_vec_zero(aux, nb);
acb_poly_zero(pol);
acb_add(c, acb_mat_entry(N, k, j), acb_mat_entry(N, j, k), prec);
if (j == k)
{
acb_mul_2exp_si(c, c, -1);
}
acb_poly_set_coeff_acb(pol, 1, c);
acb_poly_exp_series(pol, pol, (ord / 2) + 1, prec);
for (l = 0; l <= (ord / 2); l++)
{
for (i = 0; i < g; i++)
{
tup[i] = 0;
}
tup[j] += l;
tup[k] += l;
acb_poly_get_coeff_acb(&aux[acb_theta_jet_index(tup, g)], pol, l);
}
acb_theta_jet_mul(res, res, aux, ord, g, prec);
}
}

acb_mat_clear(tp);
acb_poly_clear(pol);
_acb_vec_clear(aux, nb);
_acb_vec_clear(y, g);
acb_clear(c);
flint_free(tup);
}

void
acb_theta_jet_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
slong ord, slong prec)
{
slong g = acb_mat_nrows(tau);
slong nbth = acb_theta_jet_nb(ord, g);
fmpz_mat_t mat, gamma;
acb_mat_t new_tau, c, cinv, N;
acb_ptr new_zs, aux, units;
fmpz_mat_t mat;
acb_mat_t new_tau, N, ct;
acb_ptr new_zs, aux, exps, units;
acb_t s, t;
ulong image_ab;
slong kappa, e, j;
int res;

if (nb <= 0)
{
Expand All @@ -116,48 +34,32 @@ acb_theta_jet_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,

fmpz_mat_init(mat, 2 * g, 2 * g);
acb_mat_init(new_tau, g, g);
acb_mat_init(c, g, g);
acb_mat_init(cinv, 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);
aux = _acb_vec_init(nb * nbth);
units = _acb_vec_init(8);
acb_init(s);
acb_init(t);

acb_siegel_reduce(mat, tau, prec);
acb_siegel_transform_cocycle_inv(new_tau, c, cinv, mat, tau, prec);
_acb_vec_unit_roots(units, 8, 8, prec);
res = acb_theta_reduce_tau(new_zs, new_tau, mat, N, ct, exps, zs, nb, tau, prec);

acb_mat_transpose(cinv, cinv);
for (j = 0; j < nb; j++)
{
acb_mat_vector_mul_col(new_zs + j * g, cinv, zs + j * g, prec);
}

if (acb_siegel_is_reduced(new_tau, -10, prec))
if (res)
{
/* todo: reduce z here */

sp2gz_inv(mat, mat);
_acb_vec_unit_roots(units, 8, 8, prec);
kappa = acb_theta_transform_kappa(s, mat, new_tau, prec);
image_ab = acb_theta_transform_char(&e, mat, 0);

fmpz_mat_window_init(gamma, mat, g, 0, 2 * g, g);
acb_mat_set_fmpz_mat(N, gamma);
acb_mat_mul(N, N, cinv, prec);
acb_const_pi(t, prec);
acb_mul_onei(t, t);
acb_mat_scalar_mul_acb(N, N, t, prec);
fmpz_mat_window_clear(gamma);

acb_theta_jet_one_notransform(aux, new_zs, nb, new_tau, ord, image_ab, prec);

for (j = 0; j < nb; j++)
{
acb_mul(t, s, &units[(kappa + e) % 8], prec);
_acb_vec_scalar_mul(th + j * nbth, aux + j * nbth, nbth, t, prec);
acb_theta_jet_compose(th + j * nbth, th + j * nbth, cinv, ord, prec);
acb_theta_jet_compose(th + j * nbth, th + j * nbth, ct, ord, prec);
acb_theta_jet_exp_qf(aux + j * nbth, zs + j * g, N, ord, prec);
acb_theta_jet_mul(th + j * nbth, th + j * nbth, aux + j * nbth, ord, g, prec);
}
Expand All @@ -169,9 +71,8 @@ acb_theta_jet_00(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,

fmpz_mat_clear(mat);
acb_mat_clear(new_tau);
acb_mat_clear(c);
acb_mat_clear(cinv);
acb_mat_clear(N);
acb_mat_clear(ct);
_acb_vec_clear(new_zs, nb * g);
_acb_vec_clear(aux, nb * nbth);
_acb_vec_clear(units, 8);
Expand Down
83 changes: 0 additions & 83 deletions src/acb_theta/jet_all.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,89 +13,6 @@
#include "acb_mat.h"
#include "acb_theta.h"

/* Compute jet of exp (z^T N z) */
static void
acb_theta_jet_exp_qf(acb_ptr res, acb_srcptr z, const acb_mat_t N, slong ord, slong prec)
{
slong g = acb_mat_nrows(N);
slong nb = acb_theta_jet_nb(ord, g);
acb_mat_t tp;
acb_poly_t pol;
acb_ptr aux;
acb_ptr y;
acb_t c;
slong * tup;
slong j, k, l, i;

acb_mat_init(tp, g, g);
acb_poly_init(pol);
aux = _acb_vec_init(nb);
y = _acb_vec_init(g);
acb_init(c);
tup = flint_malloc(g * sizeof(slong));

/* exp((z+h)^T N (z+h)) = exp(z^T N z) exp(z^T (N+N^T) h) exp(h^T N h) */
_acb_vec_zero(res, nb);
acb_mat_vector_mul_col(y, N, z, prec);
acb_dot(&res[0], NULL, 0, z, 1, y, 1, g, prec);
acb_exp(&res[0], &res[0], prec);

acb_mat_transpose(tp, N);
acb_mat_add(tp, tp, N, prec);
acb_mat_vector_mul_row(y, z, tp, prec);
for (j = 0; j < g; j++)
{
_acb_vec_zero(aux, nb);
acb_poly_zero(pol);
acb_poly_set_coeff_acb(pol, 1, &y[j]);
acb_poly_exp_series(pol, pol, ord + 1, prec);
for (l = 0; l <= ord; l++)
{
for (i = 0; i < g; i++)
{
tup[i] = 0;
}
tup[j] = l;
acb_poly_get_coeff_acb(&aux[acb_theta_jet_index(tup, g)], pol, l);
}
acb_theta_jet_mul(res, res, aux, ord, g, prec);
}

for (j = 0; j < g; j++)
{
for (k = j; k < g; k++)
{
_acb_vec_zero(aux, nb);
acb_poly_zero(pol);
acb_add(c, acb_mat_entry(N, k, j), acb_mat_entry(N, j, k), prec);
if (j == k)
{
acb_mul_2exp_si(c, c, -1);
}
acb_poly_set_coeff_acb(pol, 1, c);
acb_poly_exp_series(pol, pol, (ord / 2) + 1, prec);
for (l = 0; l <= (ord / 2); l++)
{
for (i = 0; i < g; i++)
{
tup[i] = 0;
}
tup[j] += l;
tup[k] += l;
acb_poly_get_coeff_acb(&aux[acb_theta_jet_index(tup, g)], pol, l);
}
acb_theta_jet_mul(res, res, aux, ord, g, prec);
}
}

acb_mat_clear(tp);
acb_poly_clear(pol);
_acb_vec_clear(aux, nb);
_acb_vec_clear(y, g);
acb_clear(c);
flint_free(tup);
}

void
acb_theta_jet_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
slong ord, slong prec)
Expand Down
Loading

0 comments on commit 3bb95c5

Please sign in to comment.