Ticket #723: ntls-fp-lll.patch

sage/libs/ntl/ntl_mat_ZZ.pyx

@@ -17 +17 @@
 include "../../ext/stdsage.pxi"
 include 'misc.pxi'
 include 'decl.pxi'
+
+cdef extern from "NTL/LLL.h":
+    cdef long mat_ZZ_LLL_FP   "LLL_FP"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_FP_U "LLL_FP"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_QP   "LLL_QP"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_QP_U "LLL_QP"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_XD   "LLL_XD"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_XD_U "LLL_XD"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_RR   "LLL_RR"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_LLL_RR_U "LLL_RR"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+
+    cdef long mat_ZZ_G_LLL_FP   "G_LLL_FP"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_FP_U "G_LLL_FP"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_QP   "G_LLL_QP"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_QP_U "G_LLL_QP"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_XD   "G_LLL_XD"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_XD_U "G_LLL_XD"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_RR   "G_LLL_RR"(mat_ZZ_c B, double delta, int deep, int check , int verbose)
+    cdef long mat_ZZ_G_LLL_RR_U "G_LLL_RR"(mat_ZZ_c B, mat_ZZ_c U, double delta, int deep, int check , int verbose)
+
 
 from sage.libs.ntl.ntl_ZZ cimport ntl_ZZ
 from sage.libs.ntl.ntl_ZZX cimport ntl_ZZX
@@ -352 +372 @@
             rank = int(mat_ZZ_LLL(&det2,&self.x,int(a),int(b),int(verbose)))
             _sig_off
             return rank,make_ZZ(det2)
+
+    def LLL_FP(self, delta=0.75 , return_U=False, verbose=False):
+        r"""
+        Performs approximate LLL reduction of self (puts self in an
+        LLL form) subject to the following conditions:
+
+        The precision is double.
+        
+        The return value is the rank of B.
+
+        Classical Gramm-Schmidt Orthogonalization is used:
+
+        This choice uses classical methods for computing the
+        Gramm-Schmidt othogonalization.  It is fast but prone to
+        stability problems.  This strategy was first proposed by
+        Schnorr and Euchner [C. P. Schnorr and M. Euchner,
+        Proc. Fundamentals of Computation Theory, LNCS 529, pp. 68-85,
+        1991].  The version implemented here is substantially
+        different, improving both stability and performance.
+        
+        If return_U is True, then also U is returned which is
+        the transition matrix: $U * self_{old} = self_{new}$
+
+        The optional argument 'delta' is the reduction parameter, and
+        may be set so that 0.50 <= delta < 1.  Setting it close to 1
+        yields shorter vectors, and also improves the stability, but
+        increases the running time.  Recommended value: delta =
+        0.99.
+
+        The optional parameter 'verbose' can be set to see all kinds
+        of fun things printed while the routine is executing.  A
+        status report is also printed every once in a while.
+
+        INPUT:
+           delta    -- as described above (0.5 <= delta < 1.0) (default: 0.75)
+           return_U -- return U as described above
+           verbose  -- if True NTL will produce some verbatim messages on
+                       what's going on internally (default: False)
+
+        OUTPUT:
+            (rank,[U]) where rank and U are as described above and U
+            is an optional return value if return_U is True.
+
+        EXAMPLE:
+            sage: M=ntl.mat_ZZ(3,3,[1,2,3,4,5,6,7,8,9])
+            sage: M.LLL_FP()
+            2
+            sage: M
+            [[0 0 0]
+            [2 1 0]
+            [-1 1 3]
+            ]
+            sage: M=ntl.mat_ZZ(4,4,[-6,9,-15,-18,4,-6,10,12,10,-16,18,35,-24,36,-46,-82]); M
+            [[-6 9 -15 -18]
+            [4 -6 10 12]
+            [10 -16 18 35]
+            [-24 36 -46 -82]
+            ]
+            sage: M.LLL_FP()
+            3
+            sage: M
+            [[0 0 0 0]
+            [0 -2 0 0]
+            [-2 1 -5 -6]
+            [0 -1 -7 5]
+            ]
+
+        WARNING: This method modifies self. So after applying this method your matrix
+        will be a vector of vectors.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_LLL_FP_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_LLL_FP(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def LLL_QP(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.LLL_FP using the same
+        calling conventions but with quad float precision.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_LLL_QP_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_LLL_QP(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def LLL_XD(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.LLL_FP using the same
+        calling conventions but with extended exponent double
+        precision.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_LLL_XD_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_LLL_XD(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def LLL_RR(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.LLL_FP using the same
+        calling conventions but with arbitrary precision floating
+        point numbers.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_LLL_RR_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_LLL_RR(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    # Givens Orthogonalization.  This is a bit slower, but generally
+    # much more stable, and is really the preferred orthogonalization
+    # strategy.  For a nice description of this, see Chapter 5 of
+    # [G. Golub and C. van Loan, Matrix Computations, 3rd edition,
+    # Johns Hopkins Univ. Press, 1996].
+
+    def G_LLL_FP(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.LLL_FP using the same
+        calling conventions but uses the Givens Orthogonalization.
+
+        Givens Orthogonalization.  This is a bit slower, but generally
+        much more stable, and is really the preferred
+        orthogonalization strategy.  For a nice description of this,
+        see Chapter 5 of [G. Golub and C. van Loan, Matrix
+        Computations, 3rd edition, Johns Hopkins Univ. Press, 1996].
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_FP_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_FP(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def G_LLL_QP(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.G_LLL_FP using the same
+        calling conventions but with quad float precision.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_QP_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_QP(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def G_LLL_XD(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.G_LLL_FP using the same
+        calling conventions but with extended exponent double
+        precision.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_XD_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_XD(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank
+
+    def G_LLL_RR(self, delta, return_U=False, verbose=False):
+        r"""
+        Peforms the same reduction as self.G_LLL_FP using the same
+        calling conventions but with aribitrary precision floating
+        point numbers.
+        """
+        cdef ntl_mat_ZZ U
+        if return_U:
+            U = PY_NEW(ntl_mat_ZZ)
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_RR_U(self.x, U.x, float(delta), 0, 0, int(verbose)))
+            _sig_off
+            return rank, U
+        else:
+            _sig_on
+            rank = int(mat_ZZ_G_LLL_RR(self.x,float(delta),0,0,int(verbose)))
+            _sig_off
+            return rank

sage/matrix/matrix_integer_dense.pyx

@@ -1404 +1404 @@
     ####################################################################################
     # LLL
     ####################################################################################    
-    def lllgram(self):
+    def LLL_gram(self):
         """
         LLL reduction of the lattice whose gram matrix is self.
 
