Ticket #4260: trac_4260-dense_ctypes_template.patch

File trac_4260-dense_ctypes_template.patch, 122.7 KB (added by malb, 9 years ago)

add templated classes for float and double representation

  • module_list.py

    # HG changeset patch
    # User Burcin Erocal <burcin@erocal.org>
    # Date 1267137620 -3600
    # Node ID 9b63c2421a7f5c12f483f1bba6c1629cfe4acedb
    # Parent  d263952c945f384fc3b411284d21a5ac9de7f575
    trac 4260: create templated classes for mod n dense matrices which store floats and doubles.
    
    diff -r d263952c945f -r 9b63c2421a7f module_list.py
    a b  
    899899              sources = ['sage/matrix/matrix_modn_dense.pyx'],
    900900              libraries = ['gmp']),
    901901
     902    Extension('sage.matrix.matrix_modn_dense_float',
     903              sources = ['sage/matrix/matrix_modn_dense_float.pyx'],
     904              language="c++",
     905              libraries = ['gmp', 'linbox', BLAS, BLAS2],
     906              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
     907
     908    Extension('sage.matrix.matrix_modn_dense_double',
     909              sources = ['sage/matrix/matrix_modn_dense_double.pyx'],
     910              language="c++",
     911              libraries = ['gmp', 'linbox', BLAS, BLAS2],
     912              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
     913
    902914    Extension('sage.matrix.matrix_modn_sparse',
    903915              sources = ['sage/matrix/matrix_modn_sparse.pyx'],
    904916              libraries = ['gmp']),
  • new file sage/libs/linbox/fflas.pxd

    diff -r d263952c945f -r 9b63c2421a7f sage/libs/linbox/fflas.pxd
    - +  
     1from modular cimport ModDoubleField, ModFloatField, ModDoubleFieldElement, ModFloatFieldElement
     2
     3ctypedef unsigned long mod_int
     4
     5cdef extern from "linbox/fflas/fflas.h" namespace "std":
     6    cdef cppclass vector[T]:
     7        cppclass iterator:
     8            T operator*()
     9            iterator operator++()
     10            bint operator==(iterator)
     11            bint operator!=(iterator)
     12        vector()
     13        void push_back(T&)
     14        T& operator[](int)
     15        T& at(int)
     16        iterator begin()
     17        iterator end()
     18        size_t size()
     19
     20    cdef cppclass list[T]:
     21        cppclass iterator:
     22            T operator*()
     23            iterator operator++()
     24            bint operator==(iterator)
     25            bint operator!=(iterator)
     26        void push_back(T&)
     27        void pop_front()
     28        T& front()
     29        iterator begin()
     30        iterator end()
     31        size_t size()
     32        void clear()
     33
     34cdef extern from "linbox/fflas/fflas.h":
     35    ctypedef enum fflas_trans_enum "LinBox::FFLAS::FFLAS_TRANSPOSE":
     36        fflas_no_trans  "LinBox::FFLAS::FflasNoTrans"
     37        fflas_trans  "LinBox::FFLAS::FflasTrans"
     38
     39    ctypedef enum fflas_side_enum "LinBox::FFLAS::FFLAS_SIDE":
     40        fflas_right  "LinBox::FFLAS::FflasRight"
     41
     42    # double
     43    void ModDouble_fgemv "LinBox::FFLAS::fgemv" \
     44            (ModDoubleField F, fflas_trans_enum transA,
     45                    size_t nrows, size_t ncols,
     46                    ModDoubleFieldElement alpha, ModDoubleFieldElement* A,
     47                    size_t lda, ModDoubleFieldElement* X, size_t incX,
     48                    ModDoubleFieldElement beta, ModDoubleFieldElement* Y,
     49                    size_t incY)
     50
     51    ModDoubleFieldElement* ModDouble_fgemm "LinBox::FFLAS::fgemm" \
     52            (ModDoubleField F,
     53                    fflas_trans_enum transA, fflas_trans_enum transB,
     54                    size_t nrowsA, size_t ncolsB, size_t ncolsA,
     55                    ModDoubleFieldElement alpha, ModDoubleFieldElement* A,
     56                    size_t A_stride, ModDoubleFieldElement* B, int B_stride,
     57                    ModDoubleFieldElement beta, ModDoubleFieldElement* C,
     58                    size_t C_stride)
     59
     60
     61    # float
     62    void ModFloat_fgemv "LinBox::FFLAS::fgemv" \
     63            (ModFloatField F, fflas_trans_enum transA,
     64                    size_t nrows, size_t ncols,
     65                    ModFloatFieldElement alpha, ModFloatFieldElement* A,
     66                    size_t lda, ModFloatFieldElement* X, size_t incX,
     67                    ModFloatFieldElement beta, ModFloatFieldElement* Y,
     68                    size_t incY)
     69
     70    ModFloatFieldElement* ModFloat_fgemm "LinBox::FFLAS::fgemm" \
     71            (ModFloatField F,
     72                    fflas_trans_enum transA, fflas_trans_enum transB,
     73                    size_t nrowsA, size_t ncolsB, size_t ncolsA,
     74                    ModFloatFieldElement alpha, ModFloatFieldElement* A,
     75                    size_t A_stride, ModFloatFieldElement* B, int B_stride,
     76                    ModFloatFieldElement beta, ModFloatFieldElement* C,
     77                    size_t C_stride)
     78
     79cdef extern from "linbox/ffpack/ffpack.h":
     80    # double
     81    bint ModDouble_is_singular "LinBox::FFPACK::IsSingular" (ModDoubleField F,
     82                                                             size_t nrows, size_t ncols, ModDoubleFieldElement* A,
     83                                                             size_t A_stride)
     84
     85    ModDoubleFieldElement* ModDouble_invert_in_place "LinBox::FFPACK::Invert"  (ModDoubleField F, size_t order,
     86                                                                                ModDoubleFieldElement* A, size_t A_stride, int nullity)
     87
     88    ModDoubleFieldElement ModDoubleDet "LinBox::FFPACK::Det" (ModDoubleField F,
     89                                                              size_t nrows, size_t ncols,
     90                                                              ModDoubleFieldElement* A, size_t A_stride)
     91
     92    int ModDoubleRank "LinBox::FFPACK::Rank" (ModDoubleField,
     93                                              size_t nrows, size_t ncols,
     94                                              ModDoubleFieldElement *A, size_t lda)
     95
     96    size_t ModDouble_echelon "LinBox::FFPACK::ReducedRowEchelonForm" (ModDoubleField F, size_t a, size_t b,
     97                                                                      ModDoubleFieldElement* matrix,
     98                                                                      size_t s, size_t* P, size_t* Q)
     99   
     100    void ModDouble_applyp "LinBox::FFPACK::applyP" (ModDoubleField F,
     101                                                    fflas_side_enum s, fflas_trans_enum tr,
     102                                                    size_t nr, size_t foo, size_t r,
     103                                                    ModDoubleFieldElement* matrix, size_t nc, size_t* Q)
     104   
     105    void ModDouble_MinPoly "LinBox::FFPACK::MinPoly" ( ModDoubleField F,
     106                                                       vector[ModDoubleFieldElement] minP, size_t N,
     107                                                       ModDoubleFieldElement* A, size_t lda,
     108                                                       ModDoubleFieldElement* X, size_t ldx, size_t* P)
     109
     110    void ModDouble_CharPoly "LinBox::FFPACK::CharPoly" ( ModDoubleField F,
     111                                                         list[vector[ModDoubleFieldElement]] charp, size_t N,
     112                                                         ModDoubleFieldElement* A, size_t lda)
     113
     114    # float
     115
     116    bint ModFloat_is_singular "LinBox::FFPACK::IsSingular" (ModFloatField F,
     117                                                            size_t nrows, size_t ncols, ModFloatFieldElement* A,
     118                                                            size_t A_stride)
     119
     120    ModFloatFieldElement* ModFloat_invert_in_place "LinBox::FFPACK::Invert" (ModFloatField F, size_t order,
     121                                                                             ModFloatFieldElement* A, size_t A_stride, int nullity)
     122
     123    ModFloatFieldElement ModFloatDet "LinBox::FFPACK::Det" (ModFloatField F,
     124                                                            size_t nrows, size_t ncols,
     125                                                            ModFloatFieldElement* A, size_t A_stride)
     126
     127    int ModFloatRank "LinBox::FFPACK::Rank" (ModFloatField,
     128                                             size_t nrows, size_t ncols,
     129                                             ModFloatFieldElement *A, size_t lda)
     130
     131    size_t ModFloat_echelon "LinBox::FFPACK::ReducedRowEchelonForm" (ModFloatField F, size_t a, size_t b,
     132                                                                     ModFloatFieldElement* matrix,
     133                                                                     size_t s, size_t* P, size_t* Q)
     134
     135    void ModFloat_applyp "LinBox::FFPACK::applyP" (ModFloatField F,
     136                                                   fflas_side_enum s, fflas_trans_enum tr,
     137                                                   size_t nr, size_t foo, size_t r,
     138                                                   ModFloatFieldElement* matrix, size_t nc, size_t* Q)
     139
     140    void ModFloat_MinPoly "LinBox::FFPACK::MinPoly" ( ModFloatField F,
     141                                                      vector[ModFloatFieldElement] minP, size_t N,
     142                                                      ModFloatFieldElement* A, size_t lda,
     143                                                      ModFloatFieldElement* X, size_t ldx, size_t* P)
     144
     145    void ModFloat_CharPoly "LinBox::FFPACK::CharPoly" ( ModFloatField F,
     146                                                        list[vector[ModFloatFieldElement]] charp, size_t N,
     147                                                        ModFloatFieldElement* A, size_t lda )
  • new file sage/libs/linbox/modular.pxd

    diff -r d263952c945f -r 9b63c2421a7f sage/libs/linbox/modular.pxd
    - +  
     1cdef extern from "linbox/field/modular.h":
     2    # double
     3
     4    cdef cppclass ModDoubleFieldElement "LinBox::Modular<double>::Element":
     5        pass
     6
     7    cdef cppclass ModDoubleField "LinBox::Modular<double>":
     8        ModDoubleField(int modulus)
     9        ModDoubleFieldElement init(ModDoubleFieldElement res, int v)
     10        ModDoubleFieldElement init(ModDoubleFieldElement res, double v)
     11        ModDoubleFieldElement inv(  ModDoubleFieldElement x, ModDoubleFieldElement y)
     12        ModDoubleFieldElement neg(  ModDoubleFieldElement x, ModDoubleFieldElement y)
     13        ModDoubleFieldElement mul(  ModDoubleFieldElement r, ModDoubleFieldElement x, ModDoubleFieldElement y)
     14        ModDoubleFieldElement mulin(ModDoubleFieldElement x, ModDoubleFieldElement y)
     15        ModDoubleFieldElement addin(ModDoubleFieldElement x, ModDoubleFieldElement y)
     16        ModDoubleFieldElement invin(ModDoubleFieldElement y)
     17        ModDoubleFieldElement negin(ModDoubleFieldElement y)
     18        int characteristic(int c)
     19        bint isZero(ModDoubleFieldElement x)
     20
     21
     22    # float
     23    cdef cppclass ModFloatFieldElement "LinBox::Modular<float>::Element":
     24        pass
     25
     26    cdef cppclass ModFloatField "LinBox::Modular<float>":
     27        ModFloatField(int modulus)
     28        ModFloatFieldElement init(ModFloatFieldElement res, int v)
     29        ModFloatFieldElement init(ModFloatFieldElement res, double v)
     30        ModFloatFieldElement inv(  ModFloatFieldElement x, ModFloatFieldElement y)
     31        ModFloatFieldElement neg(  ModFloatFieldElement x, ModFloatFieldElement y)
     32        ModFloatFieldElement mul(  ModFloatFieldElement r, ModFloatFieldElement x, ModFloatFieldElement y)
     33        ModFloatFieldElement mulin(ModFloatFieldElement x, ModFloatFieldElement y)
     34        ModFloatFieldElement addin(ModFloatFieldElement x, ModFloatFieldElement y)       
     35        ModFloatFieldElement invin(ModFloatFieldElement y)
     36        ModFloatFieldElement negin(ModFloatFieldElement y)
     37        int characteristic(int c)
     38        bint isZero(ModFloatFieldElement x)
     39       
     40
  • sage/matrix/matrix0.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix0.pyx
    a b  
    40964096       
    40974097            sage: m = matrix(Zmod(49),2,[2,1,3,3])
    40984098            sage: type(m)
    4099             <type 'sage.matrix.matrix_modn_dense.Matrix_modn_dense'>
     4099            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
    41004100            sage: ~m
    41014101            [ 1 16]
    41024102            [48 17]
  • sage/matrix/matrix2.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix2.pyx
    a b  
    10911091        # As of Sage 3.4, computing determinants directly in Z/nZ for
    10921092        # n composite is too slow, so we lift to Z and compute there.
    10931093        if is_IntegerModRing(R):
    1094             from matrix_modn_dense import Matrix_modn_dense
    1095             if not (isinstance(self, Matrix_modn_dense) and R.characteristic().is_prime()):
     1094            from matrix_modn_dense import is_Matrix_modn_dense
     1095            if not (is_Matrix_modn_dense(self) and R.characteristic().is_prime()):
    10961096                return R(self.lift().det())
    10971097       
    10981098        # N.B.  The following comment should be obsolete now that the generic
  • sage/matrix/matrix_cyclo_dense.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_cyclo_dense.pyx
    a b  
    4848from matrix cimport Matrix
    4949import matrix_dense
    5050from matrix_integer_dense import _lift_crt
    51 from matrix_modn_dense import _matrix_from_rows_of_matrices
    5251from sage.structure.element cimport Matrix as baseMatrix
    5352from misc import matrix_integer_dense_rational_reconstruction
    5453from sage.ext.multi_modular import MAX_MODULUS
     
    613612            Amodp, _ = self._reductions(p)
    614613            Bmodp, _ = right._reductions(p)
    615614            _,     S = self._reduction_matrix(p)
    616             X = _matrix_from_rows_of_matrices([Amodp[i] * Bmodp[i] for i
    617                                                in range(len(Amodp))])
     615            X = Amodp[0]._matrix_from_rows_of_matrices([Amodp[i] * Bmodp[i] for i in range(len(Amodp))])
    618616            v.append(S*X)
    619617            p = previous_prime(p)
    620618        M = matrix(ZZ, self._base_ring.degree(), self._nrows*right.ncols())
  • sage/matrix/matrix_integer_dense.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_integer_dense.pyx
    a b  
    101101from sage.structure.element import is_Vector
    102102from sage.structure.sequence import Sequence
    103103
    104 from matrix_modn_dense import Matrix_modn_dense
    105 from matrix_modn_dense cimport Matrix_modn_dense
     104from matrix_modn_dense_float cimport Matrix_modn_dense_template
     105from matrix_modn_dense_float cimport Matrix_modn_dense_float
     106from matrix_modn_dense_double cimport Matrix_modn_dense_double
     107
    106108from matrix_mod2_dense import Matrix_mod2_dense
    107109from matrix_mod2_dense cimport Matrix_mod2_dense
    108110
     111from matrix_modn_dense cimport is_Matrix_modn_dense
     112
     113
    109114from matrix2 import decomp_seq
    110115
    111116import sage.modules.free_module
     
    12211226        return res
    12221227
    12231228    cdef _mod_int_c(self, mod_int p):
     1229        from matrix_modn_dense_float import MAX_MODULUS as MAX_MODULUS_FLOAT
     1230        from matrix_modn_dense_double import MAX_MODULUS as MAX_MODULUS_DOUBLE
     1231
     1232        cdef Py_ssize_t i, j
     1233        cdef mpz_t* self_row
     1234
     1235        cdef float* res_row_f
     1236        cdef Matrix_modn_dense_float res_f
     1237
     1238        cdef double* res_row_d
     1239        cdef Matrix_modn_dense_double res_d
     1240
    12241241        if p == 2:
    12251242            return self._mod_two()
    1226         cdef Py_ssize_t i, j
    1227         cdef Matrix_modn_dense res
    1228         cdef mpz_t* self_row
    1229         cdef mod_int* res_row
    1230         res = Matrix_modn_dense.__new__(Matrix_modn_dense,
    1231                         matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
    1232         for i from 0 <= i < self._nrows:
    1233             self_row = self._matrix[i]
    1234             res_row = res._matrix[i]
    1235             for j from 0 <= j < self._ncols:
    1236                 res_row[j] = mpz_fdiv_ui(self_row[j], p)
    1237         return res
    1238        
     1243        elif p < MAX_MODULUS_FLOAT:
     1244            res_f = Matrix_modn_dense_float.__new__(Matrix_modn_dense_float,
     1245                                                    matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
     1246            for i from 0 <= i < self._nrows:
     1247                self_row = self._matrix[i]
     1248                res_row_f = res_f._matrix[i]
     1249                for j from 0 <= j < self._ncols:
     1250                    res_row_f[j] = <float>mpz_fdiv_ui(self_row[j], p)
     1251            return res_f
     1252
     1253        elif p < MAX_MODULUS_DOUBLE:
     1254            res_d = Matrix_modn_dense_double.__new__(Matrix_modn_dense_double,
     1255                                                     matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
     1256            for i from 0 <= i < self._nrows:
     1257                self_row = self._matrix[i]
     1258                res_row_d = res_d._matrix[i]
     1259                for j from 0 <= j < self._ncols:
     1260                    res_row_d[j] = <double>mpz_fdiv_ui(self_row[j], p)
     1261            return res_d
     1262        else:
     1263            raise ValueError("p to big.")
     1264
    12391265    def _reduce(self, moduli):
     1266        from matrix_modn_dense_float import MAX_MODULUS as MAX_MODULUS_FLOAT
     1267        from matrix_modn_dense_double import MAX_MODULUS as MAX_MODULUS_DOUBLE
    12401268       
    12411269        if isinstance(moduli, (int, long, Integer)):
    12421270            return self._mod_int(moduli)
     
    12461274        cdef MultiModularBasis mm
    12471275        mm = moduli               
    12481276
    1249         res = [Matrix_modn_dense.__new__(Matrix_modn_dense,
    1250                                          matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False),
    1251                                          None, None, None) for p in mm]
     1277        res = []
     1278        for p in mm:
     1279            if p < MAX_MODULUS_FLOAT:
     1280                res.append( Matrix_modn_dense_float.__new__(Matrix_modn_dense_float,
     1281                                                            matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False),
     1282                                                            None, None, None) )
     1283            elif p < MAX_MODULUS_DOUBLE:
     1284                res.append( Matrix_modn_dense_double.__new__(Matrix_modn_dense_double,
     1285                                                             matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False),
     1286                                                             None, None, None) )
     1287            else:
     1288                raise ValueError("p=%d too big."%p)
    12521289
    12531290        cdef size_t i, k, n
    12541291        cdef Py_ssize_t nr, nc
     
    12611298        row_list = <mod_int**>sage_malloc(sizeof(mod_int*) * n)
    12621299        if row_list == NULL:
    12631300            raise MemoryError("out of memory allocating multi-modular coefficient list")
     1301        for k in range(n):
     1302            row_list[k] = <mod_int*>sage_malloc(sizeof(mod_int)*nc)
     1303            if row_list[k] == NULL:
     1304                raise MemoryError("out of memory allocating multi-modular coefficient list")
    12641305
    12651306        sig_on()
    12661307        for i from 0 <= i < nr:
    12671308            for k from 0 <= k < n:
    1268                 row_list[k] = (<Matrix_modn_dense>res[k])._matrix[i]
    1269             mm.mpz_reduce_vec(self._matrix[i], row_list, nc)
     1309                mm.mpz_reduce_vec(self._matrix[i], row_list, nc)
     1310                if isinstance(res[k], Matrix_modn_dense_float):
     1311                    for j in range(nc):
     1312                        (<Matrix_modn_dense_float>res[k])._matrix[i][j] = (<float>row_list[k][j])%(<Matrix_modn_dense_float>res[k]).p
     1313                else:
     1314                    for j in range(nc):
     1315                        (<Matrix_modn_dense_double>res[k])._matrix[i][j] = (<double>row_list[k][j])%(<Matrix_modn_dense_double>res[k]).p
    12701316        sig_off()
    12711317       
     1318        for k in range(n):
     1319            sage_free(row_list[k])
    12721320        sage_free(row_list)
    12731321        return res
    12741322
     
    50035051    mm = moduli
    50045052   
    50055053    for b in residues:
    5006         if not PY_TYPE_CHECK(b, Matrix_modn_dense):
    5007             raise TypeError("Can only perform CRT on list of type Matrix_modn_dense.")
     5054        if not is_Matrix_modn_dense(b):
     5055            raise TypeError("Can only perform CRT on list of matrices mod n.")
    50085056
    50095057    cdef mod_int **row_list
    50105058    row_list = <mod_int**>sage_malloc(sizeof(mod_int*) * n)
    50115059    if row_list == NULL:
    50125060        raise MemoryError("out of memory allocating multi-modular coefficient list")
    5013    
     5061
    50145062    sig_on()
    5015     for i from 0 <= i < nr:
    5016         for k from 0 <= k < n:
    5017             row_list[k] = (<Matrix_modn_dense>residues[k])._matrix[i]
     5063    for k in range(n):
     5064        row_list[k] = <mod_int *>sage_malloc(sizeof(mod_int) * nc)
     5065        if row_list[k] == NULL:
     5066            raise MemoryError("out of memory allocating multi-modular coefficient list")
     5067
     5068    for i in range(nr):
     5069        for k in range(n):
     5070            (<Matrix_modn_dense_template>residues[k])._copy_row_to_mod_int_array(row_list[k],i)
    50185071        mm.mpz_crt_vec(M._matrix[i], row_list, nc)
     5072
     5073    for k in range(n):
     5074        sage_free(row_list[k])
     5075    sage_free(row_list)
     5076
    50195077    sig_off()
    5020    
    5021     sage_free(row_list)
    50225078    return M
    50235079
    50245080#######################################################
  • sage/matrix/matrix_modn_dense.pxd

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense.pxd
    a b  
    2424    cdef _init_linbox(self)
    2525    cpdef _export_as_string(self)
    2626       
     27
     28cpdef is_Matrix_modn_dense(self)
  • sage/matrix/matrix_modn_dense.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense.pyx
    a b  
    1919    sage: a.rank()
    2020    2
    2121    sage: type(a)
    22     <type 'sage.matrix.matrix_modn_dense.Matrix_modn_dense'>
     22    <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
    2323    sage: a[0,0] = 5
    2424    sage: a.rank()
    2525    3
     
    127127cimport matrix0
    128128cimport sage.structure.element
    129129
     130from sage.matrix.matrix_modn_dense_float cimport Matrix_modn_dense_float
     131from sage.matrix.matrix_modn_dense_double cimport Matrix_modn_dense_double
    130132
    131133from sage.libs.linbox.linbox cimport Linbox_modn_dense
    132134cdef Linbox_modn_dense linbox
     
    153155
    154156from sage.structure.proof.proof import get_flag as get_proof_flag
    155157
    156 
    157158cdef long num = 1
    158159cdef bint little_endian = (<char*>(&num))[0]
    159160
    160 
    161161##############################################################################
    162162#       Copyright (C) 2004,2005,2006 William Stein <wstein@gmail.com>
    163163#  Distributed under the terms of the GNU General Public License (GPL)
     
    176176    # x * __richcmp__    -- always the same
    177177    ########################################################################
    178178    def __cinit__(self, parent, entries, copy, coerce):
     179        sage.misc.misc.deprecation("This class is replaced by Matrix_modn_dense_float/Matrix_modn_dense_double.")       
    179180   
    180181        matrix_dense.Matrix_dense.__init__(self, parent)
    181182
     
    219220            [6 5]
    220221            [4 2]
    221222        """
    222        
    223223        cdef mod_int e
    224224        cdef Py_ssize_t i, j, k
    225225        cdef mod_int *v
     
    17051705            [3 4 5]
    17061706            [6 0 1]
    17071707            sage: type(a)
    1708             <type 'sage.matrix.matrix_modn_dense.Matrix_modn_dense'>
     1708            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
    17091709       
    17101710        We test the optional check flag.
    17111711       
     
    18921892            memcpy(M._entries, self._entries+i*n, sizeof(mod_int)*n)
    18931893            ans.append(M)
    18941894        return ans
    1895        
    18961895
    1897 def _matrix_from_rows_of_matrices(X):
     1896    _matrix_from_rows_of_matrices = staticmethod(__matrix_from_rows_of_matrices)
     1897   
     1898def __matrix_from_rows_of_matrices(X):
    18981899    """
     1900    Return a matrix whose ith row is ``X[i].list()``.
     1901   
    18991902    INPUT:
    19001903   
     1904    - ``X`` - a nonempty list of matrices of the same size mod a
     1905       single modulus ``p``
    19011906   
    1902     -  ``X`` - a nonempty list of matrices of the same size
    1903        mod a single prime p
     1907    OUTPUT: A single matrix mod p whose ith row is ``X[i].list()``.
    19041908   
    1905    
    1906     OUTPUT: A single matrix mod p whose ith row is X[i].list().
     1909
    19071910    """
    19081911    # The code below is just a fast version of the following:
    19091912    ##     from constructor import matrix
     
    19181921
    19191922    T = X[0]
    19201923    m = T._nrows * T._ncols
    1921     A = Matrix_modn_dense.__new__(Matrix_modn_dense,
    1922            MatrixSpace(X[0].base_ring(), n, m),
    1923                                   0, 0, 0)
     1924    A = Matrix_modn_dense.__new__(Matrix_modn_dense, MatrixSpace(X[0].base_ring(), n, m), 0, 0, 0)
    19241925    A.p = T.p
    19251926    A.gather = T.gather
    19261927
     
    19291930        memcpy(A._entries + i*m, T._entries, sizeof(mod_int)*m)
    19301931    return A
    19311932
     1933cpdef is_Matrix_modn_dense(self):
     1934    """
     1935    """
     1936    return isinstance(self, Matrix_modn_dense) | isinstance(self, Matrix_modn_dense_float) | isinstance(self, Matrix_modn_dense_double)
  • new file sage/matrix/matrix_modn_dense_double.pxd

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_double.pxd
    - +  
     1include "../ext/cdefs.pxi"
     2
     3ctypedef double celement
     4
     5include "matrix_modn_dense_template_header.pxi"
     6
     7cdef class Matrix_modn_dense_double(Matrix_modn_dense_template):
     8    pass
  • new file sage/matrix/matrix_modn_dense_double.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_double.pyx
    - +  
     1"""
     2Dense matrices over `\ZZ/n\ZZ` for `n < 2^{23}` using LinBox's ``Modular<float>``
     3
     4AUTHORS:
     5- Burcin Erocal
     6- Martin Albrecht
     7"""
     8###############################################################################
     9#   SAGE: Open Source Mathematical Software
     10#       Copyright (C) 2011 Burcin Erocal <burcin@erocal.org>
     11#       Copyright (C) 2011 Martin Albrecht <martinralbrecht@googlemail.com>
     12#  Distributed under the terms of the GNU General Public License (GPL),
     13#  version 2 or any later version.  The full text of the GPL is available at:
     14#                  http://www.gnu.org/licenses/
     15###############################################################################
     16
     17include "../ext/stdsage.pxi"
     18include "../ext/interrupt.pxi"
     19
     20# randstate in template needs this
     21include '../ext/random.pxi'
     22
     23from sage.libs.linbox.modular cimport ModDoubleField as ModField, ModDoubleFieldElement as ModFieldElement
     24
     25from sage.libs.linbox.fflas cimport ModDouble_fgemm as Mod_fgemm, ModDouble_fgemv as Mod_fgemv, \
     26    ModDoubleDet as ModDet, \
     27    ModDoubleRank as ModRank, ModDouble_echelon as Mod_echelon, \
     28    ModDouble_applyp as Mod_applyp, \
     29    ModDouble_MinPoly as Mod_MinPoly, \
     30    ModDouble_CharPoly as Mod_CharPoly
     31
     32MAX_MODULUS = 2**23
     33
     34cdef extern from "../rings/finite_rings/stdint.h":
     35    ctypedef int int_fast32_t
     36    ctypedef int int_fast64_t
     37    int_fast32_t INTEGER_MOD_INT32_LIMIT
     38    int_fast64_t INTEGER_MOD_INT64_LIMIT
     39
     40from sage.rings.finite_rings.integer_mod cimport IntegerMod_int64
     41
     42include "matrix_modn_dense_template.pxi"
     43
     44
     45cdef class Matrix_modn_dense_double(Matrix_modn_dense_template):
     46    """
     47    Dense matrices over `\ZZ/n\ZZ` for `n < 2^{23}` using LinBox's ``Modular<float>``
     48    """
     49    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value):
     50        """
     51        Set the (i,j) entry of self to the int value.
     52
     53        EXAMPLE::
     54
     55            sage: A = random_matrix(GF(3016963), 4, 4); A
     56            [ 220081 2824836  765701 2282256]
     57            [1795330  767112 2967421 1373921]
     58            [2757699 1142917 2720973 2877160]
     59            [1674049 1341486 2641133 2173280]
     60
     61            sage: A[0,0] = 220082r; A
     62            [ 220082 2824836  765701 2282256]
     63            [1795330  767112 2967421 1373921]
     64            [2757699 1142917 2720973 2877160]
     65            [1674049 1341486 2641133 2173280]
     66           
     67            sage: a = A[0,0]; a
     68            220082
     69
     70            sage: ~a
     71            2859358
     72
     73            sage: A = random_matrix(Integers(5099106), 4, 4); A
     74            [2629491 1237101 2033003 3788106]
     75            [4649912 1157595 4928315 4382585]
     76            [4252686  978867 2601478 1759921]
     77            [1303120 1860486 3405811 2203284]
     78
     79            sage: A[0,0] = 220082r; A
     80            [ 220082 1237101 2033003 3788106]
     81            [4649912 1157595 4928315 4382585]
     82            [4252686  978867 2601478 1759921]
     83            [1303120 1860486 3405811 2203284]
     84           
     85            sage: a = A[0,0]; a
     86            220082
     87
     88            sage: a*a
     89            4777936
     90        """
     91        self._matrix[i][j] = <double>value
     92
     93    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, x):
     94        """
     95        EXAMPLE::
     96
     97            sage: A = random_matrix(GF(3016963), 4, 4); A
     98            [ 220081 2824836  765701 2282256]
     99            [1795330  767112 2967421 1373921]
     100            [2757699 1142917 2720973 2877160]
     101            [1674049 1341486 2641133 2173280]
     102
     103            sage: K = A.base_ring()
     104            sage: A[0,0] = K(220082); A
     105            [ 220082 2824836  765701 2282256]
     106            [1795330  767112 2967421 1373921]
     107            [2757699 1142917 2720973 2877160]
     108            [1674049 1341486 2641133 2173280]
     109           
     110            sage: a = A[0,0]; a
     111            220082
     112
     113            sage: ~a
     114            2859358
     115
     116            sage: A = random_matrix(Integers(5099106), 4, 4); A
     117            [2629491 1237101 2033003 3788106]
     118            [4649912 1157595 4928315 4382585]
     119            [4252686  978867 2601478 1759921]
     120            [1303120 1860486 3405811 2203284]
     121
     122            sage: K = A.base_ring()
     123            sage: A[0,0] = K(220081)
     124           
     125            sage: a = A[0,0]; a
     126            220081
     127
     128            sage: a*a
     129            4337773
     130        """
     131        # doesn't do bounds or any other checks; assumes x is in self._base_ring
     132        self._matrix[i][j] = <double>(<IntegerMod_int>x).ivalue
     133
     134    cdef IntegerMod_abstract get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
     135        """
     136        EXAMPLE::
     137
     138            sage: A = random_matrix(GF(3016963), 4, 4); A
     139            [ 220081 2824836  765701 2282256]
     140            [1795330  767112 2967421 1373921]
     141            [2757699 1142917 2720973 2877160]
     142            [1674049 1341486 2641133 2173280]
     143
     144            sage: a = A[0,0]; a
     145            220081
     146
     147            sage: ~a
     148            697224
     149
     150            sage: K = A.base_ring()
     151            sage: ~K(220081)
     152            697224
     153
     154            sage: A = random_matrix(Integers(5099106), 4, 4); A
     155            [2629491 1237101 2033003 3788106]
     156            [4649912 1157595 4928315 4382585]
     157            [4252686  978867 2601478 1759921]
     158            [1303120 1860486 3405811 2203284]
     159
     160            sage: a = A[0,1]; a
     161            1237101
     162
     163            sage: a*a
     164            3803997
     165
     166            sage: K = A.base_ring()
     167            sage: K(1237101)^2
     168            3803997
     169        """
     170        # doesn't do checks
     171        if (<Matrix_modn_dense_template>self).p <= INTEGER_MOD_INT32_LIMIT:
     172            return IntegerMod_int(self._base_ring, <mod_int>(<Matrix_modn_dense_template>self)._matrix[i][j])
     173        else:
     174            return IntegerMod_int64(self._base_ring, <mod_int>(<Matrix_modn_dense_template>self)._matrix[i][j])
     175
  • new file sage/matrix/matrix_modn_dense_float.pxd

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_float.pxd
    - +  
     1include "../ext/cdefs.pxi"
     2
     3ctypedef float celement
     4
     5include "matrix_modn_dense_template_header.pxi"
     6
     7cdef class Matrix_modn_dense_float(Matrix_modn_dense_template):
     8    pass
     9
  • new file sage/matrix/matrix_modn_dense_float.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_float.pyx
    - +  
     1"""
     2Dense matrices over `\ZZ/n\ZZ` for `n < 2^{11}` using LinBox's ``Modular<float>``
     3
     4AUTHORS:
     5- Burcin Erocal
     6- Martin Albrecht
     7"""
     8###############################################################################
     9#   SAGE: Open Source Mathematical Software
     10#       Copyright (C) 2011 Burcin Erocal <burcin@erocal.org>
     11#       Copyright (C) 2011 Martin Albrecht <martinralbrecht@googlemail.com>
     12#  Distributed under the terms of the GNU General Public License (GPL),
     13#  version 2 or any later version.  The full text of the GPL is available at:
     14#                  http://www.gnu.org/licenses/
     15###############################################################################
     16
     17include "../ext/stdsage.pxi"
     18include "../ext/interrupt.pxi"
     19
     20# randstate in template needs this
     21include '../ext/random.pxi'
     22
     23from sage.libs.linbox.modular cimport ModFloatField as ModField, ModFloatFieldElement as ModFieldElement
     24
     25from sage.libs.linbox.fflas cimport ModFloat_fgemm as Mod_fgemm, ModFloat_fgemv as Mod_fgemv, \
     26        ModFloatDet as ModDet, \
     27        ModFloatRank as ModRank, ModFloat_echelon as Mod_echelon, \
     28        ModFloat_applyp as Mod_applyp, \
     29        ModFloat_MinPoly as Mod_MinPoly, \
     30        ModFloat_CharPoly as Mod_CharPoly
     31
     32# LinBox supports up to 2^11 using float but that's double dog slow,
     33# so we pick a smaller value for crossover
     34MAX_MODULUS = 2**8
     35
     36include "matrix_modn_dense_template.pxi"
     37
     38cdef class Matrix_modn_dense_float(Matrix_modn_dense_template):
     39    """
     40    Dense matrices over `\ZZ/n\ZZ` for `n < 2^{11}` using LinBox's ``Modular<float>``
     41    """
     42    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value):
     43        """
     44        Set the (i,j) entry of self to the int value.
     45        """
     46        self._matrix[i][j] = <float>value
     47
     48    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, x):
     49        # doesn't do bounds or any other checks; assumes x is in self._base_ring
     50        self._matrix[i][j] = <float>(<IntegerMod_int>x).ivalue
     51
     52    cdef IntegerMod_int get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
     53        # doesn't do checks
     54        return IntegerMod_int(self._base_ring, <mod_int>(<Matrix_modn_dense_template>self)._matrix[i][j])
  • new file sage/matrix/matrix_modn_dense_template.pxi

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_template.pxi
    - +  
     1"""
     2Dense matrices over `\ZZ/n\ZZ` for `n` small using the LinBox library (FFLAS/FFPACK).
     3
     4FFLAS/FFPACK are libraries to provide BLAS/LAPACK-style routines for
     5working with finite fields. Additionally, these routines reduce to
     6BLAS/LAPACK routines using floating point arithmetic.
     7"""
     8
     9###############################################################################
     10#   SAGE: Open Source Mathematical Software
     11#       Copyright (C) 2004,2005,2006 William Stein <wstein@gmail.com>
     12#       Copyright (C) 2011 Burcin Erocal <burcin@erocal.org>
     13#       Copyright (C) 2011 Martin Albrecht <martinralbrecht@googlemail.com>
     14#  Distributed under the terms of the GNU General Public License (GPL),
     15#  version 2 or any later version.  The full text of the GPL is available at:
     16#                  http://www.gnu.org/licenses/
     17###############################################################################
     18
     19from sage.libs.linbox.fflas cimport fflas_trans_enum, fflas_no_trans, fflas_trans, \
     20    fflas_right, vector, list as std_list
     21
     22import sage.ext.multi_modular
     23cimport sage.rings.fast_arith
     24cdef sage.rings.fast_arith.arith_int ArithIntObj
     25ArithIntObj  = sage.rings.fast_arith.arith_int()
     26
     27cdef extern from *:
     28    void* memcpy(void* dst, void* src, long n)
     29
     30cdef extern from "stdio.h":
     31    int snprintf(char *str, size_t size, char *format, ...)
     32
     33from sage.modules.vector_modn_dense cimport Vector_modn_dense
     34
     35from sage.rings.arith import is_prime
     36from sage.structure.element cimport ModuleElement
     37
     38cimport matrix_dense
     39cimport matrix
     40cimport matrix0
     41cimport sage.structure.element
     42
     43from sage.structure.element cimport Matrix
     44from sage.rings.finite_rings.integer_mod cimport IntegerMod_int, IntegerMod_abstract
     45from sage.misc.misc import verbose, get_verbose, cputime
     46from sage.rings.integer cimport Integer
     47from sage.structure.element cimport ModuleElement, RingElement, Element, Vector
     48from matrix_integer_dense cimport Matrix_integer_dense
     49from sage.rings.integer_ring  import ZZ
     50from sage.structure.proof.proof import get_flag as get_proof_flag
     51
     52cdef long num = 1
     53cdef bint little_endian = (<char*>(&num))[0]
     54
     55cdef inline celement_invert(celement a, celement n):
     56    """
     57    Invert the finite field element `a` modulo `n`.
     58    """
     59    # This is copied from linbox source linbox/field/modular-float.h
     60    # The extended Euclidean algoritm
     61    cdef int x_int, y_int, q, tx, ty, temp
     62    x_int = <int>n
     63    y_int = <int>a
     64    tx = 0
     65    ty = 1
     66
     67    while y_int != 0:
     68        # always: gcd (n,residue) = gcd (x_int,y_int)
     69        #         sx*n + tx*residue = x_int
     70        #         sy*n + ty*residue = y_int
     71        q = x_int / y_int # integer quotient
     72        temp = y_int
     73        y_int = x_int - q * y_int
     74        x_int = temp
     75        temp = ty
     76        ty = tx - q * ty;
     77        tx = temp
     78
     79    if tx < 0:
     80         tx += <int>n
     81
     82    # now x_int = gcd (n,residue)
     83    return <celement>tx
     84
     85cdef inline bint linbox_is_zero(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols) except -1:
     86    """
     87    Return 1 if all entries of this matrix are zero.
     88    """
     89    cdef Py_ssize_t i, j
     90    for i in range(nrows):
     91        for j in range(ncols):
     92            if (entries+i*ncols+j)[0] != 0:
     93                return 0
     94    return 1
     95   
     96cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols):
     97    """
     98    Return the reduced row echelon form of this matrix.
     99
     100    EXAMPLES::
     101
     102        sage: A = random_matrix(GF(7), 10, 20); A
     103        [3 1 6 6 4 4 2 2 3 5 4 5 6 2 2 1 2 5 0 5]
     104        [3 2 0 5 0 1 5 4 2 3 6 4 5 0 2 4 2 0 6 3]
     105        [2 2 4 2 4 5 3 4 4 4 2 5 2 5 4 5 1 1 1 1]
     106        [0 6 3 4 2 2 3 5 1 1 4 2 6 5 6 3 4 5 5 3]
     107        [5 2 4 3 6 2 3 6 2 1 3 3 5 3 4 2 2 1 6 2]
     108        [0 5 6 3 2 5 6 6 3 2 1 4 5 0 2 6 5 2 5 1]
     109        [4 0 4 2 6 3 3 5 3 0 0 1 2 5 5 1 6 0 0 3]
     110        [2 0 1 0 0 3 0 2 4 2 2 4 4 4 5 4 1 2 3 4]
     111        [2 4 1 4 3 0 6 2 2 5 2 5 3 6 4 2 2 6 4 4]
     112        [0 0 2 2 1 6 2 0 5 0 4 3 1 6 0 6 0 4 6 5]
     113
     114        sage: A.echelon_form()
     115        [1 0 0 0 0 0 0 0 0 0 6 2 6 0 1 1 2 5 6 2]
     116        [0 1 0 0 0 0 0 0 0 0 0 4 5 4 3 4 2 5 1 2]
     117        [0 0 1 0 0 0 0 0 0 0 6 3 4 6 1 0 3 6 5 6]
     118        [0 0 0 1 0 0 0 0 0 0 0 3 5 2 3 4 0 6 5 3]
     119        [0 0 0 0 1 0 0 0 0 0 0 6 3 4 5 3 0 4 3 2]
     120        [0 0 0 0 0 1 0 0 0 0 1 1 0 2 4 2 5 5 5 0]
     121        [0 0 0 0 0 0 1 0 0 0 1 0 1 3 2 0 0 0 5 3]
     122        [0 0 0 0 0 0 0 1 0 0 4 4 2 6 5 4 3 4 1 0]
     123        [0 0 0 0 0 0 0 0 1 0 1 0 4 2 3 5 4 6 4 0]
     124        [0 0 0 0 0 0 0 0 0 1 2 0 5 0 5 5 3 1 1 4]
     125
     126    ::
     127
     128        sage: A = random_matrix(GF(13), 10, 10); A
     129        [ 8  3 11 11  9  4  8  7  9  9]
     130        [ 2  9  6  5  7 12  3  4 11  5]
     131        [12  6 11 12  4  3  3  8  9  5]
     132        [ 4  2 10  5 10  1  1  1  6  9]
     133        [12  8  5  5 11  4  1  2  8 11]
     134        [ 2  6  9 11  4  7  1  0 12  2]
     135        [ 8  9  0  7  7  7 10  4  1  4]
     136        [ 0  8  2  6  7  5  7 12  2  3]
     137        [ 2 11 12  3  4  7  2  9  6  1]
     138        [ 0 11  5  9  4  5  5  8  7 10]
     139
     140        sage: MS = parent(A)
     141        sage: B = A.augment(MS(1))
     142        sage: B.echelonize()
     143        sage: A.rank()
     144        10
     145        sage: C = B.submatrix(0,10,10,10); C
     146        [ 4  9  4  4  0  4  7 11  9 11]
     147        [11  7  6  8  2  8  6 11  9  5]
     148        [ 3  9  9  2  4  8  9  2  9  4]
     149        [ 7  0 11  4  0  9  6 11  8  1]
     150        [12 12  4 12  3 12  6  1  7 12]
     151        [12  2 11  6  6  6  7  0 10  6]
     152        [ 0  7  3  4  7 11 10 12  4  6]
     153        [ 5 11  0  5  3 11  4 12  5 12]
     154        [ 6  7  3  5  1  4 11  7  4  1]
     155        [ 4  9  6  7 11  1  2 12  6  7]
     156
     157        sage: ~A == C
     158        True
     159       
     160    ::
     161       
     162        sage: A = random_matrix(Integers(10), 10, 20); A
     163        [1 8 0 7 7 8 0 5 9 1 3 6 7 1 5 9 6 1 1 6]
     164        [2 7 2 7 3 7 8 3 9 3 6 1 3 4 4 8 1 9 1 1]
     165        [8 6 2 1 0 0 9 0 9 1 1 2 6 8 5 7 6 3 7 0]
     166        [9 8 7 4 8 0 0 8 0 3 8 4 4 8 5 1 7 9 9 7]
     167        [3 5 9 2 2 0 6 7 1 7 4 0 9 4 3 6 8 3 5 1]
     168        [2 8 3 7 7 6 9 4 2 7 6 5 7 0 4 8 6 7 5 3]
     169        [1 4 1 2 5 3 4 8 1 3 7 5 8 9 7 8 7 6 8 1]
     170        [5 1 8 4 6 3 8 8 1 3 4 8 4 0 7 8 7 7 5 9]
     171        [7 9 6 7 4 3 7 1 0 3 1 7 5 0 4 2 1 7 1 7]
     172        [3 7 1 8 2 2 6 0 1 8 9 3 4 8 1 0 6 3 3 8]
     173
     174        sage: A.echelon_form()
     175        Traceback (most recent call last):
     176        ...
     177        NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 10'.
     178       
     179    ::
     180        sage: A = random_matrix(GF(16007), 10, 20); A
     181        [15455  1177 10072  4693  3887  4102 10746 15265  6684 14559  4535 13921  9757  9525  9301  8566  2460  9609  3887  6205]
     182        [ 8602 10035  1242  9776   162  7893 12619  6660 13250  1988 14263 11377  2216  1247  7261  8446 15081 14412  7371  7948]
     183        [12634  7602   905  9617 13557  2694 13039  4936 12208 15480  3787 11229   593 12462  5123 14167  6460  3649  5821  6736]
     184        [10554  2511 11685 12325 12287  6534 11636  5004  6468  3180  3607 11627 13436  5106  3138 13376  8641  9093  2297  5893]
     185        [ 1025 11376 10288   609 12330  3021   908 13012  2112 11505    56  5971   338  2317  2396  8561  5593  3782  7986 13173]
     186        [ 7607   588  6099 12749 10378   111  2852 10375  8996  7969   774 13498 12720  4378  6817  6707  5299  9406 13318  2863]
     187        [15545   538  4840  1885  8471  1303 11086 14168  1853 14263  3995 12104  1294  7184  1188 11901 15971  2899  4632   711]
     188        [  584 11745  7540 15826 15027  5953  7097 14329 10889 12532 13309 15041  6211  1749 10481  9999  2751 11068    21  2795]
     189        [  761 11453  3435 10596  2173  7752 15941 14610  1072  8012  9458  5440   612 10581 10400   101 11472 13068  7758  7898]
     190        [10658  4035  6662   655  7546  4107  6987  1877  4072  4221  7679 14579  2474  8693  8127 12999 11141   605  9404 10003]
     191        sage: A.echelon_form()
     192        [    1     0     0     0     0     0     0     0     0     0  8416  8364 10318  1782 13872  4566 14855  7678 11899  2652]
     193        [    0     1     0     0     0     0     0     0     0     0  4782 15571  3133 10964  5581 10435  9989 14303  5951  8048]
     194        [    0     0     1     0     0     0     0     0     0     0 15688  6716 13819  4144   257  5743 14865 15680  4179 10478]
     195        [    0     0     0     1     0     0     0     0     0     0  4307  9488  2992  9925 13984 15754  8185 11598 14701 10784]
     196        [    0     0     0     0     1     0     0     0     0     0   927  3404 15076  1040  2827  9317 14041 10566  5117  7452]
     197        [    0     0     0     0     0     1     0     0     0     0  1144 10861  5241  6288  9282  5748  3715 13482  7258  9401]
     198        [    0     0     0     0     0     0     1     0     0     0   769  1804  1879  4624  6170  7500 11883  9047   874   597]
     199        [    0     0     0     0     0     0     0     1     0     0 15591 13686  5729 11259 10219 13222 15177 15727  5082 11211]
     200        [    0     0     0     0     0     0     0     0     1     0  8375 14939 13471 12221  8103  4212 11744 10182  2492 11068]
     201        [    0     0     0     0     0     0     0     0     0     1  6534   396  6780 14734  1206  3848  7712  9770 10755   410]
     202
     203    ::
     204
     205        sage: A = random_matrix(Integers(10000), 10, 20); A
     206        [2286 1825 2376 1180 1432 3070 5298 1485 2315   11  486 8911  669 2420 8838 9006 5606 1815 1700 3259]
     207        [5553 3349 3416 9701  784 2614 7654   74 7642 7631 3667 7491 5417 8046  235 1350 8458 1781 7254 3230]
     208        [3899 2557 9043 1740 4704  408 5572 2540 6823 2684 2483 6333 6137 2803 5750 7792 5897 3759 6960 1399]
     209        [3576 6431 1352 9019  632 5164 6432 9946 9449  665 3808 8426 5446 1299 6781 8902 3723 9639 4571 5725]
     210        [7153 7274 4275 5090 7116  504 2323 8981 7961 6643 4822 8069 7183 2766 9492 2429 9059 2620 9328 6563]
     211        [4730 3360 8630 3700 1781 5762 1988 3441  487 6869 7931 2534 2783 4002 4040 9831 7860  925 8792 1284]
     212        [8437 3824 3088 2626 6094 9571 7214 2194 2261 7330 7443 8896 4586 6386 2736 4395 4433 7213 7654 5854]
     213        [2256 7334 3566 7224 8490 4137 2690 2999  609 7291 7241  908 4810 5162 4079 5581 5198 2563 7058  910]
     214        [ 160 5757 3867 1620 3893 5557 2491 5934 4831 9953 5663 8183 1187 9725 3137 3159 1374 5277 8498 5701]
     215        [2640 2586 7674 3417 4878 7052 9606 7388  459 9237 9381 9652  604 5572 2728 4393 4059 8574 7845 7072]
     216
     217        sage: A.echelon_form()
     218        Traceback (most recent call last):
     219        ...
     220        NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 10000'.
     221
     222    TESTS::
     223
     224        sage: A = random_matrix(GF(7),  0, 10)
     225        sage: A.echelon_form()
     226        []
     227
     228        sage: A = random_matrix(GF(7), 10,  0)
     229        sage: A.echelon_form()
     230        []
     231
     232        sage: A = random_matrix(GF(7),  0,  0)
     233        sage: A.echelon_form()
     234        []
     235
     236        sage: A = matrix(GF(7),  10,  10)
     237        sage: A.echelon_form()
     238        [0 0 0 0 0 0 0 0 0 0]
     239        [0 0 0 0 0 0 0 0 0 0]
     240        [0 0 0 0 0 0 0 0 0 0]
     241        [0 0 0 0 0 0 0 0 0 0]
     242        [0 0 0 0 0 0 0 0 0 0]
     243        [0 0 0 0 0 0 0 0 0 0]
     244        [0 0 0 0 0 0 0 0 0 0]
     245        [0 0 0 0 0 0 0 0 0 0]
     246        [0 0 0 0 0 0 0 0 0 0]
     247        [0 0 0 0 0 0 0 0 0 0]
     248
     249        sage: A = random_matrix(GF(16007),  0, 10)
     250        sage: A.echelon_form()
     251        []
     252
     253        sage: A = random_matrix(GF(16007), 10,  0)
     254        sage: A.echelon_form()
     255        []
     256
     257        sage: A = random_matrix(GF(16007),  0,  0)
     258        sage: A.echelon_form()
     259        []
     260
     261        sage: A = matrix(GF(16007),  10,  10)
     262        sage: A.echelon_form()
     263        [0 0 0 0 0 0 0 0 0 0]
     264        [0 0 0 0 0 0 0 0 0 0]
     265        [0 0 0 0 0 0 0 0 0 0]
     266        [0 0 0 0 0 0 0 0 0 0]
     267        [0 0 0 0 0 0 0 0 0 0]
     268        [0 0 0 0 0 0 0 0 0 0]
     269        [0 0 0 0 0 0 0 0 0 0]
     270        [0 0 0 0 0 0 0 0 0 0]
     271        [0 0 0 0 0 0 0 0 0 0]
     272        [0 0 0 0 0 0 0 0 0 0]
     273    """
     274    if linbox_is_zero(modulus, entries, nrows, ncols):
     275        return 0,[]
     276
     277    cdef Py_ssize_t i, j
     278    cdef ModField *F = new ModField(<long>modulus)
     279    cdef size_t* P = <size_t*>sage_malloc(sizeof(size_t)*nrows)
     280    cdef size_t* Q = <size_t*>sage_malloc(sizeof(size_t)*ncols)
     281
     282    if nrows*ncols > 1000: sig_on()
     283    cdef Py_ssize_t r = Mod_echelon(F[0], nrows, ncols, <ModFieldElement*>entries, ncols, P, Q)
     284    if nrows*ncols > 1000: sig_off()
     285
     286    for i in range(nrows):
     287        for j in range(r):
     288            (entries+i*ncols+j)[0] = 0
     289        if i<r:
     290            (entries + i*(ncols+1))[0] = 1
     291
     292    Mod_applyp(F[0], fflas_right, fflas_no_trans, nrows, 0, r, <ModFieldElement*>entries, ncols, Q)
     293
     294    cdef list pivots = [int(Q[i]) for i in range(r)]
     295       
     296    sage_free(P)
     297    sage_free(Q)
     298    del F
     299    return r, pivots
     300
     301cdef inline celement *linbox_copy(celement modulus, celement *entries,  Py_ssize_t nrows, Py_ssize_t ncols) except NULL:
     302    """
     303    Create a copy of the entries array.
     304    """
     305    cdef celement *entries_copy = <celement*>sage_malloc(sizeof(celement)*nrows*ncols)
     306    memcpy(entries_copy, entries, sizeof(celement)*nrows*ncols)
     307    return entries_copy
     308
     309cdef inline int linbox_rank(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols) except -1:
     310    """
     311    Return the rank of this matrix.
     312
     313    EXAMPLES::
     314
     315        sage: A = random_matrix(GF(3), 100, 100)
     316        sage: B = copy(A)
     317        sage: A.rank()
     318        99
     319        sage: B == A
     320        True
     321
     322        sage: A = random_matrix(GF(3), 100, 100, density=0.01)
     323        sage: A.rank()
     324        63
     325
     326        sage: A = matrix(GF(3), 100, 100)
     327        sage: A.rank()
     328        0
     329       
     330    ::
     331
     332        sage: A = random_matrix(GF(16007), 100, 100)
     333        sage: B = copy(A)
     334        sage: A.rank()
     335        100
     336        sage: B == A
     337        True
     338        sage: MS = A.parent()
     339        sage: MS(1) == ~A*A
     340        True
     341
     342    TESTS::
     343
     344        sage: A = random_matrix(GF(7), 0, 0)
     345        sage: A.rank()
     346        0
     347
     348        sage: A = random_matrix(GF(7), 1, 0)
     349        sage: A.rank()
     350        0
     351
     352        sage: A = random_matrix(GF(7), 0, 1)
     353        sage: A.rank()
     354        0
     355
     356        sage: A = random_matrix(GF(16007), 0, 0)
     357        sage: A.rank()
     358        0
     359
     360        sage: A = random_matrix(GF(16007), 1, 0)
     361        sage: A.rank()
     362        0
     363
     364        sage: A = random_matrix(GF(16007), 0, 1)
     365        sage: A.rank()
     366        0
     367    """
     368    cdef ModField *F = new ModField(<long>modulus)
     369   
     370    cdef celement *cpy = linbox_copy(modulus, entries, nrows, ncols)
     371   
     372    if nrows*ncols > 1000: sig_on()
     373    r = ModRank(F[0], nrows, ncols, <ModFieldElement*>cpy, ncols)
     374    if nrows*ncols > 1000: sig_off()
     375    sage_free(cpy)
     376    del F
     377    return r
     378
     379cdef inline celement linbox_det(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols):
     380    """
     381    Return the determinant of this matrix.
     382
     383    EXAMPLES::
     384
     385        sage: A = random_matrix(GF(7), 10, 10); A
     386        [3 1 6 6 4 4 2 2 3 5]
     387        [4 5 6 2 2 1 2 5 0 5]
     388        [3 2 0 5 0 1 5 4 2 3]
     389        [6 4 5 0 2 4 2 0 6 3]
     390        [2 2 4 2 4 5 3 4 4 4]
     391        [2 5 2 5 4 5 1 1 1 1]
     392        [0 6 3 4 2 2 3 5 1 1]
     393        [4 2 6 5 6 3 4 5 5 3]
     394        [5 2 4 3 6 2 3 6 2 1]
     395        [3 3 5 3 4 2 2 1 6 2]
     396
     397        sage: A.determinant()
     398        6
     399       
     400    ::
     401
     402        sage: A = random_matrix(GF(7), 100, 100)
     403        sage: A.determinant()
     404        2
     405       
     406        sage: A.transpose().determinant()
     407        2
     408       
     409        sage: B = random_matrix(GF(7), 100, 100)
     410        sage: B.determinant()
     411        4
     412
     413        sage: (A*B).determinant() == A.determinant() * B.determinant()
     414        True
     415       
     416    ::
     417
     418        sage: A = random_matrix(GF(16007), 10, 10); A
     419        [ 5037  2388  4150  1400   345  5945  4240 14022 10514   700]
     420        [15552  8539  1927  3870  9867  3263 11637   609 15424  2443]
     421        [ 3761 15836 12246 15577 10178 13602 13183 15918 13942  2958]
     422        [ 4526 10817  6887  6678  1764  9964  6107  1705  5619  5811]
     423        [13537 15004  8307 11846 14779   550 14113  5477  7271  7091]
     424        [13338  4927 11406 13065  5437 12431  6318  5119 14198   496]
     425        [ 1044   179 12881   353 12975 12567  1092 10433 12304   954]
     426        [10072  8821 14118 13895  6543 13484 10685 14363  2612 11070]
     427        [15113   237  2612 14127 11589  5808   117  9656 15957 14118]
     428        [15233 11080  5716  9029 11402  9380 13045 13986 14544  5771]
     429
     430        sage: A.determinant()
     431        10207
     432       
     433    ::
     434
     435        sage: A = random_matrix(GF(16007), 100, 100)
     436        sage: A.determinant()
     437        3576
     438       
     439       
     440        sage: A.transpose().determinant()
     441        3576
     442
     443        sage: B = random_matrix(GF(16007), 100, 100)
     444        sage: B.determinant()
     445        4075
     446
     447        sage: (A*B).determinant() == A.determinant() * B.determinant()
     448        True
     449
     450    TESTS::
     451
     452        sage: A = random_matrix(GF(7), 0, 0); A.det()
     453        1
     454
     455        sage: A = random_matrix(GF(7), 0, 1); A.det()
     456        Traceback (most recent call last):
     457        ...
     458        ValueError: self must be a square matrix
     459
     460        sage: A = random_matrix(GF(7), 1, 0); A.det()
     461        Traceback (most recent call last):
     462        ...
     463        ValueError: self must be a square matrix
     464
     465        sage: A = matrix(GF(7), 5, 5); A.det()
     466        0
     467
     468    ::
     469
     470        sage: A = random_matrix(GF(16007), 0, 0); A.det()
     471        1
     472
     473        sage: A = random_matrix(GF(16007), 0, 1); A.det()
     474        Traceback (most recent call last):
     475        ...
     476        ValueError: self must be a square matrix
     477       
     478        sage: A = random_matrix(GF(16007), 1, 0); A.det()
     479        Traceback (most recent call last):
     480        ...
     481        ValueError: self must be a square matrix
     482
     483        sage: A = matrix(GF(16007), 5, 5); A.det()
     484        0
     485    """
     486    cdef ModField *F = new ModField(<long>modulus)
     487    cdef celement *cpy = linbox_copy(modulus, entries, nrows, ncols)
     488    if nrows*ncols > 1000: sig_on()
     489    d =  <celement>ModDet(F[0], nrows, ncols, <ModFieldElement*>cpy, ncols)
     490    if nrows*ncols > 1000: sig_off()
     491    sage_free(cpy)
     492    del F
     493    return d
     494
     495cdef inline int linbox_matrix_matrix_multiply(celement modulus, celement* ans, celement* A, celement* B, Py_ssize_t m, Py_ssize_t n, Py_ssize_t k):
     496    """
     497    C = A*B
     498
     499    EXAMPLES::
     500
     501        sage: A = random_matrix(GF(7),2,2); A
     502        [3 1]
     503        [6 6]
     504
     505        sage: B = random_matrix(GF(7),2,2); B
     506        [4 4]
     507        [2 2]
     508
     509        sage: A*B
     510        [0 0]
     511        [1 1]
     512
     513        sage: 3*A
     514        [2 3]
     515        [4 4]
     516
     517        sage: MS = parent(A)
     518        sage: MS(3) * A
     519        [2 3]
     520        [4 4]
     521
     522    ::
     523   
     524        sage: A = random_matrix(GF(17), 201, 117)
     525        sage: B = random_matrix(GF(17), 117, 195)
     526        sage: C = random_matrix(GF(17), 201, 117)
     527        sage: D = random_matrix(GF(17), 117, 195)
     528
     529        sage: E = (A+C)*(B+D)
     530
     531        sage: F = A*B + A*D + C*B + C*D
     532
     533        sage: E == F
     534        True
     535
     536    ::
     537
     538        sage: A = random_matrix(GF(17), 200, 200)
     539        sage: MS = parent(A)
     540        sage: (MS(0) * A) == 0
     541        True
     542
     543        sage: (MS(1) * A) == A
     544        True
     545
     546    ::
     547
     548        sage: A = random_matrix(Integers(8),2,2); A
     549        [7 2]
     550        [6 1]
     551
     552        sage: B = random_matrix(Integers(8),2,2); B
     553        [4 0]
     554        [5 6]
     555
     556        sage: A*B
     557        [6 4]
     558        [5 6]
     559
     560        sage: 3*A
     561        [5 6]
     562        [2 3]
     563
     564        sage: MS = parent(A)
     565        sage: MS(3) * A
     566        [5 6]
     567        [2 3]
     568
     569    ::
     570   
     571        sage: A = random_matrix(Integers(16), 201, 117)
     572        sage: B = random_matrix(Integers(16), 117, 195)
     573        sage: C = random_matrix(Integers(16), 201, 117)
     574        sage: D = random_matrix(Integers(16), 117, 195)
     575
     576        sage: E = (A+C)*(B+D)
     577
     578        sage: F = A*B + A*D + C*B + C*D
     579
     580        sage: E == F
     581        True
     582
     583    ::
     584
     585        sage: A = random_matrix(Integers(16), 200, 200)
     586        sage: MS = parent(A)
     587        sage: (MS(0) * A) == 0
     588        True
     589
     590        sage: (MS(1) * A) == A
     591        True
     592
     593    ::
     594
     595        sage: A = random_matrix(GF(16007),2,2); A
     596        [ 7856  5786]
     597        [10134 14607]
     598
     599        sage: B = random_matrix(GF(16007),2,2); B
     600        [10839  6194]
     601        [13327  5985]
     602
     603        sage: A*B
     604        [14254  4853]
     605        [ 8754 15217]
     606
     607        sage: 3*A
     608        [ 7561  1351]
     609        [14395 11807]
     610
     611        sage: MS = parent(A)
     612        sage: MS(3) * A
     613        [ 7561  1351]
     614        [14395 11807]
     615
     616    ::
     617   
     618        sage: A = random_matrix(GF(15991), 201, 117)
     619        sage: B = random_matrix(GF(15991), 117, 195)
     620        sage: C = random_matrix(GF(15991), 201, 117)
     621        sage: D = random_matrix(GF(15991), 117, 195)
     622
     623        sage: E = (A+C)*(B+D)
     624
     625        sage: F = A*B + A*D + C*B + C*D
     626
     627        sage: E == F
     628        True
     629
     630    ::
     631
     632        sage: A = random_matrix(GF(16007), 200, 200)
     633        sage: MS = parent(A)
     634        sage: (MS(0) * A) == 0
     635        True
     636
     637        sage: (MS(1) * A) == A
     638        True
     639
     640    ::
     641
     642        sage: A = random_matrix(Integers(1008),2,2); A
     643        [354 413]
     644        [307 499]
     645
     646        sage: B = random_matrix(Integers(1008),2,2); B
     647        [952  41]
     648        [973 851]
     649
     650        sage: A*B
     651        [1001   73]
     652        [ 623  772]
     653
     654        sage: 3*A
     655        [ 54 231]
     656        [921 489]
     657
     658        sage: MS = parent(A)
     659        sage: MS(3) * A
     660        [ 54 231]
     661        [921 489]
     662
     663    ::
     664   
     665        sage: A = random_matrix(Integers(1600), 201, 117)
     666        sage: B = random_matrix(Integers(1600), 117, 195)
     667        sage: C = random_matrix(Integers(1600), 201, 117)
     668        sage: D = random_matrix(Integers(1600), 117, 195)
     669
     670        sage: E = (A+C)*(B+D)
     671
     672        sage: F = A*B + A*D + C*B + C*D
     673
     674        sage: E == F
     675        True
     676    """
     677    cdef ModField *F = new ModField(<long>modulus)
     678    cdef ModFieldElement one, mone, zero
     679    F[0].init(one, <int>1)
     680    F[0].init(zero, <int>0)
     681    if m*n*k > 100000: sig_on()
     682    Mod_fgemm(F[0], fflas_no_trans, fflas_no_trans, m, n, k,
     683              one, <ModFieldElement*>A, k, <ModFieldElement*>B, n, zero,
     684              <ModFieldElement*>ans, n)
     685    if m*n*k > 100000: sig_off()
     686    del F
     687
     688cdef inline int linbox_matrix_vector_multiply(celement modulus, celement* C, celement* A, celement* b, Py_ssize_t m, Py_ssize_t n, fflas_trans_enum trans):
     689    """
     690    C = A*v
     691    """
     692    cdef ModField *F = new ModField(<long>modulus)
     693    cdef ModFieldElement one, mone, zero
     694    F.init(one, <int>1)
     695    F.init(zero, <int>0)
     696
     697    Mod_fgemv(F[0], trans,  m, n,
     698              one, <ModFieldElement*>A, n,
     699              <ModFieldElement*>b, 1,
     700              zero, <ModFieldElement*>C, 1)
     701    del F
     702
     703cdef inline linbox_minpoly(celement modulus, Py_ssize_t nrows, celement* entries):
     704    """
     705    Compute the minimal polynomial.
     706
     707    EXAMPLE::
     708
     709        sage: A = random_matrix(GF(17), 10, 10); A
     710        [ 2 14  0 15 11 10 16  2  9  4]
     711        [10 14  1 14  3 14 12 14  3 13]
     712        [10  1 14  6  2 14 13  7  6 14]
     713        [10  3  9 15  8  1  5  8 10 11]
     714        [ 5 12  4  9 15  2  6 11  2 12]
     715        [ 6 10 12  0  6  9  7  7  3  8]
     716        [ 2  9  1  5 12 13  7 16  7 11]
     717        [11  1  0  2  0  4  7  9  8 15]
     718        [ 5  3 16  2 11 10 12 14  0  7]
     719        [16  4  6  5  2  3 14 15 16  4]
     720
     721        sage: B = copy(A)
     722        sage: min_p = A.minimal_polynomial(proof=True); min_p
     723        x^10 + 13*x^9 + 10*x^8 + 9*x^7 + 10*x^6 + 4*x^5 + 10*x^4 + 10*x^3 + 12*x^2 + 14*x + 7
     724        sage: min_p(A) == 0
     725        True
     726        sage: B == A
     727        True
     728       
     729        sage: char_p = A.characteristic_polynomial(); char_p
     730        x^10 + 13*x^9 + 10*x^8 + 9*x^7 + 10*x^6 + 4*x^5 + 10*x^4 + 10*x^3 + 12*x^2 + 14*x + 7
     731        sage: min_p.divides(char_p)
     732        True
     733
     734    ::
     735
     736        sage: A = random_matrix(GF(1214471), 10, 10); A
     737        [ 266673  745841  418200  521668  905837  160562  831940   65852  173001  515930]
     738        [ 714380  778254  844537  584888  392730  502193  959391  614352  775603  240043]
     739        [1156372  104118 1175992  612032 1049083  660489 1066446  809624   15010 1002045]
     740        [ 470722  314480 1155149 1173111   14213 1190467 1079166  786442  429883  563611]
     741        [ 625490 1015074  888047 1090092  892387    4724  244901  696350  384684  254561]
     742        [ 898612   44844   83752 1091581  349242  130212  580087  253296  472569  913613]
     743        [ 919150   38603  710029  438461  736442  943501  792110  110470  850040  713428]
     744        [ 668799 1122064  325250 1084368  520553 1179743  791517   34060 1183757 1118938]
     745        [ 642169   47513   73428 1076788  216479  626571  105273  400489 1041378 1186801]
     746        [ 158611  888598 1138220 1089631   56266 1092400  890773 1060810  211135  719636]
     747
     748        sage: B = copy(A)
     749        sage: min_p = A.minimal_polynomial(proof=True); min_p
     750        x^10 + 283013*x^9 + 252503*x^8 + 512435*x^7 + 742964*x^6 + 130817*x^5 + 581471*x^4 + 899760*x^3 + 207023*x^2 + 470831*x + 381978
     751
     752        sage: min_p(A) == 0
     753        True
     754        sage: B == A
     755        True
     756       
     757        sage: char_p = A.characteristic_polynomial(); char_p
     758        x^10 + 283013*x^9 + 252503*x^8 + 512435*x^7 + 742964*x^6 + 130817*x^5 + 581471*x^4 + 899760*x^3 + 207023*x^2 + 470831*x + 381978
     759
     760        sage: min_p.divides(char_p)
     761        True       
     762
     763    TESTS::
     764
     765        sage: A = random_matrix(GF(17), 0, 0)
     766        sage: A.minimal_polynomial()
     767        1
     768
     769        sage: A = random_matrix(GF(17), 0, 1)
     770        sage: A.minimal_polynomial()
     771        Traceback (most recent call last):
     772        ...
     773        ValueError: matrix must be square
     774
     775        sage: A = random_matrix(GF(17), 1, 0)
     776        sage: A.minimal_polynomial()
     777        Traceback (most recent call last):
     778        ...
     779        ValueError: matrix must be square
     780
     781        sage: A = matrix(GF(17), 10, 10)
     782        sage: A.minimal_polynomial()
     783        x
     784
     785    ::
     786
     787        sage: A = random_matrix(GF(2535919), 0, 0)
     788        sage: A.minimal_polynomial()
     789        1
     790
     791        sage: A = random_matrix(GF(2535919), 0, 1)
     792        sage: A.minimal_polynomial()
     793        Traceback (most recent call last):
     794        ...
     795        ValueError: matrix must be square
     796
     797        sage: A = random_matrix(GF(2535919), 1, 0)
     798        sage: A.minimal_polynomial()
     799        Traceback (most recent call last):
     800        ...
     801        ValueError: matrix must be square
     802
     803        sage: A = matrix(GF(2535919), 10, 10)
     804        sage: A.minimal_polynomial()
     805        x
     806    """
     807    cdef Py_ssize_t i
     808    cdef ModField *F = new ModField(<long>modulus)
     809    cdef vector[ModFieldElement] *minP = new vector[ModFieldElement]()
     810    cdef ModFieldElement *X = <ModFieldElement*>sage_malloc(nrows*(nrows+1)*sizeof(ModFieldElement))
     811    cdef size_t *P = <size_t*>sage_malloc(nrows*sizeof(size_t))
     812
     813    cdef celement *cpy = linbox_copy(modulus, entries, nrows, nrows)
     814
     815    if nrows*nrows > 1000: sig_on()
     816    Mod_MinPoly(F[0], minP[0], nrows, <ModFieldElement*>cpy, nrows, X, nrows, P)
     817    if nrows*nrows > 1000: sig_off()
     818
     819    sage_free(cpy)
     820
     821    l = []
     822    for i in range(minP.size()):
     823        l.append( <celement>minP.at(i) )
     824
     825    sage_free(P)
     826    sage_free(X)
     827    del F
     828    return l
     829
     830cdef inline linbox_charpoly(celement modulus, Py_ssize_t nrows, celement* entries):
     831    """
     832    Compute the minimal polynomial.
     833
     834    EXAMPLE::
     835
     836        sage: A = random_matrix(GF(19), 10, 10); A
     837        [ 3  1  8 10  5 16 18  9  6  1]
     838        [ 5 14  4  4 14 15  5 11  3  0]
     839        [ 4  1  0  7 11  6 17  8  5  6]
     840        [ 4  6  9  4  8  1 18 17  8 18]
     841        [11  2  0  6 13  7  4 11 16 10]
     842        [12  6 12  3 15 10  5 11  3  8]
     843        [15  1 16  2 18 15 14  7  2 11]
     844        [16 16 17  7 14 12  7  7  0  5]
     845        [13 15  9  2 12 16  1 15 18  7]
     846        [10  8 16 18  9 18  2 13  5 10]
     847
     848        sage: B = copy(A)
     849        sage: min_p = A.minimal_polynomial(proof=True); min_p
     850        x^10 + 2*x^9 + 18*x^8 + 4*x^7 + 13*x^6 + 11*x^5 + 2*x^4 + 5*x^3 + 7*x^2 + 16*x + 6
     851       
     852        sage: min_p(A) == 0
     853        True
     854        sage: B == A
     855        True
     856       
     857        sage: char_p = A.characteristic_polynomial(); char_p
     858        x^10 + 2*x^9 + 18*x^8 + 4*x^7 + 13*x^6 + 11*x^5 + 2*x^4 + 5*x^3 + 7*x^2 + 16*x + 6
     859
     860        sage: min_p.divides(char_p)
     861        True
     862
     863    ::
     864
     865        sage: A = random_matrix(GF(2916337), 7, 7); A
     866        [ 446196 2267054   36722 2092388 1694559  514193 1196222]
     867        [1242955 1040744   99523 2447069   40527  930282 2685786]
     868        [2892660 1347146 1126775 2131459  869381 1853546 2266414]
     869        [2897342 1342067 1054026  373002   84731 1270068 2421818]
     870        [ 569466  537440  572533  297105 1415002 2079710  355705]
     871        [2546914 2299052 2883413 1558788 1494309 1027319 1572148]
     872        [ 250822  522367 2516720  585897 2296292 1797050 2128203]
     873
     874        sage: B = copy(A)
     875        sage: min_p = A.minimal_polynomial(proof=True); min_p
     876        x^7 + 1191770*x^6 + 547840*x^5 + 215639*x^4 + 2434512*x^3 + 1039968*x^2 + 483592*x + 733817
     877
     878        sage: min_p(A) == 0
     879        True
     880        sage: B == A
     881        True
     882       
     883        sage: char_p = A.characteristic_polynomial(); char_p
     884        x^7 + 1191770*x^6 + 547840*x^5 + 215639*x^4 + 2434512*x^3 + 1039968*x^2 + 483592*x + 733817
     885
     886        sage: min_p.divides(char_p)
     887        True       
     888
     889    TESTS::
     890
     891        sage: A = random_matrix(GF(19), 0, 0)
     892        sage: A.minimal_polynomial()
     893        1
     894
     895        sage: A = random_matrix(GF(19), 0, 1)
     896        sage: A.minimal_polynomial()
     897        Traceback (most recent call last):
     898        ...
     899        ValueError: matrix must be square
     900
     901        sage: A = random_matrix(GF(19), 1, 0)
     902        sage: A.minimal_polynomial()
     903        Traceback (most recent call last):
     904        ...
     905        ValueError: matrix must be square
     906
     907        sage: A = matrix(GF(19), 10, 10)
     908        sage: A.minimal_polynomial()
     909        x
     910
     911    ::
     912
     913        sage: A = random_matrix(GF(4198973), 0, 0)
     914        sage: A.minimal_polynomial()
     915        1
     916
     917        sage: A = random_matrix(GF(4198973), 0, 1)
     918        sage: A.minimal_polynomial()
     919        Traceback (most recent call last):
     920        ...
     921        ValueError: matrix must be square
     922
     923        sage: A = random_matrix(GF(4198973), 1, 0)
     924        sage: A.minimal_polynomial()
     925        Traceback (most recent call last):
     926        ...
     927        ValueError: matrix must be square
     928
     929        sage: A = matrix(GF(4198973), 10, 10)
     930        sage: A.minimal_polynomial()
     931        x
     932    """
     933    cdef Py_ssize_t i
     934    cdef ModField *F = new ModField(<long>modulus)
     935    cdef std_list[vector[ModFieldElement]] P_list
     936    P_list.clear()
     937
     938    cdef celement *cpy = linbox_copy(modulus, entries, nrows, nrows)
     939
     940    if nrows*nrows > 1000: sig_on()
     941    Mod_CharPoly(F[0], P_list, nrows, <ModFieldElement*>cpy, nrows)
     942    if nrows*nrows > 1000: sig_off()
     943
     944    sage_free(cpy)
     945
     946    cdef vector[ModFieldElement] tmp
     947    l = []
     948    while P_list.size():
     949        l.append([])
     950        tmp = P_list.front()
     951        for i in range(tmp.size()):
     952            l[-1].append(<celement>tmp.at(i))
     953        P_list.pop_front()
     954
     955    del F
     956    return l
     957
     958cdef class Matrix_modn_dense_template(matrix_dense.Matrix_dense):
     959    def __cinit__(self, parent, entries, copy, coerce):
     960        """
     961        Create a new matrix.
     962
     963        INPUT:
     964
     965        - ``parent`` - a matrix space
     966
     967        - ``entries`` - a list of entries or a scalar
     968
     969        - ``copy`` - ignroed
     970
     971        - ``coerce`` - perform modular reduction first?
     972
     973        EXAMPLES::
     974
     975            sage: A = random_matrix(GF(3),1000,1000)
     976            sage: type(A)
     977            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     978            sage: A = random_matrix(Integers(10),1000,1000)
     979            sage: type(A)
     980            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     981            sage: A = random_matrix(Integers(2^16),1000,1000)
     982            sage: type(A)
     983            <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>           
     984        """
     985        matrix_dense.Matrix_dense.__init__(self, parent)
     986
     987        cdef long p = self._base_ring.characteristic()
     988        self.p = p
     989        if p >= MAX_MODULUS:
     990            raise OverflowError("p (=%s) must be < %s."%(p, MAX_MODULUS))
     991
     992        sig_on()
     993        self._entries = <celement *> sage_malloc(sizeof(celement)*self._nrows*self._ncols)
     994        sig_off()
     995        if self._entries == NULL:
     996           raise MemoryError("Error allocating matrix.")
     997
     998        sig_on()
     999        self._matrix = <celement **> sage_malloc(sizeof(celement*)*self._nrows)
     1000        sig_off()
     1001        if self._matrix == NULL:
     1002            sage_free(self._entries)
     1003            self._entries = NULL
     1004            raise MemoryError("Error allocating memory.")
     1005
     1006        cdef unsigned int k
     1007        cdef Py_ssize_t i
     1008        k = 0
     1009        for i in range(self._nrows):
     1010            self._matrix[i] = self._entries + k
     1011            k = k + self._ncols
     1012
     1013    def __dealloc__(self):
     1014        """
     1015        TESTS::
     1016
     1017        sage: import gc
     1018        sage: for i in range(10):
     1019        ...      A = random_matrix(GF(7),1000,1000)
     1020        ...      B = random_matrix(Integers(10),1000,1000)
     1021        ...      C = random_matrix(GF(16007),1000,1000)
     1022        ...      D = random_matrix(Integers(1000),1000,1000)
     1023        ...      del A
     1024        ...      del B
     1025        ...      del C
     1026        ...      del D
     1027        ...      _ = gc.collect()
     1028
     1029        """
     1030        if self._entries == NULL:
     1031            return
     1032        sage_free(self._entries)
     1033        sage_free(self._matrix)
     1034
     1035    def __init__(self, parent, entries, copy, coerce):
     1036        """
     1037        Create a new matrix.
     1038
     1039        INPUT:
     1040
     1041        - ``parent`` - a matrix space
     1042
     1043        - ``entries`` - a list of entries or a scalar
     1044
     1045        - ``copy`` - ignroed
     1046
     1047        - ``coerce`` - perform modular reduction first?
     1048
     1049        EXAMPLES::
     1050
     1051            sage: A = random_matrix(GF(3),1000,1000)
     1052            sage: type(A)
     1053            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     1054            sage: A = random_matrix(Integers(10),1000,1000)
     1055            sage: type(A)
     1056            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     1057            sage: A = random_matrix(Integers(2^16),1000,1000)
     1058            sage: type(A)
     1059            <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>           
     1060
     1061        TESTS::
     1062
     1063            sage: Matrix(GF(7), 2, 2, [-1, int(-2), GF(7)(-3), 1/4])
     1064            [6 5]
     1065            [4 2]
     1066
     1067            sage: Matrix(GF(6434383), 2, 2, [-1, int(-2), GF(7)(-3), 1/4])
     1068            [6434382 6434381]
     1069            [      4 1608596]
     1070
     1071            sage: Matrix(Integers(4618990), 2, 2, [-1, int(-2), GF(7)(-3), 1/7])
     1072            [4618989 4618988]
     1073            [      4 2639423]
     1074        """
     1075        cdef celement e
     1076        cdef Py_ssize_t i, j, k
     1077        cdef celement *v
     1078        cdef long p
     1079        p = self._base_ring.characteristic()
     1080
     1081        R = self.base_ring()
     1082
     1083        # scalar?
     1084        if not isinstance(entries, list) and not isinstance(entries, tuple):
     1085            sig_on()
     1086            for i in range(self._nrows*self._ncols):
     1087                self._entries[i] = 0
     1088            sig_off()
     1089            if entries is None:
     1090                # zero matrix
     1091                pass
     1092            else:
     1093                e = R(entries)
     1094                if e != 0:
     1095                    for i in range(min(self._nrows, self._ncols)):
     1096                        self._matrix[i][i] = e
     1097            return
     1098
     1099        # all entries are given as a long list
     1100        if len(entries) != self._nrows * self._ncols:
     1101            raise IndexError("The vector of entries has the wrong length.")
     1102
     1103        k = 0
     1104        cdef celement n
     1105        cdef long tmp
     1106
     1107        for i in range(self._nrows):
     1108            sig_check()
     1109            v = self._matrix[i]
     1110            for j in range(self._ncols):
     1111                x = entries[k]
     1112                if PY_TYPE_CHECK_EXACT(x, int):
     1113                    tmp = (<long>x) % p
     1114                    v[j] = tmp + (tmp<0)*p
     1115                elif PY_TYPE_CHECK_EXACT(x, IntegerMod_int) and (<IntegerMod_int>x)._parent is R:
     1116                    v[j] = <celement>(<IntegerMod_int>x).ivalue
     1117                elif PY_TYPE_CHECK_EXACT(x, Integer):
     1118                    if coerce:
     1119                        v[j] = mpz_fdiv_ui((<Integer>x).value, p)
     1120                    else:
     1121                        v[j] = mpz_get_ui((<Integer>x).value)
     1122                elif coerce:
     1123                    v[j] = R(entries[k])
     1124                else:
     1125                    v[j] = <float>(entries[k])
     1126                k = k + 1
     1127
     1128    def __richcmp__(Matrix_modn_dense_template self, right, int op):  # always need for mysterious reasons.
     1129        return self._richcmp(right, op)
     1130
     1131    def __hash__(self):
     1132        """
     1133        EXAMPLE::
     1134
     1135            sage: B = random_matrix(GF(127),3,3)
     1136            sage: B.set_immutable()
     1137            sage: {B:0} # indirect doctest
     1138            {[  9  75  94]
     1139            [  4  57 112]
     1140            [ 59  85  45]: 0}
     1141
     1142        ::
     1143
     1144            sage: M = random_matrix(GF(7), 10, 10)
     1145            sage: M.set_immutable()
     1146            sage: hash(M)
     1147            143
     1148            sage: MZ = M.change_ring(ZZ)
     1149            sage: MZ.set_immutable()
     1150            sage: hash(MZ)
     1151            143
     1152            sage: MS = M.sparse_matrix()
     1153            sage: MS.set_immutable()
     1154            sage: hash(MS)
     1155            143
     1156       
     1157        TEST::
     1158       
     1159            sage: A = matrix(GF(2),2,0)
     1160            sage: hash(A)
     1161            Traceback (most recent call last):
     1162            ...
     1163            TypeError: mutable matrices are unhashable
     1164            sage: A.set_immutable()
     1165            sage: hash(A)
     1166            0
     1167        """
     1168        if self.is_mutable():
     1169            raise TypeError("Mutable matrices are unhashable.")
     1170        x = self.fetch('hash')
     1171        if not x is None:
     1172            return x
     1173
     1174        cdef long _hash = 0
     1175        cdef celement *_matrix
     1176        cdef long n = 0
     1177        cdef Py_ssize_t i, j
     1178
     1179        if self._nrows == 0 or self._ncols == 0:
     1180            return 0
     1181
     1182        sig_on()
     1183        for i in range(self._nrows):
     1184            _matrix = self._matrix[i]
     1185            for j in range(self._ncols):
     1186                _hash ^= <long>(n * _matrix[j])
     1187                n+=1
     1188        sig_off()
     1189
     1190        if _hash == -1:
     1191            return -2
     1192
     1193        self.cache('hash', _hash)
     1194
     1195        return _hash
     1196
     1197    def _pickle(self):
     1198        """
     1199        Utility function for pickling.
     1200
     1201        If the prime is small enough to fit in a byte, then it is
     1202        stored as a contiguous string of bytes (to save
     1203        space). Otherwise, memcpy is used to copy the raw data in the
     1204        platforms native format. Any byte-swapping or word size
     1205        difference is taken care of in unpickling (optimizing for
     1206        unpickling on the same platform things were pickled on).
     1207
     1208        The upcoming buffer protocol would be useful to not have to do
     1209        any copying.
     1210
     1211        EXAMPLES::
     1212
     1213            sage: m = matrix(Integers(128), 3, 3, [ord(c) for c in "Hi there!"]); m
     1214            [ 72 105  32]
     1215            [116 104 101]
     1216            [114 101  33]
     1217            sage: m._pickle()
     1218            ((1, ..., 'Hi there!'), 10)
     1219        """
     1220        cdef Py_ssize_t i, j
     1221        cdef unsigned char* us
     1222        cdef mod_int *um
     1223        cdef unsigned char* row_us
     1224        cdef mod_int *row_um
     1225        cdef long word_size
     1226        cdef celement *row_self
     1227
     1228        if self.p <= 0xFF:
     1229            word_size = sizeof(char)
     1230        else:
     1231            word_size = sizeof(mod_int)
     1232
     1233        s = PyString_FromStringAndSize(NULL, word_size * self._nrows * self._ncols)
     1234
     1235        sig_on()
     1236        if word_size == sizeof(char):
     1237            us = <unsigned char*><char *>s
     1238            for i in range(self._nrows):
     1239                row_self = self._matrix[i]
     1240                row_us = &us[i*self._ncols]
     1241                for j in range(self._ncols):
     1242                    row_us[j] = <unsigned char><mod_int>row_self[j]
     1243        else:
     1244            um = <mod_int*><char *>s
     1245            for i in range(self._nrows):
     1246                row_self = self._matrix[i]
     1247                row_um = &um[i*self._ncols]
     1248                for j in range(self._ncols):
     1249                    row_um[j] = <mod_int>row_self[j]
     1250        sig_off()
     1251
     1252        return (word_size, little_endian, s), 10
     1253
     1254
     1255    def _unpickle(self, data, int version):
     1256        """
     1257        TESTS:
     1258
     1259        Test for char-sized modulus::
     1260
     1261            sage: A = random_matrix(GF(7), 5, 9)
     1262            sage: data, version = A._pickle()
     1263            sage: B = A.parent()(0)
     1264            sage: B._unpickle(data, version)
     1265            sage: B == A
     1266            True
     1267
     1268        And for larger modulus::
     1269
     1270            sage: A = random_matrix(GF(1009), 51, 5)
     1271            sage: data, version = A._pickle()
     1272            sage: B = A.parent()(0)
     1273            sage: B._unpickle(data, version)
     1274            sage: B == A
     1275            True
     1276
     1277        Now test all the bit-packing options::
     1278
     1279            sage: A = matrix(Integers(1000), 2, 2)
     1280            sage: A._unpickle((1, True, '\x01\x02\xFF\x00'), 10)
     1281            sage: A
     1282            [  1   2]
     1283            [255   0]
     1284
     1285        ::
     1286
     1287            sage: A = matrix(Integers(1000), 1, 2)
     1288            sage: A._unpickle((4, True, '\x02\x01\x00\x00\x01\x00\x00\x00'), 10)
     1289            sage: A
     1290            [258   1]
     1291            sage: A._unpickle((4, False, '\x00\x00\x02\x01\x00\x00\x01\x03'), 10)
     1292            sage: A
     1293            [513 259]
     1294            sage: A._unpickle((8, True, '\x03\x01\x00\x00\x00\x00\x00\x00\x05\x00\x00\x00\x00\x00\x00\x00'), 10)
     1295            sage: A
     1296            [259   5]
     1297            sage: A._unpickle((8, False, '\x00\x00\x00\x00\x00\x00\x02\x08\x00\x00\x00\x00\x00\x00\x01\x04'), 10)
     1298            sage: A
     1299            [520 260]
     1300
     1301        Now make sure it works in context::
     1302
     1303            sage: A = random_matrix(Integers(33), 31, 31)
     1304            sage: loads(dumps(A)) == A
     1305            True
     1306            sage: A = random_matrix(Integers(3333), 31, 31)
     1307            sage: loads(dumps(A)) == A
     1308            True
     1309        """
     1310        if version < 10:
     1311            return matrix_dense.Matrix_dense._unpickle(self, data, version)
     1312
     1313        cdef Py_ssize_t i, j
     1314        cdef unsigned char* us
     1315        cdef mod_int *um
     1316        cdef unsigned char* row_us
     1317        cdef mod_int *row_um
     1318        cdef long word_size
     1319        cdef celement *row_self
     1320        cdef bint little_endian_data
     1321        cdef unsigned char* udata
     1322
     1323        if version == 10:
     1324            word_size, little_endian_data, s = data
     1325
     1326            sig_on()
     1327            if word_size == 1:
     1328                us = <unsigned char*><char *>s
     1329                for i from 0 <= i < self._nrows:
     1330                    row_self = self._matrix[i]
     1331                    row_us = &us[i*self._ncols]
     1332                    for j from 0<= j < self._ncols:
     1333                        row_self[j] = <celement>row_us[j]
     1334
     1335            elif word_size == sizeof(mod_int) and little_endian == little_endian_data:
     1336                um = <mod_int*><char *>s
     1337                for i from 0 <= i < self._nrows:
     1338                    row_self = self._matrix[i]
     1339                    row_um = &um[i*self._ncols]
     1340                    for j from 0<= j < self._ncols:
     1341                        row_self[j] =  <celement>row_um[j]
     1342                   
     1343            elif little_endian_data:
     1344                us = <unsigned char*><char *>s
     1345                for i from 0 <= i < self._nrows:
     1346                    row_self = self._matrix[i]
     1347                    for j from 0<= j < self._ncols:
     1348                        udata = &us[(i*self._ncols+j)*word_size]
     1349                        row_self[j] =  <celement>((<mod_int>udata[0]) + ((<mod_int>udata[1]) << 8) + ((<mod_int>udata[2]) << 16) + ((<mod_int>udata[3]) << 24))
     1350            else:
     1351                us = <unsigned char*><char *>s
     1352                for i from 0 <= i < self._nrows:
     1353                    row_self = self._matrix[i]
     1354                    for j from 0<= j < self._ncols:
     1355                        udata = &us[(i*self._ncols+j)*word_size]
     1356                        row_self[j] =  <celement>((<mod_int>udata[word_size-1]) + ((<mod_int>udata[word_size-2]) << 8) +((<mod_int>udata[word_size-3]) << 16) + ((<mod_int>udata[word_size-4]) << 24))
     1357            sig_off()
     1358        else:
     1359            raise ValueError("Unknown matrix version.")
     1360
     1361    def __neg__(self):
     1362        """
     1363        EXAMPLES::
     1364
     1365            sage: m = matrix(GF(19), 3, 3, range(9)); m
     1366            [0 1 2]
     1367            [3 4 5]
     1368            [6 7 8]
     1369            sage: -m
     1370            [ 0 18 17]
     1371            [16 15 14]
     1372            [13 12 11]
     1373        """
     1374        cdef Py_ssize_t i, j
     1375        cdef Matrix_modn_dense_template M
     1376        cdef celement p = self.p
     1377
     1378        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
     1379
     1380        sig_on()
     1381        for i in range(self._nrows*self._ncols):
     1382            if self._entries[i]:
     1383                M._entries[i] = p - self._entries[i]
     1384            else:
     1385                M._entries[i] = 0
     1386        sig_off()
     1387        return M
     1388
     1389
     1390    cpdef ModuleElement _lmul_(self, RingElement right):
     1391        """
     1392        EXAMPLES::
     1393       
     1394            sage: a = random_matrix(Integers(60), 400, 500)
     1395            sage: 3*a + 9*a == 12*a
     1396            True
     1397        """
     1398        return self._rmul_(right)
     1399       
     1400    cpdef ModuleElement _rmul_(self, RingElement left):
     1401        """
     1402        EXAMPLES::
     1403       
     1404            sage: a = matrix(GF(101), 3, 3, range(9)); a
     1405            [0 1 2]
     1406            [3 4 5]
     1407            [6 7 8]
     1408            sage: a * 5
     1409            [ 0  5 10]
     1410            [15 20 25]
     1411            [30 35 40]
     1412            sage: a * 50
     1413            [  0  50 100]
     1414            [ 49  99  48]
     1415            [ 98  47  97]
     1416        """
     1417        cdef Py_ssize_t i,j
     1418        cdef Matrix_modn_dense_template M
     1419        cdef celement p = self.p
     1420        cdef celement a = left
     1421
     1422        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
     1423
     1424        sig_on()
     1425        for i in range(self._nrows*self._ncols):
     1426            M._entries[i] = (a*self._entries[i]) % p
     1427        sig_off()
     1428        return M
     1429
     1430    def __copy__(self):
     1431        cdef Matrix_modn_dense_template A
     1432        A = self.__class__.__new__(self.__class__, self._parent, 0, 0, 0)
     1433        memcpy(A._entries, self._entries, sizeof(celement)*self._nrows*self._ncols)
     1434        if self._subdivisions is not None:
     1435            A.subdivide(*self.subdivisions())
     1436        return A
     1437
     1438
     1439    cpdef ModuleElement _add_(self, ModuleElement right):
     1440        """
     1441        Add two dense matrices over Z/nZ
     1442       
     1443        EXAMPLES::
     1444       
     1445            sage: a = MatrixSpace(GF(19),3)(range(9))
     1446            sage: a+a
     1447            [ 0  2  4]
     1448            [ 6  8 10]
     1449            [12 14 16]
     1450            sage: b = MatrixSpace(GF(19),3)(range(9))
     1451            sage: b.swap_rows(1,2)
     1452            sage: a+b
     1453            [ 0  2  4]
     1454            [ 9 11 13]
     1455            [ 9 11 13]
     1456            sage: b+a
     1457            [ 0  2  4]
     1458            [ 9 11 13]
     1459            [ 9 11 13]
     1460        """
     1461        cdef Py_ssize_t i
     1462        cdef celement k, p
     1463        cdef Matrix_modn_dense_template M
     1464
     1465        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
     1466        p = self.p
     1467        cdef celement* other_ent = (<Matrix_modn_dense_template>right)._entries
     1468
     1469        sig_on()
     1470        for i in range(self._nrows*self._ncols):
     1471            k = self._entries[i] + other_ent[i]
     1472            M._entries[i] = k - (k >= p) * p
     1473        sig_off()
     1474        return M
     1475
     1476
     1477    cpdef ModuleElement _sub_(self, ModuleElement right):
     1478        """
     1479        Subtract two dense matrices over Z/nZ
     1480       
     1481        EXAMPLES::
     1482       
     1483            sage: a = matrix(GF(11), 3, 3, range(9)); a
     1484            [0 1 2]
     1485            [3 4 5]
     1486            [6 7 8]
     1487            sage: a - 4
     1488            [7 1 2]
     1489            [3 0 5]
     1490            [6 7 4]
     1491            sage: a - matrix(GF(11), 3, 3, range(1, 19, 2))
     1492            [10  9  8]
     1493            [ 7  6  5]
     1494            [ 4  3  2]
     1495        """
     1496        cdef Py_ssize_t i
     1497        cdef celement k, p
     1498        cdef Matrix_modn_dense_template M
     1499
     1500        M = self.__class__.__new__(self.__class__, self._parent, None, None, None)
     1501        p = self.p
     1502        cdef celement* other_ent = (<Matrix_modn_dense_template>right)._entries
     1503
     1504        sig_on()
     1505        for i in range(self._nrows*self._ncols):
     1506            k = p + self._entries[i] - other_ent[i]
     1507            M._entries[i] = k - (k >= p) * p
     1508        sig_off()
     1509        return M
     1510
     1511
     1512    cdef int _cmp_c_impl(self, Element right) except -2:
     1513        """
     1514        Compare two dense matrices over Z/nZ
     1515       
     1516        EXAMPLES::
     1517       
     1518            sage: a = matrix(GF(17), 4, range(3, 83, 5)); a
     1519            [ 3  8 13  1]
     1520            [ 6 11 16  4]
     1521            [ 9 14  2  7]
     1522            [12  0  5 10]
     1523            sage: a == a
     1524            True
     1525            sage: b = a - 3; b
     1526            [ 0  8 13  1]
     1527            [ 6  8 16  4]
     1528            [ 9 14 16  7]
     1529            [12  0  5  7]
     1530            sage: b < a
     1531            True
     1532            sage: b > a
     1533            False
     1534            sage: b == a
     1535            False
     1536            sage: b + 3 == a
     1537            True
     1538        """
     1539        cdef Py_ssize_t i
     1540        cdef celement* other_ent = (<Matrix_modn_dense_template>right)._entries
     1541        sig_on()
     1542        for i in range(self._nrows*self._ncols):
     1543            if self._entries[i] < other_ent[i]:
     1544                sig_off()
     1545                return -1
     1546            elif self._entries[i] > other_ent[i]:
     1547                sig_off()
     1548                return 1
     1549        sig_off()
     1550        return 0
     1551
     1552
     1553    cdef Matrix _matrix_times_matrix_(self, Matrix right):
     1554        """
     1555        Multiply matrices using LinBox.
     1556
     1557        INPUT:
     1558
     1559        - ``right``- Matrix
     1560
     1561        """
     1562        if get_verbose() >= 2:
     1563            verbose('mod-p multiply of %s x %s matrix by %s x %s matrix modulo %s'%(
     1564                    self._nrows, self._ncols, right._nrows, right._ncols, self.p))
     1565
     1566        cdef int e
     1567        cdef Matrix_modn_dense_template ans, B
     1568       
     1569        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     1570
     1571        B = right
     1572
     1573        linbox_matrix_matrix_multiply(self.p, ans._entries, self._entries,
     1574                                      B._entries, self._nrows, B._ncols, B._nrows)
     1575       
     1576        return ans
     1577
     1578    cdef Vector _vector_times_matrix_(self, Vector v):
     1579        """
     1580        EXAMPLES::
     1581
     1582            sage: A = random_matrix(GF(17), 10, 20)
     1583            sage: v = random_vector(GF(17), 10)
     1584            sage: matrix(v*A) == matrix(v)*A
     1585            True
     1586
     1587            sage: A = random_matrix(Integers(126), 10, 20)
     1588            sage: v = random_vector(Integers(126), 10)
     1589            sage: matrix(v*A) == matrix(v)*A
     1590            True
     1591
     1592            sage: A = random_matrix(GF(4796509), 10, 20)
     1593            sage: v = random_vector(GF(4796509), 10)
     1594            sage: matrix(v*A) == matrix(v)*A
     1595            True
     1596
     1597            sage: A = random_matrix(Integers(16337), 10, 20)
     1598            sage: v = random_vector(Integers(16337), 10)
     1599            sage: matrix(v*A) == matrix(v)*A
     1600            True
     1601        """
     1602        if not isinstance(v, Vector_modn_dense):
     1603            return (self.new_matrix(1,self._nrows, entries=v.list()) * self)[0]
     1604
     1605        cdef Py_ssize_t i
     1606        cdef Vector_modn_dense b = v
     1607
     1608        cdef celement *_b = <celement*>sage_malloc(sizeof(celement)*self._nrows)
     1609        cdef celement *_c = <celement*>sage_malloc(sizeof(celement)*self._ncols)
     1610
     1611        for i in range(self._nrows):
     1612            _b[i] = <celement>b._entries[i]
     1613
     1614        linbox_matrix_vector_multiply(self.p, _c, self._entries, _b, self._nrows, self._ncols, fflas_trans)
     1615
     1616        M = self._row_ambient_module()
     1617        cdef Vector_modn_dense c = M.zero_vector()
     1618        for i in range(self._ncols):
     1619            c._entries[i] = <mod_int>_c[i]
     1620        sage_free(_b)
     1621        sage_free(_c)
     1622        return c
     1623
     1624    cdef Vector _matrix_times_vector_(self, Vector v):
     1625        if not isinstance(v, Vector_modn_dense):
     1626            return matrix_dense.Matrix_dense._matrix_times_vector(self, v)
     1627
     1628        cdef Py_ssize_t i
     1629        cdef Vector_modn_dense b = v
     1630
     1631        cdef celement *_b = <celement*>sage_malloc(sizeof(celement)*self._ncols)
     1632        cdef celement *_c = <celement*>sage_malloc(sizeof(celement)*self._nrows)
     1633
     1634        for i in range(self._ncols):
     1635            _b[i] = <celement>b._entries[i]
     1636
     1637        linbox_matrix_vector_multiply(self.p, _c, self._entries, _b, self._nrows, self._ncols, fflas_no_trans)
     1638
     1639        M = self._column_ambient_module()
     1640        cdef Vector_modn_dense c = M.zero_vector()
     1641        for i in range(self._nrows):
     1642            c._entries[i] = <mod_int>_c[i]
     1643        sage_free(_b)
     1644        sage_free(_c)
     1645        return c
     1646   
     1647
     1648    ########################################################################
     1649    # LEVEL 3 functionality (Optional)
     1650    #  x * cdef _sub_
     1651    #    * __deepcopy__
     1652    #    * __invert__
     1653    #    * Matrix windows -- only if you need strassen for that base
     1654    #    * Other functions (list them here):
     1655    #        - all row/column operations, but optimized
     1656    # x      - echelon form in place
     1657    #        - Hessenberg forms of matrices
     1658    ########################################################################
     1659
     1660
     1661    def charpoly(self, var='x', algorithm='linbox'):
     1662        """
     1663        Returns the characteristic polynomial of self.
     1664       
     1665        INPUT:
     1666       
     1667       
     1668        -  ``var`` - a variable name
     1669       
     1670        -  ``algorithm`` - 'generic', 'linbox' or 'all' (default)
     1671       
     1672       
     1673        EXAMPLES::
     1674       
     1675            sage: A = Mat(GF(7),3,3)(range(3)*3)
     1676            sage: A.charpoly()
     1677            x^3 + 4*x^2
     1678       
     1679        ::
     1680       
     1681            sage: A = Mat(Integers(6),3,3)(range(9))
     1682            sage: A.charpoly()
     1683            x^3
     1684       
     1685        ALGORITHM: Uses LinBox if self.base_ring() is a field, otherwise
     1686        use Hessenberg form algorithm.
     1687        """
     1688        if algorithm == 'linbox' and (self.p == 2 or not self.base_ring().is_field()):
     1689            algorithm = 'generic' # LinBox only supports Z/pZ (p prime)
     1690       
     1691        if algorithm == 'linbox':
     1692            g = self._charpoly_linbox(var)
     1693        elif algorithm == 'generic':
     1694            g = matrix_dense.Matrix_dense.charpoly(self, var)
     1695        elif algorithm == 'all':
     1696            g = self._charpoly_linbox(var)
     1697            h = matrix_dense.Matrix_dense.charpoly(self, var)
     1698            if g != h:
     1699                raise ArithmeticError("Characteristic polynomials do not match.")
     1700        else:
     1701            raise ValueError, "no algorithm '%s'"%algorithm
     1702        self.cache('charpoly_%s_%s'%(algorithm, var), g)
     1703        return g
     1704
     1705    def minpoly(self, var='x', algorithm='linbox', proof=None):
     1706        """
     1707        Returns the minimal polynomial of`` self``.
     1708
     1709        INPUT:
     1710
     1711        - ``var`` - a variable name
     1712           
     1713        - ``algorithm`` - ``generic`` or ``linbox`` (default: ``linbox``)
     1714
     1715        - ``proof`` -- (default: ``True``); whether to provably return
     1716          the true minimal polynomial; if ``False``, we only guarantee
     1717          to return a divisor of the minimal polynomial.  There are
     1718          also certainly cases where the computed results is
     1719          frequently not exactly equal to the minimal polynomial (but
     1720          is instead merely a divisor of it).
     1721
     1722         .. warning::
     1723
     1724             If ``proof=True``, minpoly is insanely slow compared to
     1725             ``proof=False``. This matters since proof=True is the
     1726             default, unless you first type
     1727             ``proof.linear_algebra(False)``.
     1728
     1729        EXAMPLES::
     1730
     1731            sage: R.<x>=GF(3)[]
     1732            sage: A = matrix(GF(3),2,[0,0,1,2])
     1733            sage: A.minpoly()
     1734            x^2 + x
     1735
     1736            sage: A.minpoly(proof=False) in [x, x+1, x^2+x]
     1737            True
     1738        """
     1739        proof = get_proof_flag(proof, "linear_algebra")
     1740
     1741        if algorithm == 'linbox' and (self.p == 2 or not self.base_ring().is_field()):       
     1742            algorithm='generic' # LinBox only supports fields
     1743
     1744        if algorithm == 'linbox':
     1745            if self._nrows != self._ncols:
     1746                raise ValueError("matrix must be square")
     1747
     1748            if self._nrows <= 1:
     1749                return matrix_dense.Matrix_dense.minpoly(self, var)
     1750
     1751            R = self._base_ring[var]
     1752            v = linbox_minpoly(self.p, self._nrows, self._entries)
     1753            g = R(v)
     1754
     1755            if proof == True:
     1756                while g(self):  # insanely toy slow (!)
     1757                    g = g.lcm(R(linbox_minpoly(self.p, self._nrows, self._entries)))
     1758
     1759        elif algorithm == 'generic':
     1760            raise NotImplementedError("Minimal polynomials are not implemented for Z/nZ.")
     1761
     1762        else:
     1763            raise ValueError("no algorithm '%s'"%algorithm)
     1764
     1765        self.cache('minpoly_%s_%s'%(algorithm, var), g)
     1766        return g
     1767
     1768    def _charpoly_linbox(self, var='x'):
     1769        """
     1770        Computes the characteristic polynomial using LinBox. No checks
     1771        are performed.
     1772        """
     1773        verbose('_charpoly_linbox...')
     1774       
     1775        if self._nrows != self._ncols:
     1776            raise ValueError, "matrix must be square"
     1777        if self._nrows <= 1:
     1778            return matrix_dense.Matrix_dense.charpoly(self, var)
     1779        R = self._base_ring[var]
     1780        # call linbox for charpoly
     1781        v = linbox_charpoly(self.p, self._nrows, self._entries)
     1782        r = R(1)
     1783        for e in v:
     1784            r *= R(e)
     1785        return r
     1786
     1787    def echelonize(self, algorithm="linbox", **kwds):
     1788        """
     1789        Puts self in row echelon form.
     1790
     1791        INPUT:
     1792               
     1793        - ``self`` - a mutable matrix
     1794        - ``algorithm`` - 'linbox' - uses the C++ linbox library
     1795        - ``'gauss'`` - uses a custom slower `O(n^3)` Gauss
     1796           elimination implemented in Sage.
     1797        - ``'all'`` - compute using both algorithms and verify that
     1798           the results are the same (for the paranoid).
     1799        - ``**kwds`` - these are all ignored
     1800               
     1801        OUTPUT:
     1802               
     1803        - ``self`` is put in reduced row echelon form.
     1804        - the rank of self is computed and cached
     1805        - the pivot columns of self are computed and cached.
     1806        - the fact that self is now in echelon form is recorded and
     1807          cached so future calls to echelonize return immediately.
     1808               
     1809        EXAMPLES::
     1810       
     1811            sage: A = matrix(GF(97),3,4,range(12))
     1812            sage: A.echelonize(); A
     1813            [ 1  0 96 95]
     1814            [ 0  1  2  3]
     1815            [ 0  0  0  0]
     1816            sage: A.pivots()
     1817            (0, 1)
     1818
     1819            sage: for p in (3,17,97,127,1048573):
     1820            ...      for i in range(10):
     1821            ...          A = random_matrix(GF(3), 100, 100)
     1822            ...          A.echelonize(algorithm='all')
     1823        """
     1824        x = self.fetch('in_echelon_form')
     1825        if not x is None:
     1826            return  # already known to be in echelon form
     1827
     1828        if not self.base_ring().is_field():
     1829            raise NotImplementedError("Echelon form not implemented over '%s'."%self.base_ring())
     1830
     1831        if algorithm == 'linbox':
     1832            self._echelonize_linbox()
     1833
     1834        elif algorithm == 'gauss':
     1835            self._echelon_in_place_classical()
     1836
     1837        elif algorithm == 'all':
     1838            A = self.__copy__()
     1839            self._echelonize_linbox()
     1840            A._echelon_in_place_classical()
     1841            if A != self:
     1842                raise ArithmeticError("Bug in echelon form.")
     1843        else:
     1844            raise ValueError("Algorithm '%s' not known"%algorithm)
     1845
     1846    def _echelonize_linbox(self):
     1847        """
     1848        Puts self in row echelon form using LinBox.
     1849        """
     1850        self.check_mutability()
     1851        self.clear_cache()
     1852
     1853        t = verbose('Calling echelonize mod %d.'%self.p)
     1854        r, pivots = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols)
     1855        verbose('done with echelonize',t)
     1856
     1857        self.cache('in_echelon_form',True)
     1858        self.cache('rank', r)
     1859        self.cache('pivots', tuple(pivots))
     1860
     1861    def _echelon_in_place_classical(self):
     1862        self.check_mutability()
     1863        self.clear_cache()
     1864
     1865        cdef Py_ssize_t start_row, c, r, nr, nc, i
     1866        cdef celement p, a, s, t, b
     1867        cdef celement **m
     1868
     1869        start_row = 0
     1870        p = self.p
     1871        m = self._matrix
     1872        nr = self._nrows
     1873        nc = self._ncols
     1874        pivots = []
     1875        fifth = self._ncols / 10 + 1
     1876        do_verb = (get_verbose() >= 2)
     1877        for c from 0 <= c < nc:
     1878            if do_verb and (c % fifth == 0 and c>0):
     1879                tm = verbose('on column %s of %s'%(c, self._ncols),
     1880                             level = 2,
     1881                             caller_name = 'matrix_modn_dense echelon')
     1882            #end if
     1883            sig_check()
     1884            for r from start_row <= r < nr:
     1885                a = m[r][c]
     1886                if a:
     1887                    pivots.append(c)
     1888                    a_inverse = celement_invert(a, p)
     1889                    self.rescale_row_c(r, a_inverse, c)
     1890                    self.swap_rows_c(r, start_row)
     1891                    for i from 0 <= i < nr:
     1892                        if i != start_row:
     1893                            b = m[i][c]
     1894                            if b != 0:
     1895                                self.add_multiple_of_row_c(i, start_row, p-b, c)
     1896                    start_row = start_row + 1
     1897                    break
     1898        self.cache('pivots', tuple(pivots))
     1899        self.cache('in_echelon_form',True)
     1900
     1901    cdef xgcd_eliminate(self, celement * row1, celement* row2, Py_ssize_t start_col):
     1902        """
     1903        Reduces row1 and row2 by a unimodular transformation using the xgcd
     1904        relation between their first coefficients a and b.
     1905
     1906        INPUT:
     1907       
     1908        - ``row1, row2`` - the two rows to be transformed (within
     1909          self)
     1910       
     1911        -``start_col`` - the column of the pivots in ``row1`` and
     1912         ``row2``. It is assumed that all entries before ``start_col``
     1913         in ``row1`` and ``row2`` are zero.
     1914       
     1915       
     1916        OUTPUT:
     1917               
     1918        - g: the gcd of the first elements of row1 and
     1919          row2 at column start_col
     1920
     1921        - put row1 = s \* row1 + t \* row2 row2 = w \*
     1922          row1 + t \* row2 where g = sa + tb
     1923        """
     1924        cdef int p = <int>self.p
     1925        cdef celement * row1_p, * row2_p
     1926        cdef celement tmp
     1927        cdef int g, s, t, v, w
     1928        cdef Py_ssize_t nc, i
     1929        cdef int a = <int>row1[start_col]
     1930        cdef int b = <int>row2[start_col]
     1931        g = ArithIntObj.c_xgcd_int (a,b,<int*>&s,<int*>&t)
     1932        v = a/g
     1933        w = -<int>b/g
     1934        nc = self.ncols()
     1935
     1936        for i from start_col <= i < nc:
     1937            tmp = ( s * <int>row1[i] + t * <int>row2[i]) % p
     1938            row2[i] = (w* <int>row1[i] + v*<int>row2[i]) % p
     1939            row1[i] = tmp
     1940        return g
     1941
     1942    cdef rescale_row_c(self, Py_ssize_t row, multiple, Py_ssize_t start_col):
     1943        cdef celement r, p
     1944        cdef celement* v
     1945        cdef Py_ssize_t i
     1946        p = self.p
     1947        v = self._matrix[row]
     1948        for i from start_col <= i < self._ncols:
     1949            v[i] = (v[i]*<celement>multiple) % p
     1950
     1951    cdef rescale_col_c(self, Py_ssize_t col, multiple, Py_ssize_t start_row):
     1952        """
     1953        EXAMPLES::
     1954
     1955            sage: n=3; b = MatrixSpace(Integers(37),n,n,sparse=False)([1]*n*n)
     1956            sage: b
     1957            [1 1 1]
     1958            [1 1 1]
     1959            [1 1 1]
     1960            sage: b.rescale_col(1,5)
     1961            sage: b
     1962            [1 5 1]
     1963            [1 5 1]
     1964            [1 5 1]
     1965
     1966        Recaling need not include the entire row.::
     1967
     1968            sage: b.rescale_col(0,2,1); b
     1969            [1 5 1]
     1970            [2 5 1]
     1971            [2 5 1]           
     1972       
     1973        Bounds are checked.::
     1974       
     1975            sage: b.rescale_col(3,2)
     1976            Traceback (most recent call last):
     1977            ...
     1978            IndexError: matrix column index out of range
     1979       
     1980        Rescaling by a negative number::
     1981
     1982            sage: b.rescale_col(2,-3); b
     1983            [ 1  5 34]
     1984            [ 2  5 34]
     1985            [ 2  5 34]
     1986        """
     1987        cdef celement r, p
     1988        cdef celement* v
     1989        cdef Py_ssize_t i
     1990        p = self.p
     1991        for i from start_row <= i < self._nrows:
     1992            self._matrix[i][col] = (self._matrix[i][col]*<celement>multiple) % p
     1993
     1994    cdef add_multiple_of_row_c(self, Py_ssize_t row_to, Py_ssize_t row_from, multiple, Py_ssize_t start_col):
     1995        cdef celement p
     1996        cdef celement *v_from, *v_to
     1997
     1998        p = self.p
     1999        v_from = self._matrix[row_from]
     2000        v_to = self._matrix[row_to]
     2001
     2002        cdef Py_ssize_t i, nc
     2003        nc = self._ncols
     2004        for i from start_col <= i < nc:
     2005            v_to[i] = ((<celement>multiple) * v_from[i] +  v_to[i]) % p
     2006
     2007    cdef add_multiple_of_column_c(self, Py_ssize_t col_to, Py_ssize_t col_from, multiple, Py_ssize_t start_row):
     2008        cdef celement  p
     2009        cdef celement **m
     2010
     2011        m = self._matrix
     2012        p = self.p
     2013
     2014        cdef Py_ssize_t i, nr
     2015        nr = self._nrows
     2016        for i from start_row <= i < self._nrows:
     2017            m[i][col_to] = (m[i][col_to] + (<celement>multiple) * m[i][col_from]) %p
     2018
     2019    cdef swap_rows_c(self, Py_ssize_t row1, Py_ssize_t row2):
     2020        """
     2021        EXAMPLES::
     2022       
     2023            sage: A = matrix(Integers(8), 2,[1,2,3,4])
     2024            sage: A.swap_rows(0,1)
     2025            sage: A
     2026            [3 4]
     2027            [1 2]
     2028        """
     2029        cdef celement* r1 = self._matrix[row1]
     2030        cdef celement* r2 = self._matrix[row2]
     2031        cdef celement temp
     2032        for i in range(self._ncols):
     2033            temp = r1[i]
     2034            r1[i] = r2[i]
     2035            r2[i] = temp
     2036
     2037    cdef swap_columns_c(self, Py_ssize_t col1, Py_ssize_t col2):
     2038        cdef Py_ssize_t i, nr
     2039        cdef celement t
     2040        cdef celement **m
     2041        m = self._matrix
     2042        nr = self._nrows
     2043        for i from 0 <= i < nr:
     2044            t = m[i][col1]
     2045            m[i][col1] = m[i][col2]
     2046            m[i][col2] = t
     2047           
     2048    def hessenbergize(self):
     2049        """
     2050        Transforms self in place to its Hessenberg form.
     2051        """
     2052        self.check_mutability()
     2053        x = self.fetch('in_hessenberg_form')
     2054        if not x is None and x: return  # already known to be in Hessenberg form
     2055       
     2056        if self._nrows != self._ncols:
     2057            raise ArithmeticError, "Matrix must be square to compute Hessenberg form."
     2058
     2059        cdef Py_ssize_t n
     2060        n = self._nrows
     2061
     2062        cdef celement **h
     2063        h = self._matrix
     2064
     2065        cdef celement p, t, t_inv, u
     2066        cdef Py_ssize_t i, j, m
     2067        p = self.p
     2068
     2069        sig_on()
     2070        for m from 1 <= m < n-1:
     2071            # Search for a nonzero entry in column m-1
     2072            i = -1
     2073            for r from m+1 <= r < n:
     2074                if h[r][m-1]:
     2075                     i = r
     2076                     break
     2077                 
     2078            if i != -1:
     2079                 # Found a nonzero entry in column m-1 that is strictly
     2080                 # below row m.  Now set i to be the first nonzero position >=
     2081                 # m in column m-1.
     2082                 if h[m][m-1]:
     2083                     i = m
     2084                 t = h[i][m-1]
     2085                 t_inv = celement_invert(t,p)
     2086                 if i > m:
     2087                     self.swap_rows_c(i,m)
     2088                     self.swap_columns_c(i,m)
     2089
     2090                 # Now the nonzero entry in position (m,m-1) is t.
     2091                 # Use t to clear the entries in column m-1 below m.
     2092                 for j from m+1 <= j < n:
     2093                     if h[j][m-1]:
     2094                         u = (h[j][m-1] * t_inv) % p
     2095                         self.add_multiple_of_row_c(j, m, p - u, 0)  # h[j] -= u*h[m]
     2096                         # To maintain charpoly, do the corresponding
     2097                         # column operation, which doesn't mess up the
     2098                         # matrix, since it only changes column m, and
     2099                         # we're only worried about column m-1 right
     2100                         # now.  Add u*column_j to column_m.
     2101                         self.add_multiple_of_column_c(m, j, u, 0)
     2102                 # end for
     2103            # end if
     2104        # end for
     2105        sig_off()
     2106        self.cache('in_hessenberg_form',True)
     2107
     2108    def _charpoly_hessenberg(self, var):
     2109        """
     2110        Transforms self in place to its Hessenberg form then computes and
     2111        returns the coefficients of the characteristic polynomial of this
     2112        matrix.
     2113       
     2114        INPUT:
     2115       
     2116       
     2117        -  ``var`` - name of the indeterminate of the
     2118           charpoly.
     2119       
     2120       
     2121        The characteristic polynomial is represented as a vector of ints,
     2122        where the constant term of the characteristic polynomial is the 0th
     2123        coefficient of the vector.
     2124        """
     2125        if self._nrows != self._ncols:
     2126            raise ArithmeticError, "charpoly not defined for non-square matrix."
     2127
     2128        cdef Py_ssize_t i, m, n,
     2129        n = self._nrows
     2130
     2131        cdef celement p, t
     2132        p = self.p
     2133
     2134        # Replace self by its Hessenberg form, and set H to this form
     2135        # for notation below.
     2136        cdef Matrix_modn_dense_template H
     2137        H = self.__copy__()
     2138        H.hessenbergize()
     2139
     2140
     2141        # We represent the intermediate polynomials that come up in
     2142        # the calculations as rows of an (n+1)x(n+1) matrix, since
     2143        # we've implemented basic arithmetic with such a matrix.
     2144        # Please see the generic implementation of charpoly in
     2145        # matrix.py to see more clearly how the following algorithm
     2146        # actually works.  (The implementation is clearer (but slower)
     2147        # if one uses polynomials to represent polynomials instead of
     2148        # using the rows of a matrix.)  Also see Cohen's first GTM,
     2149        # Algorithm 2.2.9.
     2150
     2151        cdef Matrix_modn_dense_template c
     2152        c = self.new_matrix(nrows=n+1,ncols=n+1)    # the 0 matrix
     2153        c._matrix[0][0] = 1
     2154        for m from 1 <= m <= n:
     2155            # Set the m-th row of c to (x - H[m-1,m-1])*c[m-1] = x*c[m-1] - H[m-1,m-1]*c[m-1]
     2156            # We do this by hand by setting the m-th row to c[m-1]
     2157            # shifted to the right by one.  We then add
     2158            # -H[m-1,m-1]*c[m-1] to the resulting m-th row.
     2159            for i from 1 <= i <= n:
     2160                c._matrix[m][i] = c._matrix[m-1][i-1]
     2161            # the p-.. below is to keep scalar normalized between 0 and p.
     2162            c.add_multiple_of_row_c(m, m-1, p - H._matrix[m-1][m-1], 0)
     2163            t = 1
     2164            for i from 1 <= i < m:
     2165                t = (t*H._matrix[m-i][m-i-1]) % p
     2166                # Set the m-th row of c to c[m] - t*H[m-i-1,m-1]*c[m-i-1]
     2167                c.add_multiple_of_row_c(m, m-i-1, p - (t*H._matrix[m-i-1][m-1])%p, 0)
     2168               
     2169        # The answer is now the n-th row of c.
     2170        v = []
     2171        for i from 0 <= i <= n:
     2172            v.append(int(c._matrix[n][i]))
     2173        R = self._base_ring[var]    # polynomial ring over the base ring
     2174        return R(v)
     2175
     2176    def rank(self):
     2177        """
     2178        Return the rank of this matrix.
     2179       
     2180        EXAMPLES::
     2181       
     2182            sage: m = matrix(GF(7),5,range(25))
     2183            sage: m.rank()
     2184            2
     2185       
     2186        Rank is not implemented over the integers modulo a composite yet.::
     2187       
     2188            sage: m = matrix(Integers(4), 2, [2,2,2,2])
     2189            sage: m.rank()
     2190            Traceback (most recent call last):
     2191            ...
     2192            NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 4'.
     2193        """
     2194        cdef Matrix_modn_dense_template A
     2195        if self.p > 2 and is_prime(self.p):
     2196            x = self.fetch('rank')
     2197            if not x is None:
     2198                return x
     2199            r = Integer(linbox_rank(self.p, self._entries, self._nrows, self._ncols))
     2200            self.cache('rank', r)
     2201            return r
     2202        else:
     2203            # linbox is very buggy for p=2
     2204            return matrix_dense.Matrix_dense.rank(self)
     2205 
     2206    def determinant(self):
     2207        """
     2208        Return the determinant of this matrix.
     2209       
     2210        EXAMPLES::
     2211       
     2212            sage: m = matrix(GF(101),5,range(25))     
     2213            sage: m.det()                             
     2214            0
     2215       
     2216        ::
     2217       
     2218            sage: m = matrix(Integers(4), 2, [2,2,2,2])           
     2219            sage: m.det()                             
     2220            0
     2221       
     2222        TESTS::
     2223
     2224            sage: m = random_matrix(GF(3), 3, 4)
     2225            sage: m.determinant()
     2226            Traceback (most recent call last):
     2227            ...
     2228            ValueError: self must be a square matrix
     2229        """
     2230        if self._nrows != self._ncols:
     2231            raise ValueError, "self must be a square matrix"
     2232        if self._nrows == 0:
     2233            return self._coerce_element(1)
     2234
     2235        if self.p > 2 and is_prime(self.p):
     2236            x = self.fetch('det')
     2237            if not x is None:
     2238                return x
     2239            d = linbox_det(self.p, self._entries, self._nrows, self._ncols)
     2240            d2 = self._coerce_element(d)
     2241            self.cache('det', d2)
     2242            return d2
     2243        else:
     2244            return matrix_dense.Matrix_dense.determinant(self)
     2245
     2246    def randomize(self, density=1, nonzero=False):
     2247        """
     2248        Randomize ``density`` proportion of the entries of this matrix,
     2249        leaving the rest unchanged.
     2250       
     2251        INPUT:
     2252       
     2253        -  ``density`` - Integer; proportion (roughly) to be considered for
     2254           changes
     2255        -  ``nonzero`` - Bool (default: ``False``); whether the new entries
     2256           are forced to be non-zero
     2257       
     2258        OUTPUT:
     2259       
     2260        -  None, the matrix is modified in-space
     2261       
     2262        EXAMPLES::
     2263       
     2264            sage: A = matrix(GF(5), 5, 5, 0)
     2265            sage: A.randomize(0.5); A
     2266            [0 0 0 2 0]
     2267            [0 3 0 0 2]
     2268            [4 0 0 0 0]
     2269            [4 0 0 0 0]
     2270            [0 1 0 0 0]
     2271            sage: A.randomize(); A
     2272            [3 3 2 1 2]
     2273            [4 3 3 2 2]
     2274            [0 3 3 3 3]
     2275            [3 3 2 2 4]
     2276            [2 2 2 1 4]
     2277        """
     2278        density = float(density)
     2279        if density <= 0:
     2280            return
     2281        if density > 1:
     2282            density = float(1)
     2283
     2284        self.check_mutability()
     2285        self.clear_cache()
     2286
     2287        cdef randstate rstate = current_randstate()
     2288
     2289        cdef int nc, p = <int>self.p
     2290        cdef long pm1
     2291       
     2292        if not nonzero:
     2293            # Original code, before adding the ``nonzero`` option.
     2294            if density == 1:
     2295                for i from 0 <= i < self._nrows*self._ncols:
     2296                    self._entries[i] = rstate.c_random() % p
     2297            else:
     2298                nc = self._ncols
     2299                num_per_row = int(density * nc)
     2300                sig_on()
     2301                for i from 0 <= i < self._nrows:
     2302                    for j from 0 <= j < num_per_row:
     2303                        k = rstate.c_random() % nc
     2304                        self._matrix[i][k] = rstate.c_random() % p
     2305                sig_off()
     2306        else:
     2307            # New code, to implement the ``nonzero`` option.
     2308            pm1 = p - 1
     2309            if density == 1:
     2310                for i from 0 <= i < self._nrows*self._ncols:
     2311                    self._entries[i] = (rstate.c_random() % pm1) + 1
     2312            else:
     2313                nc = self._ncols
     2314                num_per_row = int(density * nc)
     2315                sig_on()
     2316                for i from 0 <= i < self._nrows:
     2317                    for j from 0 <= j < num_per_row:
     2318                        k = rstate.c_random() % nc
     2319                        self._matrix[i][k] = (rstate.c_random() % pm1) + 1
     2320                sig_off()
     2321
     2322    def _magma_init_(self, magma):
     2323        """
     2324        Returns a string representation of self in Magma form.
     2325       
     2326        INPUT:
     2327       
     2328       
     2329        -  ``magma`` - a Magma session
     2330       
     2331       
     2332        OUTPUT: string
     2333       
     2334        EXAMPLES::
     2335       
     2336            sage: a = matrix(GF(389),2,2,[1..4])
     2337            sage: magma(a)                                         # optional - magma
     2338            [  1   2]
     2339            [  3   4]
     2340            sage: a._magma_init_(magma)                            # optional - magma
     2341            'Matrix(GF(389),2,2,StringToIntegerSequence("1 2 3 4"))'
     2342       
     2343        A consistency check::
     2344       
     2345            sage: a = random_matrix(GF(13),50); b = random_matrix(GF(13),50)
     2346            sage: magma(a*b) == magma(a)*magma(b)                  # optional - magma
     2347            True
     2348        """
     2349        s = self.base_ring()._magma_init_(magma)
     2350        return 'Matrix(%s,%s,%s,StringToIntegerSequence("%s"))'%(
     2351            s, self._nrows, self._ncols, self._export_as_string())
     2352
     2353    cpdef _export_as_string(self):
     2354        """
     2355        Return space separated string of the entries in this matrix.
     2356
     2357        EXAMPLES::
     2358
     2359            sage: w = matrix(GF(997),2,3,[1,2,5,-3,8,2]); w
     2360            [  1   2   5]
     2361            [994   8   2]
     2362            sage: w._export_as_string()
     2363            '1 2 5 994 8 2'
     2364        """
     2365        cdef int ndigits = len(str(self.p))
     2366
     2367        cdef Py_ssize_t i, n
     2368        cdef char *s, *t
     2369
     2370        if self._nrows == 0 or self._ncols == 0:
     2371            data = ''
     2372        else:
     2373            n = self._nrows*self._ncols*(ndigits + 1) + 2  # spaces between each number plus trailing null
     2374            s = <char*> sage_malloc(n * sizeof(char))
     2375            t = s
     2376            sig_on()
     2377            for i in range(self._nrows * self._ncols):
     2378                t += snprintf(t, ndigits+2, "%ld ", <long>self._entries[i])
     2379               
     2380            sig_off()
     2381            data = str(s)[:-1]
     2382            sage_free(s)
     2383        return data
     2384
     2385    def _list(self):
     2386        """
     2387        Return list of elements of self.  This method is called by ``self.list()``.
     2388
     2389        EXAMPLES::
     2390
     2391            sage: w = matrix(GF(19), 2, 3, [1..6])
     2392            sage: w.list()
     2393            [1, 2, 3, 4, 5, 6]
     2394            sage: w._list()
     2395            [1, 2, 3, 4, 5, 6]
     2396            sage: w.list()[0].parent()
     2397            Finite Field of size 19
     2398
     2399        TESTS::
     2400
     2401            sage: w = random_matrix(GF(3),100)
     2402            sage: w.parent()(w.list()) == w
     2403            True
     2404        """
     2405        cdef Py_ssize_t i
     2406        F = self.base_ring()
     2407        return [F(<int>self._entries[i]) for i in range(self._nrows*self._ncols)]
     2408
     2409    def lift(self):
     2410        """
     2411        Return the lift of this matrix to the integers.
     2412
     2413        EXAMPLES::
     2414
     2415            sage: A = matrix(GF(7),2,3,[1..6])
     2416            sage: A.lift()
     2417            [1 2 3]
     2418            [4 5 6]
     2419            sage: A.lift().parent()
     2420            Full MatrixSpace of 2 by 3 dense matrices over Integer Ring
     2421
     2422
     2423            sage: A = matrix(GF(16007),2,3,[1..6])
     2424            sage: A.lift()
     2425            [1 2 3]
     2426            [4 5 6]
     2427            sage: A.lift().parent()
     2428            Full MatrixSpace of 2 by 3 dense matrices over Integer Ring
     2429
     2430        Subdivisions are preserved when lifting::
     2431
     2432            sage: A.subdivide([], [1,1]); A
     2433            [1||2 3]
     2434            [4||5 6]
     2435            sage: A.lift()
     2436            [1||2 3]
     2437            [4||5 6]
     2438        """
     2439        cdef Py_ssize_t i, j
     2440
     2441        cdef Matrix_integer_dense L
     2442        L = Matrix_integer_dense.__new__(Matrix_integer_dense,
     2443                                         self.parent().change_ring(ZZ),
     2444                                         0, 0, 0)
     2445        cdef mpz_t* L_row
     2446        cdef celement* A_row
     2447        for i from 0 <= i < self._nrows:
     2448            L_row = L._matrix[i]
     2449            A_row = self._matrix[i]
     2450            for j from 0 <= j < self._ncols:
     2451                mpz_init_set_d(L_row[j], A_row[j])
     2452        L._initialized = 1
     2453        L.subdivide(self.subdivisions())
     2454        return L
     2455
     2456
     2457    def _matrices_from_rows(self, Py_ssize_t nrows, Py_ssize_t ncols):
     2458        """
     2459        Make a list of matrix from the rows of this matrix.  This is a
     2460        fairly technical function which is used internally, e.g., by
     2461        the cyclotomic field linear algebra code.
     2462       
     2463        INPUT:
     2464       
     2465        - ``nrows, ncols`` - integers
     2466           
     2467        OUTPUT:
     2468           
     2469        - ``list`` - a list of matrices
     2470        """
     2471        if nrows * ncols != self._ncols:
     2472            raise ValueError, "nrows * ncols must equal self's number of columns"
     2473
     2474        from matrix_space import MatrixSpace
     2475        F = self.base_ring()
     2476        MS = MatrixSpace(F, nrows, ncols)
     2477
     2478        cdef Matrix_modn_dense_template M
     2479        cdef Py_ssize_t i
     2480        cdef Py_ssize_t n = nrows * ncols
     2481        ans = []
     2482        for i from 0 <= i < self._nrows:
     2483            # Quickly construct a new mod-p matrix
     2484            M = self.__class__.__new__(self.__class__, MS, 0,0,0)
     2485            M.p = self.p
     2486            # Set the entries
     2487            memcpy(M._entries, self._entries+i*n, sizeof(celement)*n)
     2488            ans.append(M)
     2489        return ans
     2490
     2491    def __nonzero__(self):
     2492        return not linbox_is_zero(self.p, self._entries, self._nrows, self._ncols)
     2493
     2494    _matrix_from_rows_of_matrices = staticmethod(__matrix_from_rows_of_matrices)
     2495
     2496    cdef int _copy_row_to_mod_int_array(self, mod_int *to, Py_ssize_t i):
     2497        cdef Py_ssize_t j
     2498        cdef celement *_from = self._entries+(i*self._ncols)
     2499        for j in range(self._ncols):
     2500            to[j] = <mod_int>_from[j]
     2501
     2502    # def new_matrix(self, nrows=None, ncols=None, entries=0, coerce=True, copy=True, sparse=None):
     2503    #     """
     2504    #     Create a matrix in the parent of this matrix with the given number
     2505    #     of rows, columns, etc. The default parameters are the same as for
     2506    #     self.
     2507
     2508    #     INPUT:
     2509
     2510    #     These three variables get sent to :func:`matrix_space`:
     2511
     2512    #     - ``nrows``, ``ncols`` - number of rows and columns in returned
     2513    #       matrix. If not specified, defaults to ``None`` and will give a
     2514    #       matrix of the same size as self.
     2515    #     - ``sparse`` - whether returned matrix is sparse or
     2516    #       not. Defaults to same value as self.
     2517
     2518    #     The remaining three variables (``coerce``, ``entries``, and
     2519    #     ``copy``) are used by
     2520    #     :func:`sage.matrix.matrix_space.MatrixSpace` to construct the
     2521    #     new matrix.
     2522
     2523    #     .. warning::
     2524
     2525    #        This function called with no arguments returns the zero
     2526    #        matrix of the same dimension and sparseness of self.
     2527
     2528    #     EXAMPLE::
     2529
     2530    #         sage: from sage.matrix.matrix_modn_dense_float import Matrix_modn_dense_float
     2531    #         sage: from sage.matrix.matrix_modn_dense_double import Matrix_modn_dense_double
     2532
     2533    #         sage: K = GF(7)
     2534    #         sage: MS = MatrixSpace(K,100,100)
     2535    #         sage: A = Matrix_modn_dense_float(MS, None, False, False)
     2536    #         sage: B = Matrix_modn_dense_double(MS, None, False, False)
     2537    #         sage: type(A)
     2538    #         <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     2539    #         sage: type(A.new_matrix())
     2540    #         <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
     2541    #         sage: type(B)
     2542    #         <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
     2543    #         sage: type(B.new_matrix())
     2544    #         <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
     2545    #     """
     2546    #     MS = self.matrix_space(nrows, ncols, sparse=sparse)
     2547    #     if sparse:
     2548    #         return MS(entries=entries, coerce=coerce, copy=copy)
     2549    #     else:
     2550    #         return self.__class__(MS, entries, coerce, copy)
     2551
     2552cpdef __matrix_from_rows_of_matrices(X):
     2553    """
     2554    INPUT:
     2555       
     2556    - ``X`` - a nonempty list of matrices of the same size mod a
     2557       single modulus `n`
     2558
     2559    EXAMPLES::
     2560
     2561        sage: X = [random_matrix(GF(17), 4, 4) for _ in range(10)]; X
     2562        [
     2563        [ 2 14  0 15]  [12 14  3 13]  [ 9 15  8  1]  [ 2 12  6 10]
     2564        [11 10 16  2]  [10  1 14  6]  [ 5  8 10 11]  [12  0  6  9]
     2565        [ 9  4 10 14]  [ 2 14 13  7]  [ 5 12  4  9]  [ 7  7  3  8]
     2566        [ 1 14  3 14], [ 6 14 10  3], [15  2  6 11], [ 2  9  1  5],
     2567        <BLANKLINE>
     2568        [12 13  7 16]  [ 5  3 16  2]  [14 15 16  4]  [ 1 15 11  0]
     2569        [ 7 11 11  1]  [11 10 12 14]  [14  1 12 13]  [16 13  8 14]
     2570        [ 0  2  0  4]  [ 0  7 16  4]  [ 5  5 16 13]  [13 14 16  4]
     2571        [ 7  9  8 15], [ 6  5  2  3], [10 12  1  7], [15  6  6  6],
     2572        <BLANKLINE>
     2573        [ 4 10 11 15]  [13 12  5  1]
     2574        [11  2  9 14]  [16 13 16  7]
     2575        [12  5  4  4]  [12  2  0 11]
     2576        [ 2  0 12  8], [13 11  6 15]
     2577        ]
     2578        sage: X[0]._matrix_from_rows_of_matrices(X) # indirect doctest
     2579        [ 2 14  0 15 11 10 16  2  9  4 10 14  1 14  3 14]
     2580        [12 14  3 13 10  1 14  6  2 14 13  7  6 14 10  3]
     2581        [ 9 15  8  1  5  8 10 11  5 12  4  9 15  2  6 11]
     2582        [ 2 12  6 10 12  0  6  9  7  7  3  8  2  9  1  5]
     2583        [12 13  7 16  7 11 11  1  0  2  0  4  7  9  8 15]
     2584        [ 5  3 16  2 11 10 12 14  0  7 16  4  6  5  2  3]
     2585        [14 15 16  4 14  1 12 13  5  5 16 13 10 12  1  7]
     2586        [ 1 15 11  0 16 13  8 14 13 14 16  4 15  6  6  6]
     2587        [ 4 10 11 15 11  2  9 14 12  5  4  4  2  0 12  8]
     2588        [13 12  5  1 16 13 16  7 12  2  0 11 13 11  6 15]           
     2589   
     2590    OUTPUT: A single matrix mod p whose ith row is X[i].list().
     2591    """
     2592    # The code below is just a fast version of the following:
     2593    ##     from constructor import matrix
     2594    ##     K = X[0].base_ring()
     2595    ##     v = sum([y.list() for y in X],[])
     2596    ##     return matrix(K, len(X), X[0].nrows()*X[0].ncols(), v)
     2597
     2598    from matrix_space import MatrixSpace
     2599    cdef Matrix_modn_dense_template A, T
     2600    cdef Py_ssize_t i, n, m
     2601    n = len(X)
     2602
     2603    T = X[0]
     2604    m = T._nrows * T._ncols
     2605    A = T.__class__.__new__(T.__class__, MatrixSpace(X[0].base_ring(), n, m), 0, 0, 0)
     2606    A.p = T.p
     2607
     2608    for i from 0 <= i < n:
     2609        T = X[i]
     2610        memcpy(A._entries + i*m, T._entries, sizeof(celement)*m)
     2611    return A
     2612
  • new file sage/matrix/matrix_modn_dense_template_header.pxi

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_dense_template_header.pxi
    - +  
     1"""
     2Dense Matrix Template for C/C++ Library Interfaces
     3"""
     4
     5cimport matrix_dense
     6
     7cdef extern from "../ext/multi_modular.h":
     8    ctypedef unsigned long mod_int
     9    mod_int MOD_INT_OVERFLOW
     10
     11cdef class Matrix_modn_dense_template(matrix_dense.Matrix_dense):
     12    cdef celement **_matrix
     13    cdef celement *_entries
     14    cdef mod_int p
     15    cdef xgcd_eliminate (self, celement * row1, celement* row2, Py_ssize_t start_col)
     16    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value)
     17    cpdef _export_as_string(self)
     18    cdef int _copy_row_to_mod_int_array(self, mod_int *to, Py_ssize_t i)
  • sage/matrix/matrix_modn_sparse.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_modn_sparse.pyx
    a b  
    389389            [ 9 12 15]
    390390            [19 26 33]
    391391            sage: type(c)
    392             <type 'sage.matrix.matrix_modn_dense.Matrix_modn_dense'>
     392            <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
    393393
    394394            sage: a = matrix(GF(2), 20, 20, sparse=True)
    395395            sage: a*a == a._matrix_times_matrix_dense(a)
  • sage/matrix/matrix_space.py

    diff -r d263952c945f -r 9b63c2421a7f sage/matrix/matrix_space.py
    a b  
    3434import matrix_generic_dense
    3535import matrix_generic_sparse
    3636
     37import matrix_modn_dense
    3738import matrix_modn_sparse
    3839
    3940import matrix_mod2_dense
     
    848849                import matrix_modn_dense_double, matrix_modn_dense_float
    849850                if R.order() == 2:
    850851                    return matrix_mod2_dense.Matrix_mod2_dense
     852                elif R.order() < matrix_modn_dense_float.MAX_MODULUS:
     853                    return matrix_modn_dense_float.Matrix_modn_dense_float
    851854                elif R.order() < matrix_modn_dense_double.MAX_MODULUS:
    852855                    return matrix_modn_dense_double.Matrix_modn_dense_double
    853                 elif R.order() < matrix_modn_dense_float.MAX_MODULUS:
    854                     return matrix_modn_dense_float.Matrix_modn_dense_float
     856                # elif R.order() < matrix_modn_dense.MAX_MODULUS:
     857                #     return matrix_modn_dense.Matrix_modn_dense
    855858                return matrix_generic_dense.Matrix_generic_dense
    856859            elif sage.rings.polynomial.multi_polynomial_ring_generic.is_MPolynomialRing(R) and R.base_ring().is_field():
    857860                return matrix_mpolynomial_dense.Matrix_mpolynomial_dense
  • sage/structure/sage_object.pyx

    diff -r d263952c945f -r 9b63c2421a7f sage/structure/sage_object.pyx
    a b  
    10751075    all "standard pickles" unpickle::
    10761076
    10771077        sage: sage.structure.sage_object.unpickle_all()  # (4s on sage.math, 2011)
     1078        doctest:... DeprecationWarning: This class is replaced by Matrix_modn_dense_float/Matrix_modn_dense_double.
    10781079        Successfully unpickled ... objects.
    10791080        Failed to unpickle 0 objects.
    10801081