include/flint/acb_theta.h

0001 /*
0002     Copyright (C) 2023 Jean Kieffer
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef ACB_THETA_H
0013 #define ACB_THETA_H
0014
0015 #ifdef ACB_THETA_INLINES_C
0016 #define ACB_THETA_INLINE
0017 #else
0018 #define ACB_THETA_INLINE static inline
0019 #endif
0020
0021 #include "fmpz_types.h"
0022 #include "acb_types.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 #define ACB_THETA_LOW_PREC 32
0029
0030 /* The Siegel modular group */
0031
0032 FLINT_FORCE_INLINE slong
0033 sp2gz_dim(const fmpz_mat_t mat)
0034 {
0035     return (mat->r) / 2;
0036 }
0037
0038 void sp2gz_set_blocks(fmpz_mat_t mat, const fmpz_mat_t alpha, const fmpz_mat_t beta,
0039     const fmpz_mat_t gamma, const fmpz_mat_t delta);
0040 void sp2gz_j(fmpz_mat_t mat);
0041 void sp2gz_block_diag(fmpz_mat_t mat, const fmpz_mat_t U);
0042 void sp2gz_trig(fmpz_mat_t mat, const fmpz_mat_t S);
0043 void sp2gz_embed(fmpz_mat_t res, const fmpz_mat_t mat);
0044 void sp2gz_restrict(fmpz_mat_t res, const fmpz_mat_t mat);
0045
0046 FLINT_FORCE_INLINE slong
0047 sp2gz_nb_fundamental(slong g)
0048 {
0049     if (g == 1)
0050         return 1;
0051     if (g == 2)
0052         return 19;
0053     else
0054         return 19 * ((g * (g - 1)) / 2) + (1 << g);
0055 }
0056
0057 void sp2gz_fundamental(fmpz_mat_t mat, slong j);
0058
0059 int sp2gz_is_correct(const fmpz_mat_t mat);
0060 int sp2gz_is_j(const fmpz_mat_t mat);
0061 int sp2gz_is_block_diag(const fmpz_mat_t mat);
0062 int sp2gz_is_trig(const fmpz_mat_t mat);
0063 int sp2gz_is_embedded(fmpz_mat_t res, const fmpz_mat_t mat);
0064
0065 void sp2gz_inv(fmpz_mat_t inv, const fmpz_mat_t mat);
0066 fmpz_mat_struct * sp2gz_decompose(slong * nb, const fmpz_mat_t mat);
0067
0068 void sp2gz_randtest(fmpz_mat_t mat, flint_rand_t state, slong bits);
0069
0070 /* The Siegel half space */
0071
0072 void acb_siegel_cocycle(acb_mat_t c, const fmpz_mat_t mat, const acb_mat_t tau, slong prec);
0073 void acb_siegel_transform_cocycle_inv(acb_mat_t w, acb_mat_t c, acb_mat_t cinv,
0074     const fmpz_mat_t mat, const acb_mat_t tau, slong prec);
0075 void acb_siegel_transform(acb_mat_t w, const fmpz_mat_t mat, const acb_mat_t tau, slong prec);
0076 void acb_siegel_cho_yinv(arb_mat_t cho, arb_mat_t yinv, const acb_mat_t tau, slong prec);
0077
0078 void acb_siegel_reduce(fmpz_mat_t mat, const acb_mat_t tau, slong prec);
0079 int acb_siegel_is_reduced(const acb_mat_t tau, slong tol_exp, slong prec);
0080
0081 slong acb_siegel_kappa(acb_t sqrtdet, const fmpz_mat_t mat, const acb_mat_t tau, int sqr, slong prec);
0082 slong acb_siegel_kappa2(const fmpz_mat_t mat);
0083
0084 void acb_siegel_randtest(acb_mat_t tau, flint_rand_t state, slong prec, slong mag_bits);
0085 void acb_siegel_randtest_reduced(acb_mat_t tau, flint_rand_t state, slong prec, slong mag_bits);
0086 void acb_siegel_randtest_compact(acb_mat_t tau, flint_rand_t state, int exact, slong prec);
0087 void acb_siegel_randtest_vec(acb_ptr z, flint_rand_t state, slong g, slong prec);
0088 void acb_siegel_randtest_vec_reduced(acb_ptr zs, flint_rand_t state,
0089     slong nb, const acb_mat_t tau, int exact, slong prec);
0090
0091 /* Theta characteristics */
0092
0093 FLINT_FORCE_INLINE int
0094 acb_theta_char_bit(ulong ch, slong j, slong n)
0095 {
0096     return (ch >> (n - 1 - j)) & 1;
0097 }
0098
0099 void acb_theta_char_get_arb(arb_ptr v, ulong a, slong g);
0100 void acb_theta_char_get_acb(acb_ptr v, ulong a, slong g);
0101 ulong acb_theta_char_set_slong_vec(const slong * vec, slong len);
0102
0103 slong acb_theta_char_dot(ulong a, ulong b, slong g);
0104 slong acb_theta_char_dot_slong(ulong a, const slong * n, slong g);
0105
0106 FLINT_FORCE_INLINE int
0107 acb_theta_char_is_even(ulong ab, slong g)
0108 {
0109     return acb_theta_char_dot(ab >> g, ab, g) % 2 == 0;
0110 }
0111
0112 void acb_theta_char_table(ulong * ch, slong * e, const fmpz_mat_t mat, ulong ab, int all);
0113 void acb_theta_char_shuffle(acb_ptr res, const fmpz_mat_t mat, acb_srcptr th,
0114     int sqr, slong prec);
0115
0116 /* Toolbox for derivatives */
0117
0118 slong acb_theta_jet_nb(slong ord, slong g);
0119 slong acb_theta_jet_total_order(const slong * tup, slong g);
0120 void acb_theta_jet_tuples(slong * tups, slong ord, slong g);
0121 slong acb_theta_jet_index(const slong * tup, slong g);
0122
0123 void acb_theta_jet_mul(acb_ptr res, acb_srcptr v1, acb_srcptr v2, slong ord,
0124     slong g, slong prec);
0125 void acb_theta_jet_compose(acb_ptr res, acb_srcptr v, const acb_mat_t N,
0126     slong ord, slong prec);
0127 void acb_theta_jet_exp_pi_i(acb_ptr res, arb_srcptr a, slong ord, slong g, slong prec);
0128 void acb_theta_jet_exp_qf(acb_ptr res, acb_srcptr z, const acb_mat_t N, slong ord, slong prec);
0129
0130 /* Ellipsoids */
0131
0132 struct acb_theta_eld_struct
0133 {
0134     slong dim, ambient_dim;
0135     slong * last_coords;
0136     slong min, mid, max, nr, nl;
0137     struct acb_theta_eld_struct * rchildren;
0138     struct acb_theta_eld_struct * lchildren;
0139     slong nb_pts, nb_border;
0140     slong * box;
0141 };
0142
0143 typedef struct acb_theta_eld_struct acb_theta_eld_t[1];
0144
0145 void acb_theta_eld_init(acb_theta_eld_t E, slong d, slong g);
0146 void acb_theta_eld_clear(acb_theta_eld_t E);
0147
0148 int acb_theta_eld_set(acb_theta_eld_t E, const arb_mat_t C, const arf_t R2, arb_srcptr v);
0149
0150 FLINT_FORCE_INLINE slong
0151 acb_theta_eld_nb_pts(const acb_theta_eld_t E)
0152 {
0153     return E->nb_pts;
0154 }
0155
0156 void acb_theta_eld_points(slong * pts, const acb_theta_eld_t E);
0157
0158 FLINT_FORCE_INLINE slong
0159 acb_theta_eld_box(const acb_theta_eld_t E, slong j)
0160 {
0161     return E->box[j];
0162 }
0163
0164 FLINT_FORCE_INLINE slong
0165 acb_theta_eld_nb_border(const acb_theta_eld_t E)
0166 {
0167     return E->nb_border;
0168 }
0169
0170 void acb_theta_eld_border(slong * pts, const acb_theta_eld_t E);
0171 int acb_theta_eld_contains(const acb_theta_eld_t E, const slong * pt);
0172 void acb_theta_eld_print(const acb_theta_eld_t E);
0173
0174 void acb_theta_eld_distances(arb_ptr ds, acb_srcptr zs, slong nb,
0175     const acb_mat_t tau, slong prec);
0176
0177 /* Error bounds in summation algorithms */
0178
0179 void acb_theta_sum_radius(arf_t R2, arf_t eps, const arb_mat_t cho, slong ord, slong prec);
0180 void acb_theta_sum_jet_radius(arf_t R2, arf_t eps, const arb_mat_t cho, arb_srcptr v,
0181     slong ord, slong prec);
0182 void acb_theta_sum_term(acb_t res, acb_srcptr z, const acb_mat_t tau, const slong * tup,
0183     const slong * n, slong prec);
0184 slong acb_theta_sum_addprec(const arb_t d);
0185
0186 /* Context structures in summation algorithms */
0187
0188 struct acb_theta_ctx_tau_struct
0189 {
0190     slong g;
0191     int allow_shift;
0192     arb_mat_struct yinv;
0193     arb_mat_struct cho;
0194
0195     acb_mat_t exp_tau_div_4;
0196     acb_mat_t exp_tau_div_2;
0197     acb_mat_t exp_tau;
0198     acb_mat_t exp_tau_div_4_inv;
0199     acb_mat_t exp_tau_div_2_inv;
0200     acb_mat_t exp_tau_inv;
0201
0202     /* allow_shift only */
0203     acb_ptr exp_tau_a;
0204     acb_ptr exp_tau_a_inv;
0205     acb_ptr exp_a_tau_a_div_4;
0206 };
0207
0208 typedef struct acb_theta_ctx_tau_struct acb_theta_ctx_tau_t[1];
0209
0210 typedef struct
0211 {
0212     slong g;
0213     acb_ptr exp_z;
0214     acb_ptr exp_2z;
0215     acb_ptr exp_z_inv;
0216     acb_ptr exp_2z_inv;
0217     arb_ptr v;
0218     arb_struct u;
0219     arb_struct uinv;
0220     int is_real;
0221 }
0222 acb_theta_ctx_z_struct;
0223
0224 typedef acb_theta_ctx_z_struct acb_theta_ctx_z_t[1];
0225
0226 void acb_theta_ctx_tau_init(acb_theta_ctx_tau_t ctx, int allow_shift, slong g);
0227 void acb_theta_ctx_tau_clear(acb_theta_ctx_tau_t ctx);
0228 void acb_theta_ctx_z_init(acb_theta_ctx_z_t ctx, slong g);
0229 void acb_theta_ctx_z_clear(acb_theta_ctx_z_t ctx);
0230 acb_theta_ctx_z_struct * acb_theta_ctx_z_vec_init(slong nb, slong g);
0231 void acb_theta_ctx_z_vec_clear(acb_theta_ctx_z_struct * vec, slong nb);
0232
0233 void acb_theta_ctx_exp_inv(acb_t exp_inv, const acb_t exp, const acb_t x, int is_real, slong prec);
0234 void acb_theta_ctx_sqr_inv(acb_t sqr_inv, const acb_t inv, const acb_t sqr, int is_real, slong prec);
0235
0236 void acb_theta_ctx_tau_set(acb_theta_ctx_tau_t ctx, const acb_mat_t tau, slong prec);
0237 void acb_theta_ctx_tau_dupl(acb_theta_ctx_tau_t ctx, slong prec);
0238 int acb_theta_ctx_tau_overlaps(const acb_theta_ctx_tau_t ctx1, const acb_theta_ctx_tau_t ctx2);
0239
0240 void acb_theta_ctx_z_set(acb_theta_ctx_z_t ctx, acb_srcptr z, const acb_theta_ctx_tau_t ctx_tau, slong prec);
0241 void acb_theta_ctx_z_dupl(acb_theta_ctx_z_t ctx, slong prec);
0242 void acb_theta_ctx_z_add_real(acb_theta_ctx_z_t res, const acb_theta_ctx_z_t ctx,
0243     const acb_theta_ctx_z_t ctx_real, slong prec);
0244 void acb_theta_ctx_z_common_v(arb_ptr v, const acb_theta_ctx_z_struct * vec, slong nb, slong prec);
0245 int acb_theta_ctx_z_overlaps(const acb_theta_ctx_z_t ctx1, const acb_theta_ctx_z_t ctx2);
0246
0247 /* Summation algorithms */
0248
0249 typedef void (*acb_theta_sum_worker_t)(acb_ptr, acb_srcptr, acb_srcptr, const slong *,
0250     slong, const acb_t, const slong *, slong, slong, slong, slong);
0251
0252 void acb_theta_sum_sqr_pow(acb_ptr * sqr_pow, const acb_mat_t exp_tau, const acb_theta_eld_t E, slong prec);
0253 void acb_theta_sum_work(acb_ptr th, slong len, acb_srcptr exp_z, acb_srcptr exp_z_inv,
0254     const acb_mat_t exp_tau, const acb_mat_t exp_tau_inv, const acb_ptr * sqr_pow,
0255     const acb_theta_eld_t E, slong ord, slong prec, acb_theta_sum_worker_t worker);
0256 void acb_theta_sum(acb_ptr th, const acb_theta_ctx_z_struct * vec, slong nb,
0257     const acb_theta_ctx_tau_t ctx_tau, arb_srcptr distances, int all_a,
0258     int all_b, int tilde, slong prec);
0259 void acb_theta_sum_jet(acb_ptr th, const acb_theta_ctx_z_struct * vec, slong nb,
0260     const acb_theta_ctx_tau_t ctx_tau, slong ord, int all_a, int all_b, slong prec);
0261
0262 /* AGM steps */
0263
0264 void acb_theta_agm_sqrt(acb_ptr res, acb_srcptr a, acb_srcptr roots, slong nb, slong prec);
0265 void acb_theta_agm_mul(acb_ptr res, acb_srcptr a1, acb_srcptr a2, slong g, int all, slong prec);
0266 void acb_theta_agm_mul_tight(acb_ptr res, acb_srcptr a0, acb_srcptr a,
0267     arb_srcptr d0, arb_srcptr d, slong g, int all, slong prec);
0268
0269 /* Quasilinear algorithms on reduced input */
0270
0271 int acb_theta_ql_nb_steps(slong * pattern, const acb_mat_t tau, int cst, slong prec);
0272
0273 int acb_theta_ql_lower_dim(acb_ptr * new_zs, acb_ptr * cofactors, slong ** pts,
0274     slong * nb, arf_t err, slong * fullprec, acb_srcptr z, const acb_mat_t tau,
0275     arb_srcptr distances, slong s, ulong a, slong prec);
0276 void acb_theta_ql_recombine(acb_ptr th, acb_srcptr th0, acb_srcptr cofactors,
0277     const slong * pts, slong nb, const arf_t err, slong fullprec,
0278     slong s, ulong a, int all, slong g, slong prec);
0279 int acb_theta_ql_setup(acb_ptr rts, acb_ptr rts_all, acb_ptr t, slong * guard, slong * easy_steps,
0280     acb_srcptr zs, slong nb, const acb_mat_t tau, arb_srcptr distances,
0281     slong nb_steps, int all, slong prec);
0282 void acb_theta_ql_exact(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
0283     const slong * pattern, int all, int shifted_prec, slong prec);
0284
0285 void acb_theta_ql_local_bound(arb_t c, arb_t rho, acb_srcptr z, const acb_mat_t tau, slong ord);
0286 void acb_theta_ql_jet_error(arb_ptr err, acb_srcptr z, const acb_mat_t tau,
0287     acb_srcptr dth, slong ord, slong prec);
0288
0289 void acb_theta_ql_jet_fd(acb_ptr th, acb_srcptr zs, slong nb,
0290     const acb_mat_t tau, slong ord, int all, slong prec);
0291 void acb_theta_ql_jet(acb_ptr th, acb_srcptr zs, slong nb,
0292     const acb_mat_t tau, slong ord, int all, slong prec);
0293
0294 /* Reduction and main functions */
0295
0296 void acb_theta_jet_notransform(acb_ptr th, acb_srcptr zs, slong nb,
0297     const acb_mat_t tau, slong ord, ulong ab, int all, int sqr, slong prec);
0298
0299 int acb_theta_reduce_tau(acb_ptr new_zs, acb_mat_t new_tau, fmpz_mat_t mat, acb_mat_t N,
0300     acb_mat_t ct, acb_ptr exps, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec);
0301 int acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
0302     slong nb, const acb_mat_t tau, slong prec);
0303
0304 void acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb,
0305     const acb_mat_t tau, slong ord, ulong ab, int all, int sqr, slong prec);
0306
0307 ACB_THETA_INLINE void
0308 acb_theta_one(acb_ptr th, acb_srcptr z, const acb_mat_t tau, ulong ab, slong prec)
0309 {
0310     acb_theta_jet(th, z, 1, tau, 0, ab, 0, 0, prec);
0311 }
0312
0313 ACB_THETA_INLINE void
0314 acb_theta_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)
0315 {
0316     acb_theta_jet(th, z, 1, tau, 0, 0, 1, sqr, prec);
0317 }
0318
0319 /* Dimension 2 specifics */
0320
0321 void acb_theta_g2_detk_symj(acb_poly_t res, const acb_mat_t m, const acb_poly_t f,
0322     slong k, slong j, slong prec);
0323 void acb_theta_g2_transvectant(acb_poly_t res, const acb_poly_t g, const acb_poly_t h,
0324     slong m, slong n, slong k, int lead, slong prec);
0325 slong acb_theta_g2_character(const fmpz_mat_t mat);
0326
0327 void acb_theta_g2_even_weight(acb_t psi4, acb_t psi6, acb_t chi10, acb_t chi12,
0328     acb_srcptr th2, slong prec);
0329 void acb_theta_g2_chi5(acb_t res, acb_srcptr th, slong prec);
0330 void acb_theta_g2_chi35(acb_t res, acb_srcptr th, slong prec);
0331 void acb_theta_g2_chi3_6(acb_poly_t res, acb_srcptr dth, slong prec);
0332
0333 void acb_theta_g2_sextic_chi5(acb_poly_t f, acb_t chi5, const acb_mat_t tau, slong prec);
0334 void acb_theta_g2_covariants(acb_poly_struct * res, const acb_poly_t f, int lead, slong prec);
0335
0336 #ifdef __cplusplus
0337 }
0338 #endif
0339
0340 #endif