Write t-reduce_z; further comments on code coverage

src/acb_theta/eld_distances.c

@@ -15,6 +15,7 @@
 #include "acb_mat.h"
 #include "acb_theta.h"
 
+/* Should not be used in tests */
 static void
 acb_theta_dist_unif(arb_t d, const arb_mat_t cho, slong prec)
 {

src/acb_theta/g2_sextic_chi5.c

@@ -62,6 +62,7 @@ acb_theta_g2_sextic_chi5(acb_poly_t f, acb_t chi5, const acb_mat_t tau, slong pr
     }
     else
     {
+        /* Should not happen in tests */
         acb_indeterminate(chi5);
         for (k = 0; k < 6; k++)
         {

src/acb_theta/ql_exact.c

@@ -543,8 +543,10 @@ acb_theta_ql_exact_steps(acb_ptr th, acb_srcptr zs, slong nb,
         acb_theta_ql_perform_steps(th, th_init, rts, rts_all, nb, nb_steps,
             distances, easy_steps, all, g, prec + guard);
     }
-    else /* setup or intput did not succeed, fall back to summation */
+    else
     {
+        /* Setup or intput did not succeed, fall back to summation */
+        /* Should not happen in tests */
         acb_theta_ql_exact_sum(th, zs, nb, tau, distances, all, 1, prec);
     }
 

src/acb_theta/ql_nb_steps.c

@@ -97,6 +97,7 @@ acb_theta_ql_nb_steps(slong * pattern, const acb_mat_t tau, int cst, slong prec)
       dimensions will depend on whether or not duplications/dimension-lowerings
       have already been applied. */
     /* See /path/to/flint/build/acb_theta/profile/p-acb_theta_ql_exact */
+    /* Some of these branches will not show up in tests. */
     for (s = g - 1; s >= 0; s--)
     {
         /* Find out how many duplication steps have been performed so far
@@ -134,7 +135,6 @@ acb_theta_ql_nb_steps(slong * pattern, const acb_mat_t tau, int cst, slong prec)
            the dimension-lowering strategy. */
         if (s == 0)
         {
-            /* Some of these branches will not show up in tests. */
             nb = pattern[s] - 5;
             if (nb >= 10)
             {

src/acb_theta/reduce_z.c

@@ -46,6 +46,7 @@ acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
     arb_mat_t cho, yinv;
     arb_ptr y;
     acb_ptr t, x;
+    arf_t e;
     slong j, k;
     int res = 1;
 
@@ -56,6 +57,7 @@ acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
     y = _arb_vec_init(g);
     t = _acb_vec_init(g);
     x = _acb_vec_init(g);
+    arf_init(e);
 
     acb_siegel_cho_yinv(cho, yinv, tau, prec);
 
@@ -98,14 +100,14 @@ acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
                 acb_sub_arb(&x[k], &x[k], &y[k], prec);
             }
         }
-        else
+
+        /* Set real part to [-1,1] if error is too large */
+        for (k = 0; k < g; k++)
         {
-            /* Still OK; set real part to [-1,1] */
-            for (k = 0; k < g; k++)
+            if (mag_cmp_2exp_si(arb_radref(acb_realref(&x[k])), 0) > 0)
             {
                 arb_zero_pm_one(acb_realref(&x[k]));
             }
-            res = 1;
         }
 
         /* Set new_z */
@@ -117,5 +119,6 @@ acb_theta_reduce_z(acb_ptr new_zs, arb_ptr rs, acb_ptr cs, acb_srcptr zs,
     _arb_vec_clear(y, g);
     _acb_vec_clear(t, g);
     _acb_vec_clear(x, g);
+    arf_clear(e);
     return res;
 }

src/acb_theta/siegel_randtest_vec.c

@@ -34,6 +34,9 @@ acb_siegel_randtest_vec(acb_ptr z, flint_rand_t state, slong g, slong prec)
         case 3:
             acb_randtest(&z[k], state, prec, 20);
             break;
+        case 4:
+            arb_urandom(acb_imagref(&z[k]), state, prec);
+            arb_randtest(acb_realref(&z[k]), state, prec, mag_bits);
         default:
             acb_urandom(&z[k], state, prec);
         }

src/acb_theta/test/main.c

@@ -46,6 +46,7 @@
 #include "t-ql_local_bound.c"
 #include "t-ql_lower_dim.c"
 #include "t-ql_setup.c"
+#include "t-reduce_z.c"
 #include "t-siegel_cocycle.c"
 #include "t-siegel_is_reduced.c"
 #include "t-siegel_kappa.c"
@@ -99,6 +100,7 @@ test_struct tests[] =
     TEST_FUNCTION(acb_theta_ql_local_bound),
     TEST_FUNCTION(acb_theta_ql_lower_dim),
     TEST_FUNCTION(acb_theta_ql_setup),
+    TEST_FUNCTION(acb_theta_reduce_z),
     TEST_FUNCTION(acb_theta_siegel_cocycle),
     TEST_FUNCTION(acb_theta_siegel_is_reduced),
     TEST_FUNCTION(acb_theta_siegel_kappa),

src/acb_theta/test/t-reduce_z.c

@@ -0,0 +1,95 @@
+/*
+    Copyright (C) 2025 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 <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "fmpz.h"
+#include "arb.h"
+#include "acb.h"
+#include "arb_mat.h"
+#include "acb_mat.h"
+#include "acb_theta.h"
+
+TEST_FUNCTION_START(acb_theta_reduce_z, state)
+{
+    slong iter;
+
+    /* Test: entries of r are even integers, im(new_zs) is im(zs - tau * r) */
+    for (iter = 0; iter < 100 * flint_test_multiplier(); iter++)
+    {
+        slong g = 1 + n_randint(state, 3);
+        slong nb = n_randint(state, 3);
+        slong prec = 100 + n_randint(state, 200);
+        slong bits = n_randint(state, 4);
+        acb_ptr zs, new_zs, cs;
+        acb_mat_t tau;
+        arb_mat_t im;
+        arb_ptr rs, test, y;
+        fmpz_t x;
+        slong j;
+        int res;
+
+        acb_mat_init(tau, g, g);
+        arb_mat_init(im, g, g);
+        zs = _acb_vec_init(nb * g);
+        new_zs = _acb_vec_init(nb * g);
+        cs = _acb_vec_init(nb);
+        test = _arb_vec_init(g);
+        y = _arb_vec_init(g);
+        rs = _arb_vec_init(nb * g);
+        fmpz_init(x);
+
+        acb_siegel_randtest_reduced(tau, state, prec, bits);
+        acb_siegel_randtest_vec(zs, state, nb * g, prec);
+        acb_mat_get_imag(im, tau);
+
+        res = acb_theta_reduce_z(new_zs, rs, cs, zs, nb, tau, prec);
+
+        if (res)
+        {
+            for (j = 0; j < g * nb; j++)
+            {
+                res = arb_get_unique_fmpz(x, &rs[j]);
+                if (!res || !arb_is_exact(&rs[j]) || fmpz_mod_ui(x, x, 2) != 0)
+                {
+                    flint_printf("FAIL (integers)\n");
+                    flint_abort();
+                }
+            }
+
+            for (j = 0; j < nb; j++)
+            {
+                arb_mat_vector_mul_col(test, im, rs + j * g, prec);
+                _acb_vec_get_imag(y, zs + j * g, g);
+                _arb_vec_sub(test, y, test, g, prec);
+                _acb_vec_get_imag(y, new_zs + j * g, g);
+                if (!_arb_vec_overlaps(test, y, g))
+                {
+                    flint_printf("FAIL (overlap)\n");
+                    _arb_vec_printd(y, g, 5);
+                    _arb_vec_printd(test, g, 5);
+                    flint_abort();
+                }
+            }
+        }
+
+        acb_mat_clear(tau);
+        arb_mat_clear(im);
+        _acb_vec_clear(zs, nb * g);
+        _acb_vec_clear(new_zs, nb * g);
+        _acb_vec_clear(cs, nb);
+        _arb_vec_clear(test, g);
+        _arb_vec_clear(y, g);
+        _arb_vec_clear(rs, nb * g);
+        fmpz_clear(x);
+    }
+
+    TEST_FUNCTION_END(state);
+}