Packed basis matrices in LLL

doc/source/fmpz_lll.rst

@@ -3,6 +3,93 @@
 **fmpz_lll.h** -- LLL reduction
 ==================================================================================================
 
+Packed basis matrix
+--------------------------------------------------------------------------------
+
+The floating-point LLL variants operate internally on a packed copy of the
+basis matrix, stored as a :type:`fmpz_lll_packed_t`, in which every entry
+occupies a fixed number of limbs in two's complement and the entries of a
+row are contiguous. Row operations are performed slot-wise modulo
+a power of two, which is exact by virtue of an upper bound on the bit length
+of the entries of each row which is maintained by all operations; the number
+of limbs is grown or shrunk as needed. This representation is much cheaper
+to update than a matrix of ``fmpz`` entries for the small (one to a few limbs)
+entries typical of lattice reduction, and Gram matrix entries can be computed
+exactly with fixed-size dot product kernels.
+
+.. type:: fmpz_lll_packed_struct
+
+.. type:: fmpz_lll_packed_t
+
+.. function:: slong fmpz_lll_packed_limbs(const fmpz_mat_t B)
+
+    Returns the number of limbs per entry to use for packing ``B``,
+    including some headroom. This is determined by the largest entry of
+    ``B``; since all entries use the same number of limbs, packing is
+    only used by the LLL functions when this number is small
+    (``FMPZ_LLL_PACKED_MAX_LIMBS``, respectively
+    ``FMPZ_LLL_PACKED_MAX_LIMBS_MPF`` in the multiprecision LLL where
+    the Gram matrix entries are computed exactly) and the packed matrix
+    is not too large (``FMPZ_LLL_PACKED_MAX_SIZE`` limbs). Slots are
+    shrunk during the reduction once all entries have become smaller.
+
+.. function:: void fmpz_lll_packed_init(fmpz_lll_packed_t P, slong d, slong n, slong m)
+              void fmpz_lll_packed_clear(fmpz_lll_packed_t P)
+
+    Initialises resp. clears a packed matrix with `d` rows, `n` columns
+    and `m` limbs per entry.
+
+.. function:: void fmpz_lll_packed_set_fmpz_mat(fmpz_lll_packed_t P, const fmpz_mat_t B)
+              void fmpz_lll_packed_get_fmpz_mat(fmpz_mat_t B, const fmpz_lll_packed_t P)
+
+    Packs ``B`` (whose entries must fit) into ``P``, resp. unpacks ``P``
+    into ``B``.
+
+.. function:: void fmpz_lll_packed_tighten(fmpz_lll_packed_t P, slong i)
+              void fmpz_lll_packed_grow(fmpz_lll_packed_t P)
+              void fmpz_lll_packed_maybe_shrink(fmpz_lll_packed_t P)
+
+    Recomputes the exact bit bound of row `i`; increases the number of
+    limbs per entry by one; decreases the number of limbs per entry by one
+    if all rows comfortably fit.
+
+.. function:: void fmpz_lll_packed_move_row(fmpz_lll_packed_t P, slong i, slong j)
+
+    Moves row `i` to position `j`, shifting the rows in between.
+
+.. function:: void fmpz_lll_packed_row_sub(fmpz_lll_packed_t P, slong i, slong j)
+              void fmpz_lll_packed_row_add(fmpz_lll_packed_t P, slong i, slong j)
+              void fmpz_lll_packed_row_submul_si(fmpz_lll_packed_t P, slong i, slong j, slong x)
+              void fmpz_lll_packed_row_submul_fmpz(fmpz_lll_packed_t P, slong i, slong j, const fmpz_t x)
+              void fmpz_lll_packed_row_submul_si_2exp(fmpz_lll_packed_t P, slong i, slong j, slong x, ulong e)
+
+    Sets row `i` to row `i` minus (or plus) row `j`, respectively minus
+    `x` (times `2^e`) times row `j`. The results are exact.
+
+.. function:: slong fmpz_lll_packed_get_d_vec_2exp(double * appv, fmpz_lll_packed_t P, slong i)
+
+    Equivalent to :func:`_fmpz_vec_get_d_vec_2exp` applied to row `i`.
+    Also sets the bit bound of the row to its exact value.
+
+.. function:: void fmpz_lll_packed_dot(fmpz_t res, const fmpz_lll_packed_t P, slong i, slong j, slong len)
+
+    Sets ``res`` to the exact dot product of the first ``len`` entries of
+    rows `i` and `j`.
+
+.. function:: double fmpz_lll_packed_heuristic_dot(const double * vec1, const double * vec2, slong len2, const fmpz_lll_packed_t P, slong k, slong j, slong exp_adj)
+
+    Equivalent to :func:`fmpz_lll_heuristic_dot` for a packed matrix.
+
+.. function:: int fmpz_lll_check_babai_packed(int kappa, fmpz_lll_packed_t P, fmpz_mat_t U, d_mat_t mu, d_mat_t r, double *s, d_mat_t appB, int *expo, fmpz_gram_t A, int a, int zeros, int kappamax, int n, const fmpz_lll_t fl, int heuristic)
+              int fmpz_lll_advance_check_babai_packed(int cur_kappa, int kappa, fmpz_lll_packed_t P, fmpz_mat_t U, d_mat_t mu, d_mat_t r, double *s, d_mat_t appB, int *expo, fmpz_gram_t A, int a, int zeros, int kappamax, int n, const fmpz_lll_t fl, int heuristic)
+              int fmpz_lll_check_babai_heuristic_packed(int kappa, fmpz_lll_packed_t P, fmpz_mat_t U, gr_mat_t mu, gr_mat_t r, gr_ptr s, fmpz_gram_t A, int a, int zeros, int kappamax, int n, gr_ptr tmp, gr_ptr rtmp, gr_ctx_t ctx, const fmpz_lll_t fl)
+
+    Versions of :func:`fmpz_lll_check_babai`, :func:`fmpz_lll_advance_check_babai`
+    and :func:`fmpz_lll_check_babai_heuristic` operating on a packed basis
+    matrix (``heuristic`` selects the heuristic variant). In the
+    multiprecision version the Gram matrix entries are computed exactly.
+
+
 Parameter manipulation
 --------------------------------------------------------------------------------
 

