# HG changeset patch # User Martin Albrecht # Date 1190417056 -3600 # Node ID 2bb5a85fe302ae44f6147b176ef3f45e987d2fbd # Parent 5f0e731bb66c2c4ad46d97bdfab24a55c689d43e NTL's floating point LLL exposed diff -r 5f0e731bb66c -r 2bb5a85fe302 sage/libs/ntl/ntl_mat_ZZ.pyx --- a/sage/libs/ntl/ntl_mat_ZZ.pyx Fri Sep 21 10:29:03 2007 +0100 +++ b/sage/libs/ntl/ntl_mat_ZZ.pyx Sat Sep 22 00:24:16 2007 +0100 @@ -17,6 +17,26 @@ include "../../ext/stdsage.pxi" 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,3 +372,227 @@ cdef class ntl_mat_ZZ: 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 diff -r 5f0e731bb66c -r 2bb5a85fe302 sage/matrix/matrix_integer_dense.pyx --- a/sage/matrix/matrix_integer_dense.pyx Fri Sep 21 10:29:03 2007 +0100 +++ b/sage/matrix/matrix_integer_dense.pyx Sat Sep 22 00:24:16 2007 +0100 @@ -1404,7 +1404,7 @@ cdef class Matrix_integer_dense(matrix_d #################################################################################### # LLL #################################################################################### - def lllgram(self): + def LLL_gram(self): """ LLL reduction of the lattice whose gram matrix is self. @@ -1425,7 +1425,7 @@ cdef class Matrix_integer_dense(matrix_d 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,11 +1434,11 @@ cdef class Matrix_integer_dense(matrix_d 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,20 +1464,21 @@ cdef class Matrix_integer_dense(matrix_d 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,50 +1487,116 @@ cdef class Matrix_integer_dense(matrix_d 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 = 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):