gh-38259: Fixed and improvements in `is_LLL_reduced` and `approximate_closest_vector`
    
This fixes #38115. Along with the goal stated there I have fixed some
stuff I wasn't happy with from the first iteration of
`approximate_closest_vector`. These changes could be breaking but seeing
as there has been no release since its introduction this shouldn't be a
problem.

Changes:
 - Added the `algorithm` parameter to `is_LLL_reduced` along with the
new algorithm `fpLLL` which uses fpylll's `LLL.is_reduced` method, **and
made it the default**, as I see no advantages with the old
implementation.
 - Added the `algorithm` parameter to `approximate_closest_vector` along
with the new algorithms `embedding` and `rounding_off`, **and made
`embedding` the default.**
 - ~Changed the alias of `approximate_closest_vector` from `babai` to
`cvp` to be more accurate but still short.~

I'm getting weird errors when building documentation locally but I don't
think it's related to my changes, so I haven't been able to check if the
doc build looks alright yet. (and it won't build on my own repo b/c my
GitHub username has uppercase letters in it :)

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.
    
URL: #38259
Reported by: TheBlupper
Reviewer(s): Giacomo Pope, Matthias Köppe, TheBlupper

Author: Release Manager

src/sage/matrix/matrix_integer_dense.pyx

@@ -3297,7 +3297,7 @@ cdef class Matrix_integer_dense(Matrix_dense):
         else:
             return R
 
-    def is_LLL_reduced(self, delta=None, eta=None):
+    def is_LLL_reduced(self, delta=None, eta=None, algorithm='fpLLL'):
         r"""
         Return ``True`` if this lattice is `(\delta, \eta)`-LLL reduced.
         See ``self.LLL`` for a definition of LLL reduction.
@@ -3308,6 +3308,8 @@ cdef class Matrix_integer_dense(Matrix_dense):
 
         - ``eta`` -- (default: `0.501`) parameter `\eta` as described above
 
+        - ``algorithm`` -- either ``'fpLLL'`` (default) or ``'sage'``
+
         EXAMPLES::
 
             sage: A = random_matrix(ZZ, 10, 10)
@@ -3316,6 +3318,15 @@ cdef class Matrix_integer_dense(Matrix_dense):
             False
             sage: L.is_LLL_reduced()
             True
+
+        The ``'sage'`` algorithm currently does not work for matrices with
+        linearly dependent rows::
+
+            sage: A = matrix(ZZ, [[1, 2, 3], [2, 4, 6]])
+            sage: A.is_LLL_reduced(algorithm='sage')
+            Traceback (most recent call last):
+            ...
+            ValueError: linearly dependent input for module version of Gram-Schmidt
         """
         if eta is None:
             eta = 0.501
@@ -3330,20 +3341,27 @@ cdef class Matrix_integer_dense(Matrix_dense):
         if eta < 0.5:
             raise TypeError("eta must be >= 0.5")
 
-        # this is pretty slow
-        import sage.modules.misc
-        G, mu = sage.modules.misc.gram_schmidt(self.rows())
-        #For any $i>j$, we have $|mu_{i, j}| <= \eta$
-        for e in mu.list():
-            if e.abs() > eta:
-                return False
-
-        #For any $i<d$, we have $\delta |b_i^*|^2 <= |b_{i+1}^* + mu_{i+1, i} b_i^* |^2$
-        norms = [G[i].norm()**2 for i in range(len(G))]
-        for i in range(1,self.nrows()):
-            if norms[i] < (delta - mu[i,i-1]**2) * norms[i-1]:
-                return False
-        return True
+        if algorithm == 'fpLLL':
+            from fpylll import LLL, IntegerMatrix
+            A = IntegerMatrix.from_matrix(self)
+            return LLL.is_reduced(A, delta=delta, eta=eta)
+        elif algorithm == 'sage':
+            # This is pretty slow
+            import sage.modules.misc
+            G, mu = sage.modules.misc.gram_schmidt(self.rows())
+            # For any $i>j$, we have $|mu_{i, j}| <= \eta$
+            for e in mu.list():
+                if e.abs() > eta:
+                    return False
+
+            # For any $i<d$, we have $\delta |b_i^*|^2 <= |b_{i+1}^* + mu_{i+1, i} b_i^* |^2$
+            norms = [G[i].norm()**2 for i in range(len(G))]
+            for i in range(1,self.nrows()):
+                if norms[i] < (delta - mu[i,i-1]**2) * norms[i-1]:
+                    return False
+            return True
+        else:
+            raise ValueError("algorithm must be one of 'fpLLL' or 'sage'")
 
     def prod_of_row_sums(self, cols):
         """

src/sage/modules/free_module_integer.py

@@ -32,6 +32,7 @@
 ##############################################################################
 
 from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
 from sage.matrix.constructor import matrix
 from sage.misc.cachefunc import cached_method
 from sage.modules.free_module import FreeModule_submodule_with_basis_pid, FreeModule_ambient_pid
@@ -784,16 +785,11 @@ def CVPP_2V(t, V, voronoi_cell):
             i -= 1
         return t - t_new
 
-    def approximate_closest_vector(self, t, delta=None, *args, **kwargs):
+    def approximate_closest_vector(self, t, delta=None, algorithm='embedding', *args, **kwargs):
         r"""