src/fmpz_lll.h

@@ -69,6 +69,61 @@ typedef union
 
 typedef fmpz_gram_union fmpz_gram_t[1];
 
+/* Packed basis matrix  ******************************************************/
+
+typedef struct
+{
+    nn_ptr entries;     /* d * n * m limbs */
+    nn_ptr * rows;      /* row pointers, permuted by row moves */
+    slong * bits;       /* upper bound on the bit length of the entries of each row */
+    slong d;            /* rows */
+    slong n;            /* columns */
+    slong m;            /* limbs per entry (two's complement) */
+}
+fmpz_lll_packed_struct;
+
+typedef fmpz_lll_packed_struct fmpz_lll_packed_t[1];
+
+/*
+    Only pack if the entries fit in this many limbs. In the multiprecision
+    LLL the Gram matrix is computed exactly from the packed rows, which is
+    only competitive with approximate dot products for small entries.
+*/
+#define FMPZ_LLL_PACKED_MAX_LIMBS 32
+#define FMPZ_LLL_PACKED_MAX_LIMBS_MPF FLINT_MPN_DOT_TAB_N
+
+/* do not pack if the packed matrix would need more limbs than this */
+#define FMPZ_LLL_PACKED_MAX_SIZE (WORD(1) << 26)
+
+slong fmpz_lll_packed_limbs(const fmpz_mat_t B);
+void fmpz_lll_packed_init(fmpz_lll_packed_t P, slong d, slong n, slong m);
+void fmpz_lll_packed_clear(fmpz_lll_packed_t P);
+void fmpz_lll_packed_set_fmpz_mat(fmpz_lll_packed_t P, const fmpz_mat_t B);
+void fmpz_lll_packed_get_fmpz_mat(fmpz_mat_t B, const fmpz_lll_packed_t P);
+void fmpz_lll_packed_tighten(fmpz_lll_packed_t P, slong i);
+void fmpz_lll_packed_grow(fmpz_lll_packed_t P);
+void fmpz_lll_packed_maybe_shrink(fmpz_lll_packed_t P);
+void fmpz_lll_packed_move_row(fmpz_lll_packed_t P, slong i, slong j);
+void fmpz_lll_packed_row_sub(fmpz_lll_packed_t P, slong i, slong j);
+void fmpz_lll_packed_row_add(fmpz_lll_packed_t P, slong i, slong j);
+void fmpz_lll_packed_row_submul_si(fmpz_lll_packed_t P, slong i, slong j, slong x);
+void fmpz_lll_packed_row_submul_fmpz(fmpz_lll_packed_t P, slong i, slong j, const fmpz_t x);
+void fmpz_lll_packed_row_submul_si_2exp(fmpz_lll_packed_t P, slong i, slong j, slong x, ulong e);
+slong fmpz_lll_packed_get_d_vec_2exp(double * appv, fmpz_lll_packed_t P, slong i);
+void fmpz_lll_packed_dot(fmpz_t res, const fmpz_lll_packed_t P, slong i, slong j, slong len);
+double fmpz_lll_packed_heuristic_dot(const double * vec1, const double * vec2, slong len2,
+       const fmpz_lll_packed_t P, slong k, slong j, slong exp_adj);
+
+int fmpz_lll_check_babai_packed(int kappa, fmpz_lll_packed_t P, fmpz_mat_t U, d_mat_t mu, d_mat_t r, double *s,
+       d_mat_t appB, int *expo, fmpz_gram_t A,
+       int a, int zeros, int kappamax, int n, const fmpz_lll_t fl, int heuristic);
+int fmpz_lll_advance_check_babai_packed(int cur_kappa, int kappa, fmpz_lll_packed_t P, fmpz_mat_t U, d_mat_t mu, d_mat_t r, double *s,
+       d_mat_t appB, int *expo, fmpz_gram_t A,
+       int a, int zeros, int kappamax, int n, const fmpz_lll_t fl, int heuristic);
+int fmpz_lll_check_babai_heuristic_packed(int kappa, fmpz_lll_packed_t P, fmpz_mat_t U,
+       gr_mat_t mu, gr_mat_t r, gr_ptr s, fmpz_gram_t A, int a, int zeros,
+       int kappamax, int n, gr_ptr tmp, gr_ptr rtmp, gr_ctx_t ctx, const fmpz_lll_t fl);
+
 /* Parameter manipulation  ***************************************************/
 
 void fmpz_lll_context_init_default(fmpz_lll_t fl);