@@ -1425 +1425 @@
             sage: M = Matrix(ZZ, 2, 2, [5,3,3,2]) ; M
             [5 3]
             [3 2]
-            sage: U = M.lllgram(); U
+            sage: U = M.LLL_gram(); U
             [-1  1]
             [ 1 -2]
             sage: U.transpose() * M * U
@@ -1434 +1434 @@
             
         Semidefinite and indefinite forms raise a ValueError:
 
-            sage: Matrix(ZZ,2,2,[2,6,6,3]).lllgram()
+            sage: Matrix(ZZ,2,2,[2,6,6,3]).LLL_gram()
             Traceback (most recent call last):
             ...
             ValueError: not a definite matrix
-            sage: Matrix(ZZ,2,2,[1,0,0,-1]).lllgram()
+            sage: Matrix(ZZ,2,2,[1,0,0,-1]).LLL_gram()
             Traceback (most recent call last):
             ...
             ValueError: not a definite matrix
@@ -1464 +1464 @@
                 U[i,n-1] = - U[i,n-1]
         return U
 
-    def lll(self, delta=None):
+    def LLL(self, delta=None, algorithm=None, **kwargs):
         r"""
-        Returns LLL reduced lattice R for self.
+        Returns LLL reduced or approximated LLL reduced lattice R for
+        self.
 
         The lattice is returned as a matrix. Also the rank (and the
         determinant) of self are cached.
 
         More specifically, elementary row transformations are
-        performed on a copy of self so that the non-zero rows of
-        R form an LLL-reduced basis for the lattice spanned by
-        the rows of self. The default reduction parameter is
-        $\delta=3/4$, which means that the squared length of the first
-        non-zero basis vector is no more than $2^{r-1}$ times that of
-        the shortest vector in the lattice.
+        performed on a copy of self so that the non-zero rows of R
+        form an LLL-reduced basis for the lattice spanned by the rows
+        of self. The default reduction parameter is $\delta=3/4$,
+        which means that the squared length of the first non-zero
+        basis vector is no more than $2^{r-1}$ times that of the
+        shortest vector in the lattice.
 
         For a basis reduced with parameter $\delta$, the squared
         length of the first non-zero basis vector is no more than
@@ -1486 +1487 @@
 
         If we can compute the determinant of self using this method,
         we also cache it. Note that in general this only happens when
-        self.rank() == self.ncols().
+        self.rank() == self.ncols() and the exact algorithm is used.
 
         INPUT:
-           delta -- arameter a as described above (default: 3/4)
+            delta -- arameter a as described above (default: 3/4)
+            algorithm -- string, one of the algorithms mentioned above
+                        or None (default: None)
+            use_givens -- use Givens orthogonalization (default: False)
+                          only applicable to approximate reductions and NTL.
+                          This is more stable but slower.
+
+        Also, if the verbose level is >= 2, some more verbose output
+        is printed during the calculation if NTL is used.
+
+        AVAILABLE ALGORITHMS:
+            NTL:LLL -- default, exact reduction
+            NTL:LLL_FP -- approximate reduction over double precision
+                          floating point numbers.
+            NTL:LLL_QP -- approximate reduction over quad precision
+                          floating point numbers.
+            NTL:LLL_XD -- approximate reduction over extended exponent
+                          double precision floating point numbers.
+            NTL:LLL_RR -- approximate reduction over arbitrary precision
+                          floating point numbers.
 
         OUTPUT:
             a matrix over the integers
 
         EXAMPLE:
             sage: A = Matrix(ZZ,3,3,range(1,10))
-            sage: A.lll()
+            sage: A.LLL()
             [ 0  0  0]
             [ 2  1  0]
             [-1  1  3]
             
         ALGORITHM: Uses NTL.
+
+        REFERENCES:
+            ntl.mat_ZZ for details on the used algorithms.
         """
 
         import sage.libs.ntl.all
         ntl_ZZ = sage.libs.ntl.all.ZZ
 