-        Compute a vector `w` such that
-
-        .. MATH::
-
-            |t-w|<(\frac{1}{\delta-\frac{1}{4}})^{d/2}|t-u|
-
-        where `u` is the closest lattice point to `t`, `\delta` is the LLL
-        reduction parameter, and `d` is the dimension of the lattice.
+        Compute a vector `w` in this lattice which is close to the target vector `t`.
+        The ratio `\frac{|t-w|}{|t-u|}`, where `u` is the closest lattice vector to `t`,
+        is exponential in the dimension of the lattice.
 
         This will check whether the basis is already `\delta`-LLL-reduced
         and otherwise it will run LLL to make sure that it is. For more
@@ -805,6 +801,18 @@ def approximate_closest_vector(self, t, delta=None, *args, **kwargs):
 
         - ``delta`` -- (default: ``0.99``) the LLL reduction parameter
 
+        - ``algorithm`` -- string (default: 'embedding'):
+
+          - ``'embedding'`` -- embeds the lattice in a d+1 dimensional space
+            and seeks short vectors using LLL. This calls LLL twice but is
+            usually still quick.
+
+          - ``'nearest_plane'`` -- uses the "NEAREST PLANE" algorithm from [Bab86]_
+
+          - ``'rounding_off'`` -- uses the "ROUNDING OFF" algorithm from [Bab86]_.
+            This yields slightly worse results than the other algorithms but is
+            at least faster than ``'nearest_plane'``.
+
         - ``*args`` -- passed through to :meth:`LLL`
 
         - ``**kwds`` -- passed through to :meth:`LLL`
@@ -825,30 +833,62 @@ def approximate_closest_vector(self, t, delta=None, *args, **kwargs):
             ....:                     [0, 0, 101, 0], [-28, 39, 45, 1]], lll_reduce=False)
             sage: t = vector([1337]*4)
             sage: L.approximate_closest_vector(t, delta=0.26)
-            (1348, 1340, 1383, 1337)
+            (1331, 1324, 1349, 1334)
             sage: L.approximate_closest_vector(t, delta=0.99)
             (1326, 1349, 1339, 1345)
             sage: L.closest_vector(t)
             (1326, 1349, 1339, 1345)
 
-        ALGORITHM:
-
-        Uses the algorithm from [Bab86]_.
+            sage: # Checking that the other algorithms work
+            sage: L.approximate_closest_vector(t, algorithm='nearest_plane')
+            (1326, 1349, 1339, 1345)
+            sage: L.approximate_closest_vector(t, algorithm='rounding_off')
+            (1331, 1324, 1349, 1334)
         """
         if delta is None:
             delta = ZZ(99)/ZZ(100)
 
-        # bound checks on delta are performed in is_LLL_reduced
+        # Bound checks on delta are performed in is_LLL_reduced
         if not self._reduced_basis.is_LLL_reduced(delta=delta):
             self.LLL(*args, delta=delta, **kwargs)
 
         B = self._reduced_basis
-        G = B.gram_schmidt()[0]
         t = vector(t)
 
-        b = t
-        for i in reversed(range(G.nrows())):
-            b -= B[i] * ((b * G[i]) / (G[i] * G[i])).round("even")
-        return (t - b).change_ring(ZZ)
+        if algorithm == 'embedding':
+            L = matrix(QQ, B.nrows()+1, B.ncols()+1)
+            L.set_block(0, 0, B)
+            L.set_block(B.nrows(), 0, matrix(t))
+            weight = (B[-1]*B[-1]).isqrt()+1  # Norm of the largest vector
+            L[-1, -1] = weight
+
+            # The vector should be the last row but we iterate just in case
+            for v in reversed(L.LLL(delta=delta, *args, **kwargs).rows()):
+                if abs(v[-1]) == weight:
+                    return t - v[:-1]*v[-1].sign()
+            raise ValueError('No suitable vector found in basis.'
+                             'This is a bug, please report it.')
+
+        elif algorithm == 'nearest_plane':
+            G = B.gram_schmidt()[0]
+
+            b = t
+            for i in reversed(range(G.nrows())):
+                b -= B[i] * ((b * G[i]) / (G[i] * G[i])).round("even")
+            return (t - b).change_ring(ZZ)
+
+        elif algorithm == 'rounding_off':
+            # t = x*B might not have a solution over QQ so we instead solve
+            # the system x*B*B^T = t*B^T which will be the "closest" solution
+            # if it does not exist, same effect as using the psuedo-inverse
+            sol = (B*B.T).solve_left(t*B.T)
+            return vector(ZZ, [QQ(x).round('even') for x in sol])*B
 
-    babai = approximate_closest_vector
+        else:
+            raise ValueError("algorithm must be one of 'embedding', 'nearest_plane' or 'rounding_off'")
+
+    def babai(self, *args, **kwargs):
+        """
+        Alias for :meth:`approximate_closest_vector`.
+        """
+        return self.approximate_closest_vector(*args, **kwargs)