src/fmpz_lll/check_babai_heuristic.c

@@ -70,8 +70,8 @@ static int _gr_vec_norm2(gr_ptr res, gr_srcptr vec, slong len, gr_ctx_t ctx)
 
 /* XXX: dubious use of DBL_MIN */
 
-int
-fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
+static int
+_fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_lll_packed_struct * P, fmpz_mat_t U,
                                gr_mat_t mu, gr_mat_t r, gr_ptr s,
                                gr_mat_t appB, fmpz_gram_t A, int a, int zeros,
                                int kappamax, int n, gr_ptr tmp, gr_ptr rtmp,
@@ -112,7 +112,13 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
             {
                 if (_gr_cmp_d(ENTRY(A->appSP2, kappa, j), DBL_MIN, ctx) == 0)
                 {
-                    status |= _gr_vec_dot(ENTRY(A->appSP2, kappa, j), NULL, 0, ROW(appB, kappa), ROW(appB, j), n, ctx);
+                    if (P != NULL)
+                    {
+                        fmpz_lll_packed_dot(ztmp, P, kappa, j, n);
+                        status |= gr_set_fmpz(ENTRY(A->appSP2, kappa, j), ztmp, ctx);
+                    }
+                    else
+                        status |= _gr_vec_dot(ENTRY(A->appSP2, kappa, j), NULL, 0, ROW(appB, kappa), ROW(appB, j), n, ctx);
 
 #if 0
                     /* If a heuristic told us that some cancellation probably happened,
@@ -167,14 +173,20 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
                         if (sgn >= 0)   /* in this case, X is 1 */
                         {
                             status |= _gr_vec_sub(ENTRY(mu, kappa, zeros + 1), ENTRY(mu, kappa, zeros + 1), ENTRY(mu, j, zeros + 1), j - (zeros + 1), ctx);
-                            _fmpz_vec_sub(fmpz_mat_row(B, kappa), fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n);
+                            if (P != NULL)
+                                fmpz_lll_packed_row_sub(P, kappa, j);
+                            else
+                                _fmpz_vec_sub(fmpz_mat_row(B, kappa), fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n);
                             if (U != NULL)
                                 _fmpz_vec_sub(fmpz_mat_row(U, kappa), fmpz_mat_row(U, kappa), fmpz_mat_row(U, j), U->c);
                         }
                         else    /* otherwise X is -1 */
                         {
                             status |= _gr_vec_add(ENTRY(mu, kappa, zeros + 1), ENTRY(mu, kappa, zeros + 1), ENTRY(mu, j, zeros + 1), j - (zeros + 1), ctx);
-                            _fmpz_vec_add(fmpz_mat_row(B, kappa), fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n);
+                            if (P != NULL)
+                                fmpz_lll_packed_row_add(P, kappa, j);
+                            else
+                                _fmpz_vec_add(fmpz_mat_row(B, kappa), fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n);
                             if (U != NULL)
                                 _fmpz_vec_add(fmpz_mat_row(U, kappa), fmpz_mat_row(U, kappa), fmpz_mat_row(U, j), U->c);
                         }
@@ -200,7 +212,10 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
                         status |= _gr_vec_submul_scalar(ENTRY(mu, kappa, zeros + 1), ENTRY(mu, j, zeros + 1), j - (zeros + 1), tmp, ctx);
                         status |= gr_get_fmpz(ztmp, tmp, ctx);
 
-                        _fmpz_vec_scalar_submul_fmpz(fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n, ztmp);
+                        if (P != NULL)
+                            fmpz_lll_packed_row_submul_fmpz(P, kappa, j, ztmp);
+                        else
+                            _fmpz_vec_scalar_submul_fmpz(fmpz_mat_row(B, kappa), fmpz_mat_row(B, j), n, ztmp);
                         if (U != NULL)
                             _fmpz_vec_scalar_submul_fmpz(fmpz_mat_row(U, kappa), fmpz_mat_row(U, j), U->c, ztmp);
                     }
@@ -209,7 +224,10 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
 
             if (test)           /* Anything happened? */
             {
-                status |= _gr_vec_set_fmpz_vec(ROW(appB, kappa), fmpz_mat_row(B, kappa), n, ctx);
+                if (P != NULL)
+                    fmpz_lll_packed_tighten(P, kappa);
+                else
+                    status |= _gr_vec_set_fmpz_vec(ROW(appB, kappa), fmpz_mat_row(B, kappa), n, ctx);
                 aa = zeros + 1;
 
                 for (i = zeros + 1; i <= kappa; i++)
@@ -222,7 +240,15 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
         } while (test);
 
         if (_gr_cmp_d(ENTRY(A->appSP2, kappa, kappa), DBL_MIN, ctx) == 0)
-            status |= _gr_vec_norm2(ENTRY(A->appSP2, kappa, kappa), ROW(appB, kappa), n, ctx);
+        {
+            if (P != NULL)
+            {
+                fmpz_lll_packed_dot(ztmp, P, kappa, kappa, n);
+                status |= gr_set_fmpz(ENTRY(A->appSP2, kappa, kappa), ztmp, ctx);
+            }
+            else
+                status |= _gr_vec_norm2(ENTRY(A->appSP2, kappa, kappa), ROW(appB, kappa), n, ctx);
+        }
 
         status |= gr_set(GR_ENTRY(s, zeros + 1, sz), ENTRY(A->appSP2, kappa, kappa), ctx);
 
@@ -462,3 +488,22 @@ fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
 }
 
 #undef GM