-        if delta is None:
-            delta = ZZ(3)/ZZ(4)
-        elif delta <= ZZ(1)/ZZ(4):
-            raise TypeError, "delta must be > 1/4"
-        elif delta > 1:
-            raise TypeError, "delta must be <= 1/4"
-
-        delta = delta/ZZ(1) # QQ(delta)
-
-        a = delta.numer()
-        b = delta.denom()
+        if get_verbose() >= 2: verb = True
+        else: verb = False
+
+        use_givens = kwargs.get("use_givens",False)
+
+        if algorithm is None:
+            algorithm = "NTL:LLL"
+
+        if algorithm == "NTL:LLL":
+            if delta is None:
+                delta = ZZ(3)/ZZ(4)
+            elif delta <= ZZ(1)/ZZ(4):
+                raise TypeError, "delta must be > 1/4"
+            elif delta > 1:
+                raise TypeError, "delta must be <= 1"
+
+            delta = QQ(delta) # QQ(delta)
+            a = delta.numer()
+            b = delta.denom()
+
+        elif algorithm in ("NTL:LLL_FP","NTL:LLL_QP","NTL:LLL_XD","NTL:LLL_RR"):
+            if delta is None:
+                delta = 0.99
+            elif delta < 0.5:
+                raise TypeError, "delta must be >= 0.5"
+            elif delta > 1:
+                raise TypeError, "delta must be <= 1"
+            delta = float(delta)
 
         A = sage.libs.ntl.all.mat_ZZ(self.nrows(),self.ncols(),map(ntl_ZZ,self.list()))
-        r, det2 = A.LLL(a,b)
-        r,det2 = ZZ(r), ZZ(det2)
+
+        if algorithm == "NTL:LLL":
+            r, det2 = A.LLL(a,b, verbose=verb)
+            det2 = ZZ(det2)
+            try:
+                det = ZZ(det2.sqrt_approx())
+                self.cache("det", det)
+            except TypeError:
+                pass
+        elif algorithm == "NTL:LLL_FP":
+            if use_givens:
+                r = A.G_LLL_FP(delta, verbose=verb)
+            else:
+                r = A.LLL_FP(delta, verbose=verb)
+        elif algorithm == "NTL:LLL_QP":
+            if use_givens:
+                r = A.G_LLL_QP(delta, verbose=verb)
+            else:
+                r = A.LLL_QP(delta, verbose=verb)
+        elif algorithm == "NTL:LLL_XD":
+            if use_givens:
+                r = A.G_LLL_XD(delta, verbose=verb)
+            else:
+                r = A.LLL_XD(delta, verbose=verb)
+        elif algorithm == "NTL:LLL_RR":
+            if use_givens:
+                r = A.G_LLL_RR(delta, verbose=verb)
+            else:
+                r = A.LLL_XD(delta, verbose=verb)
+        else:
+            raise TypeError, "algorithm %s not supported"%algorithm
+
+        r = ZZ(r)
 
         cdef Matrix_integer_dense R = <Matrix_integer_dense>self.new_matrix(entries=map(ZZ,A.list()))
         self.cache("rank",r)
-        try:
-            det = ZZ(det2.sqrt_approx())
-            self.cache("det", det)
-        except TypeError:
-            pass
+
         return R
 
     def prod_of_row_sums(self, cols):