Fix bug in t-spd_lll_reduce

doc/source/arb_mat.rst

@@ -796,19 +796,18 @@ LLL reduction
     is not symmetric or the result of rounding is not a positive definite
     matrix. The warnings of :func:`arf_get_fmpz` apply.
 
-.. function:: void arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, slong prec)
+.. function:: void arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, double delta, double eta, slong prec)
 
     Given a symmetric positive definite matrix *A*, sets *U* to an invertible
-    matrix such that `U^T A U` is close to being LLL-reduced. If
-    :func:`arb_mat_spd_get_fmpz_mat` succeeds at the chosen precision, we call
-    :func:`fmpz_lll`, and otherwise set *U* to the identity matrix. The
-    warnings of :func:`arf_get_fmpz` apply.
+    matrix such that `U^T A U` is close to being LLL-reduced for the given
+    parameters `\delta,\eta`. If :func:`arb_mat_spd_get_fmpz_mat` succeeds at
+    the chosen precision, we call :func:`fmpz_lll`, and otherwise set *U* to
+    the identity matrix. The warnings of :func:`arf_get_fmpz` apply.
 
-.. function:: int arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec)
+.. function:: int arb_mat_spd_is_lll_reduced(const arb_mat_t A, double delta, double eta, slong prec)
 
     Given a symmetric positive definite matrix *A*, returns nonzero iff *A* is
-    certainly LLL-reduced with a tolerance of `\varepsilon = 2^{tol\_exp}`,
-    meaning that it satisfies the inequalities `|\mu_{j,k}|\leq \eta +
-    \varepsilon` and `(\delta - \varepsilon) \lVert b_{k-1}^*\rVert^2 \leq
-    \lVert b_k^*\rVert^2 + \mu_{k,k-1}^2 \lVert b_{k-1}^*\rVert^2` (with the
-    usual notation) for the default parameters `\eta = 0.51`, `\delta = 0.99`.
+    certainly LLL-reduced for the given LLL parameters `\eta, \delta`, meaning
+    that it satisfies the inequalities `|\mu_{j,k}|\leq \eta` and `\delta
+    \lVert b_{k-1}^*\rVert^2 \leq \lVert b_k^*\rVert^2 + \mu_{k,k-1}^2 \lVert
+    b_{k-1}^*\rVert^2` (with the usual notation).

src/acb_theta/siegel_is_reduced.c

@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "math.h"
 #include "fmpz_mat.h"
 #include "arb_mat.h"
 #include "acb.h"
@@ -24,6 +25,7 @@ acb_siegel_is_reduced(const acb_mat_t tau, slong tol_exp, slong prec)
     arb_mat_t im;
     acb_t det;
     arb_t abs, t, u;
+    double v;
     slong j, k;
     int res = 1;
 
@@ -55,8 +57,10 @@ acb_siegel_is_reduced(const acb_mat_t tau, slong tol_exp, slong prec)
 
     if (res)
     {
+        v = pow((double) 2, (double) tol_exp);
         acb_mat_get_imag(im, tau);
-        res = arb_mat_spd_is_lll_reduced(im, tol_exp, prec);
+        /* cf. choice of LLL parameters in siegel_reduce.c */
+        res = arb_mat_spd_is_lll_reduced(im, 0.99 - v, 0.51 + v, prec);
     }
 
     arb_one(t);

src/acb_theta/siegel_reduce.c

@@ -105,7 +105,7 @@ acb_siegel_reduce_imag(fmpz_mat_t mat, const acb_mat_t tau, slong prec)
     fmpz_mat_init(U, g, g);
 
     acb_mat_get_imag(im, tau);
-    arb_mat_spd_lll_reduce(U, im, prec);
+    arb_mat_spd_lll_reduce(U, im, 0.99, 0.51, prec);
     sp2gz_block_diag(mat, U);
 
     arb_mat_clear(im);
@@ -164,7 +164,7 @@ acb_siegel_reduce(fmpz_mat_t mat, const acb_mat_t tau, slong prec)
         lp = acb_siegel_reduce_real_lowprec(ntau, nmat, g, prec);
         acb_siegel_transform(w, m, w, lp);
         acb_mat_get_imag(im, w);
-        if (!arb_mat_spd_is_lll_reduced(im, -10, lp))
+        if (!arb_mat_spd_is_lll_reduced(im, 0.989, 0.511, lp))
         {
             stop = 1;
             break;

src/arb_mat.h

@@ -317,9 +317,9 @@ slong arb_mat_allocated_bytes(const arb_mat_t x);
 
 int arb_mat_spd_get_fmpz_mat(fmpz_mat_t B, const arb_mat_t A, slong prec);
 
-void arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, slong prec);
+void arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, double delta, double eta, slong prec);
 
-int arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec);
+int arb_mat_spd_is_lll_reduced(const arb_mat_t A, double delta, double eta, slong prec);
 
 #ifdef __cplusplus
 }

src/arb_mat/spd_is_lll_reduced.c

@@ -16,12 +16,12 @@
 
 /* Adapted from fmpz_mat_is_reduced_gram */
 int
-arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec)
+arb_mat_spd_is_lll_reduced(const arb_mat_t A, double delta, double eta, slong prec)
 {
     slong d = arb_mat_nrows(A);
     arb_mat_t r, mu;
     arb_ptr s;
-    arb_t delta, eta, t;
+    arb_t delta_arb, eta_arb, t;
     slong i, j, k;
     int res = 1;
 
@@ -33,18 +33,12 @@ arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec)
     arb_mat_init(r, d, d);
     arb_mat_init(mu, d, d);
     s = _arb_vec_init(d);
-    arb_init(delta);
-    arb_init(eta);
+    arb_init(delta_arb);
+    arb_init(eta_arb);
     arb_init(t);
 
-    arb_one(t);
-    arb_mul_2exp_si(t, t, tol_exp);
-    arb_set_si(delta, 99);
-    arb_div_si(delta, delta, 100, prec);
-    arb_sub(delta, delta, t, prec);
-    arb_set_si(eta, 51);
-    arb_div_si(eta, eta, 100, prec);
-    arb_add(eta, eta, t, prec);
+    arb_set_d(delta_arb, delta);
+    arb_set_d(eta_arb, eta);
 
     arb_set(arb_mat_entry(r, 0, 0), arb_mat_entry(A, 0, 0));
 
@@ -62,15 +56,15 @@ arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec)
             arb_div(arb_mat_entry(mu, i, j), arb_mat_entry(r, i, j),
                 arb_mat_entry(r, j, j), prec);
             arb_abs(t, arb_mat_entry(mu, i, j));
-            if (!arb_le(t, eta))
+            if (!arb_le(t, eta_arb))
             {
                 res = 0;
             }
             arb_set(&s[j + 1], &s[j]);
             arb_submul(&s[j + 1], arb_mat_entry(mu, i, j), arb_mat_entry(r, i, j), prec);
         }
         arb_set(arb_mat_entry(r, i, i), &s[i]);
-        arb_mul(t, delta, arb_mat_entry(r, i - 1, i - 1), prec);
+        arb_mul(t, delta_arb, arb_mat_entry(r, i - 1, i - 1), prec);
         if (!arb_le(t, &s[i - 1]))
         {
             res = 0;
@@ -80,8 +74,8 @@ arb_mat_spd_is_lll_reduced(const arb_mat_t A, slong tol_exp, slong prec)
     arb_mat_clear(r);
     arb_mat_clear(mu);
     _arb_vec_clear(s, d);
-    arb_clear(delta);
-    arb_clear(eta);
+    arb_clear(delta_arb);
+    arb_clear(eta_arb);
     arb_clear(t);
     return res;
 }

src/arb_mat/spd_lll_reduce.c

@@ -14,7 +14,7 @@
 #include "arb_mat.h"
 
 void
-arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, slong prec)
+arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, double delta, double eta, slong prec)
 {
     fmpz_lll_t fl;
     fmpz_mat_t N;
@@ -33,7 +33,7 @@ arb_mat_spd_lll_reduce(fmpz_mat_t U, const arb_mat_t A, slong prec)
     if (r)
     {
         /* Default Flint LLL values, except Gram */
-        fmpz_lll_context_init(fl, 0.99, 0.51, GRAM, EXACT);
+        fmpz_lll_context_init(fl, delta, eta, GRAM, EXACT);
         fmpz_lll(N, U, fl);
     }
 

src/arb_mat/test/t-spd_lll_reduce.c

@@ -24,7 +24,6 @@ TEST_FUNCTION_START(arb_mat_spd_lll_reduce, state)
         slong g = 1 + n_randint(state, 4);
         slong prec = 200 + n_randint(state, 500);
         slong mag_bits = 1 + n_randint(state, 5);
-        slong tol_exp = -10 - n_randint(state, 20);
         arb_mat_t M;
         arb_mat_t R;
         arb_mat_t T;
@@ -36,14 +35,14 @@ TEST_FUNCTION_START(arb_mat_spd_lll_reduce, state)
         fmpz_mat_init(U, g, g);
 
         arb_mat_randtest_spd(M, state, prec, mag_bits);
-        arb_mat_spd_lll_reduce(U, M, prec);
+        arb_mat_spd_lll_reduce(U, M, 0.99, 0.51, prec);
 
         arb_mat_set_fmpz_mat(T, U);
         arb_mat_mul(R, T, M, prec);
         arb_mat_transpose(T, T);
         arb_mat_mul(R, R, T, prec);
 
-        if (!arb_mat_spd_is_lll_reduced(R, tol_exp, prec))
+        if (!arb_mat_spd_is_lll_reduced(R, 0.989, 0.511, prec))
         {
             flint_printf("FAIL (reduction)\n");
             arb_mat_printd(M, 10);
@@ -53,13 +52,14 @@ TEST_FUNCTION_START(arb_mat_spd_lll_reduce, state)
             flint_abort();
         }
 
-        if (arb_mat_spd_is_lll_reduced(M, tol_exp, prec)
+        if (arb_mat_spd_is_lll_reduced(M, 1.01, 0.499, prec)
             && !fmpz_mat_is_one(U))
         {
             flint_printf("FAIL (identity)\n");
             arb_mat_printd(M, 10);
             fmpz_mat_print_pretty(U);
             flint_printf("\n");
+            flint_abort();
         }
 
         arb_mat_clear(M);