diff --git a/src/module_list.py b/src/module_list.py index 9bbc05b91c4..679f8b1dd27 100755 --- a/src/module_list.py +++ b/src/module_list.py @@ -129,6 +129,7 @@ def uname_specific(name, value, alternative): Extension('sage.algebras.quatalg.quaternion_algebra_cython', sources = ['sage/algebras/quatalg/quaternion_algebra_cython.pyx'], language='c++', + depends = flint_depends, libraries = ["flint", "gmp", "gmpxx", "m", "stdc++", "ntl"]), ################################ @@ -699,7 +700,7 @@ def uname_specific(name, value, alternative): language="c++", include_dirs = [SAGE_INC + '/fplll'], extra_compile_args=["-DFPLLL_V3_COMPAT"], - depends = [SAGE_INC + "/fplll/fplll.h"]), + depends = [SAGE_INC + "/fplll/fplll.h"] + flint_depends), Extension('sage.libs.linbox.linbox', sources = ['sage/libs/linbox/linbox.pyx'], @@ -750,7 +751,11 @@ def uname_specific(name, value, alternative): Extension('sage.libs.pari.pari_instance', sources = ["sage/libs/pari/pari_instance.pyx"], - libraries = ['pari', 'gmp']), + language="c", + extra_compile_args = ["-std=c99", "-D_XPG6"], + libraries = [ 'pari', 'gmp','flint'], + depends = flint_depends), + Extension('sage.libs.ppl', sources = ['sage/libs/ppl.pyx', 'sage/libs/ppl_shim.cc'], @@ -1052,13 +1057,12 @@ def uname_specific(name, value, alternative): sources = ['sage/matrix/matrix_integer_2x2.pyx'], libraries = ['gmp']), - # TODO -- change to use BLAS at some point. Extension('sage.matrix.matrix_integer_dense', sources = ['sage/matrix/matrix_integer_dense.pyx'], - extra_compile_args = m4ri_extra_compile_args, - depends = [SAGE_INC + '/m4ri/m4ri.h'], + extra_compile_args = ['-std=c99'] + m4ri_extra_compile_args, # order matters for cygwin!! - libraries = ['iml', 'pari', 'gmp', 'm', BLAS, BLAS2]), + libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2], + depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends), Extension('sage.matrix.matrix_integer_sparse', sources = ['sage/matrix/matrix_integer_sparse.pyx'], @@ -1080,7 +1084,11 @@ def uname_specific(name, value, alternative): Extension('sage.matrix.matrix_modn_dense', sources = ['sage/matrix/matrix_modn_dense.pyx'], - libraries = ['gmp']), + language="c++", + extra_compile_args = ["-D_XPG6"]+ m4ri_extra_compile_args, + # order matters for cygwin!! + libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2], + depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends), Extension('sage.matrix.matrix_modn_dense_float', sources = ['sage/matrix/matrix_modn_dense_float.pyx'], @@ -1092,7 +1100,7 @@ def uname_specific(name, value, alternative): sources = ['sage/matrix/matrix_modn_dense_double.pyx'], language="c++", libraries = ['linbox', 'givaro', 'mpfr', 'gmpxx', 'gmp', BLAS, BLAS2], - extra_compile_args = ['-DDISABLE_COMMENTATOR'] + givaro_extra_compile_args), + extra_compile_args = ["-D_XPG6" '-DDISABLE_COMMENTATOR']+ m4ri_extra_compile_args + givaro_extra_compile_args), Extension('sage.matrix.matrix_modn_sparse', sources = ['sage/matrix/matrix_modn_sparse.pyx'], @@ -1107,7 +1115,11 @@ def uname_specific(name, value, alternative): Extension('sage.matrix.matrix_rational_dense', sources = ['sage/matrix/matrix_rational_dense.pyx'], - libraries = ['pari', 'gmp']), + language="c", + extra_compile_args = ["-std=c99", "-D_XPG6"]+ m4ri_extra_compile_args, + # order matters for cygwin!! + libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2], + depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends), Extension('sage.matrix.matrix_rational_sparse', sources = ['sage/matrix/matrix_rational_sparse.pyx'], @@ -1918,7 +1930,6 @@ def uname_specific(name, value, alternative): flint_depends, libraries = ['flint', 'gmp', 'ratpoints']), - Extension('sage.schemes.elliptic_curves.period_lattice_region', sources = ['sage/schemes/elliptic_curves/period_lattice_region.pyx'], include_dirs = numpy_include_dirs), @@ -2095,13 +2106,13 @@ def uname_specific(name, value, alternative): Extension('sage.sat.solvers.satsolver', sources = ['sage/sat/solvers/satsolver.pyx']), - + ################################ - ## + ## ## sage.schemes ## ################################ - + Extension('sage.schemes.projective.projective_morphism_helper', sources = ['sage/schemes/projective/projective_morphism_helper.pyx']), ] diff --git a/src/sage/algebras/quatalg/quaternion_algebra_cython.pyx b/src/sage/algebras/quatalg/quaternion_algebra_cython.pyx index 913f698128b..124ffa4441b 100644 --- a/src/sage/algebras/quatalg/quaternion_algebra_cython.pyx +++ b/src/sage/algebras/quatalg/quaternion_algebra_cython.pyx @@ -88,21 +88,32 @@ def integral_matrix_and_denom_from_rational_quaternions(v, reverse=False): # Now fill in each row x of A, multiplying it by q = d/denom(x) cdef mpz_t q cdef mpz_t* row + cdef mpz_t tmp mpz_init(q) + mpz_init(tmp) for i in range(n): x = v[i] mpz_fdiv_q(q, d.value, x.d) if reverse: - mpz_mul(A._matrix[n-i-1][3], q, x.x) - mpz_mul(A._matrix[n-i-1][2], q, x.y) - mpz_mul(A._matrix[n-i-1][1], q, x.z) - mpz_mul(A._matrix[n-i-1][0], q, x.w) + mpz_mul(tmp, q, x.x) + A.set_unsafe_mpz(n-i-1,3,tmp) + mpz_mul(tmp, q, x.y) + A.set_unsafe_mpz(n-i-1,2,tmp) + mpz_mul(tmp, q, x.z) + A.set_unsafe_mpz(n-i-1,1,tmp) + mpz_mul(tmp, q, x.w) + A.set_unsafe_mpz(n-i-1,0,tmp) else: - mpz_mul(A._matrix[i][0], q, x.x) - mpz_mul(A._matrix[i][1], q, x.y) - mpz_mul(A._matrix[i][2], q, x.z) - mpz_mul(A._matrix[i][3], q, x.w) + mpz_mul(tmp, q, x.x) + A.set_unsafe_mpz(i,0,tmp) + mpz_mul(tmp, q, x.y) + A.set_unsafe_mpz(i,1,tmp) + mpz_mul(tmp, q, x.z) + A.set_unsafe_mpz(i,2,tmp) + mpz_mul(tmp, q, x.w) + A.set_unsafe_mpz(i,3,tmp) mpz_clear(q) + mpz_clear(tmp) return A, d def rational_matrix_from_rational_quaternions(v, reverse=False): @@ -195,6 +206,8 @@ def rational_quaternions_from_integral_matrix_and_denom(A, Matrix_integer_dense a = Integer(A.invariants()[0]) b = Integer(A.invariants()[1]) cdef Py_ssize_t i, j + cdef mpz_t tmp + mpz_init(tmp) if reverse: rng = range(H.nrows()-1,-1,-1) @@ -207,19 +220,28 @@ def rational_quaternions_from_integral_matrix_and_denom(A, Matrix_integer_dense mpz_set(x.a, a.value) mpz_set(x.b, b.value) if reverse: - mpz_init_set(x.x, H._matrix[i][3]) - mpz_init_set(x.y, H._matrix[i][2]) - mpz_init_set(x.z, H._matrix[i][1]) - mpz_init_set(x.w, H._matrix[i][0]) + H.get_unsafe_mpz(i,3,tmp) + mpz_init_set(x.x, tmp) + H.get_unsafe_mpz(i,2,tmp) + mpz_init_set(x.y, tmp) + H.get_unsafe_mpz(i,1,tmp) + mpz_init_set(x.z, tmp) + H.get_unsafe_mpz(i,0,tmp) + mpz_init_set(x.w, tmp) else: - mpz_init_set(x.x, H._matrix[i][0]) - mpz_init_set(x.y, H._matrix[i][1]) - mpz_init_set(x.z, H._matrix[i][2]) - mpz_init_set(x.w, H._matrix[i][3]) + H.get_unsafe_mpz(i,0,tmp) + mpz_init_set(x.x, tmp) + H.get_unsafe_mpz(i,1,tmp) + mpz_init_set(x.y, tmp) + H.get_unsafe_mpz(i,2,tmp) + mpz_init_set(x.z, tmp) + H.get_unsafe_mpz(i,3,tmp) + mpz_init_set(x.w, tmp) mpz_init_set(x.d, d.value) # WARNING -- we do *not* canonicalize the entries in the quaternion. This is # I think _not_ needed for quaternion_element.pyx v.append(x) + mpz_clear(tmp) return v diff --git a/src/sage/libs/flint/fmpz_mat.pxi b/src/sage/libs/flint/fmpz_mat.pxi new file mode 100644 index 00000000000..666060f5656 --- /dev/null +++ b/src/sage/libs/flint/fmpz_mat.pxi @@ -0,0 +1,32 @@ +include "sage/libs/flint/fmpz.pxi" +include "sage/libs/flint/fmpz_poly.pxi" + +cdef extern from "flint/fmpz_mat.h": + ctypedef void* fmpz_mat_t + void fmpz_mat_init(fmpz_mat_t mat, unsigned long rows, unsigned long cols) + void fmpz_mat_init_set(fmpz_mat_t mat, const fmpz_mat_t src) + void fmpz_mat_set(fmpz_mat_t result, fmpz_mat_t mat) + void fmpz_mat_clear(fmpz_mat_t mat) + int fmpz_mat_print_pretty(const fmpz_mat_t mat) + fmpz_t fmpz_mat_entry(fmpz_mat_t mat, long i, long j) + void fmpz_mat_one(fmpz_mat_t mat) + void fmpz_mat_zero(fmpz_mat_t mat) + void fmpz_mat_neg(fmpz_mat_t f, fmpz_mat_t g) + void fmpz_mat_scalar_mul_si(fmpz_mat_t B, const fmpz_mat_t A, long c) + void fmpz_mat_scalar_mul_fmpz(fmpz_mat_t B, const fmpz_mat_t A, const fmpz_t c) + void fmpz_mat_mul(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B) + void fmpz_mat_sqr(fmpz_mat_t B, const fmpz_mat_t A) + void fmpz_mat_add(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B) + void fmpz_mat_sub(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B) + void fmpz_mat_pow(fmpz_mat_t C, const fmpz_mat_t A, unsigned long n) + int fmpz_mat_is_zero(const fmpz_mat_t mat) + void fmpz_mat_charpoly(fmpz_poly_t cp, const fmpz_mat_t mat) + void fmpz_mat_det(fmpz_t det, const fmpz_mat_t A) + int fmpz_mat_inv(fmpz_mat_t Ainv, fmpz_t den, const fmpz_mat_t A) + void fmpz_mat_transpose(fmpz_mat_t B, const fmpz_mat_t A) + long fmpz_mat_rref(fmpz_mat_t B, fmpz_t den, const fmpz_mat_t A) + long fmpz_mat_fflu(fmpz_mat_t B, fmpz_poly_t den, long *perm, const fmpz_mat_t A, int rank_check) + long fmpz_mat_rref_fraction_free(long * perm, fmpz_mat_t B, fmpz_t den, const fmpz_mat_t A) + long fmpz_mat_rank(const fmpz_mat_t A) + int fmpz_mat_solve(fmpz_mat_t X, fmpz_t den, const fmpz_mat_t A, const fmpz_mat_t B) + long fmpz_mat_nullspace(fmpz_mat_t B, const fmpz_mat_t A) diff --git a/src/sage/libs/fplll/fplll.pxd b/src/sage/libs/fplll/fplll.pxd index 924855f614b..4f21421c90f 100644 --- a/src/sage/libs/fplll/fplll.pxd +++ b/src/sage/libs/fplll/fplll.pxd @@ -4,3 +4,12 @@ include "fplll.pxi" cdef class FP_LLL: cdef object fp_map cdef ZZ_mat[mpz_t] *_lattice + +cdef extern from "flint/fmpz.h": + ctypedef void* fmpz_t + void fmpz_get_mpz(mpz_t x, const fmpz_t f) + void fmpz_set_mpz(fmpz_t f, const mpz_t val) + +cdef extern from "flint/fmpz_mat.h": + ctypedef void* fmpz_mat_t + fmpz_t fmpz_mat_entry(fmpz_mat_t mat ,long i ,long j) diff --git a/src/sage/libs/fplll/fplll.pyx b/src/sage/libs/fplll/fplll.pyx index 2fd70912c2a..e7952e48779 100644 --- a/src/sage/libs/fplll/fplll.pyx +++ b/src/sage/libs/fplll/fplll.pyx @@ -139,19 +139,20 @@ cdef class FP_LLL: ValueError: fpLLL cannot handle matrices with zero rows. """ - cdef int i,j - - cdef Z_NR[mpz_t] t + cdef unsigned long i,j if A._nrows == 0: raise ValueError('fpLLL cannot handle matrices with zero rows.') self._lattice = new ZZ_mat[mpz_t](A._nrows,A._ncols) - - for i in range(A._nrows): - for j in range(A._ncols): - t.set(A._matrix[i][j]) - self._lattice[0][i][j] = t + cdef mpz_t tmp + mpz_init(tmp) + for i from 0 <= i < A._nrows: + for j from 0 <= j < A._ncols: + A.get_unsafe_mpz(i,j,tmp) + # mpz_set(self._lattice[0][i][j],tmp) + self._lattice[0][i][j].set(tmp) + mpz_clear(tmp) def __dealloc__(self): """ @@ -1461,5 +1462,5 @@ cdef to_sage(ZZ_mat[mpz_t] *A): for i from 0 <= i < A.getRows(): for j from 0 <= j < A.getCols(): - mpz_set(B._matrix[i][j], A[0][i][j].getData()) + B.set_unsafe_mpz(i,j,A[0][i][j].getData()) return B diff --git a/src/sage/libs/linkages/padics/mpz.pxi b/src/sage/libs/linkages/padics/mpz.pxi index 749b62902b7..0f24d8cd43a 100644 --- a/src/sage/libs/linkages/padics/mpz.pxi +++ b/src/sage/libs/linkages/padics/mpz.pxi @@ -77,7 +77,7 @@ cdef inline int ccmp(mpz_t a, mpz_t b, long prec, bint reduce_a, bint reduce_b, cdef int ans if reduce_a or reduce_b: mpz_sub(holder.value, a, b) - mpz_mod(holder.value, holder.value, prime_pow.pow_mpz_t_tmp(prec)) + mpz_mod(holder.value, holder.value, &prime_pow.pow_mpz_t_tmp(prec)[0]) return mpz_sgn(holder.value) else: ans = mpz_cmp(a,b) diff --git a/src/sage/libs/pari/pari_instance.pxd b/src/sage/libs/pari/pari_instance.pxd index 95db20db13d..bc4ad1f90cb 100644 --- a/src/sage/libs/pari/pari_instance.pxd +++ b/src/sage/libs/pari/pari_instance.pxd @@ -1,4 +1,6 @@ include 'decl.pxi' +include 'sage/libs/flint/fmpz.pxi' +include 'sage/libs/flint/fmpz_mat.pxi' cimport sage.structure.parent_base cimport cython @@ -27,8 +29,8 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase): cdef gen new_ref(self, GEN g, gen parent) cdef gen _empty_vector(self, long n) cdef long get_var(self, v) - cdef GEN _new_GEN_from_mpz_t_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc) - cdef GEN _new_GEN_from_mpz_t_matrix_rotate90(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc) - cdef gen integer_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf) + cdef GEN _new_GEN_from_fmpz_mat_t(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc) + cdef GEN _new_GEN_from_fmpz_mat_t_rotate90(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc) + cdef gen integer_matrix(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf) cdef GEN _new_GEN_from_mpq_t_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc) cdef gen rational_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc) diff --git a/src/sage/libs/pari/pari_instance.pyx b/src/sage/libs/pari/pari_instance.pyx index 21e3b367b15..aaf1e85f1b5 100644 --- a/src/sage/libs/pari/pari_instance.pyx +++ b/src/sage/libs/pari/pari_instance.pyx @@ -889,7 +889,7 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase): """ return objtogen(s) - cdef GEN _new_GEN_from_mpz_t_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc): + cdef GEN _new_GEN_from_fmpz_mat_t(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc): r""" Create a new PARI ``t_MAT`` with ``nr`` rows and ``nc`` columns from a ``mpz_t**``. @@ -900,13 +900,16 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase): cdef GEN x cdef GEN A = zeromatcopy(nr, nc) cdef Py_ssize_t i, j + cdef mpz_t tmp + mpz_init(tmp) for i in range(nr): for j in range(nc): - x = self._new_GEN_from_mpz_t(B[i][j]) + fmpz_get_mpz(tmp,fmpz_mat_entry(B,i,j)) + x = self._new_GEN_from_mpz_t(tmp) set_gcoeff(A, i+1, j+1, x) # A[i+1, j+1] = x (using 1-based indexing) return A - cdef GEN _new_GEN_from_mpz_t_matrix_rotate90(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc): + cdef GEN _new_GEN_from_fmpz_mat_t_rotate90(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc): r""" Create a new PARI ``t_MAT`` with ``nr`` rows and ``nc`` columns from a ``mpz_t**`` and rotate the matrix 90 degrees @@ -920,13 +923,16 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase): cdef GEN x cdef GEN A = zeromatcopy(nc, nr) cdef Py_ssize_t i, j + cdef mpz_t tmp + mpz_init(tmp) for i in range(nr): for j in range(nc): - x = self._new_GEN_from_mpz_t(B[i][nc-j-1]) + fmpz_get_mpz(tmp,fmpz_mat_entry(B,i,nc-j-1)) + x = self._new_GEN_from_mpz_t(tmp) set_gcoeff(A, j+1, i+1, x) # A[j+1, i+1] = x (using 1-based indexing) return A - cdef gen integer_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf): + cdef gen integer_matrix(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf): """ EXAMPLES:: @@ -936,9 +942,9 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase): pari_catch_sig_on() cdef GEN g if permute_for_hnf: - g = self._new_GEN_from_mpz_t_matrix_rotate90(B, nr, nc) + g = self._new_GEN_from_fmpz_mat_t_rotate90(B, nr, nc) else: - g = self._new_GEN_from_mpz_t_matrix(B, nr, nc) + g = self._new_GEN_from_fmpz_mat_t(B, nr, nc) return self.new_gen(g) cdef GEN _new_GEN_from_mpq_t_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc): diff --git a/src/sage/matrix/change_ring.pyx b/src/sage/matrix/change_ring.pyx index 64a76045b74..3c7957babdb 100644 --- a/src/sage/matrix/change_ring.pyx +++ b/src/sage/matrix/change_ring.pyx @@ -39,7 +39,7 @@ def integer_to_real_double_dense(Matrix_integer_dense A): S, None, None, None) for i from 0 <= i < A._nrows: for j from 0 <= j < A._ncols: - M.set_unsafe_double(i, j, mpz_get_d(A._matrix[i][j])) + M.set_unsafe_double(i,j,A.get_unsafe_double(i,j)) return M diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index ef3893dd0d2..21560cef21e 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -2635,6 +2635,8 @@ cdef class Matrix(matrix1.Matrix): - 'default' - allows the algorithm to be chosen automatically - 'generic' - naive algorithm usable for matrices over any field + - 'flint' - FLINT library code for matrices over the rationals + or the integers - 'pari' - PARI library code for matrices over number fields or the integers - 'padic' - padic algorithm from IML library for matrices @@ -2739,20 +2741,14 @@ cdef class Matrix(matrix1.Matrix): sage: B = copy(A).sparse_matrix() sage: set_verbose(1) sage: D = A.right_kernel(); D - verbose ... - verbose 1 () computing right kernel matrix over the rationals for 4x5 matrix - ... - verbose 1 () done computing right kernel matrix over the rationals for 4x5 matrix + verbose 1 () computing a right kernel for 4x5 matrix over Rational Field ... Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: [ 1 0 1 1/2 -1/2] [ 0 1 -1/2 -1/4 -1/4] sage: S = B.right_kernel(); S - verbose ... - verbose 1 () computing right kernel matrix over the rationals for 4x5 matrix - ... - verbose 1 () done computing right kernel matrix over the rationals for 4x5 matrix + verbose 1 () computing a right kernel for 4x5 matrix over Rational Field ... Vector space of degree 5 and dimension 2 over Rational Field Basis matrix: @@ -3002,8 +2998,9 @@ cdef class Matrix(matrix1.Matrix): sage: B = copy(A).sparse_matrix() sage: set_verbose(1) sage: D = A.right_kernel(); D - verbose ... + verbose 1 () computing a right kernel for 4x7 matrix over Integer Ring verbose 1 () computing right kernel matrix over the integers for 4x7 matrix + ... verbose 1 () done computing right kernel matrix over the integers for 4x7 matrix ... Free module of degree 7 and rank 3 over Integer Ring @@ -3012,8 +3009,9 @@ cdef class Matrix(matrix1.Matrix): [ 0 35 0 25 -1 -31 17] [ 0 0 7 12 -3 -1 -8] sage: S = B.right_kernel(); S - verbose ... + verbose 1 () computing a right kernel for 4x7 matrix over Integer Ring verbose 1 () computing right kernel matrix over the integers for 4x7 matrix + ... verbose 1 () done computing right kernel matrix over the integers for 4x7 matrix ... Free module of degree 7 and rank 3 over Integer Ring @@ -3154,10 +3152,12 @@ cdef class Matrix(matrix1.Matrix): algorithm = kwds.pop('algorithm', None) if algorithm is None: algorithm = 'default' - elif not algorithm in ['default', 'generic', 'pari', 'padic', 'pluq']: + elif not algorithm in ['default', 'generic', 'flint', 'pari', 'padic', 'pluq']: raise ValueError("matrix kernel algorithm '%s' not recognized" % algorithm ) elif algorithm == 'padic' and not (is_IntegerRing(R) or is_RationalField(R)): raise ValueError("'padic' matrix kernel algorithm only available over the rationals and the integers, not over %s" % R) + elif algorithm == 'flint' and not (is_IntegerRing(R) or is_RationalField(R)): + raise ValueError("'flint' matrix kernel algorithm only available over the rationals and the integers, not over %s" % R) elif algorithm == 'pari' and not (is_IntegerRing(R) or (is_NumberField(R) and not is_RationalField(R))): raise ValueError("'pari' matrix kernel algorithm only available over non-trivial number fields and the integers, not over %s" % R) elif algorithm == 'generic' and not R.is_field(): @@ -3278,6 +3278,8 @@ cdef class Matrix(matrix1.Matrix): - 'default' - allows the algorithm to be chosen automatically - 'generic' - naive algorithm usable for matrices over any field + - 'flint' - FLINT library code for matrices over the rationals + or the integers - 'pari' - PARI library code for matrices over number fields or the integers - 'padic' - padic algorithm from IML library for matrices @@ -3620,6 +3622,8 @@ cdef class Matrix(matrix1.Matrix): - 'default' - allows the algorithm to be chosen automatically - 'generic' - naive algorithm usable for matrices over any field + - 'flint' - FLINT library code for matrices over the rationals + or the integers - 'pari' - PARI library code for matrices over number fields or the integers - 'padic' - padic algorithm from IML library for matrices @@ -6640,7 +6644,7 @@ cdef class Matrix(matrix1.Matrix): going along each row. INPUT: - + - ``check`` -- (default: ``False``) If ``True`` return a tuple of the maximal matrix and the permutations taking taking ``self`` to the maximal matrix. @@ -6671,8 +6675,8 @@ cdef class Matrix(matrix1.Matrix): sage: M.permutation_normal_form(check=True) ( - [ 5 -1] - [ 4 2] + [ 5 -1] + [ 4 2] [ 3 -1], ((1,2,3), (1,2)) ) @@ -6754,7 +6758,7 @@ cdef class Matrix(matrix1.Matrix): n = [j for j in range(nrows - l) if SN[j] == nb] # Now compare to our previous max if b < nb: - # Bigger so save line + # Bigger so save line b = nb m = [n] # and delete all previous attempts diff --git a/src/sage/matrix/matrix_integer_dense.pxd b/src/sage/matrix/matrix_integer_dense.pxd index 5de47386749..eee7dabc25b 100644 --- a/src/sage/matrix/matrix_integer_dense.pxd +++ b/src/sage/matrix/matrix_integer_dense.pxd @@ -1,41 +1,48 @@ include "sage/ext/cdefs.pxi" +include "sage/libs/ntl/decl.pxi" +include "sage/libs/flint/fmpz.pxi" +include "sage/libs/flint/fmpz_mat.pxi" +include "sage/libs/flint/fmpz_poly.pxi" -from sage.ext.mod_int cimport * cimport matrix_dense cimport sage.rings.integer from sage.rings.integer cimport Integer +from sage.ext.mod_int cimport * ctypedef long* GEN cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): cdef char _initialized - cdef mpz_t *_entries - cdef mpz_t **_matrix + cdef char _initialized_linbox + cdef mpz_t * _entries + cdef mpz_t ** _rows + cdef fmpz_mat_t _matrix cdef object _pivots cdef int mpz_height(self, mpz_t height) except -1 cpdef double _log_avg_sq1(self) except -1.0 cdef _mod_int_c(self, mod_int modulus) cdef _mod_two(self) - cdef _zero_out_matrix(self) - cdef _new_unitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols) cdef _pickle_version0(self) cdef _unpickle_version0(self, data) cpdef _export_as_string(self, int base=?) - cdef _init_linbox(self) + cdef _dealloc_linbox(self) cdef void reduce_entry_unsafe(self, Py_ssize_t i, Py_ssize_t j, Integer modulus) + cdef set_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, const mpz_t value) + cdef set_unsafe_si(self, Py_ssize_t i, Py_ssize_t j, const long long value) + cdef get_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, mpz_t value) + cdef double get_unsafe_double(self, Py_ssize_t i, Py_ssize_t j) # HNF Modn cdef int _hnf_modn(Matrix_integer_dense self, Matrix_integer_dense res, mod_int det) except -1 cdef long long* _hnf_modn_impl(Matrix_integer_dense self, mod_int det, Py_ssize_t nrows, Py_ssize_t ncols) except NULL - cdef _new_uninitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols) + cdef Matrix_integer_dense _new_uninitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols) cdef extract_hnf_from_pari_matrix(self, GEN H, int flag, bint include_zero_rows) -cdef int four_dim_det(mpz_t, mpz_t *) except -1 cpdef _lift_crt(Matrix_integer_dense M, residues, moduli=*) ################################################################ diff --git a/src/sage/matrix/matrix_integer_dense.pyx b/src/sage/matrix/matrix_integer_dense.pyx index 2cec8c7f833..bb10b440e70 100644 --- a/src/sage/matrix/matrix_integer_dense.pyx +++ b/src/sage/matrix/matrix_integer_dense.pyx @@ -7,6 +7,8 @@ AUTHORS: - Robert Bradshaw +- Marc Masdeu (August 2014). Implemented using FLINT. + EXAMPLES:: sage: a = matrix(ZZ, 3,3, range(9)); a @@ -31,13 +33,8 @@ TESTS:: sage: TestSuite(a).run() sage: Matrix(ZZ,0,0).inverse() [] -""" -########## *** IMPORTANT *** -# If you're working on this code, we *always* assume that -# self._matrix[i] = self._entries[i*self._ncols] -# !!!!!!!! Do not break this! -# This is assumed in the _rational_kernel_iml +""" ###################################################################### # Copyright (C) 2006,2007 William Stein @@ -82,7 +79,6 @@ include "sage/ext/random.pxi" cdef extern from "math.h": double log(double x) - from sage.ext.multi_modular import MultiModularBasis from sage.ext.multi_modular cimport MultiModularBasis @@ -93,6 +89,7 @@ from sage.rings.integer_ring import ZZ, IntegerRing_class from sage.rings.integer_ring cimport IntegerRing_class from sage.rings.finite_rings.integer_mod_ring import IntegerModRing from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing +from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint from sage.structure.element cimport ModuleElement, RingElement, Element, Vector from sage.structure.element import is_Vector from sage.structure.sequence import Sequence @@ -118,6 +115,8 @@ cimport sage.structure.element import sage.matrix.matrix_space as matrix_space +from sage.ext.mod_int cimport * + ################ # Used for modular HNF from sage.rings.fast_arith cimport arith_int @@ -156,7 +155,7 @@ fplll_fp_map = {None: None, cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse r""" - Matrix over the integers. + Matrix over the integers, implemented using FLINT. On a 32-bit machine, they can have at most `2^{32}-1` rows or columns. On a 64-bit machine, matrices can have at most @@ -215,38 +214,32 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse self._nrows = parent.nrows() self._ncols = parent.ncols() self._pivots = None - - # Allocate an array where all the entries of the matrix are stored. - sig_on() - self._entries = sage_malloc(sizeof(mpz_t) * (self._nrows * self._ncols)) - sig_off() - if self._entries == NULL: - raise MemoryError("out of memory allocating a matrix") - - # Allocate an array of pointers to the rows, which is useful for - # certain algorithms. - ################################## - # *IMPORTANT*: FOR MATRICES OVER ZZ, WE ALWAYS ASSUME THAT - # THIS ARRAY IS *not* PERMUTED. This should be OK, since all - # algorithms are multi-modular. - ################################## - self._matrix = sage_malloc(sizeof(mpz_t*)*self._nrows) - if self._matrix == NULL: - sage_free(self._entries) - self._entries = NULL - raise MemoryError("out of memory allocating a matrix") - - # Set each of the pointers in the array self._matrix to point - # at the memory for the corresponding row. - cdef Py_ssize_t i, k - k = 0 - for i from 0 <= i < self._nrows: - self._matrix[i] = self._entries + k - k = k + self._ncols + self._initialized_linbox = False + self._entries = NULL + self._rows = NULL cdef _init_linbox(self): + cdef Py_ssize_t i, j, k + if not self._initialized_linbox: + sig_on() + self._rows = sage_malloc(sizeof(mpz_t*)*self._nrows) + self._entries = sage_malloc(sizeof(mpz_t) * self._nrows * self._ncols) + k = 0 + for i from 0 <= i < self._nrows: + self._rows[i] = self._entries + k + for j from 0 <= j < self._ncols: + mpz_init(self._entries[k]) + fmpz_get_mpz(self._entries[k],fmpz_mat_entry(self._matrix,i,j)) + k += 1 + linbox.set(self._rows, self._nrows, self._ncols) + sig_off() + self._initialized_linbox = True + + cdef _dealloc_linbox(self): + cdef Py_ssize_t k sig_on() - linbox.set(self._matrix, self._nrows, self._ncols) + for k from 0 <= k < self._nrows * self._ncols: + mpz_clear(self._entries[k]) sig_off() def __copy__(self): @@ -265,14 +258,15 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse True """ cdef Matrix_integer_dense A - A = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, - 0, 0, 0) + A = self._new_uninitialized_matrix(self._nrows,self._ncols) + cdef Py_ssize_t i sig_on() - for i from 0 <= i < self._nrows * self._ncols: - mpz_init_set(A._entries[i], self._entries[i]) + fmpz_mat_set(A._matrix,self._matrix) sig_off() A._initialized = True + if self._initialized_linbox: + A._init_linbox() if self._subdivisions is not None: A.subdivide(*self.subdivisions()) return A @@ -301,17 +295,18 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse def __dealloc__(self): """ - Frees all the memory allocated for this matrix. EXAMPLE:: + Frees all the memory allocated for this matrix. + + EXAMPLE:: sage: a = Matrix(ZZ,2,[1,2,3,4]) sage: del a """ - cdef Py_ssize_t i - if self._initialized: - for i from 0 <= i < (self._nrows * self._ncols): - mpz_clear(self._entries[i]) - sage_free(self._entries) - sage_free(self._matrix) + fmpz_mat_clear(self._matrix) + if self._initialized_linbox: + self._dealloc_linbox() + sage_free(self._rows) + sage_free(self._entries) def __init__(self, parent, entries, copy, coerce): r""" @@ -387,7 +382,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: v.parent() Full MatrixSpace of 1 by 100000 dense matrices over Integer Ring """ - cdef Py_ssize_t i, j + cdef Py_ssize_t i, j, k cdef bint is_list cdef Integer x @@ -404,46 +399,43 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse else: entries = list(entries) is_list = 1 - + fmpz_mat_init(self._matrix,self._nrows,self._ncols) if is_list: - # Create the matrix whose entries are in the given entry list. if len(entries) != self._nrows * self._ncols: - sage_free(self._entries) - sage_free(self._matrix) - self._entries = NULL raise TypeError("entries has the wrong length") if coerce: - for i from 0 <= i < self._nrows * self._ncols: - x = ZZ(entries[i]) - # todo -- see integer.pyx and the TODO there; perhaps this could be - # sped up by creating a mpz_init_set_sage function. - mpz_init_set(self._entries[i], x.value) + k = 0 + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + x = ZZ(entries[k]) + k += 1 + # todo -- see integer.pyx and the TODO there; perhaps this could be + # sped up by creating a mpz_init_set_sage function. + fmpz_set_mpz(fmpz_mat_entry(self._matrix, i, j),(x).value) self._initialized = True else: - for i from 0 <= i < self._nrows * self._ncols: - mpz_init_set(self._entries[i], ( entries[i]).value) + k = 0 + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + fmpz_set_mpz(fmpz_mat_entry(self._matrix, i,j),( entries[k]).value) + k += 1 self._initialized = True else: - # If x is zero, make the zero matrix and be done. if mpz_sgn(x.value) == 0: self._zero_out_matrix() + self._initialized = True return # the matrix must be square: if self._nrows != self._ncols: - sage_free(self._entries) - sage_free(self._matrix) - self._entries = NULL raise TypeError("nonzero scalar matrix must be square") # Now we set all the diagonal entries to x and all other entries to 0. self._zero_out_matrix() - j = 0 for i from 0 <= i < self._nrows: - mpz_set(self._entries[j], x.value) - j = j + self._nrows + 1 + fmpz_set_mpz(fmpz_mat_entry(self._matrix,i,i), x.value) self._initialized = True cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value): @@ -473,11 +465,78 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [10 1 2] [ 3 4 5] """ - mpz_set(self._matrix[i][j], (value).value) + fmpz_set_mpz(fmpz_mat_entry(self._matrix,i,j),(value).value) + self._initialized_linbox = False + # if self._initialized_linbox: + # mpz_set(self._entries[i*self._ncols + j],(value).value) + + cdef set_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, const mpz_t value): + """ + Set position i,j of this matrix to x. + + (VERY UNSAFE - value *must* be of type Integer). + + INPUT: + + + - ``i`` - row + + - ``j`` - column + + - ``value`` - The value to set self[i,j] to. value + MUST be of type Integer + + + EXAMPLES:: + + sage: a = matrix(ZZ,2,3, range(6)); a + [0 1 2] + [3 4 5] + sage: a[0,0] = 10 + sage: a + [10 1 2] + [ 3 4 5] + """ + fmpz_set_mpz(fmpz_mat_entry(self._matrix,i,j),value) + self._initialized_linbox = False + # if self._initialized_linbox: + # mpz_set(self._entries[i*self._ncols + j],value) + + cdef set_unsafe_si(self, Py_ssize_t i, Py_ssize_t j, const long long value): + """ + Set position i,j of this matrix to x. + + (VERY UNSAFE - value *must* be of type Integer). + + INPUT: + + + - ``i`` - row + + - ``j`` - column + + - ``value`` - The value to set self[i,j] to. value + MUST be of type Integer + + + EXAMPLES:: + + sage: a = matrix(ZZ,2,3, range(6)); a + [0 1 2] + [3 4 5] + sage: a[0,0] = 10 + sage: a + [10 1 2] + [ 3 4 5] + """ + fmpz_set_si(fmpz_mat_entry(self._matrix,i,j),value) + self._initialized_linbox = False + # if self._initialized_linbox: + # mpz_set_si(self._entries[i*self._ncols + j],value) cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): """ - Returns (j, i) entry of self as a new Integer. + Returns (i, j) entry of self as a new Integer. .. warning:: @@ -501,9 +560,62 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef Integer z z = PY_NEW(Integer) - mpz_set(z.value, self._matrix[i][j]) + fmpz_get_mpz(z.value,fmpz_mat_entry(self._matrix, i, j)) return z + cdef get_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, mpz_t value): + """ + Returns (i, j) entry of self as a new Integer. + + .. warning:: + + This is very unsafe; it assumes i and j are in the right + range. + + EXAMPLES:: + + sage: a = MatrixSpace(ZZ,3)(range(9)); a + [0 1 2] + [3 4 5] + [6 7 8] + sage: a[1,2] + 5 + sage: a[4,7] + Traceback (most recent call last): + ... + IndexError: matrix index out of range + sage: a[-1,0] + 6 + """ + # mpz_init(value) + fmpz_get_mpz(value,fmpz_mat_entry(self._matrix, i, j)) + + cdef double get_unsafe_double(self, Py_ssize_t i, Py_ssize_t j): + """ + Returns (j, i) entry of self as a new Integer. + + .. warning:: + + This is very unsafe; it assumes i and j are in the right + range. + + EXAMPLES:: + + sage: a = MatrixSpace(ZZ,3)(range(9)); a + [0 1 2] + [3 4 5] + [6 7 8] + sage: a[1,2] + 5 + sage: a[4,7] + Traceback (most recent call last): + ... + IndexError: matrix index out of range + sage: a[-1,0] + 6 + """ + return fmpz_get_d(fmpz_mat_entry(self._matrix, i, j)) + def _pickle(self): """ EXAMPLES:: @@ -547,7 +659,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # from sec 5.14 of the GMP manual. ?? cdef int i, j, len_so_far, m, n cdef char *a - cdef char *s, *t, *tmp + cdef char *s + cdef char *t + cdef char *tmp if self._nrows == 0 or self._ncols == 0: data = '' @@ -558,24 +672,26 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse len_so_far = 0 sig_on() - for i from 0 <= i < self._nrows * self._ncols: - m = mpz_sizeinbase (self._entries[i], base) - if len_so_far + m + 2 >= n: - # copy to new string with double the size - n = 2*n + m + 1 - tmp = sage_malloc(n * sizeof(char)) - strcpy(tmp, s) - sage_free(s) - s = tmp - t = s + len_so_far - #endif - mpz_get_str(t, base, self._entries[i]) - m = strlen(t) - len_so_far = len_so_far + m + 1 - t = t + m - t[0] = 32 - t[1] = 0 - t = t + 1 + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + # mat_entry = fmpz_mat_entry(self._matrix,i,j) + m = fmpz_sizeinbase(fmpz_mat_entry(self._matrix,i,j), base) + if len_so_far + m + 2 >= n: + # copy to new string with double the size + n = 2*n + m + 1 + tmp = sage_malloc(n * sizeof(char)) + strcpy(tmp, s) + sage_free(s) + s = tmp + t = s + len_so_far + #endif + fmpz_get_str(t, base, fmpz_mat_entry(self._matrix,i,j)) + m = strlen(t) + len_so_far = len_so_far + m + 1 + t = t + m + t[0] = 32 + t[1] = 0 + t = t + 1 sig_off() data = str(s)[:-1] sage_free(s) @@ -588,15 +704,20 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse raise RuntimeError("unknown matrix version (=%s)"%version) cdef _unpickle_version0(self, data): - cdef Py_ssize_t i, n + cdef Py_ssize_t i, j, n, k + fmpz_mat_init(self._matrix,self._nrows,self._ncols) data = data.split() n = self._nrows * self._ncols if len(data) != n: raise RuntimeError("invalid pickle data.") - for i from 0 <= i < n: - s = data[i] - if mpz_init_set_str(self._entries[i], s, 32): - raise RuntimeError("invalid pickle data") + k = 0 + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + s = data[k] + k += 1 + # fmpz_init(fmpz_mat_entry(self._matrix,i,j)) + if fmpz_set_str(fmpz_mat_entry(self._matrix,i,j), s, 32): + raise RuntimeError("invalid pickle data") self._initialized = True @@ -607,32 +728,30 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # LEVEL 1 helpers: # These function support the implementation of the level 1 functionality. ######################################################################## + cdef Matrix_integer_dense _new_uninitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols): + """ + Return a new matrix over the integers from given parent + All memory is allocated for this matrix, but its + entries have not yet been filled in. + """ + cdef object P = matrix_space.MatrixSpace(ZZ, nrows, ncols, sparse=False) + cdef Matrix_integer_dense ans = Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None) + fmpz_mat_init(ans._matrix,nrows,ncols) + if ans._matrix == NULL: + raise MemoryError("out of memory allocating a matrix") + return ans + cdef _zero_out_matrix(self): """ Set this matrix to be the zero matrix. This is only for internal use. (Note: this matrix must NOT already have initialised entries.) """ - # TODO: This is about 6-10 slower than MAGMA doing what seems to be the same thing. - # Moreover, NTL can also do this quickly. Why? I think both have specialized - # small integer classes. (dmharvey: yes, NTL does not allocate any memory when - # initializing a ZZ to zero.) sig_on() - cdef Py_ssize_t i - for i from 0 <= i < self._nrows * self._ncols: - mpz_init(self._entries[i]) + fmpz_mat_zero(self._matrix) sig_off() self._initialized = True - - cdef _new_unitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols): - """ - Return a new matrix over the integers with the given number of rows - and columns. All memory is allocated for this matrix, but its - entries have not yet been filled in. - """ - P = self._parent.matrix_space(nrows, ncols) - return Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None) - + self._initialized_linbox = False ######################################################################## # LEVEL 2 functionality @@ -682,15 +801,10 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: a.__nonzero__() False """ - cdef mpz_t *a, *b - cdef Py_ssize_t i, j - cdef int k - for i from 0 <= i < self._nrows * self._ncols: - if mpz_sgn(self._entries[i]): - return True - return False + return not fmpz_mat_is_zero(self._matrix) - def _multiply_linbox(self, Matrix right): + + def _multiply_linbox(self, Matrix_integer_dense right): """ Multiply matrices over ZZ using linbox. @@ -709,16 +823,32 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [14 50] """ cdef int e - cdef Matrix_integer_dense ans, B - B = right - ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols()) - self._init_linbox() + cdef long int i,j + cdef Matrix_integer_dense ans + cdef Matrix_integer_dense left = self + + if self._nrows == right._nrows: + # self acts on the space of right + parent = right.parent() + if self._ncols == right._ncols: + # right acts on the space of self + parent = self.parent() + else: + parent = self.matrix_space(left._nrows, right._ncols) + + right._init_linbox() + ans = self._new_uninitialized_matrix(parent.nrows(),parent.ncols()) + ans._init_linbox() + left._init_linbox() sig_on() - linbox.matrix_matrix_multiply(ans._matrix, B._matrix, B._nrows, B._ncols) + linbox.matrix_matrix_multiply(ans._rows, right._rows, right._nrows, right._ncols) + for i from 0 <= i < ans._nrows: + for j from 0 <= j < ans._ncols: + fmpz_set_mpz(fmpz_mat_entry(ans._matrix,i,j),ans._rows[i][j]) sig_off() return ans - def _multiply_classical(self, Matrix right): + def _multiply_classical(self, Matrix_integer_dense right): """ EXAMPLE:: @@ -733,8 +863,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse if self._ncols != right._nrows: raise IndexError("Number of columns of self must equal number of rows of right.") - cdef Py_ssize_t i, j, k, l, nr, nc, snc - cdef mpz_t *v + cdef Py_ssize_t i, j, k, nr, nc, snc cdef object parent nr = self._nrows @@ -753,46 +882,50 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef Matrix_integer_dense M, _right _right = right - M = Matrix_integer_dense.__new__(Matrix_integer_dense, parent, None, None, None) - Matrix.__init__(M, parent) - - # M has memory allocated but entries are not initialized - cdef mpz_t *entries - entries = M._entries - - cdef mpz_t s, z - mpz_init(s) - mpz_init(z) + M = self._new_uninitialized_matrix(parent.nrows(),parent.ncols()) + cdef fmpz_t s + fmpz_init(s) sig_on() - l = 0 for i from 0 <= i < nr: for j from 0 <= j < nc: - mpz_set_si(s,0) # set s = 0 - v = self._matrix[i] + fmpz_set_si(s,0) # set s = 0 for k from 0 <= k < snc: - mpz_mul(z, v[k], _right._matrix[k][j]) - mpz_add(s, s, z) - mpz_init_set(entries[l], s) - l += 1 + fmpz_addmul(s, fmpz_mat_entry(self._matrix,i,k), fmpz_mat_entry(_right._matrix,k,j)) + fmpz_set(fmpz_mat_entry(M._matrix,i,j),s) sig_off() - mpz_clear(s) - mpz_clear(z) + fmpz_clear(s) M._initialized = True return M cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right): + cdef Py_ssize_t i, j, k, l, nr, nc, snc + cdef object parent + + if self._ncols != right._nrows: + raise IndexError("Number of columns of self must equal number of rows of right.") + + nr = self._nrows + nc = right._ncols + snc = self._ncols - ############# - # see the tune_multiplication function below. - n = max(self._nrows, self._ncols, right._nrows, right._ncols) - if n <= 20: - return self._multiply_classical(right) - a = self.height(); b = right.height() - if float(max(a,b)) / float(n) >= 0.70: - return self._multiply_classical(right) + if self._nrows == right._nrows: + # self acts on the space of right + parent = right.parent() + if self._ncols == right._ncols: + # right acts on the space of self + parent = self.parent() else: - return self._multiply_multi_modular(right) + parent = self.matrix_space(nr, nc) + + cdef Matrix_integer_dense M = self._new_uninitialized_matrix(parent.nrows(),parent.ncols()) + + sig_on() + fmpz_mat_mul(M._matrix,self._matrix,(right)._matrix) + sig_off() + M._initialized = True + return M + cpdef ModuleElement _lmul_(self, RingElement right): """ @@ -805,12 +938,14 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef Py_ssize_t i cdef Integer _x + cdef fmpz_t z _x = Integer(right) cdef Matrix_integer_dense M - M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None) - for i from 0 <= i < self._nrows * self._ncols: - mpz_init(M._entries[i]) - mpz_mul(M._entries[i], self._entries[i], _x.value) + M = self._new_uninitialized_matrix(self._nrows,self._ncols) + sig_on() + fmpz_set_mpz(z,_x.value) + fmpz_mat_scalar_mul_fmpz(M._matrix,self._matrix,z) + sig_off() M._initialized = True return M @@ -835,17 +970,10 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef Py_ssize_t i, j cdef Matrix_integer_dense M - M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None) + M = self._new_uninitialized_matrix(self._nrows,self._ncols) sig_on() - cdef mpz_t *row_self, *row_right, *row_ans - for i from 0 <= i < self._nrows: - row_self = self._matrix[i] - row_right = ( right)._matrix[i] - row_ans = M._matrix[i] - for j from 0 <= j < self._ncols: - mpz_init(row_ans[j]) - mpz_add(row_ans[j], row_self[j], row_right[j]) + fmpz_mat_add(M._matrix,self._matrix,( right)._matrix) sig_off() M._initialized = True return M @@ -866,17 +994,10 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef Py_ssize_t i, j cdef Matrix_integer_dense M - M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None) + M = self._new_uninitialized_matrix(self._nrows,self._ncols) sig_on() - cdef mpz_t *row_self, *row_right, *row_ans - for i from 0 <= i < self._nrows: - row_self = self._matrix[i] - row_right = ( right)._matrix[i] - row_ans = M._matrix[i] - for j from 0 <= j < self._ncols: - mpz_init(row_ans[j]) - mpz_sub(row_ans[j], row_self[j], row_right[j]) + fmpz_mat_sub(M._matrix,self._matrix,( right)._matrix) sig_off() M._initialized = True return M @@ -898,23 +1019,23 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: Matrix(ZZ, [[5, 10], [20, 30]]).__cmp__(Matrix(ZZ, [[0, 10], [25, 30]])) 1 """ - cdef mpz_t *a, *b cdef Py_ssize_t i, j cdef int k + sig_on() for i from 0 <= i < self._nrows: - a = self._matrix[i] - b = (right)._matrix[i] for j from 0 <= j < self._ncols: - k = mpz_cmp(a[j], b[j]) + k = fmpz_cmp(fmpz_mat_entry(self._matrix,i,j),fmpz_mat_entry((right)._matrix,i,j)) if k: + sig_off() if k < 0: return -1 else: return 1 + sig_off() return 0 - + # TODO: Implement better cdef Vector _vector_times_matrix_(self, Vector v): """ Returns the vector times matrix product. @@ -937,19 +1058,25 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef Vector_integer_dense w, ans cdef Py_ssize_t i, j - cdef mpz_t x + cdef fmpz_t x + cdef fmpz_t z M = self._row_ambient_module() w = v ans = M.zero_vector() - mpz_init(x) + sig_on() + fmpz_init(x) + fmpz_init(z) for i from 0 <= i < self._ncols: - mpz_set_si(x, 0) + fmpz_set_si(x, 0) for j from 0 <= j < self._nrows: - mpz_addmul(x, w._entries[j], self._matrix[j][i]) - mpz_set(ans._entries[i], x) - mpz_clear(x) + fmpz_set_mpz(z,w._entries[j]) + fmpz_addmul(x, z, fmpz_mat_entry(self._matrix,j,i)) + fmpz_get_mpz(ans._entries[i], x) + fmpz_clear(x) + fmpz_clear(z) + sig_off() return ans @@ -981,14 +1108,14 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ return self, ZZ(1) - def charpoly(self, var='x', algorithm='linbox'): + def charpoly(self, var='x', algorithm='generic'): """ INPUT: - ``var`` - a variable name - - ``algorithm`` - 'linbox' (default) 'generic' + - ``algorithm`` - 'generic' (default), 'flint' or 'linbox' .. note:: @@ -1028,25 +1155,31 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse x^2 - 3*x - 2 """ + cdef long i,n + cdef fmpz_t x + cdef Integer z + cdef Polynomial_integer_dense_flint g + if algorithm == 'generic': + algorithm = 'linbox' cache_key = 'charpoly_%s' % algorithm g = self.fetch(cache_key) if g is not None: return g.change_variable_name(var) - if algorithm == 'linbox' and not USE_LINBOX_POLY: - algorithm = 'generic' - if algorithm == 'linbox': + if algorithm == 'flint' or (algorithm == 'linbox' and not USE_LINBOX_POLY): + g = PolynomialRing(ZZ,names = var).gen() + sig_on() + fmpz_mat_charpoly(g.__poly,self._matrix) + sig_off() + elif algorithm == 'linbox': g = self._charpoly_linbox(var) - elif algorithm == 'generic': - g = matrix_dense.Matrix_dense.charpoly(self, var) else: raise ValueError("no algorithm '%s'"%algorithm) - self.cache(cache_key, g) return g - def minpoly(self, var='x', algorithm='linbox'): + def minpoly(self, var='x', algorithm = 'linbox'): """ INPUT: @@ -1074,16 +1207,15 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse x = self.fetch(key) if x: return x + if algorithm == 'linbox' and not USE_LINBOX_POLY: algorithm = 'generic' - if algorithm == 'linbox': g = self._minpoly_linbox(var) elif algorithm == 'generic': g = matrix_dense.Matrix_dense.minpoly(self, var) else: raise ValueError("no algorithm '%s'"%algorithm) - self.cache(key, g) return g @@ -1104,6 +1236,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``var`` - 'x' - ``typ`` - 'minpoly' or 'charpoly' + """ time = verbose('computing %s of %s x %s matrix using linbox'%(typ, self._nrows, self._ncols)) if self._nrows != self._ncols: @@ -1145,7 +1278,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef mpz_t h cdef Integer x + sig_on() self.mpz_height(h) + sig_off() x = Integer() x.set_from_mpz(h) mpz_clear(h) @@ -1166,23 +1301,22 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse OUTPUT: sets the value of height to the height of this matrix, i.e., the max absolute value of the entries of the matrix. """ - cdef mpz_t x, h - cdef Py_ssize_t i + cdef fmpz_t x,h + cdef Py_ssize_t i,j sig_on() - mpz_init_set_si(h, 0) - mpz_init(x) - - for i from 0 <= i < self._nrows * self._ncols: - mpz_abs(x, self._entries[i]) - if mpz_cmp(h, x) < 0: - mpz_set(h, x) - - mpz_init_set(height, h) - mpz_clear(h) - mpz_clear(x) + fmpz_init_set_ui(h, 0) + fmpz_init(x) + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + fmpz_abs(x, fmpz_mat_entry(self._matrix,i,j)) + if fmpz_cmp(h, x) < 0: + fmpz_set(h, x) + mpz_init(height) + fmpz_get_mpz(height,h) + fmpz_clear(h) + fmpz_clear(x) sig_off() - return 0 # no error occurred. cpdef double _log_avg_sq1(self) except -1.0: @@ -1202,32 +1336,32 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: matrix(ZZ,10)._log_avg_sq1() 0.0 """ - cdef Py_ssize_t i - cdef Py_ssize_t N = self._nrows * self._ncols - + cdef unsigned long i, j + cdef unsigned long N = self._nrows * self._ncols # wsq * 4^wsq_e = sum of entries squared plus number of entries cdef double wsq = N cdef long wsq_e = 0 cdef double d - cdef long e + cdef long e = 0 sig_on() - for i in range(N): - d = mpz_get_d_2exp(&e, self._entries[i]) - if (e > wsq_e): - wsq = ldexp(wsq, 2*(wsq_e - e)) - wsq_e = e - elif (e < wsq_e): - d = ldexp(d, e - wsq_e) - wsq += d*d + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + d = fmpz_get_d_2exp(&e,fmpz_mat_entry(self._matrix,i,j)) + if (e > wsq_e): + wsq = ldexp(wsq, 2*(wsq_e - e)) + wsq_e = e + elif (e < wsq_e): + d = ldexp(d, e - wsq_e) + wsq += d*d sig_off() # Compute log(wsq * 4^wsq_e / N) = # log(wsq * 4^wsq_e / N) = log(wsq/N) + wsq_e * log(4) return log(wsq/N) + wsq_e * log(4.0) - def _multiply_multi_modular(left, Matrix_integer_dense right): + def _multiply_multi_modular(self, Matrix_integer_dense right): """ Multiply this matrix by ``left`` using a multi modular algorithm. @@ -1243,8 +1377,18 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [-463 -490] """ cdef Integer h + cdef Matrix_integer_dense left = self cdef mod_int *moduli cdef int i, n, k + cdef object parent + + nr = left._nrows + nc = right._ncols + snc = left._ncols + + + cdef Matrix_integer_dense result + h = left.height() * right.height() * left.ncols() verbose('multiplying matrices of height %s and %s'%(left.height(),right.height())) @@ -1256,14 +1400,27 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse t = cputime() res[i] *= res_right[i] verbose('multiplied matrices modulo a prime (%s/%s)'%(i+1,k), t) - result = left.new_matrix(left.nrows(), right.ncols()) + result = left.new_matrix(nr,nc) _lift_crt(result, res, mm) # changes result return result cdef void reduce_entry_unsafe(self, Py_ssize_t i, Py_ssize_t j, Integer modulus): # Used for p-adic matrices. - if mpz_cmp(self._matrix[i][j], modulus.value) >= 0 or mpz_cmp_ui(self._matrix[i][j], 0) < 0: - mpz_mod(self._matrix[i][j], self._matrix[i][j], modulus.value) + cdef fmpz_t z + cdef fmpz_t mod + cdef fmpz_t *aij = fmpz_mat_entry(self._matrix,i,j) + sig_on() + fmpz_init(mod) + fmpz_set_mpz(mod,modulus.value) + fmpz_init(z) + if fmpz_cmp(aij[0], mod) >= 0 or fmpz_cmp_ui(aij[0], 0) < 0: + fmpz_mod(z, aij[0], mod) + fmpz_set(aij[0],z) + if self._initialized_linbox: + fmpz_get_mpz(self._entries[i*self._ncols + j],z) + fmpz_clear(z) + fmpz_clear(mod) + sig_off() def _mod_int(self, modulus): """ @@ -1306,20 +1463,18 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse res_f = Matrix_modn_dense_float.__new__(Matrix_modn_dense_float, matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None) for i from 0 <= i < self._nrows: - self_row = self._matrix[i] res_row_f = res_f._matrix[i] for j from 0 <= j < self._ncols: - res_row_f[j] = mpz_fdiv_ui(self_row[j], p) + res_row_f[j] = fmpz_fdiv_ui(fmpz_mat_entry(self._matrix,i,j), p) return res_f elif p < MAX_MODULUS_DOUBLE: res_d = Matrix_modn_dense_double.__new__(Matrix_modn_dense_double, matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None) for i from 0 <= i < self._nrows: - self_row = self._matrix[i] res_row_d = res_d._matrix[i] for j from 0 <= j < self._ncols: - res_row_d[j] = mpz_fdiv_ui(self_row[j], p) + res_row_d[j] = fmpz_fdiv_ui(fmpz_mat_entry(self._matrix,i,j), p) return res_d else: raise ValueError("p to big.") @@ -1351,48 +1506,39 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef size_t i, k, n cdef Py_ssize_t nr, nc - + cdef mpz_t tmp + mpz_init(tmp) n = len(mm) nr = self._nrows nc = self._ncols - cdef mod_int **row_list - row_list = sage_malloc(sizeof(mod_int*) * n) - if row_list == NULL: + cdef mod_int *entry_list + entry_list = sage_malloc(sizeof(mod_int) * n) + if entry_list == NULL: raise MemoryError("out of memory allocating multi-modular coefficient list") - for k in range(n): - row_list[k] = sage_malloc(sizeof(mod_int)*nc) - if row_list[k] == NULL: - for i in range(k): - sage_free(row_list[i]) - sage_free(row_list) - raise MemoryError("out of memory allocating multi-modular coefficient list") sig_on() for i from 0 <= i < nr: - for k from 0 <= k < n: - mm.mpz_reduce_vec(self._matrix[i], row_list, nc) - if isinstance(res[k], Matrix_modn_dense_float): - for j in range(nc): - (res[k])._matrix[i][j] = (row_list[k][j])%(res[k]).p - else: - for j in range(nc): - (res[k])._matrix[i][j] = (row_list[k][j])%(res[k]).p + for j from 0 <= j < nc: + self.get_unsafe_mpz(i,j,tmp) + mm.mpz_reduce(tmp, entry_list) + for k from 0 <= k < n: + if isinstance(res[k], Matrix_modn_dense_float): + (res[k])._matrix[i][j] = (entry_list[k])%(res[k]).p + else: + (res[k])._matrix[i][j] = (entry_list[k])%(res[k]).p sig_off() - - for k in range(n): - sage_free(row_list[k]) - sage_free(row_list) + mpz_clear(tmp) + sage_free(entry_list) return res def _echelon_in_place_classical(self): cdef Matrix_integer_dense E + sig_on() E = self.echelon_form() - - cdef int i - for i from 0 <= i < self._ncols * self._nrows: - mpz_set(self._entries[i], E._entries[i]) - + sig_off() + fmpz_mat_set(self._matrix,E._matrix) + self._initialized_linbox = False self.clear_cache() def _echelon_strassen(self): @@ -1481,7 +1627,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse transformation=False, D=None): r""" Return the echelon form of this matrix over the integers, also known - as the hermit normal form (HNF). + as the hermite normal form (HNF). INPUT: @@ -1686,6 +1832,8 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse ans = self.fetch(key) if ans is not None: return ans + cdef Matrix_integer_dense H_m,w + cdef fmpz_t den cdef Py_ssize_t nr, nc, n, i, j nr = self._nrows nc = self._ncols @@ -1702,14 +1850,13 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse if self.height().ndigits() > 10000 or n >= 50: pari_big = 1 - cdef Matrix_integer_dense H_m - proof = get_proof_flag(proof, "linear_algebra") pivots = None rank = None if algorithm == "padic": import matrix_integer_dense_hnf + self._init_linbox() if transformation: H_m, U, pivots = matrix_integer_dense_hnf.hnf_with_transformation(self, proof=proof) if not include_zero_rows: @@ -1761,27 +1908,49 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse v[i,j] = self.get_unsafe(nr-i-1,nc-j-1) try: - w = v.HNF(D=D) + w1 = v.HNF(D=D) except RuntimeError: # HNF may fail if a nxm matrix has rank < m raise ValueError("ntl only computes HNF for square matrices of full rank.") - rank = w.nrows() + rank = w1.nrows() if include_zero_rows: H_m = self.new_matrix() else: - H_m = self.new_matrix(nrows=w.nrows()) + H_m = self.new_matrix(nrows=w1.nrows()) - nr = w.nrows() - nc = w.ncols() + nr = w1.nrows() + nc = w1.ncols() - for i from 0 <= i < w.nrows(): - for j from 0 <= j < w.ncols(): - H_m[i,j] = w[nr-i-1,nc-j-1] + for i from 0 <= i < w1.nrows(): + for j from 0 <= j < w1.ncols(): + H_m[i,j] = w1[nr-i-1,nc-j-1] + elif algorithm == 'flint': + raise NotImplementedError,'Not yet implemented' + # This is how it should like, but fmpz_mat_rref_fraction_free() + # looks like is missing! -- Marc Masdeu + if transformation: + raise ValueError("transformation matrix only available with p-adic algorithm") + w = self.new_matrix() + sig_on() + fmpz_init(den) + rank = fmpz_mat_rref_fraction_free(NULL , w._matrix,den,self._matrix) + fmpz_clear(den) + sig_off() + if include_zero_rows: + H_m = self.new_matrix() + else: + H_m = self.new_matrix(nrows=rank) + nr = H_m.nrows() + nc = H_m.ncols() + for i from 0 <= i < nr: + for j from 0 <= j < nc: + fmpz_set(fmpz_mat_entry(H_m._matrix,i,j),fmpz_mat_entry(w._matrix,i,j)) else: raise TypeError("algorithm '%s' not understood"%(algorithm)) + H_m._initialized = True H_m.set_immutable() if pivots is None: from matrix_integer_dense_hnf import pivots_of_hnf_matrix @@ -1959,13 +2128,11 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # Now we determine the pivots from the matrix E as quickly as we can. # For each row, we find the first nonzero position in that row -- it is the pivot. cdef Py_ssize_t i, j, k - cdef mpz_t *row p = [] k = 0 for i from 0 <= i < E._nrows: - row = E._matrix[i] # pointer to ith row for j from k <= j < E._ncols: - if mpz_cmp_si(row[j], 0) != 0: # nonzero position + if fmpz_cmp_si(fmpz_mat_entry(E._matrix,i,j), 0) != 0: # nonzero position p.append(j) k = j+1 # so start at next position next time break @@ -2250,9 +2417,10 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``algorithm`` - determines which algorithm to use, options are: + - 'flint' - use the algorithm from the FLINT library - 'pari' - use the ``matkerint()`` function from the PARI library - 'padic' - use the p-adic algorithm from the IML library - - 'default' - use a heuristic to decide which of the two above + - 'default' - use a heuristic to decide which of the three above routines is fastest. This is the default value. - ``proof`` - this is passed to the p-adic IML algorithm. @@ -2261,11 +2429,11 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse OUTPUT: Returns a pair. First item is the string is either - 'computed-pari-int' or 'computed-iml-int', which identifies + 'computed-flint-int', 'computed-pari-int', 'computed-flint-int', which identifies the nature of the basis vectors. Second item is a matrix whose rows are a basis for the right kernel, - over the integers, as computed by either the IML or PARI libraries. + over the integers, as computed by either the FLINT, IML or PARI libraries. EXAMPLES:: @@ -2294,20 +2462,20 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: result = A._right_kernel_matrix(algorithm='default') sage: result[0] - 'computed-pari-int' + 'computed-flint-int' sage: result[1] - [-26 31 -30 21 2 -10] - [-47 -13 48 -14 -11 18] + [ 469 -214 30 -119 37 0] + [-370 165 -18 91 -30 2] - sage: result = A._right_kernel_matrix() + sage: result = A._right_kernel_matrix(algorithm='flint') sage: result[0] - 'computed-pari-int' + 'computed-flint-int' sage: result[1] - [-26 31 -30 21 2 -10] - [-47 -13 48 -14 -11 18] + [ 469 -214 30 -119 37 0] + [-370 165 -18 91 -30 2] With the 'default' given as the algorithm, several heuristics are - used to determine if PARI or IML ('padic') is used. The code has + used to determine if FLINT, PARI or IML ('padic') is used. The code has exact details, but roughly speaking, relevant factors are: the absolute size of the matrix, or the relative dimensions, or the magnitude of the entries. :: @@ -2327,13 +2495,12 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse TESTS: - We test three trivial cases. PARI is used for small matrices, + We test three trivial cases. FLINT is used for small matrices, but we let the heuristic decide that. :: sage: A = matrix(ZZ, 0, 2) sage: A._right_kernel_matrix()[1] - [1 0] - [0 1] + [] sage: A = matrix(ZZ, 2, 0) sage: A._right_kernel_matrix()[1].parent() Full MatrixSpace of 0 by 0 dense matrices over Integer Ring @@ -2356,10 +2523,12 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # *dramatically* better than doing absolutely nothing # (i.e., always choosing 'padic'), but is of course # far from optimal. -- William Stein + + # I sometimes favor FLINT over PARI, but this should be better tuned. -- Marc Masdeu if max(self._nrows, self._ncols) <= 10: - # pari much better for very small matrices, as long as entries aren't huge. - algorithm = 'pari' - if max(self._nrows, self._ncols) <= 50: + # Use FLINT for very small matrices, as long as entries aren't huge. + algorithm = 'flint' + elif max(self._nrows, self._ncols) <= 50: # when entries are huge, padic relatively good. h = self.height().ndigits() if h < 100: @@ -2374,17 +2543,25 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse else: algorithm = 'pari' - if algorithm == 'pari': + if algorithm == 'flint': + proof = kwds.pop('proof', None) + proof = get_proof_flag(proof, "linear_algebra") + K = self._rational_kernel_flint().transpose().saturation(proof=proof) + format = 'computed-flint-int' + elif algorithm == 'pari': K = self._pari_().matkerint().mattranspose().python() format = 'computed-pari-int' - if algorithm == 'padic': + elif algorithm == 'padic': proof = kwds.pop('proof', None) proof = get_proof_flag(proof, "linear_algebra") K = self._rational_kernel_iml().transpose().saturation(proof=proof) format = 'computed-iml-int' + else: + raise ValueError('unknown algorithm: %s'%algorithm) tm = verbose("done computing right kernel matrix over the integers for %sx%s matrix" % (self.nrows(), self.ncols()),level=1, t=tm) return (format, K) + # TODO: implement using flint function def _adjoint(self): """ Return the adjoint of this matrix. @@ -2722,7 +2899,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``use_givens`` -- (default: ``False``) use Givens orthogonalization only applicable to approximate reductions and NTL; this is more - stable but slower + stable but slower - ``use_siegel`` -- (default: ``False``) use Siegel's condition instead of Lovasz's condition, ignored by NTL @@ -2989,23 +3166,23 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef Py_ssize_t c, row - cdef mpz_t s, pr - mpz_init(s) - mpz_init(pr) + cdef fmpz_t s,pr + fmpz_init(s) + fmpz_init(pr) - mpz_set_si(pr, 1) + fmpz_set_si(pr, 1) for row from 0 <= row < self._nrows: - mpz_set_si(s, 0) + fmpz_set_si(s, 0) for c in cols: if c<0 or c >= self._ncols: raise IndexError("matrix column index out of range") - mpz_add(s, s, self._matrix[row][c]) - mpz_mul(pr, pr, s) + fmpz_add(s, s, fmpz_mat_entry(self._matrix,row,c)) + fmpz_mul(pr, pr, s) cdef Integer z z = PY_NEW(Integer) - mpz_set(z.value, pr) - mpz_clear(s) - mpz_clear(pr) + fmpz_get_mpz(z.value, pr) + fmpz_clear(s) + fmpz_clear(pr) return z def _linbox_sparse(self): @@ -3013,8 +3190,8 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse v = ['%s %s M'%(self._nrows, self._ncols)] for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: - if mpz_cmp_si(self._matrix[i][j], 0): - v.append('%s %s %s'%(i+1,j+1,self.get_unsafe(i,j))) + if fmpz_cmp_si(fmpz_mat_entry(self._matrix,i,j), 0): + v.append('%s %s %s'%(i+1,j+1,self.get_unsafe_double(i,j))) v.append('0 0 0\n') return '\n'.join(v) @@ -3023,7 +3200,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse v = ['%s %s x'%(self._nrows, self._ncols)] for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: - v.append(str(self.get_unsafe(i,j))) + v.append(str(self.get_unsafe_double(i,j))) return ' '.join(v) def rational_reconstruction(self, N): @@ -3136,6 +3313,8 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef randstate rstate = current_randstate() + cdef mpz_t tmp + mpz_init(tmp) cdef Py_ssize_t i, j, k, nc, num_per_row global state, ZZ @@ -3145,17 +3324,20 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # Original code, before adding the ``nonzero`` option. sig_on() if density == 1: - for i from 0 <= i < self._nrows*self._ncols: - the_integer_ring._randomize_mpz(self._entries[i], x, y, \ + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + the_integer_ring._randomize_mpz(tmp, x, y, \ distribution) + self.set_unsafe_mpz(i,j,tmp) else: nc = self._ncols num_per_row = int(density * nc) for i from 0 <= i < self._nrows: for j from 0 <= j < num_per_row: k = rstate.c_random()%nc - the_integer_ring._randomize_mpz(self._matrix[i][k], \ + the_integer_ring._randomize_mpz(tmp, \ x, y, distribution) + self.set_unsafe_mpz(i,k,tmp) sig_off() else: # New code, to implement the ``nonzero`` option. Note that this @@ -3163,24 +3345,29 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # each entry is set until it's non-zero. sig_on() if density == 1: - for i from 0 <= i < self._nrows*self._ncols: - while mpz_sgn(self._entries[i]) == 0: - the_integer_ring._randomize_mpz(self._entries[i], \ - x, y, distribution) + for i from 0 <= i < self._nrows: + for j from 0 <= j < self._ncols: + while fmpz_sgn(fmpz_mat_entry(self._matrix,i,j)) == 0: + the_integer_ring._randomize_mpz(tmp, \ + x, y, distribution) + self.set_unsafe_mpz(i,j,tmp) else: nc = self._ncols num_per_row = int(density * nc) for i from 0 <= i < self._nrows: for j from 0 <= j < num_per_row: k = rstate.c_random() % nc - while mpz_sgn(self._matrix[i][k]) == 0: - the_integer_ring._randomize_mpz(self._matrix[i][k],\ + while fmpz_sgn(fmpz_mat_entry(self._matrix,i,k)) == 0: + the_integer_ring._randomize_mpz(tmp,\ x, y, distribution) + self.set_unsafe_mpz(i,k,tmp) + sig_off() + mpz_clear(tmp) #### Rank - def rank(self): + def rank(self, algorithm = 'modp'): """ Return the rank of this matrix. @@ -3189,12 +3376,13 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``nonnegative integer`` - the rank + - ``algorithm`` - either ``'modp'`` (default) or ``'flint'`` or ``'linbox'`` .. note:: The rank is cached. - ALGORITHM: First check if the matrix has maxim possible rank by + ALGORITHM: If set to ``'modp'``, first check if the matrix has maximum possible rank by working modulo one random prime. If not call LinBox's rank function. @@ -3220,14 +3408,17 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse r = self.fetch('rank') if not r is None: return r - if self._nrows <= 6 and self._ncols <= 6 and self.height().ndigits() <= 10000: - r = self._rank_pari() - # Can very quickly detect full rank by working modulo p. - r = self._rank_modp() - if r == self._nrows or r == self._ncols: + if algorithm == 'flint' or (self._nrows <= 6 and self._ncols <= 6 and self.height().ndigits() <= 100): + r = fmpz_mat_rank(self._matrix) self.cache('rank', r) return r - # Detecting full rank didn't work -- use LinBox's general algorithm. + elif algorithm == 'modp': + # Can very quickly detect full rank by working modulo p. + r = self._rank_modp() + if r == self._nrows or r == self._ncols: + self.cache('rank', r) + return r + # Algorithm is 'linbox' or detecting full rank didn't work -- use LinBox's general algorithm. r = self._rank_linbox() self.cache('rank', r) return r @@ -3260,6 +3451,8 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``'default'`` -- automatically determine which algorithm to use depending on the matrix. + - ``'flint'`` -- let flint do the determinant + - ``'padic'`` - uses a p-adic / multimodular algorithm that relies on code in IML and linbox @@ -3349,7 +3542,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse 80. The check that the determinant of a squared matrix is a square is a sanity check that the result is probably correct:: - sage: for s in [1..80]: # long time (6s on sage.math, 2013) + sage: for s in [1..200]: # long time (6s on sage.math, 2013) ... M = random_matrix(ZZ, s) ... d = (M*M).determinant() ... assert d.is_square() @@ -3361,21 +3554,21 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse raise ValueError("self must be a square matrix") cdef Py_ssize_t n = self.nrows() + cdef Integer det = Integer() + cdef fmpz_t e + cdef double difficulty - cdef Integer det4x4 - if n <= 4: - if n <= 3: - # Use generic special cased code. - return matrix_dense.Matrix_dense.determinant(self) - else: - det4x4 = Integer() - four_dim_det(det4x4.value, self._entries) - return det4x4 + if n <=4 or algorithm == 'flint': + fmpz_init(e) + fmpz_mat_det(e, self._matrix) + fmpz_get_mpz(det.value,e) + fmpz_clear(e) + return det proof = get_proof_flag(proof, "linear_algebra") - cdef double difficulty if algorithm == 'default': + algorithm = 'flint' # These heuristics are by Jeroen Demeyer (#14007). There # is no mathematics behind this, it was experimentally # observed to work well. I tried various estimates for @@ -3387,13 +3580,22 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # mostly small entries, but it is never much faster than # padic (and it only works with proof=False), so we never # default to using linbox. - difficulty = (self._log_avg_sq1() + 2.0) * (n * n) - if difficulty <= 800: - algorithm = 'pari' - elif n <= 48 or (proof and n <= 72) or (proof and n <= 400 and difficulty <= 600000): - algorithm = 'ntl' - else: - algorithm = 'padic' + # I favor FLINT instead of PARI. -- Marc Masdeu. + # (the thresholds have to be recalculated) + #difficulty = (self._log_avg_sq1() + 2.0) * (n * n) + #if difficulty <= 800: + # algorithm = 'flint' + #elif n <= 48 or (proof and n <= 72) or (proof and n <= 400 and difficulty <= 600000): + # algorithm = 'ntl' + #else: + # algorithm = 'padic' + + if algorithm == 'flint': + fmpz_init(e) + fmpz_mat_det(e, self._matrix) + fmpz_get_mpz(det.value,e) + fmpz_clear(e) + return det if algorithm == 'padic': import matrix_integer_dense_hnf @@ -3412,6 +3614,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse self.cache('det', d) return d + def _det_linbox(self): """ Compute the determinant of this matrix using Linbox. @@ -3422,6 +3625,31 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sig_off() return Integer(d) + def _det_pari(self, int flag=0): + """ + Determinant of this matrix using Gauss-Bareiss. If (optional) + flag is set to 1, use classical Gaussian elimination. + + For efficiency purposes, this det is computed entirely on the + PARI stack then the PARI stack is cleared. This function is + most useful for very small matrices. + + EXAMPLES:: + + sage: matrix(ZZ,3,[1..9])._det_pari() + 0 + sage: matrix(ZZ,3,[1..9])._det_pari(1) + 0 + """ + pari_catch_sig_on() + cdef GEN d = det0(pari_GEN(self), flag) + # now convert d to a Sage integer e + cdef Integer e = Integer() + t_INT_to_ZZ(e.value, d) + pari.clear_stack() + return e + + def _det_ntl(self): """ Compute the determinant of this matrix using NTL. @@ -3447,25 +3675,58 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse return self.matrix_space(self._ncols, 0).zero_matrix() cdef long dim + cdef unsigned long i,j,k cdef mpz_t *mp_N time = verbose('computing null space of %s x %s matrix using IML'%(self._nrows, self._ncols)) + self._init_linbox() sig_on() - dim = nullspaceMP (self._nrows, self._ncols, self._entries, &mp_N) + dim = nullspaceMP(self._nrows, self._ncols, self._entries, &mp_N) sig_off() - P = self.matrix_space(self._ncols, dim) - # Now read the answer as a matrix. cdef Matrix_integer_dense M - M = Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None) - for i from 0 <= i < dim*self._ncols: - mpz_init_set(M._entries[i], mp_N[i]) - mpz_clear(mp_N[i]) + M = self._new_uninitialized_matrix(self._ncols,dim) + k = 0 + for i from 0 <= i < self._ncols: + for j from 0 <= j < dim: + fmpz_set_mpz(fmpz_mat_entry(M._matrix,i,j), mp_N[k]) + mpz_clear(mp_N[k]) + k += 1 sage_free(mp_N) - verbose("finished computing null space", time) M._initialized = True return M + #### Rational kernel, via flint + def _rational_kernel_flint(self): + """ + Return the rational kernel of this matrix (acting from the + left), considered as a matrix over QQ. I.e., returns a matrix K + such that self\*K = 0, and the number of columns of K equals the + nullity of self. + + AUTHORS: + + - Marc Masdeu + """ + if self._nrows == 0 or self._ncols == 0: + return self.matrix_space(self._ncols, 0).zero_matrix() + + cdef long dim + cdef fmpz_mat_t M0 + sig_on() + fmpz_mat_init(M0,self._ncols,self._ncols) + dim = fmpz_mat_nullspace(M0, self._matrix) + sig_off() + # Now read the answer as a matrix. + cdef Matrix_integer_dense M + M = self._new_uninitialized_matrix(self._ncols,dim) + for i from 0 <= i < M._nrows: + for j from 0 <= j < M._ncols: + fmpz_set(fmpz_mat_entry(M._matrix,i,j),fmpz_mat_entry(M0,i,j)) + M._initialized = True + fmpz_mat_clear(M0) + return M + def _invert_iml(self, use_nullspace=False, check_invertible=True): """ Invert this matrix using IML. The output matrix is an integer @@ -3531,7 +3792,72 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse raise ZeroDivisionError("input matrix must be nonsingular") return self._solve_iml(P.identity_matrix(), right=True) - def _solve_right_nonsingular_square(self, B, check_rank=True): + def _invert_flint(self, check_invertible=True): + """ + Invert this matrix using FLINT. The output matrix is an integer + matrix and a denominator. + + INPUT: + + + - ``self`` - an invertible matrix + + - ``check_invertible`` - (default: True) whether to + check that the matrix is invertible. + + + OUTPUT: A, d such that A\*self = d + + + - ``A`` - a matrix over ZZ + + - ``d`` - an integer + + + EXAMPLES:: + + sage: a = matrix(ZZ,3,[1,2,5, 3,7,8, 2,2,1]) + sage: b, d = a._invert_flint(); b,d + ( + [ 9 -8 19] + [-13 9 -7] + [ 8 -2 -1], 23 + ) + sage: a*b + [23 0 0] + [ 0 23 0] + [ 0 0 23] + + AUTHORS: + + - William Stein + + - Marc Masdeu -- (08/2014) Use FLINT + """ + if self._nrows != self._ncols: + raise TypeError("self must be a square matrix.") + + cdef Matrix_integer_dense M + cdef int res + cdef Integer den = Integer(0) + cdef fmpz_t fden + fmpz_init(fden) + M = self._new_uninitialized_matrix(self._nrows,self._ncols) + verbose('computing inverse of %s x %s matrix using FLINT'%(self._nrows, self._ncols)) + sig_on() + res = fmpz_mat_inv(M._matrix,fden,self._matrix) + fmpz_get_mpz(den.value,fden) + sig_off() + M._initialized = True + fmpz_clear(fden) + if check_invertible and res == 0: + raise ZeroDivisionError('Matrix is singular') + if fmpz_cmp_si(fden,0) < 0: + return -M,-den + else: + return M,den + + def _solve_right_nonsingular_square(self, B, check_rank=True, algorithm = 'iml'): r""" If self is a matrix `A` of full rank, then this function returns a vector or matrix `X` such that `A X = B`. @@ -3558,6 +3884,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``check_rank`` - bool (default: True); if True verify that in fact the rank is full. + - ``algorithm`` - ``'iml'`` (default) or ``'flint'`` OUTPUT: a matrix or vector over `\QQ` @@ -3618,8 +3945,17 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse sage: x = a \ v sage: a * x == v True + + sage: n = 100 + sage: a = random_matrix(ZZ,n) + sage: v = vector(ZZ,n,range(n)) + sage: x = a._solve_right_nonsingular_square(v,algorithm = 'flint') + sage: a * x == v + True + """ - t = verbose('starting IML solve_right...') + t = verbose('starting %s solve_right...'%algorithm) + # It would probably be much better to rewrite linbox so it # throws an error instead of ** going into an infinite loop ** # in the non-full rank case. In any case, we do this for now, @@ -3654,13 +3990,18 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse verbose('finished solve_right (via inverse)', t) return X - X, d = self._solve_iml(C, right=True) + if algorithm == 'flint': + X, d = self._solve_flint(C, right=True) + elif algorithm == 'iml': # iml + X, d = self._solve_iml(C, right = True) + else: + raise ValueError("Unknown algorithm '%s'"%algorithm) if d != 1: X = (1/d) * X if not matrix: # Convert back to a vector X = (X.base_ring() ** X.nrows())(X.list()) - verbose('finished solve_right via IML', t) + verbose('finished solve_right via %s'%algorithm, t) return X def _solve_iml(self, Matrix_integer_dense B, right=True): @@ -3753,7 +4094,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - Martin Albrecht """ - cdef int i + cdef unsigned long i, j, k cdef mpz_t *mp_N, mp_D cdef Matrix_integer_dense M cdef Integer D @@ -3762,7 +4103,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # This is *required* by the IML function we call below. raise ValueError("self must be a square matrix") - if self.nrows() == 1: + if self._nrows == 1: return B, self[0,0] cdef SOLU_POS solu_pos @@ -3800,25 +4141,28 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse solu_pos = LeftSolu sig_check() - + verbose("Initializing mp_N and mp_D") mp_N = sage_malloc( n * m * sizeof(mpz_t) ) for i from 0 <= i < n * m: mpz_init(mp_N[i]) mpz_init(mp_D) - + verbose("Done with initializing mp_N and mp_D") + self._init_linbox() + B._init_linbox() try: + verbose('Calling solver n = %s, m = %s'%(n,m)) sig_on() nonsingSolvLlhsMM(solu_pos, n, m, self._entries, B._entries, mp_N, mp_D) sig_off() - - M = Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None) - for i from 0 <= i < n*m: - mpz_init_set(M._entries[i], mp_N[i]) + M = self._new_uninitialized_matrix(P.nrows(),P.ncols()) + k = 0 + for i from 0 <= i < M._nrows: + for j from 0 <= j < M._ncols: + fmpz_set_mpz(fmpz_mat_entry(M._matrix,i,j), mp_N[k]) + k += 1 M._initialized = True - D = PY_NEW(Integer) mpz_set(D.value, mp_D) - return M, D finally: mpz_clear(mp_D) @@ -3826,7 +4170,140 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse mpz_clear(mp_N[i]) sage_free(mp_N) - def _rational_echelon_via_solve(self): + def _solve_flint(self, Matrix_integer_dense B, right=True): + """ + Let A equal self be a square matrix. Given B return an integer + matrix C and an integer d such that self C\*A == d\*B if right is + False or A\*C == d\*B if right is True. + + OUTPUT: + + + - ``C`` - integer matrix + + - ``d`` - integer denominator + + + EXAMPLES:: + + sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) + sage: B = matrix(ZZ,3,4, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) + sage: C,d = A._solve_flint(B,right=False); C + [ 6 -18 -15 27] + [ 0 24 24 -36] + [ 4 -12 -6 -2] + + :: + + sage: d + 12 + + :: + + sage: C*A == d*B + True + + :: + + sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) + sage: B = matrix(ZZ,4,3, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) + sage: C,d = A._solve_flint(B) + sage: C + [ 12 40 28] + [-12 -4 -4] + [ -6 -25 -16] + [ 12 34 16] + + :: + + sage: d + 12 + + :: + + sage: A*C == d*B + True + + Test wrong dimensions:: + + sage: A = random_matrix(ZZ, 4, 4) + sage: B = random_matrix(ZZ, 2, 3) + sage: B._solve_flint(A) + Traceback (most recent call last): + ... + ValueError: self must be a square matrix + sage: A._solve_flint(B, right=False) + Traceback (most recent call last): + ... + ArithmeticError: B's number of columns must match self's number of rows + sage: A._solve_flint(B, right=True) + Traceback (most recent call last): + ... + ArithmeticError: B's number of rows must match self's number of columns + + Check that this can be interrupted properly (:trac:`15453`):: + + sage: A = random_matrix(ZZ, 2000, 2000) + sage: B = random_matrix(ZZ, 2000, 2000) + sage: t0 = walltime() + sage: alarm(2); A._solve_flint(B) # long time + Traceback (most recent call last): + ... + AlarmInterrupt + sage: t = walltime(t0) + sage: t < 10 or t + True + + AUTHORS: + + - Marc Masdeu (08/2014) following _solve_iml implementation of Martin Albrecht + """ + cdef Matrix_integer_dense M + cdef fmpz_t tmp + cdef Integer den = Integer(0) + if self._nrows != self._ncols: + # This is *required* by the FLINT function we call below. + raise ValueError("self must be a square matrix") + + if self.nrows() == 1: + return B, self[0,0] + + if right: + fmpz_init(tmp) + if self._ncols != B._nrows: + raise ArithmeticError("B's number of rows must match self's number of columns") + + n = self._ncols + m = B._ncols + + if m == 0 or n == 0: + return self.new_matrix(nrows = n, ncols = m), Integer(1) + M = self._new_uninitialized_matrix(self._ncols,B._ncols) + sig_on() + fmpz_mat_solve(M._matrix,tmp,self._matrix,B._matrix) + fmpz_get_mpz(den.value,tmp) + fmpz_clear(tmp) + if mpz_cmp_si(den.value,0) < 0: + mpz_neg(den.value,den.value) + fmpz_mat_neg(M._matrix,M._matrix) + sig_off() + M._initialized = True + return M,den + + else: # left + if self._nrows != B._ncols: + raise ArithmeticError("B's number of columns must match self's number of rows") + + n = self._nrows + m = B._nrows + + if m == 0 or n == 0: + return self.new_matrix(nrows = n, ncols = m), Integer(1) + + M,d = self.transpose()._solve_flint(B.transpose(), right=True) + return M.transpose(),d + + def _rational_echelon_via_solve(self, solver = 'iml'): r""" Computes information that gives the reduced row echelon form (over QQ!) of a matrix with integer entries. @@ -3836,7 +4313,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``self`` - a matrix over the integers. - + - ``solver`` - either ``'iml'`` (default) or ``'flint'`` OUTPUT: @@ -3980,9 +4457,12 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse pivots_ = set(pivots) non_pivots = [i for i in range(B.ncols()) if not i in pivots_] D = B.matrix_from_columns(non_pivots) - t = verbose('calling IML solver', level=2, caller_name='p-adic echelon') - X, d = C._solve_iml(D, right=True) - t = verbose('finished IML solver', level=2, caller_name='p-adic echelon', t=t) + t = verbose('calling %s solver'%solver, level=2, caller_name='p-adic echelon') + if solver == 'iml': + X, d = C._solve_iml(D, right=True) + else: + X, d = C._solve_flint(D, right=True) + t = verbose('finished %s solver'%solver, level=2, caller_name='p-adic echelon', t=t) # Step 6: Verify that we had chosen the correct pivot columns. pivots_are_right = True @@ -4222,15 +4702,6 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse # Hermite form modulo D # This code below is by E. Burcin. Thanks! ##################################################################################### - cdef _new_uninitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols): - """ - Return a new matrix over the integers with the given number of rows - and columns. All memory is allocated for this matrix, but its - entries have not yet been filled in. - """ - P = self._parent.matrix_space(nrows, ncols) - return Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None) - def _hnf_mod(self, D): """ INPUT: @@ -4246,7 +4717,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse - ``matrix`` - the Hermite normal form of self. """ t = verbose('hermite mod %s'%D, caller_name='matrix_integer_dense') - cdef Matrix_integer_dense res = self._new_uninitialized_matrix(self._nrows, self._ncols) + cdef Matrix_integer_dense res = self._new_uninitialized_matrix(self._nrows,self._ncols) self._hnf_modn(res, D) verbose('finished hnf mod', t, caller_name='matrix_integer_dense') return res @@ -4262,7 +4733,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse k = 0 for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: - mpz_init_set_si(res._matrix[i][j], res_l[k]) + res.set_unsafe_si(i,j,res_l[k]) k += 1 res._initialized = True sage_free(res_l) @@ -4270,11 +4741,15 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef long long* _hnf_modn_impl(Matrix_integer_dense self, mod_int det, Py_ssize_t nrows, Py_ssize_t ncols) except NULL: - cdef long long *res, *T_ent, **res_rows, **T_rows, *B + cdef long long *res + cdef long long *T_ent + cdef long long **res_rows + cdef long long **T_rows + cdef long long *B cdef Py_ssize_t i, j, k cdef long long R, mod, T_i_i, T_j_i, c1, c2, q, t cdef int u, v, d - cdef mpz_t m + cdef mpz_t m,tmp # allocate memory for result matrix res = sage_malloc(sizeof(long long)*ncols*nrows) @@ -4316,15 +4791,17 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse mpz_init(m) + mpz_init(tmp) # copy entries from self to temporary matrix k = 0 for i from 0 <= i < nrows: for j from 0 <= j < ncols: - mpz_mod_ui(m, self._matrix[i][j], det) + self.get_unsafe_mpz(i,j,tmp) + mpz_mod_ui(m, tmp, det) T_ent[k] = mpz_get_si(m) k += 1 mpz_clear(m) - + mpz_clear(tmp) # initialize variables i = 0 @@ -4456,12 +4933,12 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef Matrix_integer_dense other = bottom.dense_matrix() cdef Matrix_integer_dense M M = self.new_matrix(nrows = self._nrows + other._nrows, ncols = self.ncols()) - cdef Py_ssize_t i, k - k = self._nrows * self._ncols - for i from 0 <= i < k: - mpz_set(M._entries[i], self._entries[i]) - for i from 0 <= i < other._nrows * other._ncols: - mpz_set(M._entries[k + i], other._entries[i]) + cdef Py_ssize_t i, j, k + for j from 0 <= j < self._ncols: + for i from 0 <= i < self._nrows: + fmpz_set(fmpz_mat_entry(M._matrix,i,j),fmpz_mat_entry(self._matrix,i,j)) + for i from 0 <= i < other._nrows: + fmpz_set(fmpz_mat_entry(M._matrix,self._nrows + i,j),fmpz_mat_entry(other._matrix,i,j)) if subdivide: M._subdivide_on_stack(self, other) return M @@ -4544,13 +5021,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse p, qs, qa = 0, 0, 0 for i from 0 <= i < m: for j from 0 <= j < ns: - mpz_set(Z._entries[p], self._entries[qs]) - p = p + 1 - qs = qs + 1 + fmpz_set(fmpz_mat_entry(Z._matrix,i,j),fmpz_mat_entry(self._matrix,i,j)) for j from 0 <= j < na: - mpz_set(Z._entries[p], other._entries[qa]) - p = p + 1 - qa = qa + 1 + fmpz_set(fmpz_mat_entry(Z._matrix,i,j + ns),fmpz_mat_entry(other._matrix,i,j)) if subdivide: Z._subdivide_on_augment(self, other) return Z @@ -4588,25 +5061,28 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [ 1 5 -10] """ cdef Matrix_integer_dense res = self._new_uninitialized_matrix(self._nrows + 1, self._ncols) - cdef Py_ssize_t j, k + cdef Py_ssize_t j cdef Integer z + cdef fmpz_t zflint if index < 0: raise ValueError("index must be nonnegative") if index > self._nrows: raise ValueError("index must be less than number of rows") - for j from 0 <= j < self._ncols * index: - mpz_init_set(res._entries[j], self._entries[j]) + fmpz_init(zflint) - k = 0 - for j from self._ncols * index <= j < self._ncols * (index+1): - z = row[k] - mpz_init_set(res._entries[j], z.value) - k += 1 + for j from 0 <= j < self._ncols: + for i from 0 <= i < index: + fmpz_init_set(fmpz_mat_entry(res._matrix,i,j), fmpz_mat_entry(self._matrix,i,j)) + + z = row[j] + fmpz_set_mpz(zflint,z.value) + fmpz_init_set(fmpz_mat_entry(res._matrix,index,j), zflint) - for j from self._ncols * (index+1) <= j < (self._nrows + 1)*self._ncols: - mpz_init_set(res._entries[j], self._entries[j - self._ncols]) + for i from index <= i < self._nrows: + fmpz_init_set(fmpz_mat_entry(res._matrix,i+1,j), fmpz_mat_entry(self._matrix,i,j)) res._initialized = True + fmpz_clear(zflint) return res def _delete_zero_columns(self): @@ -4670,7 +5146,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse i = 0 for c in cols_ins: # The following does this quickly: A[r,c] = self[r,i] - mpz_set(A._matrix[r][c], self._matrix[r][i]) + fmpz_set(fmpz_mat_entry(A._matrix,r,c),fmpz_mat_entry(self._matrix,r,i)) i += 1 return A @@ -4694,23 +5170,27 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ self.check_mutability() - cdef mpz_t g - mpz_init(g) - cdef Py_ssize_t i, j - cdef mpz_t* row + cdef fmpz_t g + cdef fmpz_t tmp + fmpz_init(g) + fmpz_init(tmp) + cdef long i, j for i from 0 <= i < self._nrows: - mpz_set_ui(g, 0) - row = self._matrix[i] + fmpz_set_ui(g, 0) for j from 0 <= j < self._ncols: - mpz_gcd(g, g, row[j]) - if mpz_cmp_ui(g, 1) == 0: + fmpz_gcd(g, g, fmpz_mat_entry(self._matrix,i,j)) + if fmpz_cmp_ui(g, 1) == 0: break - if mpz_cmp_ui(g, 1) != 0: + if fmpz_cmp_ui(g, 1) != 0: # divide through row for j from 0 <= j < self._ncols: - mpz_divexact(row[j], row[j], g) - mpz_clear(g) + fmpz_set(tmp,fmpz_mat_entry(self._matrix,i,j)) + fmpz_divexact(tmp, tmp, g) + fmpz_set(fmpz_mat_entry(self._matrix,i,j),tmp) + fmpz_clear(g) + fmpz_clear(tmp) + return def gcd(self): """ @@ -4724,14 +5204,16 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef Integer g = Integer(0) cdef Py_ssize_t i, j - cdef mpz_t* row - + cdef fmpz_t tmpgcd + fmpz_init(tmpgcd) for i from 0 <= i < self._nrows: - row = self._matrix[i] for j from 0 <= j < self._ncols: - mpz_gcd(g.value, g.value, row[j]) - if mpz_cmp_ui(g.value, 1) == 0: + fmpz_gcd(tmpgcd,tmpgcd,fmpz_mat_entry(self._matrix,i,j)) + if fmpz_cmp_ui(tmpgcd, 1) == 0: + fmpz_get_mpz(g.value,tmpgcd) return g + fmpz_get_mpz(g.value,tmpgcd) + fmpz_clear(tmpgcd) return g def _change_ring(self, ring): @@ -4826,14 +5308,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [2 5] """ cdef Matrix_integer_dense A - A = Matrix_integer_dense.__new__(Matrix_integer_dense, - self._parent.matrix_space(self._ncols,self._nrows), - 0,False,False) - cdef Py_ssize_t i,j + A = self._new_uninitialized_matrix(self._ncols,self._nrows) sig_on() - for i from 0 <= i < self._nrows: - for j from 0 <= j < self._ncols: - mpz_init_set(A._matrix[j][i], self._matrix[i][j]) + fmpz_mat_transpose(A._matrix,self._matrix) sig_off() A._initialized = True if self._subdivisions is not None: @@ -4871,10 +5348,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse nr , nc = (self._nrows, self._ncols) cdef Matrix_integer_dense A - A = Matrix_integer_dense.__new__(Matrix_integer_dense, - self._parent.matrix_space(nc,nr), - 0,False,False) - + A = self._new_uninitialized_matrix(nc,nr) cdef Py_ssize_t i,j cdef Py_ssize_t ri,rj # reversed i and j sig_on() @@ -4884,7 +5358,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse ri = ri-1 for j from 0 <= j < nc: rj = rj-1 - mpz_init_set(A._matrix[rj][ri], self._matrix[i][j]) + fmpz_init_set(fmpz_mat_entry(A._matrix,rj,ri),fmpz_mat_entry(self._matrix,i,j)) sig_off() A._initialized = True @@ -5013,7 +5487,7 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse """ cdef GEN A pari_catch_sig_on() - A = pari._new_GEN_from_mpz_t_matrix_rotate90(self._matrix, self._nrows, self._ncols) + A = pari._new_GEN_from_fmpz_mat_t_rotate90(self._matrix, self._nrows, self._ncols) cdef GEN H = mathnf0(A, flag) B = self.extract_hnf_from_pari_matrix(H, flag, include_zero_rows) pari.clear_stack() # This calls pari_catch_sig_off() @@ -5083,6 +5557,8 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse cdef extract_hnf_from_pari_matrix(self, GEN H, int flag, bint include_zero_rows): # Throw away the transformation matrix (yes, we should later # code this to keep track of it). + cdef mpz_t tmp + mpz_init(tmp) if flag > 0: H = gel(H,1) @@ -5096,7 +5572,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse B = self.new_matrix(nrows=H_nc) for i in range(self._ncols): for j in range(H_nc): - t_INT_to_ZZ(B._matrix[j][self._ncols-i-1], gcoeff(H, i+1, H_nc-j)) + t_INT_to_ZZ(tmp, gcoeff(H, i+1, H_nc-j)) + fmpz_set_mpz(fmpz_mat_entry(B._matrix,j,self._ncols-i-1),tmp) + mpz_clear(tmp) return B cdef inline GEN pari_GEN(Matrix_integer_dense B): @@ -5108,7 +5586,7 @@ cdef inline GEN pari_GEN(Matrix_integer_dense B): For internal use only; this directly uses the PARI stack. One should call ``sig_on()`` before and ``sig_off()`` after. """ - cdef GEN A = pari._new_GEN_from_mpz_t_matrix(B._matrix, B._nrows, B._ncols) + cdef GEN A = pari._new_GEN_from_fmpz_mat_t(B._matrix, B._nrows, B._ncols) return A @@ -5117,23 +5595,23 @@ cdef inline GEN pari_GEN(Matrix_integer_dense B): cdef _clear_columns(Matrix_integer_dense A, pivots, Py_ssize_t n): # Clear all columns cdef Py_ssize_t i, k, p, l, m = A._ncols - cdef mpz_t** matrix = A._matrix - cdef mpz_t c, t + cdef fmpz_t c,t sig_on() - mpz_init(c) - mpz_init(t) + fmpz_init(c) + fmpz_init(t) for i from 0 <= i < len(pivots): p = pivots[i] for k from 0 <= k < n: if k != i: - if mpz_cmp_si(matrix[k][p],0): - mpz_fdiv_q(c, matrix[k][p], matrix[i][p]) + if fmpz_cmp_si(fmpz_mat_entry(A._matrix,k,p),0): + fmpz_fdiv_q(c, fmpz_mat_entry(A._matrix,k,p), fmpz_mat_entry(A._matrix,i,p)) # subtract off c*v from row k; resulting A[k,i] entry will be < b, hence in Echelon form. for l from 0 <= l < m: - mpz_mul(t, c, matrix[i][l]) - mpz_sub(matrix[k][l], matrix[k][l], t) - mpz_clear(c) - mpz_clear(t) + fmpz_mul(t, c, fmpz_mat_entry(A._matrix,i,l)) + fmpz_sub(fmpz_mat_entry(A._matrix,k,l), fmpz_mat_entry(A._matrix,k,l), t) + + fmpz_clear(c) + fmpz_clear(t) sig_off() ############################################################### @@ -5193,9 +5671,9 @@ cpdef _lift_crt(Matrix_integer_dense M, residues, moduli=None): [-1 0 1 0 0 1 1 0 0 0 -1 -2] """ - cdef size_t n, i, j, k - cdef Py_ssize_t nr, nc - + cdef size_t i, j, k + cdef Py_ssize_t nr, n + cdef mpz_t *tmp = sage_malloc(sizeof(mpz_t) * M._ncols) n = len(residues) if n == 0: # special case: obviously residues[0] wouldn't make sense here. return M @@ -5226,191 +5704,22 @@ cpdef _lift_crt(Matrix_integer_dense M, residues, moduli=None): if row_list[k] == NULL: raise MemoryError("out of memory allocating multi-modular coefficient list") + for j in range(M._ncols): + mpz_init(tmp[j]) + for i in range(nr): for k in range(n): (residues[k])._copy_row_to_mod_int_array(row_list[k],i) - mm.mpz_crt_vec(M._matrix[i], row_list, nc) + mm.mpz_crt_vec(tmp, row_list, nc) + for j in range(nc): + M.set_unsafe_mpz(i,j,tmp[j]) for k in range(n): sage_free(row_list[k]) + for j in range(M._ncols): + mpz_clear(tmp[j]) sage_free(row_list) - + sage_free(tmp) sig_off() return M -####################################################### - -# Conclusions: -# On OS X Intel, at least: -# - if log_2(height) >= 0.70 * nrows, use classical - -def tune_multiplication(k, nmin=10, nmax=200, bitmin=2,bitmax=64): - """ - Compare various multiplication algorithms. - - INPUT: - - - ``k`` - integer; affects numbers of trials - - - ``nmin`` - integer; smallest matrix to use - - - ``nmax`` - integer; largest matrix to use - - - ``bitmin`` - integer; smallest bitsize - - - ``bitmax`` - integer; largest bitsize - - OUTPUT: - - - ``prints what doing then who wins`` - multimodular - or classical - - EXAMPLES:: - - sage: from sage.matrix.matrix_integer_dense import tune_multiplication - sage: tune_multiplication(2, nmin=10, nmax=60, bitmin=2,bitmax=8) - 10 2 0.2 - ... - """ - from constructor import random_matrix - from sage.rings.integer_ring import ZZ - for n in range(nmin,nmax,10): - for i in range(bitmin, bitmax, 4): - A = random_matrix(ZZ, n, n, x = 2**i) - B = random_matrix(ZZ, n, n, x = 2**i) - t0 = cputime() - for j in range(k//n + 1): - C = A._multiply_classical(B) - t0 = cputime(t0) - t1 = cputime() - for j in range(k//n+1): - C = A._multiply_multi_modular(B) - t1 = cputime(t1) - print n, i, float(i)/float(n) - if t0 < t1: - print 'classical' - else: - print 'multimod' - - -############################################################## -# Some people really really really want to make sure their -# 4x4 determinant is really really really fast. -############################################################## - -cdef int four_dim_det(mpz_t r,mpz_t *x) except -1: - """ - Internal function used in computing determinants of 4x4 matrices. - - TESTS:: - - sage: A = matrix(ZZ,4,[1,0,3,0,4,3,2,1,0,5,0,0,9,1,2,3]) - sage: A.determinant() # indirect doctest - 25 - """ - cdef mpz_t a,b - - sig_on() - mpz_init(a) - mpz_init(b) - - mpz_mul(a,x[3], x[6] ); mpz_submul(a,x[2], x[7] ) - mpz_mul(b,x[9], x[12]); mpz_submul(b,x[8], x[13]) - mpz_mul(r,a,b) - mpz_mul(a,x[1], x[7] ); mpz_submul(a,x[3], x[5] ) - mpz_mul(b,x[10],x[12]); mpz_submul(b,x[8], x[14]) - mpz_addmul(r,a,b) - mpz_mul(a,x[2], x[5] ); mpz_submul(a,x[1], x[6] ) - mpz_mul(b,x[11],x[12]); mpz_submul(b,x[8], x[15]) - mpz_addmul(r,a,b) - mpz_mul(a,x[3], x[4] ); mpz_submul(a,x[0], x[7] ) - mpz_mul(b,x[10],x[13]); mpz_submul(b,x[9], x[14]) - mpz_addmul(r,a,b) - mpz_mul(a,x[0], x[6] ); mpz_submul(a,x[2], x[4] ) - mpz_mul(b,x[11],x[13]); mpz_submul(b,x[9], x[15]) - mpz_addmul(r,a,b) - mpz_mul(a,x[1], x[4] ); mpz_submul(a,x[0], x[5] ) - mpz_mul(b,x[11],x[14]); mpz_submul(b,x[10],x[15]) - mpz_addmul(r,a,b) - - mpz_clear(a) - mpz_clear(b) - sig_off() - return 0 - -#The above was generated by the following Sage code: -#def idx(x): -# return [i for i in range(16) if x[i]] -# -#def Det4Maker(): -# Q = PolynomialRing(QQ,['x%s'%ZZ(i).str(16) for i in range(16)]) -# M = Matrix(Q,4,4,Q.gens()) -# D = M.det() -# Dm = [(m.exponents()[0],D[m]) for m in D.monomials()] -# Dm.sort(reverse=True) -# MM = [[Dm[i],Dm[i+1]] for i in range(0,len(Dm),2)] -# -# S = [] -# for j,((t,s),(r,o)) in enumerate(MM): -# a,b,c,d = [ZZ(i).str(16) for i in range(16) if t[i]] -# e,f,g,h = [ZZ(i).str(16) for i in range(16) if r[i]] -# S.append((c+d+g+h,j)) -# S.sort(lambda x,y: cmp(x[0],y[0])) -# MM = [MM[s[1]] for s in S] -# MMM = [[MM[i],MM[i+1]] for i in range(0,len(MM),2)] -# -# N = 0 -# Cstrux = "" -# Pstrux = "" -# for ((t1,s1),(t2,s2)),((t3,s3),(t4,s4)) in MMM: -# a,b,c,d = idx(t1) -# e,f = idx(t2)[2:] -# h,i = idx(t3)[:2] -# -# if s1 == 1: -# a,b,h,i = h,i,a,b -# -# Pstrux+= 'a = x[%s]*x[%s]; '%(a,b) -# Cstrux+= 'mpz_mul(a,x[%s],x[%s]); '%(a,b) -# -# Pstrux+= 'a-=x[%s]*x[%s]; '%(h,i) -# Cstrux+= 'mpz_submul(a,x[%s],x[%s])\n'%(h,i) -# -# if s4*s1 == 1: -# c,d,e,f = e,f,c,d -# -# Pstrux+= 'b = x[%s]*x[%s]; '%(c,d) -# Cstrux+= 'mpz_mul(b,x[%s],x[%s]); '%(c,d) -# -# Pstrux+= 'b-=x[%s]*x[%s]; '%(e,f) -# Cstrux+= 'mpz_submul(b,x[%s],x[%s])\n'%(e,f) -# -# -# if N: -# Pstrux+= 'r+= a*b\n' -# Cstrux+= 'mpz_addmul(r,a,b)\n' -# else: -# Pstrux+= 'r = a*b\n' -# Cstrux+= 'mpz_mul(r,a,b)\n' -# N = 1 -# return Cstrux,Pstrux -# -#print Det4Maker()[0] - -#Note, one can prove correctness of the above by setting -# sage: Q = PolynomialRing(QQ,['x%s'%ZZ(i).str(16) for i in range(16)]) -# sage: M = Matrix(Q,4,4,Q.gens()) -# sage: D = M.det() -# sage: x = Q.gens() -#and then evaluating the contents of Det4Maker()[1]... -# sage: a = x[3]*x[6]; a-=x[2]*x[7]; b = x[9]*x[12]; b-=x[8]*x[13]; r = a*b -# sage: a = x[1]*x[7]; a-=x[3]*x[5]; b = x[10]*x[12]; b-=x[8]*x[14]; r+= a*b -# sage: a = x[2]*x[5]; a-=x[1]*x[6]; b = x[11]*x[12]; b-=x[8]*x[15]; r+= a*b -# sage: a = x[3]*x[4]; a-=x[0]*x[7]; b = x[10]*x[13]; b-=x[9]*x[14]; r+= a*b -# sage: a = x[0]*x[6]; a-=x[2]*x[4]; b = x[11]*x[13]; b-=x[9]*x[15]; r+= a*b -# sage: a = x[1]*x[4]; a-=x[0]*x[5]; b = x[11]*x[14]; b-=x[10]*x[15]; r+= a*b -#and comparing results: -# sage: r-D -# 0 -#bingo! - diff --git a/src/sage/matrix/matrix_integer_dense_hnf.py b/src/sage/matrix/matrix_integer_dense_hnf.py index 8c589739ff3..f0ca7db39d3 100644 --- a/src/sage/matrix/matrix_integer_dense_hnf.py +++ b/src/sage/matrix/matrix_integer_dense_hnf.py @@ -851,7 +851,7 @@ def probable_hnf(A, include_zero_rows, proof): H = hnf_square(C, proof=proof) except NotImplementedError: # raise - # this signals that we must fallback to pari + # this signals that we must fallback to PARI verbose("generic random modular HNF algorithm failed -- we fall back to PARI") H = A.hermite_form(algorithm='pari', include_zero_rows=include_zero_rows, proof=proof) return H, H.pivots() @@ -1180,7 +1180,7 @@ def __do_check(v): print 'a = matrix(ZZ, %s, %s, %s)'%(a.nrows(), a.ncols(), a.list()) return else: - if hnf(a)[0] != a.echelon_form('pari'): + if hnf(a)[0] != a.echelon_form(algorithm = 'pari'): print "bug computing hnf of a matrix" print 'a = matrix(ZZ, %s, %s, %s)'%(a.nrows(), a.ncols(), a.list()) return diff --git a/src/sage/matrix/matrix_integer_dense_saturation.py b/src/sage/matrix/matrix_integer_dense_saturation.py index 34db058c223..0c3c4af8f00 100644 --- a/src/sage/matrix/matrix_integer_dense_saturation.py +++ b/src/sage/matrix/matrix_integer_dense_saturation.py @@ -134,7 +134,7 @@ def solve_system_with_difficult_last_row(B, A): else: break D = B.matrix_from_rows(range(C.nrows()-1)) - N = D._rational_kernel_iml() + N = D._rational_kernel_flint() if N.ncols() != 1: verbose("Difficult solve quickly failed. Using direct approach.") return B.solve_right(A) diff --git a/src/sage/matrix/matrix_integer_sparse.pyx b/src/sage/matrix/matrix_integer_sparse.pyx index b9d64673382..7aca71a4c91 100644 --- a/src/sage/matrix/matrix_integer_sparse.pyx +++ b/src/sage/matrix/matrix_integer_sparse.pyx @@ -409,7 +409,7 @@ cdef class Matrix_integer_sparse(matrix_sparse.Matrix_sparse): OUTPUT: Returns a pair. First item is the string is either - 'computed-pari-int' or 'computed-iml-int', which identifies + 'computed-pari-int', 'computed-iml-int' or 'computed-flint-int', which identifies the nature of the basis vectors. Second item is a matrix whose rows are a basis for the right kernel, @@ -438,22 +438,23 @@ cdef class Matrix_integer_sparse(matrix_sparse.Matrix_sparse): sage: X = result[1]; X [-469 214 -30 119 -37 0] [ 370 -165 18 -91 30 -2] + sage: A*X.transpose() == zero_matrix(ZZ, 4, 2) True sage: result = A._right_kernel_matrix(algorithm='default') sage: result[0] - 'computed-pari-int' + 'computed-flint-int' sage: result[1] - [-26 31 -30 21 2 -10] - [-47 -13 48 -14 -11 18] + [ 469 -214 30 -119 37 0] + [-370 165 -18 91 -30 2] sage: result = A._right_kernel_matrix() sage: result[0] - 'computed-pari-int' + 'computed-flint-int' sage: result[1] - [-26 31 -30 21 2 -10] - [-47 -13 48 -14 -11 18] + [ 469 -214 30 -119 37 0] + [-370 165 -18 91 -30 2] With the 'default' given as the algorithm, several heuristics are used to determine if PARI or IML ('padic') is used. The code has @@ -481,8 +482,7 @@ cdef class Matrix_integer_sparse(matrix_sparse.Matrix_sparse): sage: A = matrix(ZZ, 0, 2, sparse=True) sage: A._right_kernel_matrix()[1] - [1 0] - [0 1] + [] sage: A = matrix(ZZ, 2, 0, sparse=True) sage: A._right_kernel_matrix()[1].parent() Full MatrixSpace of 0 by 0 dense matrices over Integer Ring diff --git a/src/sage/matrix/matrix_modn_dense.pxd b/src/sage/matrix/matrix_modn_dense.pxd index 1541f1f44d1..7733ad541db 100644 --- a/src/sage/matrix/matrix_modn_dense.pxd +++ b/src/sage/matrix/matrix_modn_dense.pxd @@ -1,6 +1,13 @@ -from sage.ext.mod_int cimport * +include "sage/ext/cdefs.pxi" +include "sage/libs/ntl/decl.pxi" +include "sage/libs/flint/fmpz.pxi" +include "sage/libs/flint/fmpz_mat.pxi" cimport matrix_dense +cimport matrix_integer_dense +cimport sage.rings.integer +from sage.rings.integer cimport Integer +from sage.ext.mod_int cimport mod_int, MOD_INT_OVERFLOW cdef class Matrix_modn_dense(matrix_dense.Matrix_dense): cdef mod_int **_matrix @@ -8,9 +15,6 @@ cdef class Matrix_modn_dense(matrix_dense.Matrix_dense): cdef mod_int p cdef mod_int gather cdef xgcd_eliminate (self, mod_int * row1, mod_int* row2, Py_ssize_t start_col) - #cdef set_matrix(Matrix_modn_dense self, mod_int **m) - #cdef mod_int **get_matrix(Matrix_modn_dense self) - #cdef mod_int entry(self, mod_int i, mod_int j) cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value) cdef _rescale_row_c(self, Py_ssize_t row, mod_int multiple, Py_ssize_t start_col) cdef _rescale_col_c(self, Py_ssize_t col, mod_int multiple, Py_ssize_t start_row) @@ -21,5 +25,4 @@ cdef class Matrix_modn_dense(matrix_dense.Matrix_dense): cpdef _export_as_string(self) - cpdef is_Matrix_modn_dense(self) diff --git a/src/sage/matrix/matrix_modn_dense.pyx b/src/sage/matrix/matrix_modn_dense.pyx index cfe9cee5a7f..553a1b0a78d 100644 --- a/src/sage/matrix/matrix_modn_dense.pyx +++ b/src/sage/matrix/matrix_modn_dense.pyx @@ -1717,16 +1717,12 @@ cdef class Matrix_modn_dense(matrix_dense.Matrix_dense): cdef Py_ssize_t i, j cdef Matrix_integer_dense L - L = Matrix_integer_dense.__new__(Matrix_integer_dense, - self.parent().change_ring(ZZ), - 0, 0, 0) - cdef mpz_t* L_row + L = Matrix_integer_dense._new_uninitialized_matrix(Matrix_integer_dense,self._nrows,self._ncols) cdef mod_int* A_row for i from 0 <= i < self._nrows: - L_row = L._matrix[i] A_row = self._matrix[i] for j from 0 <= j < self._ncols: - mpz_init_set_si(L_row[j], A_row[j]) + fmpz_set_si(fmpz_mat_entry(L._matrix,i,j),A_row[j]) L._initialized = 1 L.subdivide(self.subdivisions()) return L diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index ef9b4e08ae2..38913805758 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -115,6 +115,7 @@ from matrix_integer_dense cimport Matrix_integer_dense from sage.rings.integer_ring import ZZ from sage.structure.proof.proof import get_flag as get_proof_flag from sage.misc.randstate cimport randstate, current_randstate +import sage.matrix.matrix_space as matrix_space cdef long num = 1 cdef bint little_endian = ((&num))[0] @@ -2898,18 +2899,20 @@ cdef class Matrix_modn_dense_template(matrix_dense.Matrix_dense): cdef Py_ssize_t i, j cdef Matrix_integer_dense L - L = Matrix_integer_dense.__new__(Matrix_integer_dense, - self.parent().change_ring(ZZ), - 0, 0, 0) + cdef object P = matrix_space.MatrixSpace(ZZ, self._nrows, self._ncols, sparse=False) + L = Matrix_integer_dense(P,ZZ(0),False,False) cdef mpz_t* L_row + cdef mpz_t tmp cdef celement* A_row + mpz_init(tmp) for i from 0 <= i < self._nrows: - L_row = L._matrix[i] A_row = self._matrix[i] for j from 0 <= j < self._ncols: - mpz_init_set_d(L_row[j], A_row[j]) + mpz_set_si(tmp,(A_row[j])) + L.set_unsafe_mpz(i,j,tmp) L._initialized = 1 L.subdivide(self.subdivisions()) + mpz_clear(tmp) return L diff --git a/src/sage/matrix/matrix_rational_dense.pxd b/src/sage/matrix/matrix_rational_dense.pxd index 3aa9ce62477..ce3fae47415 100644 --- a/src/sage/matrix/matrix_rational_dense.pxd +++ b/src/sage/matrix/matrix_rational_dense.pxd @@ -1,6 +1,14 @@ include "sage/ext/cdefs.pxi" +include "sage/libs/ntl/decl.pxi" +include "sage/libs/flint/fmpz.pxi" +include "sage/libs/flint/fmpz_poly.pxi" +include "sage/libs/flint/fmpz_mat.pxi" cimport matrix_dense +cimport matrix_integer_dense +cimport sage.rings.integer +from sage.rings.integer cimport Integer +from sage.ext.mod_int cimport * cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): diff --git a/src/sage/matrix/matrix_rational_dense.pyx b/src/sage/matrix/matrix_rational_dense.pyx index 781db98e6f8..f48bb86d2c8 100644 --- a/src/sage/matrix/matrix_rational_dense.pyx +++ b/src/sage/matrix/matrix_rational_dense.pyx @@ -641,6 +641,7 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): - "iml" -- use an iml-based algorithm + - "flint" -- use FLINT OUTPUT: the inverse of self @@ -736,6 +737,10 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): AZ, denom = self._clear_denom() B, d = AZ._invert_iml(check_invertible=check_invertible) return (denom/d)*B + elif algorithm == "flint": + AZ, denom = self._clear_denom() + B, d = AZ._invert_flint(check_invertible=check_invertible) + return (denom/d)*B else: raise ValueError("unknown algorithm '%s'"%algorithm) @@ -860,25 +865,31 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): cdef Py_ssize_t i, j cdef Matrix_integer_dense A cdef mpq_t *self_row - cdef mpz_t *A_row + cdef mpz_t tmp + cdef fmpz_t aij + fmpz_init(aij) + mpz_init(tmp) D = PY_NEW(Integer) self.mpz_denom(D.value) from sage.matrix.matrix_space import MatrixSpace MZ = MatrixSpace(ZZ, self._nrows, self._ncols, sparse=self.is_sparse()) - A = Matrix_integer_dense.__new__(Matrix_integer_dense, MZ, 0, 0, 0) + A = Matrix_integer_dense.__new__(Matrix_integer_dense, MZ, None, None, None) + fmpz_mat_init(A._matrix,self._nrows,self._ncols) sig_on() for i from 0 <= i < self._nrows: - A_row = A._matrix[i] self_row = self._matrix[i] for j from 0 <= j < self._ncols: - mpz_init(A_row[0]) - mpz_divexact(A_row[0], D.value, mpq_denref(self_row[0])) - mpz_mul(A_row[0], A_row[0], mpq_numref(self_row[0])) - A_row += 1 + fmpz_init(fmpz_mat_entry(A._matrix,i,j)) + mpz_divexact(tmp, D.value, mpq_denref(self_row[0])) + mpz_mul(tmp, tmp, mpq_numref(self_row[0])) + fmpz_set_mpz(fmpz_mat_entry(A._matrix,i,j),tmp) self_row += 1 sig_off() A._initialized = 1 self.cache('clear_denom', (A,D)) + fmpz_clear(aij) + mpz_clear(tmp) + return A, D def charpoly(self, var='x', algorithm='linbox'): @@ -1064,15 +1075,12 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): cdef Matrix_integer_dense A, B, AB cdef Matrix_rational_dense res cdef Integer D - cdef mpz_t* AB_row, cdef mpq_t* res_row sig_on() A, A_denom = self._clear_denom() B, B_denom = right._clear_denom() - if algorithm == 'default': + if algorithm == 'default' or algorithm == 'multimodular': AB = A*B - elif algorithm == 'multimodular': - AB = A._multiply_multi_modular(B) else: sig_off() raise ValueError("unknown algorithm '%s'"%algorithm) @@ -1086,13 +1094,11 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): else: res = Matrix_rational_dense.__new__(Matrix_rational_dense, self.matrix_space(AB._nrows, AB._ncols), 0, 0, 0) for i from 0 <= i < res._nrows: - AB_row = AB._matrix[i] res_row = res._matrix[i] for j from 0 <= j < res._ncols: - mpz_set(mpq_numref(res_row[0]), AB_row[0]) + fmpz_get_mpz(mpq_numref(res_row[0]), fmpz_mat_entry(AB._matrix,i,j)) mpz_set(mpq_denref(res_row[0]), D.value) mpq_canonicalize(res_row[0]) - AB_row = AB_row + 1 res_row = res_row + 1 sig_off() return res @@ -1269,7 +1275,7 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): [0 0 1] """ tm = verbose("computing right kernel matrix over the rationals for %sx%s matrix" % (self.nrows(), self.ncols()),level=1) - # _rational_kernel_iml() gets the zero-row case wrong, fix it there + # _rational_kernel_flint() gets the zero-row case wrong, fix it there if self.nrows()==0: import constructor K = constructor.identity_matrix(QQ, self.ncols()) @@ -1533,7 +1539,6 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): cdef Matrix_rational_dense E cdef Integer d cdef mpq_t* E_row - cdef mpz_t* X_row t = verbose('Computing echelon form of %s x %s matrix over QQ using p-adic nullspace algorithm.'%( self.nrows(), self.ncols())) @@ -1554,9 +1559,8 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): # Fill in the non-pivot part of the matrix for i from 0 <= i < X.nrows(): E_row = E._matrix[i] - X_row = X._matrix[i] for j from 0 <= j < X.ncols(): - mpz_set(mpq_numref(E_row[nonpivots[j]]), X_row[j]) + fmpz_get_mpz(mpq_numref(E_row[nonpivots[j]]), fmpz_mat_entry(X._matrix,i,j)) mpz_set(mpq_denref(E_row[nonpivots[j]]), d.value) mpq_canonicalize(E_row[nonpivots[j]]) @@ -1595,9 +1599,8 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): # Fill in the non-pivot part of self. for i from 0 <= i < X.nrows(): E_row = self._matrix[i] - X_row = X._matrix[i] for j from 0 <= j < X.ncols(): - mpz_set(mpq_numref(E_row[nonpivots[j]]), X_row[j]) + fmpz_get_mpz(mpq_numref(E_row[nonpivots[j]]), fmpz_mat_entry(X._matrix,i,j)) mpz_set(mpq_denref(E_row[nonpivots[j]]), d.value) mpq_canonicalize(E_row[nonpivots[j]]) @@ -2088,91 +2091,20 @@ cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): cdef Py_ssize_t i, j, nr, nc cdef mpz_t* Z_row cdef mpq_t* Q_row - + cdef mpz_t tmp + mpz_init(tmp) ZA = _lift_crt(res, mm) nr = ZA._nrows nc = ZA._ncols QA = Matrix_rational_dense.__new__(Matrix_rational_dense, self.parent(), None, None, None) m = mm.prod() for i from 0 <= i < nr: - Z_row = ZA._matrix[i] + #Z_row = ZA._matrix[i] Q_row = QA._matrix[i] for j from 0 <= j < nc: - mpq_rational_reconstruction(Q_row[j], Z_row[j], m.value) - return QA - - def _lift_crt_rr_with_lcm(self, res, mm): - """ - Optimizations: When doing the rational_recon lift of a (mod m) - first see if |a| < sqrt(m/2) in which case it lifts to an integer - (often a=0 or 1). - - If that fails, keep track of the lcm d of denominators found so - far, and check to see if z = a\*d lifts to an integer with |z| <= - sqrt(m/2). If so, no need to do rational recon. This should be the - case for most a after a while, and should saves substantial time! - """ - cdef Integer m - cdef Matrix_integer_dense ZA - cdef Matrix_rational_dense QA - cdef Py_ssize_t i, j, nr, nc - cdef mpz_t* Z_row - cdef mpq_t* Q_row - cdef mpz_t lcm_denom, sqrt_m, neg_sqrt_m, z - - mpz_init(z) - mpz_init(sqrt_m) - mpz_init(neg_sqrt_m) - mpz_init_set_ui(lcm_denom, 1) - - m = mm.prod() - mpz_fdiv_q_2exp(sqrt_m, m.value, 1) - mpz_sqrt(sqrt_m, sqrt_m) - mpz_sub(neg_sqrt_m, m.value, sqrt_m) - - t = verbose("Starting crt", level=2) - ZA = _lift_crt(res, mm) - t = verbose("crt finished", t, level=2) - nr = ZA._nrows - nc = ZA._ncols - QA = Matrix_rational_dense.__new__(Matrix_rational_dense, self.parent(), None, None, None) - - cdef bint is_integral, lcm_trick - is_integral = 0 - lcm_trick = 0 - - t = verbose("Starting rational reconstruction", level=2) - for i from 0 <= i < nr: - Z_row = ZA._matrix[i] - Q_row = QA._matrix[i] - for j from 0 <= j < nc: - if mpz_cmp(Z_row[j], sqrt_m) < 0: - mpz_set(mpq_numref(Q_row[j]), Z_row[j]) - is_integral += 1 - elif mpz_cmp(Z_row[j], neg_sqrt_m) > 0: - mpz_sub(mpq_numref(Q_row[j]), Z_row[j], m.value) - is_integral += 1 - else: - mpz_mul(z, Z_row[j], lcm_denom) - mpz_fdiv_r(z, z, m.value) - if mpz_cmp(z, sqrt_m) < 0: - mpz_set(mpq_numref(Q_row[j]), z) - mpz_set(mpq_denref(Q_row[j]), lcm_denom) - mpq_canonicalize(Q_row[j]) - lcm_trick += 1 - elif mpz_cmp(z, neg_sqrt_m) > 0: - mpz_sub(mpq_numref(Q_row[j]), z, m.value) - mpz_set(mpq_denref(Q_row[j]), lcm_denom) - mpq_canonicalize(Q_row[j]) - lcm_trick += 1 - else: - mpq_rational_reconstruction(Q_row[j], Z_row[j], m.value) - mpz_lcm(lcm_denom, lcm_denom, mpq_denref(Q_row[j])) - mpz_clear(z) - mpz_clear(sqrt_m) - mpz_clear(neg_sqrt_m) - mpz_clear(lcm_denom) - t = verbose("rr finished. integral entries: %s, lcm trick: %s, other: %s"%(is_integral, lcm_trick, nr*nc - is_integral - lcm_trick), t, level=2) + fmpz_get_mpz(tmp,fmpz_mat_entry(ZA._matrix,i,j)) + mpq_rational_reconstruction(Q_row[j], tmp, m.value) + mpz_clear(tmp) return QA def randomize(self, density=1, num_bound=2, den_bound=2, \ diff --git a/src/sage/matrix/matrix_space.py b/src/sage/matrix/matrix_space.py index 2e93fccee6c..9f20bf99d79 100644 --- a/src/sage/matrix/matrix_space.py +++ b/src/sage/matrix/matrix_space.py @@ -143,7 +143,7 @@ class MatrixSpace(UniqueRepresentation, parent_gens.ParentWithGens): _no_generic_basering_coercion = True @staticmethod - def __classcall__(cls, base_ring, nrows, ncols=None, sparse=False): + def __classcall__(cls, base_ring, nrows, ncols=None, sparse=False, implementation = 'flint'): """ Create with the command @@ -160,7 +160,7 @@ def __classcall__(cls, base_ring, nrows, ncols=None, sparse=False): columns - ``sparse`` - (default false) whether or not matrices are given a sparse representation - + - ``implementation`` - (default 'flint') use 'sage' or 'flint' OUTPUT: @@ -216,12 +216,13 @@ def __classcall__(cls, base_ring, nrows, ncols=None, sparse=False): raise TypeError("base_ring (=%s) must be a ring"%base_ring) if ncols is None: ncols = nrows nrows = int(nrows); ncols = int(ncols); sparse=bool(sparse) - return super(MatrixSpace, cls).__classcall__(cls, base_ring, nrows, ncols, sparse) + return super(MatrixSpace, cls).__classcall__(cls, base_ring, nrows, ncols, sparse,implementation) def __init__(self, base_ring, nrows, ncols=None, - sparse=False): + sparse=False, + implementation = 'flint'): """ TEST: @@ -258,6 +259,8 @@ def __init__(self, base_ring, sage: B = MatrixSpace(RDF,1000,1000).random_element() sage: C = A * B """ + self._implementation = implementation + if ncols is None: ncols = nrows from sage.categories.all import Modules, Algebras parent_gens.ParentWithGens.__init__(self, base_ring) # category = Modules(base_ring) @@ -293,6 +296,7 @@ def __init__(self, base_ring, category = Algebras(base_ring.category()) else: category = Modules(base_ring.category()) + sage.structure.parent.Parent.__init__(self, category=category) #sage.structure.category_object.CategoryObject._init_category_(self, category) @@ -371,7 +375,7 @@ def _copy_zero(self): else: return True - def __call__(self, entries=None, coerce=True, copy=True): + def __call__(self, entries=None, coerce=True, copy=True, sparse = False): """ EXAMPLES:: @@ -930,6 +934,7 @@ def _get_matrix_class(self): R = self.base_ring() if self.is_dense(): if sage.rings.integer_ring.is_IntegerRing(R): + assert self._implementation == 'flint' return matrix_integer_dense.Matrix_integer_dense elif sage.rings.rational_field.is_RationalField(R): return matrix_rational_dense.Matrix_rational_dense diff --git a/src/sage/matrix/misc.pyx b/src/sage/matrix/misc.pyx index d402ca5b563..54720b54001 100644 --- a/src/sage/matrix/misc.pyx +++ b/src/sage/matrix/misc.pyx @@ -77,7 +77,7 @@ def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer R = Matrix_rational_dense.__new__(Matrix_rational_dense, A.parent().change_ring(QQ), 0,0,0) - cdef mpz_t a, bnd, other_bnd, one, denom + cdef mpz_t a, bnd, other_bnd, one, denom, tmp cdef Integer _bnd cdef Py_ssize_t i, j cdef int do_it @@ -87,6 +87,7 @@ def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer try: mpz_init_set_si(denom, 1) mpz_init(a) + mpz_init(tmp) mpz_init_set_si(one, 1) mpz_init(other_bnd) @@ -96,7 +97,7 @@ def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer for i from 0 <= i < A._nrows: for j from 0 <= j < A._ncols: - mpz_set(a, A._matrix[i][j]) + A.get_unsafe_mpz(i,j,a) if mpz_cmp(denom, one) != 0: mpz_mul(a, a, denom) mpz_fdiv_r(a, a, N.value) @@ -115,11 +116,13 @@ def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer mpz_set_si(mpq_denref(R._matrix[i][j]), 1) else: # Otherwise have to do it the hard way - mpq_rational_reconstruction(R._matrix[i][j], A._matrix[i][j], N.value) + A.get_unsafe_mpz(i,j,tmp) + mpq_rational_reconstruction(R._matrix[i][j], tmp, N.value) mpz_lcm(denom, denom, mpq_denref(R._matrix[i][j])) mpz_clear(denom) mpz_clear(a) + mpz_clear(tmp) mpz_clear(one) mpz_clear(other_bnd) mpz_clear(bnd) diff --git a/src/sage/rings/padics/padic_ZZ_pX_FM_element.pyx b/src/sage/rings/padics/padic_ZZ_pX_FM_element.pyx index 13cb60d8d99..2bf108e6a06 100644 --- a/src/sage/rings/padics/padic_ZZ_pX_FM_element.pyx +++ b/src/sage/rings/padics/padic_ZZ_pX_FM_element.pyx @@ -310,7 +310,7 @@ cdef class pAdicZZpXFMElement(pAdicZZpXElement): mpz_init(tmp_m) mpz_invert(tmp_m, mpq_denref(x), self.prime_pow.pow_mpz_t_top()) mpz_mul(tmp_m, tmp_m, mpq_numref(x)) - mpz_mod(tmp_m, tmp_m, self.prime_pow.pow_mpz_t_top()) + mpz_mod(tmp_m, tmp_m, &self.prime_pow.pow_mpz_t_top()[0]) mpz_to_ZZ(&tmp_z, tmp_m) ZZ_pX_SetCoeff(self.value, 0, ZZ_to_ZZ_p(tmp_z)) mpz_clear(tmp_m) diff --git a/src/sage/rings/padics/pow_computer.pyx b/src/sage/rings/padics/pow_computer.pyx index 329c5fb93a3..8a37861ad8d 100644 --- a/src/sage/rings/padics/pow_computer.pyx +++ b/src/sage/rings/padics/pow_computer.pyx @@ -489,7 +489,7 @@ cdef class PowComputer_base(PowComputer_class): sage: PC._pow_mpz_t_top_test() #indirect doctest 59049 """ - return self.top_power + return &(self.top_power[0]) cdef mpz_srcptr pow_mpz_t_tmp(self, long n): """ @@ -502,11 +502,11 @@ cdef class PowComputer_base(PowComputer_class): 81 """ if n <= self.cache_limit: - return self.small_powers[n] + return &(self.small_powers[n][0]) if n == self.prec_cap: - return self.top_power + return &(self.top_power[0]) mpz_pow_ui(self.temp_m, self.prime.value, n) - return self.temp_m + return &(self.temp_m[0]) pow_comp_cache = {} cdef PowComputer_base PowComputer_c(Integer m, Integer cache_limit, Integer prec_cap, in_field): diff --git a/src/sage/rings/padics/pow_computer_ext.pyx b/src/sage/rings/padics/pow_computer_ext.pyx index ae4983cd39a..7cb55f30eb0 100644 --- a/src/sage/rings/padics/pow_computer_ext.pyx +++ b/src/sage/rings/padics/pow_computer_ext.pyx @@ -618,7 +618,22 @@ cdef class PowComputer_ext(PowComputer_class): ZZ_to_mpz(self.temp_m, &self.top_power) else: mpz_pow_ui(self.temp_m, self.prime.value, n) - return self.temp_m + return address_of_mpz(self.temp_m) + + #def _pow_mpz_t_tmp_test(self, n): + # """ + # Test for the pow_mpz_t_tmp function. See that function's documentation for important warnings. + # + # EXAMPLES: + # sage: PC = PowComputer_ext_maker(5, 10, 10, 20, False, ntl.ZZ_pX([-5, 0, 1], 5^10), 'small', 'e',ntl.ZZ_pX([1],5^10)) + # sage: PC._pow_mpz_t_tmp_test(4) #indirect doctest + # 625 + # """ + # cdef Integer _n = Integer(n) + # if _n < 0: raise ValueError + # cdef Integer ans = PY_NEW(Integer) + # mpz_set(ans.value, self.pow_mpz_t_tmp(mpz_get_si(_n.value))[0]) + # return ans cdef ZZ_c* pow_ZZ_tmp(self, long n): """ @@ -716,7 +731,20 @@ cdef class PowComputer_ext(PowComputer_class): 15625 """ ZZ_to_mpz(self.temp_m, &self.top_power) - return self.temp_m + return address_of_mpz(self.temp_m) + + #def _pow_mpz_t_top_test(self): + # """ + # Tests the pow_mpz_t_top function + # + # EXAMPLES: + # sage: PC = PowComputer_ext_maker(5, 6, 6, 12, False, ntl.ZZ_pX([-5,0,1],5^6),'small', 'e',ntl.ZZ_pX([1],5^6)) + # sage: PC._pow_mpz_t_top_test() + # 15625 + # """ + # cdef Integer ans = PY_NEW(Integer) + # mpz_set(ans.value, self.pow_mpz_t_top()[0]) + # return ans cdef ZZ_c* pow_ZZ_top(self): """ diff --git a/src/sage/rings/polynomial/polynomial_integer_dense_flint.pxd b/src/sage/rings/polynomial/polynomial_integer_dense_flint.pxd index 13fea9c91bb..db78bd22ad3 100644 --- a/src/sage/rings/polynomial/polynomial_integer_dense_flint.pxd +++ b/src/sage/rings/polynomial/polynomial_integer_dense_flint.pxd @@ -2,8 +2,6 @@ include "sage/ext/cdefs.pxi" include "sage/libs/flint/fmpz.pxi" include "sage/libs/flint/fmpz_poly.pxi" -from sage.libs.flint.ntl_interface cimport * - from sage.rings.polynomial.polynomial_element cimport Polynomial from sage.rings.integer cimport Integer from sage.structure.parent cimport Parent diff --git a/src/sage/rings/polynomial/polynomial_integer_dense_flint.pyx b/src/sage/rings/polynomial/polynomial_integer_dense_flint.pyx index 5266b735876..e917f5738b6 100644 --- a/src/sage/rings/polynomial/polynomial_integer_dense_flint.pyx +++ b/src/sage/rings/polynomial/polynomial_integer_dense_flint.pyx @@ -43,7 +43,7 @@ from sage.rings.arith import lcm from sage.libs.flint.fmpz_poly cimport fmpz_poly_reverse, fmpz_poly_revert_series -from sage.libs.flint.ntl_interface cimport fmpz_poly_set_ZZX, fmpz_poly_get_ZZX +from sage.libs.flint.ntl_interface cimport fmpz_set_ZZ, fmpz_poly_set_ZZX, fmpz_poly_get_ZZX from sage.libs.ntl.ntl_ZZX_decl cimport *, vec_pair_ZZX_long_c cdef extern from "limits.h": diff --git a/src/sage/rings/polynomial/real_roots.pyx b/src/sage/rings/polynomial/real_roots.pyx index 48d67cdb902..edf1c14c1a8 100644 --- a/src/sage/rings/polynomial/real_roots.pyx +++ b/src/sage/rings/polynomial/real_roots.pyx @@ -4521,11 +4521,12 @@ def dprod_imatrow_vec(Matrix_integer_dense m, Vector_integer_dense v, int k): cdef int vsize = len(v) cdef int ra cdef int a - + mpz_init(tmp) for a from 0 <= a < msize: ra = subsample_vec(a, msize, vsize) - mpz_addmul(sum.value, m._matrix[k][a], v._entries[ra]) - + m.get_unsafe_mpz(k,a,tmp) + mpz_addmul(sum.value, tmp, v._entries[ra]) + mpz_clear(tmp) return sum def min_max_delta_intvec(Vector_integer_dense a, Vector_integer_dense b): @@ -4622,4 +4623,3 @@ def min_max_diff_doublevec(Vector_real_double_dense c): max_diff = diff return (min_diff, max_diff) - diff --git a/src/sage/schemes/toric/ideal.py b/src/sage/schemes/toric/ideal.py index 3a207a06db5..92ba7b5b97a 100644 --- a/src/sage/schemes/toric/ideal.py +++ b/src/sage/schemes/toric/ideal.py @@ -24,7 +24,7 @@ sage: IA.ker() Free module of degree 3 and rank 1 over Integer Ring User basis matrix: - [ 1 -2 1] + [-1 2 -1] sage: IA Ideal (-z1^2 + z0*z2) of Multivariate Polynomial Ring in z0, z1, z2 over Rational Field @@ -39,8 +39,8 @@ sage: IA.ker() Free module of degree 4 and rank 2 over Integer Ring User basis matrix: - [ 1 -1 -1 1] - [ 1 -2 1 0] + [-1 1 1 -1] + [-1 2 -1 0] sage: IA Ideal (-z1*z2 + z0*z3, -z1^2 + z0*z2, z2^2 - z1*z3) of Multivariate Polynomial Ring in z0, z1, z2, z3 over Rational Field @@ -300,7 +300,7 @@ def ker(self): sage: IA.ker() Free module of degree 3 and rank 1 over Integer Ring User basis matrix: - [ 1 -2 1] + [-1 2 -1] """ if '_ker' in self.__dict__: @@ -381,9 +381,9 @@ def _naive_ideal(self, ring): sage: IA.ker() Free module of degree 3 and rank 1 over Integer Ring User basis matrix: - [ 1 -2 1] + [-1 2 -1] sage: IA._naive_ideal(IA.ring()) - Ideal (-z1^2 + z0*z2) of Multivariate Polynomial Ring in z0, z1, z2 over Rational Field + Ideal (z1^2 - z0*z2) of Multivariate Polynomial Ring in z0, z1, z2 over Rational Field """ x = ring.gens() binomials = [] diff --git a/src/sage/structure/element.pyx b/src/sage/structure/element.pyx index 2511476155c..dc0a0c10be5 100644 --- a/src/sage/structure/element.pyx +++ b/src/sage/structure/element.pyx @@ -2736,11 +2736,54 @@ cdef class Matrix(ModuleElement): sage: a = matrix(ZZ, 2, range(4)) sage: a / 5 - [ 0 1/5] + [ 0 1/5] [2/5 3/5] - """ + sage: a = matrix(ZZ, 2, range(4)) + sage: b = matrix(ZZ, 2, [1,1,0,5]) + sage: a / b + [ 0 1/5] + [ 2 1/5] + sage: c = matrix(QQ, 2, [3,2,5,7]) + sage: c / a + [-5/2 3/2] + [-1/2 5/2] + sage: a / c + [-5/11 3/11] + [-1/11 5/11] + sage: a / 7 + [ 0 1/7] + [2/7 3/7] + + Other rings work just as well:: + + sage: a = matrix(GF(3),2,2,[0,1,2,0]) + sage: b = matrix(ZZ,2,2,[4,6,1,2]) + sage: a / b + [1 2] + [2 0] + sage: c = matrix(GF(3),2,2,[1,2,1,1]) + sage: a / c + [1 2] + [1 1] + sage: a = matrix(RDF,2,2,[.1,-.4,1.2,-.6]) + sage: b = matrix(RDF,2,2,[.3,.1,-.5,1.3]) + sage: a / b # rel tol 1e-10 + [-0.15909090909090906 -0.29545454545454547] + [ 2.863636363636364 -0.6818181818181817] + sage: R. = ZZ['t'] + sage: a = matrix(R,2,2,[t^2,t+1,-t,t+2]) + sage: b = matrix(R,2,2,[t^3-1,t,-t+3,t^2]) + sage: a / b + [ (t^4 + t^2 - 2*t - 3)/(t^5 - 3*t) (t^4 - t - 1)/(t^5 - 3*t)] + [ (-t^3 + t^2 - t - 6)/(t^5 - 3*t) (t^4 + 2*t^3 + t^2 - t - 2)/(t^5 - 3*t)] + """ + cdef Matrix rightinv if have_same_parent(left, right): - return (left)._matrix_times_matrix_(~right) + rightinv = ~right + if have_same_parent(left,rightinv): + return (left)._matrix_times_matrix_(rightinv) + else: + return ((left.change_ring(rightinv.parent().base_ring())))._matrix_times_matrix_(rightinv) else: global coercion_model return coercion_model.bin_op(left, right, div)