+
+int
+fmpz_lll_check_babai_heuristic(int kappa, fmpz_mat_t B, fmpz_mat_t U,
+                               gr_mat_t mu, gr_mat_t r, gr_ptr s,
+                               gr_mat_t appB, fmpz_gram_t A, int a, int zeros,
+                               int kappamax, int n, gr_ptr tmp, gr_ptr rtmp,
+                               gr_ctx_t ctx, const fmpz_lll_t fl)
+{
+    return _fmpz_lll_check_babai_heuristic(kappa, B, NULL, U, mu, r, s, appB, A, a, zeros, kappamax, n, tmp, rtmp, ctx, fl);
+}
+
+int
+fmpz_lll_check_babai_heuristic_packed(int kappa, fmpz_lll_packed_t P, fmpz_mat_t U,
+                               gr_mat_t mu, gr_mat_t r, gr_ptr s, fmpz_gram_t A, int a, int zeros,
+                               int kappamax, int n, gr_ptr tmp, gr_ptr rtmp,
+                               gr_ctx_t ctx, const fmpz_lll_t fl)
+{
+    return _fmpz_lll_check_babai_heuristic(kappa, NULL, P, U, mu, r, s, NULL, A, a, zeros, kappamax, n, tmp, rtmp, ctx, fl);
+}

src/fmpz_lll/heuristic_dot.c

@@ -36,7 +36,7 @@ fmpz_lll_heuristic_dot(const double *vec1, const double *vec2, slong len2,
         fmpz_init(sp);
         _fmpz_vec_dot(sp, fmpz_mat_row(B, k), fmpz_mat_row(B, j), len2);
         sum = fmpz_get_d_2exp(&exp, sp);
-        sum = ldexp(sum, sum - exp_adj);
+        sum = ldexp(sum, exp - exp_adj);
         fmpz_clear(sp);
     }