Ticket #12103: trac12103_meataxe_rel11900.patch

module_list.py

@@ -946 +946 @@
               extra_compile_args = m4ri_extra_compile_args,
               language="c++"),
 
+    Extension('sage.matrix.matrix_modpn_dense',
+              sources = ['sage/matrix/matrix_modpn_dense.pyx'],
+              libraries = ['meataxe', 'csage'],
+              depends = [SAGE_INC + "meataxe.h"]),
+
     Extension('sage.matrix.matrix_modn_dense',
               sources = ['sage/matrix/matrix_modn_dense.pyx'],
               libraries = ['gmp']),

sage/matrix/constructor.py

@@ -1062 +1062 @@
 
         sage: K.<a>=FiniteField(3^2)
         sage: random_matrix(K, 2, 5)
-        [      1     2*a       1   a + 2       2]
-        [      a 2*a + 1       0       2       1]
-        
+        [      a       0       2       2   a + 1]
+        [      a   a + 2   a + 1 2*a + 1     2*a]
+
         sage: random_matrix(RR, 3, 4, density=0.66)
-        [ 0.000000000000000  0.566500636438206 0.0870635178173962  0.000000000000000]
-        [-0.662290145671671  0.000000000000000  0.475667133865666  0.000000000000000]
-        [-0.276405104068647  0.000000000000000  0.000000000000000 -0.636689607643642]
-        
+        [ 0.000000000000000 -0.615392177288427 -0.636689607643642  0.000000000000000]
+        [ 0.575979457284225  0.000000000000000 0.0870635178173962  0.000000000000000]
+        [-0.662290145671671  0.000000000000000  0.000000000000000 0.0362592516947566]
+
         sage: A = random_matrix(ComplexField(32), 3, density=0.8, sparse=True); A
-        [                             0   -0.193242896 + 0.460619713*I                              0]
-        [   0.909948633 + 0.611092515*I 0.00128001347 + 0.0659103142*I   0.199796580 + 0.0955426861*I]
-        [                             0                              0    0.354729790 - 0.184624095*I]
+        [                             0    0.909948633 + 0.611092515*I                              0]
+        [  -0.728572036 - 0.455087671*I    0.354729790 - 0.184624095*I 0.00128001347 + 0.0659103142*I]
+        [                             0                              0   -0.360952012 - 0.698699617*I]
         sage: A.is_sparse()
         True
 

sage/matrix/matrix0.pyx

@@ -3335 +3335 @@
         x = self.fetch('rank')
         if not x is None: return x
         if self._nrows == 0 or self._ncols == 0:
-            return 0        
+            return 0
         r = len(self.pivots())
         self.cache('rank', r)
         return r

sage/matrix/matrix2.pyx

@@ -2628 +2628 @@
         We test that the generic code is called for matrices over fields,
         lacking any more specific routine. ::
 
-            sage: F.<a> = FiniteField(5^2)
+            sage: F.<a> = FiniteField(5^4)
             sage: A = matrix(F, 3, 4, [[  1,   a,     1+a,  a^3+a^5],
             ...                        [  a, a^4,   a+a^4,  a^4+a^8],
             ...                        [a^2, a^6, a^2+a^6, a^5+a^10]])
             sage: set_verbose(1)
             sage: A.right_kernel(algorithm='default')
-            verbose ...
+            verbose 1 ...
             verbose 1 (<module>) computing right kernel matrix over an arbitrary field for 3x4 matrix
-            ...
-            verbose 1 (<module>) done computing right kernel matrix over an arbitrary field for 3x4 matrix
-            ...
-            Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^2
+            verbose 1 (<module>) generic in-place Gauss elimination on 3 x 4 matrix
+            verbose 1 (<module>) done with gauss echelon form (time = ...)
+            verbose 1 (<module>) done computing right kernel matrix over an arbitrary field for 3x4 matrix (time = ...)
+            verbose 1 (<module>) generic in-place Gauss elimination on 2 x 4 matrix
+            verbose 1 (<module>) done with gauss echelon form (time = ...)
+            verbose 1 ...
+            Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^4
             Basis matrix:
-            [      1       0 3*a + 4 2*a + 2]
-            [      0       1     2*a 3*a + 3]
+            [                1                 0     a^3 + a^2 + 3   2*a^3 + a^2 + 1]
+            [                0                 1 4*a^3 + 4*a^2 + 1 3*a^3 + 4*a^2 + 4]
             sage: set_verbose(0)
 
         Over the Integers:

new file sage/matrix/matrix_modpn_dense.pxd

@@ +1 @@
+#*****************************************************************************
+#
+#    Dense matrices over GF(p^n), p^n<255, using LibMeataxe as backend
+#
+#    Copyright (C) 2011 Simon A. King <simon.king@uni-jena.de>
+#
+#  Distributed under the terms of the GNU General Public License (GPL),
+#  version 2 or later (at your choice)
+#
+#    This code is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    General Public License for more details.
+#
+#  The full text of the GPL is available at:
+#
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include "../ext/interrupt.pxi"
+
+#
+# First, something from Python.h and string.h
+#
+
+cdef extern from "Python.h":
+    object PyString_FromStringAndSize(char *s, Py_ssize_t len)
+    char* PyString_AsString(object string)
+
+#
+# Next, import SOME types from meataxe.h
+# (most types are not needed, but listed here 
+# in the comments, for completeness)
+#
+
+cdef extern from "stdlib.h":
+    void* calloc ( size_t num, size_t size )
+    void free ( void * ptr )
+    void * memset ( void * ptr, int value, size_t num )
+
+cdef extern from "meataxe.h":
+    # general ctype emulations
+    ctypedef int size_t
+    ctypedef unsigned long Ulong
+    ctypedef unsigned short Ushort
+    ctypedef unsigned char Uchar
+    ctypedef struct FILE
+    
+    # assume small fields (default) 
+    ctypedef unsigned char FEL
+    ctypedef FEL *PTR
+    cdef extern FEL tmult[256][256]
+    cdef extern FEL tadd[256][256]
+    cdef extern FEL taddinv[256], tmultinv[256]
+    cdef extern FEL tinsert[8][256], textract[8][256]
+
+    # pivots
+    ctypedef struct piv_t:
+        long i,m
+
+    #
+    # cdef extern long zfl, zchar       #/* Field order and characteristic */
+    # cdef extern FEL zgen              #/* Generator */
+    # cdef extern long znoc             #/* Number of columns for row ops */
+    cdef extern size_t zrowsize         #/* Row size in memory */
+    cdef extern size_t zrowsize_io      #/* Row size in files */
+    cdef extern size_t MPB              # number of FEL stored in one byte
+    # cdef extern char zzzversion[]
+    # cdef extern char zzz_cc[]
+        
+#########################################
+# function prototypes
+    # FEL zadd(FEL a,FEL b)
+    # FEL zneg(FEL a)
+    # FEL zmul(FEL a,FEL b)
+    # FEL zinv(FEL a)
+    #
+    # FEL zembed(FEL a, long subfield)  #/* Embed from subfield */
+    # FEL zrestrict(FEL a, long subfield)       #/* Restrict to subfield */
+    void *copyrows(PTR dest, PTR src, size_t nor)
+
+    FEL zitof(long l)
+    long zftoi(FEL f)
+
+    int zsetlen(long ncols)             #/* Set row size */
+    int zsetfield(long field)
+    void zinsert(PTR row, long col, FEL mark)
+    FEL zextract(PTR row, long col)
+    void zinsert_step(PTR *pos, int *idx, FEL mark)
+    FEL zextract_step(PTR *pos, int *idx)
+
+    piv_t _zfindpiv_(PTR row)
+    # PTR zextractcol(PTR mat, long nor, long col, PTR result)
+    # PTR zaddrow(PTR dest, PTR src)
+    # PTR zsubrow(PTR dest, PTR src)
+    # PTR znegrow(PTR row)
+    PTR zmulrow(PTR row, FEL mark)
+    PTR zaddmulrow(PTR dest, PTR src, FEL f)
+    # PTR zmaprow(PTR row, PTR matrix, long nor, PTR result)
+    # PTR zpermrow(PTR row, long *perm, PTR result)
+
+    void zmoverow(PTR dest, PTR src)
+    void zcpyrow(PTR dest, PTR src, long sz)
+    # void zswaprow(PTR dest, PTR src)
+    # int zcmprow(PTR dest, PTR src)
+    size_t zsize(long nrows)
+    # PTR zalloc(long NROWS)
+    void zadvance(PTR *ptr, long nrows)
+    void zsteprow(PTR *ptr)
+    void zbackrow(PTR *ptr)
+    PTR ptrPlus(PTR base, long offset)
+    # int mtxinit()
+
+    # skip gap-conversion
+    # bitstring operations
+    # ctypedef struct bitstring_t:
+    #     Ushort id     #/* always = T_BITS */
+    #     char dummy[sizeof(long)-sizeof(Ushort)]
+    #     unsigned char buf[1]
+    # cdef extern size_t bs_size
+    # void bs_print(bitstring_t *b)
+    # char *bs_desc(bitstring_t *b)
+    # void bs_setlen(int l)
+    # void bs_reset(bitstring_t *x)
+    # bitstring_t *bs_alloc()
+    # bitstring_t *bs_dup(bitstring_t *s)
+    # void bs_free(bitstring_t *b)
+    # void bs_set(bitstring_t *b, int i)
+    # void bs_clear(bitstring_t *b, int i)
+    # int bs_test(bitstring_t *b, int i)
+    # bitstring_t *bs_read(FILE *f)
+    # int bs_write(FILE *f, bitstring_t *b)
+    # void bs_and(bitstring_t *dest, bitstring_t *src)
+    # void bs_minus(bitstring_t *dest, bitstring_t *src)
+    # void bs_or(bitstring_t *dest, bitstring_t *src)
+    # int bs_match(bitstring_t *x, bitstring_t *y)
+    # int bs_issub(bitstring_t *a, bitstring_t *b)
+    # int bs_cmp(bitstring_t *a, bitstring_t *b)
+    # void bs_cpy(bitstring_t *a, bitstring_t *b)
+
+    # file i/o
+    # size_t zreadlong(FILE *f, long *buf, size_t n)
+    # size_t zwritelong(FILE *f, long *buf, size_t n)
+    # size_t zreadvec(FILE *f, PTR buf, size_t n)
+    # size_t zwritevec(FILE *f, PTR buf, size_t n)
+    # int zseek(FILE *f, long pos)
+    # FILE *zreadhdr(char *name, long *fld, long *nor, long *noc)
+    # FILE *zwritehdr(char *name, long fld, long nor, long noc)
+
+    # Gauss elimination and related functions
+    # int zmkpivot(PTR matrix, long nor, long *piv)
+    long zmkechelon(PTR matrix, long nor, long *piv)
+    # int zcleanrow(PTR row, PTR matrix, long nor, long *piv)
+    # int zcleanrow2(PTR row,PTR matrix,long nor,long *piv,PTR row2)
+    # long znullsp(PTR matrix, long nor, long *piv, PTR nsp)
+    # int zmatinv(PTR mat, PTR result)
+
+    #  Projection on the quotient
+    # int zquotinit(PTR subspace, long dim, long *piv)
+    # int zquot(PTR space, long dim, PTR quot)
+    # int zquotop(PTR matrix, PTR quot)
+
+    # Seed vector generation
+    # cdef extern PTR zpseed_vec                #/* The seed vector */
+    # void zpseed_free()
+    # int zpseed_init(long sdim, PTR sbasis)
+    # long zpseed_make(long num)
+    #long zpseed_next()
+
+    # Command line parsing
+    # ctypedef struct proginfo_t:
+    #     char *name
+    #     char *shortdesc
+    #     char *rcsrev
+    #     char **helptext
+    # cdef extern char opt_char         #/* Current option */
+    # cdef extern char opt_text[50]     #/* Option text */
+    # cdef extern char *opt_text_ptr    #/* Current position in option text */
+    # cdef extern int opt_ind           #/* Index in argv[] list */
+    # cdef extern char MeatAxeBinDir[250]       #/* MeatAxe program directory */
+
+    # void initargs(int argc, char **argv, proginfo_t *pi)
+    # int zgetopt(char *pattern)
+    # long getint()
+
+    # Messages
+    # cdef extern int mtxerrno
+    # cdef extern int mtxerraction
+    # void mtxerror(char *text)
+    # void errexit(int code, char *text)
+    # int errhandler(int errno, char *errfile, int errline)
+    # int Message(FILE *f, char *fmt, ...)
+    # int ErrorExit(char *fmt, ...)
+    # cdef extern int msg_level
+
+    # CPU time
+    # void prtimes()
+
+    # Random numbers
+    # void RandInit(unsigned seed)
+    # long int Random()
+
+    # Miscellaneous
+    # long gcd(long a, long b)
+    # long lcm(long a, long b)
+
+    # Matrices
+    ctypedef struct matrix_t:
+        Ushort id                       #/* Always = T_MATRIX */
+        long fl, nor, noc               #/* Field, #rows, #columns */
+        PTR d                           #/* Pointer to data area */
+
+    # void matprint(char *name, matrix_t *p)
+    # char *matdesc(matrix_t *s)
+    matrix_t *matalloc(long fl, long nor, long noc)
+    matrix_t *matmalloc(long fl, long nor, long noc)
+    matrix_t *MTXmatid(long fl, long nor)
+    void matfree(matrix_t *m)
+
+    long matcmp(matrix_t *m1, matrix_t *m2)
+    matrix_t *matdup(matrix_t *src)
+    # int  matmove(matrix_t *dest, matrix_t *src)
+    matrix_t *_matextract(matrix_t *src, long first, long n)
+    matrix_t *matread(FILE *f)
+    int matwrite(FILE *f, matrix_t *mat)
+    matrix_t *matload(char *fn)
+    int matsave(matrix_t *mat, char *fn)
+
+    matrix_t *matadd(matrix_t *dest, matrix_t *src)
+    matrix_t *matmul(matrix_t *dest, matrix_t *src)
+    int matmul_result(matrix_t *dest, matrix_t *src, matrix_t *res)
+    matrix_t *matmulF(matrix_t *src, FEL mark)
+    matrix_t *matsub(matrix_t *dest, matrix_t *src)
+    matrix_t *matneg(matrix_t *src)
+    matrix_t *mattr(matrix_t *src)
+    matrix_t *matinv(matrix_t *src)
+    # int chkechelon(matrix_t *mat)
+    matrix_t *matpower(matrix_t *mat, long n)
+    long matorder(matrix_t *mat)
+    
+    # Higher-level gauss functions
+    matrix_t *echelon(matrix_t *mat)
+    long nullity(matrix_t *mat)
+    matrix_t *nullspace(matrix_t *mat)
+
+    # ycomp.c
+    # int spccomp(matrix_t *m1, matrix_t *m2, long n)
+    # int spcequal(matrix_t *m1, matrix_t *m2)
+    # int spccontains(matrix_t *m1, matrix_t *m2)
+
+    # yspin.c
+    # matrix_t *matspin(matrix_t *seed, int ngen, matrix_t *gen[])
+    # matrix_t *matspin_f(matrix_t *seed, int ngen, matrix_t *gen[], long *sdim)
+
+    # Spin-up, split, and standard basis
+    # long zspinup(PTR space, long nseed, long *piv, int ngen, PTR gen[],       int gentype)
+    # int zsbasis(PTR seed, long nseed, int ngen, PTR gen[], PTR space, long *piv, PTR basis)
+    # cdef extern long *split_pivot
+
+    # ... various options in meataxe.h
+    # int spinup(matrix_t *seed,int ngen,matrix_t *gen[],int options, matrix_t **subspace)
+    # int split(matrix_t *subspace, int ngen, matrix_t *gen[], matrix_t **sub, matrix_t **quot)
+    # matrix_t *quotproj(matrix_t *subspace, matrix_t *vectors)
+    # matrix_t *sbasis(matrix_t *seed, int ngen, matrix_t **gen)
+    # int chbasis(matrix_t *basis, int ngen, matrix_t **gen,matrix_t **newgen)
+
+    # Polynomials
+    # ctypedef struct poly_t:
+    #     Ushort id                     #/* Always = T_POLY */
+    #     long fl
+    #     long deg
+    #     size_t size  #/* deg <= size-1 */
+    #     FEL *buf
+
+    # char *poldesc(poly_t *s)
+    # poly_t *polalloc(long fl, long degree)
+    # void polfree(poly_t *x)
+    # poly_t *poldup(poly_t *x)
+    # poly_t *poladd(poly_t *dest, poly_t *src)
+    # poly_t *polmul(poly_t *dest, poly_t *src)
+    # poly_t *poldivmod(poly_t *a, poly_t *b)
+    # poly_t *polmod(poly_t *a, poly_t *b)
+    # void polprint(char *name, poly_t *p)
+    # poly_t *polread(FILE *f)
+    # int polwrite(FILE *f, poly_t *p)
+    # int polcmp(poly_t *a, poly_t *b)
+    # poly_t *polgcd(poly_t *a, poly_t *b)
+    # poly_t *polderive(poly_t *p)
+    # int vec2pol(PTR vec, poly_t *pol)
+    # int pol2vec(poly_t *pol, PTR vec)
+    # int polpack(poly_t *pol, PTR vec)
+    # poly_t *polshiftmod(poly_t *p, long n, poly_t *q)
+
+    # The word generator
+    # ctypedef struct wgdata_t:
+    #     int ngen
+    #     size_t len, max
+    #     matrix_t **b
+    # wgdata_t *WGInit(int ngen, matrix_t *gen[])
+    # int WGFree(wgdata_t *b)
+    # matrix_t *MakeWord(wgdata_t *b, long n)
+    # void makefp(wgdata_t *b, long fp[])
+    # char *SymbolicName(wgdata_t *b, long n)
+
+    # Permutations
+    # ctypedef struct perm_t:
+    #     Ushort id                     #/* Always = T_PERM */
+    #     long deg
+    #     PTR d
+    # char *permdesc(perm_t *s)
+    # void permprint(char *name, perm_t *x)
+    # perm_t *permalloc(long deg)
+    # void permfree(perm_t *m)
+    # perm_t *permdup(perm_t *src)
+    # perm_t *permmove(perm_t *dest, perm_t *src)
+    # perm_t *permread(FILE *f)
+    # int permwrite(FILE *f, perm_t *perm)
+    # perm_t *permload(char *fn)
+    # int permsave(perm_t *perm, char *fn)
+    # perm_t *permmul(perm_t *dest, perm_t *src)
+    # long permorder(perm_t *perm)
+    # perm_t *permpower(perm_t *p, long n)
+
+    # Insertion
+    # matrix_t *matinsert(matrix_t *mat, poly_t *pol)
+    # matrix_t *matinsert_(matrix_t *mat, poly_t *pol)
+
+    # Factored polynomials
+    # ctypedef struct fpoly_t:
+    #     Ushort id                     #/* Always T_FPOLY */
+    #     size_t len, max
+    #     poly_t **p                    #/* Irreducible factors */
+    #     long *e                       #/* Multiplicities */
+    # char *fpoldesc(fpoly_t *p)
+    # void fpolprint(char *name, fpoly_t *p)
+    # fpoly_t *fpolalloc()
+    # fpoly_t *fpoldup(fpoly_t *s)
+    # void fpolfree(fpoly_t *p)
+    # fpoly_t *fpolread(FILE *f)
+    # int fpolwrite(FILE *f, fpoly_t *p)
+    # fpoly_t *fpolmulp(fpoly_t *dest, poly_t *src, long pwr)
+    # fpoly_t *fpolmul(fpoly_t *dest, fpoly_t *src)
+
+    # Characteristic and minimal polynomials
+    # cdef extern long CharPolSeed
+    # poly_t *charpolfactor(matrix_t *mat)
+    # fpoly_t *charpol(matrix_t *mat)
+    # poly_t *minpolfactor(matrix_t *mat)
+    # fpoly_t *minpol(matrix_t *mat)
+
+    # Berlekamp factorization
+    # fpoly_t *factorization(poly_t *pol)
+
+    # Sets
+    # ctypedef struct set_t:
+    #     Ushort id                     #/* Always = T_SET */
+    #     size_t len, max
+    #     long *buf
+    # int set_allocstrategy(size_t first,size_t blocksize)
+    # char *set_desc(set_t *s)
+    # void set_print(char *name, set_t *x)
+    # set_t *set_alloc()
+    # set_t *set_dup(set_t *s)
+    # void set_free(set_t *x)
+    # int set_insert(set_t *set, long elem)
+    # int set_contains(set_t *set, long elem)
+    # set_t *set_read(FILE *f)
+    # int set_write(FILE *f, set_t *set)
+
+    # Sequences
+    # ctypedef struct sequence_t:
+    #     Ushort id                     #/* Always = T_SEQUENCE */
+    #     size_t len, max
+    #     void **buf
+    # char *seq_desc(sequence_t *s)
+    # void seq_print(char *name, sequence_t *s)
+    # sequence_t *seq_alloc(size_t size)
+    # sequence_t *seq_dup(sequence_t *s)
+    # void seq_free(sequence_t *x)
+    # int seq_insert(sequence_t *s, int pos, void *x)
+    # int seq_remove(sequence_t *s, int pos)
+    # sequence_t *seq_read(FILE *f)
+    # int seq_write(FILE *f, sequence_t *x)
+
+    # Integers and strings
+    # ctypedef struct integer_t:
+    #     Ushort id                     #/* Always = T_INT */
+    #     long l
+    # integer_t *intalloc(long l)
+    # integer_t *intdup(integer_t *s)
+    # char *intdesc(integer_t *s)
+    # void intprint(char *name, integer_t *x)
+    # void intfree(integer_t *x)
+    # integer_t *intread(FILE *f)
+    # int intwrite(FILE *f, integer_t *x)
+    # integer_t *intadd(integer_t *d, integer_t *s)
+    # integer_t *intsub(integer_t *d, integer_t *s)
+    # integer_t *intmul(integer_t *d, integer_t *s)
+    # integer_t *intdiv(integer_t *d, integer_t *s)
+    # integer_t *intneg(integer_t *x)
+
+    # ctypedef struct string_t:
+    #     Ushort id                     #/* Always = T_STR */
+    #     char *s
+    # string_t *stringalloc(char *c)
+    # string_t *stringdup(string_t *s)
+    # char *stringdesc(string_t *s)
+    # void stringprint(char *name, string_t *x)
+    # void stringfree(string_t *x)
+    # string_t *stringread(FILE *f)
+    # int stringwrite(FILE *f, string_t *x)
+
+    # MeatAxe object methods
+    # void mtxfree(void *x)
+    # void *mtxdup(void *x)
+    # char *mtxgetname(void *x)
+    # char *mtxgetdesc(void *x)
+    # void *mtxread(FILE *f)
+    # void mtxprint(void *x)
+    # int mtxwrite(FILE *f, void *x)
+    # void *mtxadd(void *l, void *r)
+    # void *mtxsub(void *l, void *r)
+    # void *mtxmul(void *l, void *r)
+    # void *mtxdiv(void *l, void *r)
+    # void *mtxpwr(void *l, void *r)
+    # void *mtxneg(void *x)
+    # void *mtxorder(void *x)
+
+    # Profiling
+    # cdef extern long ProfFileIO               #/* File I/O */
+    # cdef extern long ProfNullSpace    #/* Matrix null space */
+    # cdef extern long ProfSpinUp               #/* Spin-up */
+    # cdef extern long ProfMakeWord     #/* Make word */
+    # cdef extern long ProfMatInsert    #/* Matrix insertion */
+    # cdef extern long ProfPolFactor    #/* Polynomial factorization */
+    # cdef extern long ProfCharPol      #/* Char. and minimal polynomial */
+
+    # Added by Simon King:
+    cdef char* mat2str(matrix_t* m)           # Convert m->d into a char*
+    cdef void str2mat(matrix_t* m, long l, char* x) # *COPY* the data from x to m->d
+    cdef char* zrow2str(PTR d)       # Convert d into a char*
+    cdef void zstr2row(PTR d, char* x) # *COPY* the data from x to d, assuming
+                                       # that d has been assigned to at least
+                                       # zrowsize_io byte.
+
+    cdef int multiply_strassen(matrix_t* A, matrix_t* B, matrix_t *C)
+    cdef void set_cutoff(long rows, long size)
+
+# thats's all of meataxe.h !
+
+##############################################################
+
+#
+# In the .pxd-declaration file, we merely need to define the 
+# cdef-parts of the extension class. Everything else is 
+# in the .pyx-definition file.
+
+from sage.matrix.matrix_dense cimport Matrix_dense
+from sage.matrix.matrix0 cimport Matrix
+from sage.structure.element cimport Element, ModuleElement, RingElement
+cimport sage.structure.element
+
+cdef class FieldConverter_class:
+    cdef object field
+    cpdef inline object int_to_field(self, int x)
+    cpdef inline int field_to_int(self, x)
+
+cdef class Matrix_modpn_dense(Matrix_dense):
+    cdef matrix_t *Data
+    cdef FieldConverter_class _converter
+    cpdef Matrix_modpn_dense _mulInt_(self, long right)
+    cpdef int base(self)
+    cpdef int characteristic(self)
+    cpdef list _matlist_(self)
+    cpdef Matrix_modpn_dense transpose(Matrix_modpn_dense self)
+    cpdef int order(Matrix_modpn_dense self) except -1
+    cpdef Matrix_modpn_dense normalized(Matrix_modpn_dense self)
+    cpdef Matrix_modpn_dense semi_echelon(Matrix_modpn_dense self)
+    cpdef int nullity(Matrix_modpn_dense self)
+    cpdef tuple lead(self)
+    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value)
+    cdef int get_unsafe_int(self, Py_ssize_t i, Py_ssize_t j)

new file sage/matrix/matrix_modpn_dense.pyx

@@ +1 @@
+#*****************************************************************************
+#
+#    Dense matrices over GF(p^n), p^n<255, using LibMeataxe as backend
+#
+#    Copyright (C) 2011 Simon A. King <simon.king@uni-jena.de>
+#
+#  Distributed under the terms of the GNU General Public License (GPL),
+#  version 2 or later (at your choice)
+#
+#    This code is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#    General Public License for more details.
+#
+#  The full text of the GPL is available at:
+#
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+r"""
+Dense Matrices over `\mathbb F_q`, with `q<255` odd and not prime
+
+This module is a wrapper for a major modification of version 2.2.4 of the
+Aachen `C-MeatAxe <http://www.math.rwth-aachen.de/homes/MTX/download.html>`_
+and provides matrices over the finite field `\mathbb F_q`, where
+`q\le 255` is odd and not prime.
+
+The modification strips down MeatAxe, which is a collection of executables,
+to a C-library for linear algebra over fields of order at most 255. In
+addition to some smaller improvements, the modifications include an
+implementation of Winograd-Strassen multiplication.
+
+AUTHORS:
+
+- Simon King <simon.king@uni-jena.de>
+
+"""
+
+####################
+# 
+# import sage types
+#
+####################
+include "../ext/python.pxi"
+include "../ext/python_object.pxi"
+include "../ext/interrupt.pxi"
+include "../ext/stdsage.pxi"
+include "../ext/python_slice.pxi"
+
+import sage
+import sage.all
+from warnings import showwarning
+import sys
+from sage.structure.mutability import Mutability
+from sage.all import Integer
+from sage.all import cputime
+from sage.all import GF
+from sage.all import is_prime_power
+from sage.all import factor
+from sage.all import save, copy
+from sage.all import cached_method, cached_function
+from sage.all import SAGE_ROOT
+from sage.all import current_randstate
+
+## The tables for MeatAxe computations will be found in the package data folder
+import os
+os.environ['MTXLIB'] = SAGE_ROOT + '/local/share/meataxe/'
+
+#########################################################################
+#########################################################################
+# 
+# Definition of the Matrix_modpn_dense extension type
+# 
+#########################################################################
+#########################################################################
+from sage.rings.finite_rings.integer_mod import IntegerMod_int
+
+# Fast conversion from field to int and int to field
+cdef class FieldConverter_class:
+    """
+    An auxiliary class, used to convert between <int> and finite field element
+
+    This class is for non-prime fields only. The method
+    :meth:`int_to_field` exists for speed. The method
+    :meth:`field_to_int` exists in order to have a common interface
+    for elements of prime and non-prime fields; see
+    :class:`PrimeFieldConverter_class`.
+
+    EXAMPLE::
+
+        sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
+        sage: F.<y> = GF(125)
+        sage: C = FieldConverter_class(F)
+        sage: C.int_to_field(15)
+        3*y
+        sage: F.fetch_int(15)
+        3*y
+        sage: %timeit C.int_to_field(15)    #not tested
+        625 loops, best of 3: 1.04 µs per loop
+        sage: %timeit F.fetch_int(15)       #not tested
+        625 loops, best of 3: 3.97 µs per loop
+        sage: C.field_to_int(y)
+        5
+        sage: y.integer_representation()
+        5
+
+    """
+    def __init__(self, field):
+        """
+        INPUT:
+
+        A finite *non-prime* field. This assumption is not tested.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
+            sage: F.<y> = GF(125)
+            sage: C = FieldConverter_class(F)
+            sage: C.int_to_field(15)
+            3*y
+            sage: F.fetch_int(15)
+            3*y
+            sage: C.field_to_int(y)
+            5
+            sage: y.integer_representation()
+            5
+
+        """
+        self.field = field._cache.fetch_int
+    cpdef inline object int_to_field(self, int x):
+        """
+        Fetch a python int into the field.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
+            sage: F.<y> = GF(125)
+            sage: C = FieldConverter_class(F)
+            sage: C.int_to_field(15)
+            3*y
+            sage: F.fetch_int(15)
+            3*y
+
+        """
+        return self.field(x)
+    cpdef inline int field_to_int(self, x):
+        """
+        Represent a field element by a python int.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
+            sage: F.<y> = GF(125)
+            sage: C = FieldConverter_class(F)
+            sage: C.field_to_int(y)
+            5
+            sage: y.integer_representation()
+            5
+
+        """
+        return x.integer_representation()
+
+cdef class PrimeFieldConverter_class(FieldConverter_class):
+    """
+    An auxiliary class, used to convert between <int> and finite field element
+
+    This class is for prime fields only. The methods
+    :meth:`int_to_field` and :meth:`field_to_int` exist in order to
+    have a common interface for elements of prime and non-prime fields;
+    see :class:`FieldConverter_class`.
+
+    EXAMPLE::
+
+        sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
+        sage: F = GF(5)
+        sage: C = PrimeFieldConverter_class(F)
+        sage: C.int_to_field(int(2))
+        2
+        sage: F(2)
+        2
+        sage: C.field_to_int(F(2))
+        2
+        sage: int(F(2))
+        2
+
+    """
+    def __init__(self, field):
+        """
+        INPUT:
+
+        A finite *prime* field. This assumption is not tested.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
+            sage: F = GF(5)
+            sage: C = PrimeFieldConverter_class(F)
+            sage: C.int_to_field(int(2))
+            2
+            sage: F(2)
+            2
+            sage: C.field_to_int(F(2))
+            2
+            sage: int(F(2))
+            2
+
+        """
+        self.field = field
+    cpdef inline object int_to_field(self, int x):
+        """
+        Fetch a python int into the field.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
+            sage: F = GF(5)
+            sage: C = PrimeFieldConverter_class(F)
+            sage: C.int_to_field(int(2))
+            2
+            sage: F(2)
+            2
+
+        """
+        return IntegerMod_int(self.field, x)
+    cpdef inline int field_to_int(self, x):
+        """
+        Represent a field element by a python int.
+
+        EXAMPLE::
+
+            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
+            sage: F = GF(5)
+            sage: C = PrimeFieldConverter_class(F)
+            sage: C.field_to_int(F(2))
+            2
+            sage: int(F(2))
+            2
+
+        """
+        return int(x)
+
+cdef public dict _converter_cache = {}
+cdef FieldConverter_class FieldConverter(field):
+    """
+    Return a :class:`FieldConverter_class` or :class:`PrimeFieldConverter_class` instance,
+    depending whether the field is prime or not.
+
+    EXAMPLE::
+
+        sage: MS = MatrixSpace(GF(5^3,'y'),2)
+        sage: A = MS.random_element()
+        sage: A*2 == A+A    # indirect doctest
+        True
+        sage: MS = MatrixSpace(GF(5),2)
+        sage: A = MS.random_element()
+        sage: A*2 == A+A
+        True
+
+    """
+    try:
+        return _converter_cache[field]
+    except KeyError:
+        if field.is_prime_field():
+            return _converter_cache.setdefault(field, PrimeFieldConverter_class(field))
+        return _converter_cache.setdefault(field, FieldConverter_class(field))
+
+class MTX_unpickle_class:
+    """
+    Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
+
+    EXAMPLES::
+
+        sage: M = MatrixSpace(GF(25,'z'),2)([1,2,3,4])
+        sage: M == loads(dumps(M))   # indirect doctest
+        True
+
+    """
+    def __init__(self):
+        """
+        Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
+
+        TESTS::
+
+            sage: from sage.matrix.matrix_modpn_dense import MTX_unpickle_class
+            sage: unpickle = MTX_unpickle_class()   # indirect doctest
+            sage: M = MatrixSpace(GF(25,'z'),2)([1,2,3,4])
+            sage: L = M.__reduce__()[1]
+            sage: N = unpickle(*L)
+            sage: print N
+            [1 2]
+            [3 4]
+
+        """
+        self.__safe_for_unpickling__=True
+
+    def __call__(self,f,nr,nc,object Data,m):
+        """
+        Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
+
+        TESTS::
+
+            sage: from sage.matrix.matrix_modpn_dense import MTX_unpickle_class
+            sage: unpickle = MTX_unpickle_class()
+            sage: M = MatrixSpace(GF(125,'z'),1,4)([1,2,3,4])
+            sage: L = M.__reduce__()[1]
+            sage: N = unpickle(*L)   # indirect doctest
+            sage: print N
+            [1 2 3 4]
+
+        """
+        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        cdef char* x
+        #cdef list L
+        cdef int i
+        cdef long ff
+        cdef int NR
+        cdef PTR p
+        if isinstance(f, (int,long)):
+            ff = f
+            B = GF(f,'z')
+            # TODO: MatrixSpace with "implementation"
+            from sage.all import MatrixSpace
+            parent = MatrixSpace(B, nr, nc)
+        else:
+            parent = f
+            B = parent.base_ring()
+            ff = B.order()
+        OUT._parent = parent
+        OUT._base_ring = B
+        OUT._converter = FieldConverter(B)
+        OUT._ncols = nc
+        OUT._nrows = nr
+        OUT._cache = {}
+        if Data:
+            OUT.Data = matalloc(ff, nr, nc)
+            if not OUT.Data:
+                raise MemoryError, "Not enough memory to allocate %dx%d MTY matrix"%(nr,nc)
+            x=PyString_AsString(Data)
+            str2mat(OUT.Data,len(Data),x)
+        OUT._mutability = Mutability(m)
+        return OUT
+
+mtx_unpickle = MTX_unpickle_class()
+
+
+cdef class Matrix_modpn_dense(Matrix_dense):
+    r"""
+    Dense matrices over `\mathbb F_q`, `q<255` odd and not prime.
+
+    NOTE:
+
+    This class uses a major modification of the Aachen C-MeatAxe
+    as backend. In principle, it would also work for prime fields
+    and in characteristic two. However, other matrices in Sage,
+    relying on linbox, m4ri or m4rie, are more efficient in these
+    cases.
+
+    EXAMPLES::
+
+        sage: M = MatrixSpace(GF(25,'z'),2,3)([1,2,3,4,5,6])
+        sage: print M
+        [1 2 3]
+        [4 0 1]
+        sage: type(M)
+        <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+
+    The documentation of the ``__init__`` methods shows further
+    ways of creating a :class:`Matrix_modpn_dense` instance.
+    However, these should only be of internal use.
+
+    """
+
+###################
+## Init, Dealloc, Copy
+    def __cinit__(self, parent=None, entries=None, **kwds):
+        """
+        TESTS::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: Matrix_modpn_dense.__new__(Matrix_modpn_dense)   # indirect doctest
+            []
+            sage: Matrix_modpn_dense(MatrixSpace(GF(64,'z'),4), None)
+            [0 0 0 0]
+            [0 0 0 0]
+            [0 0 0 0]
+            [0 0 0 0]
+
+        """
+        if parent is None:  # this makes Matrix_modpn_dense.__new__(Matrix_modpn_dense) work,
+                            # returning a non-initialised matrix
+            return
+        if isinstance(parent,basestring): # this allows to provide a file when initialising a matrix
+            return
+        cdef long f = parent.base_ring().order()
+        cdef long nrows = parent.nrows()
+        cdef long ncols = parent.ncols()
+        if entries is None or entries==0:
+            self.Data = matalloc(f, nrows, ncols)
+            return
+        if isinstance(entries, list):
+            self.Data = matalloc(f, nrows, ncols)
+            return
+        if nrows!=ncols:
+            raise TypeError, "Matrix ought to be a square matrix"
+        self.Data = MTXmatid(f, <long>nrows)
+
+    def __init__(self, parent, entries=None, mutable=True, copy=False, coerce=False):
+        """
+        Matrix extension class using libmeataxe as backend
+
+        INPUT:
+
+        Instances of this class can be created by providing one of
+        the following input data, where ``q<255`` is a prime power,
+        ``m,n`` are non-negative integers, and `a_{11},...,a_{mn}`
+        can be coerced into ``GF(q)``. Note that a user should
+        create these instances via the matrix constructors; what
+        we explain here is for internal use only!
+
+        - None => empty matrix over an unspecified field (used for unpickling)
+        - a string ``f`` ==> load matrix from the file named ``f``
+        - A matrix space of `m\\times n` matrices over GF(q) and either
+
+          - a list `[a_{11},a_{12},...,a_{1n},a_{21},...,a_{m1},...,a_{mn}]`,
+            which results in a matrix with the given marks
+          - ``None``, which is the fastest way to creata a zero matrix.
+          - an element of GF(q), which results in a diagonal matrix with the
+            given element on the diagonal.
+
+        If the optional parameter ``mutable`` is ``False`` (by default,
+        it is ``True``), the resulting matrix can not be changed, and
+        it can be used as key in a Python dictionary.
+
+        NOTE:
+
+        If ``q`` is not a prime power, a ``MemoryError`` will be raised. 
+
+        EXAMPLES::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+
+        1. Creating an empty matrix::
+
+            sage: Matrix_modpn_dense(None)
+            []
+
+        2. Creating a zero (3x2)-matrix::
+
+            sage: Matrix_modpn_dense(MatrixSpace(GF(4,'z'),3,2))
+            [0 0]
+            [0 0]
+            [0 0]
+
+        3. Creating a matrix from a list or list of lists::
+
+            sage: Matrix_modpn_dense(MatrixSpace(GF(5),2,3),[1,2,3,4,5,6])
+            [1 2 3]
+            [4 0 1]
+            sage: Matrix_modpn_dense(MatrixSpace(GF(5),2,3),[[1,2,3],[4,5,6]])  # indirect doctest
+            [1 2 3]
+            [4 0 1]
+
+        4. Creating a diagonal matrix::
+
+            sage: M = Matrix_modpn_dense(MatrixSpace(GF(7),5),2); M
+            [2 0 0 0 0]
+            [0 2 0 0 0]
+            [0 0 2 0 0]
+            [0 0 0 2 0]
+            [0 0 0 0 2]
+
+        5. Creating a matrix from a file in MeatAxe format.
+
+           First, we have to create that file; we use a temporary file,
+           that will be removed when leaving Sage. Note that the method
+           :meth:`msave` must be used, which does not use Python pickling
+           but relies on the intrinsic C--MeatAxe way of saving.
+           ::
+
+            sage: f = tmp_filename()
+            sage: M.msave(f)
+            sage: Matrix_modpn_dense(f)
+            [2 0 0 0 0]
+            [0 2 0 0 0]
+            [0 0 2 0 0]
+            [0 0 0 2 0]
+            [0 0 0 0 2]
+
+        TESTS::
+
+            sage: MS = MatrixSpace(GF(125,'y'),2)
+            sage: A = MS(0)
+            sage: A.left_kernel()
+            Vector space of degree 2 and dimension 2 over Finite Field in y of size 5^3
+            Basis matrix:
+            [1 0]
+            [0 1]
+            sage: A.right_kernel()
+            Vector space of degree 2 and dimension 2 over Finite Field in y of size 5^3
+            Basis matrix:
+            [1 0]
+            [0 1]
+
+        """
+        if parent is None:
+            self._mutability = Mutability(False)
+            self._ncols = 0
+            self._nrows = 0
+            self._cache = {}
+            return
+        if isinstance(parent, basestring): # load from file
+            FILE = os.path.realpath(parent)
+            try:
+                fsock = open(FILE,"rb",0)
+                fsock.close()
+            except (OSError,IOError):
+                return
+            self.Data = matload(FILE)
+            if zsetfield(self.Data.fl):
+                raise ValueError, "Invalid saved date"
+            B = GF(self.Data.fl, 'z')
+            # TODO: MatrixSpace with "implementation"
+            from sage.all import MatrixSpace
+            parent = MatrixSpace(B, self.Data.nor, self.Data.noc)
+            self._mutability = Mutability(False)
+            self._parent = parent
+            self._base_ring = B
+            self._converter = FieldConverter(B)
+            self._ncols = self.Data.noc
+            self._nrows = self.Data.nor
+            self._cache = {}
+            return
+
+        if not self.Data: # should be initialised by __cinit__
+            raise MemoryError, "Error allocating memory for MeatAxe matrix"
+        Matrix_dense.__init__(self, parent)
+        self._mutability=Mutability(not mutable)
+        B = self._base_ring
+        self._converter = FieldConverter(B)
+        if entries is None:
+            return
+
+        cdef int i,j
+        cdef int f
+        if not isinstance(entries,list): # __cinit__ initialized the unit matrix
+            f = self._converter.field_to_int(self._coerce_element(entries))
+            if matmulF(self.Data,zitof(f)):
+                return
+            else:
+                raise ArithmeticError, "Matrix sizes or fields not compatible"
+        cdef PTR x
+        cdef PTR y
+        x = self.Data.d
+        cdef int nr = self.Data.nor
+        cdef int nc = self.Data.noc
+        if nr==0 or nc==0:
+            return
+        if len(entries)<nr:
+            raise ValueError, "Expected a list of size at least the number of rows"
+        cdef int idx
+        cdef list dt, dt_i
+        if isinstance(entries[0],list):
+            dt = entries
+            for i from 0 <= i < nr:
+                y = x
+                idx = 0
+                dt_i = dt[i]
+                #dtnext = dt[i].__iter__().next
+                for j from 0 <= j < nc:
+                    y[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(dt_i[j])))]
+                    idx += 1
+                    if idx == MPB:
+                        y+=1
+                        idx=0
+                    # zinsert_step is not the fastest, since tnull[*idx] and tinsert[*idx] are looked up over and over again.
+                zsteprow(&(x))
+        else:
+            dtnext = entries.__iter__().next
+            for i from 0 <= i < nr:
+                y = x
+                idx = 0
+                for j from 0 <= j < nc:
+                    bla = dtnext()
+
+                    y[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(bla)))] #dtnext())))]
+                    idx += 1
+                    if idx == MPB:
+                        y+=1
+                        idx=0
+                zsteprow(&(x))
+                            
+    def __dealloc__(self):
+        """
+        Deallocate this matrix
+
+        TESTS::
+
+            sage: M = MatrixSpace(GF(25,'z'),2,4)([1,2,3,4,4,3,2,1])
+            sage: del M    # indirect doctest
+            sage: M
+            Traceback (most recent call last):
+            ...
+            NameError: name 'M' is not defined
+        """
+        if self.Data:
+            matfree(self.Data)
+        #cdef matrix_t *p
+        self.Data = NULL
+
+    def __copy__(self):
+        """
+        Return a mutable copy of this matrix
+
+        EXAMPLES:
+
+        The following example caused a problem in a very
+        early stage of wrapping C-MeatAxe::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: M = Matrix_modpn_dense(MatrixSpace(GF(2),3,20), [20*[0],20*[0],[1]+19*[0]], mutable=False)
+            sage: M.is_mutable()
+            False
+            sage: N = copy(M)    # indirect doctest
+            sage: N.is_mutable()
+            True
+            sage: print N
+            [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
+            [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
+            [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
+            sage: N==M
+            True
+            sage: N is M
+            False
+
+        """        
+        cdef Matrix_modpn_dense retval = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        # Do the initialisation "manually"
+        retval._mutability = Mutability(False)
+        retval._parent = self._parent
+        retval._base_ring = self._base_ring
+        retval._converter = self._converter
+        retval._ncols = self._ncols
+        retval._nrows = self._nrows
+        retval._cache = {}
+        if self.Data:
+            retval.Data = matdup(self.Data)
+            if not retval.Data:
+                raise MemoryError, "Error copying a %s instance"%repr(type(self))
+        return(retval)
+    
+    def __reduce__(self):
+        r"""
+        Pickle/save this matrix.
+
+        TESTS::
+
+            sage: M = MatrixSpace(GF(25,'z'),1,4)([1,2,3,4])
+            sage: M.__reduce__()[0] 
+            <sage.matrix.matrix_modpn_dense.MTX_unpickle_class instance at ...>
+            sage: M.__reduce__()[1]
+            (Full MatrixSpace of 1 by 4 dense matrices over Finite Field in z of size 5^2, 1, 4, '\x01\x02\x03\x04', False)
+            sage: loads(dumps(M))
+            [1 2 3 4]
+            sage: loads(dumps(M)) == M  # indirect doctest
+            True
+
+        """
+        cdef char* d
+        cdef L = []
+        Lnext = L.append
+        cdef int i,NR
+        cdef PTR p
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+
+            NR = self.Data.nor
+            p = self.Data.d
+            for i from 0 <= i < NR:
+                d = zrow2str(p)
+                Lnext(PyString_FromStringAndSize(d,zrowsize_io))
+                zsteprow(&(p))
+            return mtx_unpickle, (self._parent,self.Data.nor,self.Data.noc,''.join(L),self._mutability._is_immutable) # or ...,L,self._mutable)
+        else:
+            return mtx_unpickle, (0,0,0,'',self._mutability._is_immutable)
+
+
+#######################
+## String representation is taken care of by implementing get_unsafe
+
+
+##################
+## ==,<,>, dictionary
+    def __richcmp__(left, right, int op):
+        return (<Element>left)._richcmp(right, op)
+    
+    def __cmp__(left, right):
+        return (<Element>left)._cmp(right)
+
+    cdef int _cmp_c_impl(self, Element right) except -2:
+        """
+        Compare two MeatAxe matrices
+        
+        Of course, '<' and '>' doesn't make much sense for matrices.
+
+        EXAMPLES::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: M = Matrix_modpn_dense(MatrixSpace(GF(2), 3,20), [20*[0],20*[0],[1]+19*[0]])
+            sage: N = copy(M)
+            sage: M == N   # indirect doctest
+            True
+            sage: M != N
+            False
+            sage: M<N
+            False
+            sage: N[2,19] = 1
+            sage: M == N
+            False
+            sage: M != N
+            True
+        """
+        return matcmp((<Matrix_modpn_dense>self).Data, (<Matrix_modpn_dense>right).Data)
+
+    def __hash__(self):
+        """
+        Return a hash value for M, provided M is immutable
+
+        EXAMPLES::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: M = Matrix_modpn_dense(MatrixSpace(GF(5),1,4),[[4,3,2,1]],mutable=False)
+            sage: if sys.byteorder == 'little':
+            ...     print hash(M) == 7606091044269354279  # indirect doctest
+            ... else:
+            ...     print hash(M) == 1060097699   # indirect doctest
+            True
+
+        """
+        cdef char* d
+        if not self.Data:
+            return 0
+        else:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+            d = mat2str(self.Data)
+            return hash((self.Data.fl,self.Data.nor,self.Data.noc,PyString_FromStringAndSize(d,zsize(self.Data.nor))))
+
+##########################
+## Saving should be done via pickling
+## However, we kept a method that relies on MeatAxe matsave:
+    def msave(self,f):
+        """
+        Save self into file, using the original C--MeatAxe format
+     
+        INPUT:
+
+        ``f`` - a file name (string)
+
+        NOTE: 
+
+        A matrix ``M`` should usually be saved by ``save(M,'filename')``
+        and reloaded with ``load('filename')``. Saving with
+        ``M.msave('filename')`` uses a totally different format, so that
+        the data can not be reloaded with ``load`` but with
+        ``Matrix_modpn_dense('filename')``.
+
+        EXAMPLES:
+
+        First, we create a temporary file name, that will be removed
+        when leaving Sage. Of course, in an application a permanent
+        file would usually be chosen.
+        ::
+        
+            sage: f = tmp_filename()+'.sobj'
+
+        Next, we create a matrix and use msave::
+        
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: M = MatrixSpace(GF(25,'z'),2,4)([1,2,3,4,4,3,2,1])
+            sage: M.msave(f)
+
+        It can not be reloaded with 'load'::
+        
+            sage: N = load(f)
+            Traceback (most recent call last):
+            ...
+            UnpicklingError: invalid load key, ...
+
+        But it can be reloaded with the class constructor::
+        
+            sage: N = Matrix_modpn_dense(f)
+            sage: M == N
+            True
+        """
+        matsave(self.Data,f)
+
+##########################
+## Structural parts
+
+    cpdef int base(self):
+        """
+        Return the order of the field over which the matrix is defined
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(25,'z'),1,4)([4,3,2,1])
+            sage: M.base()
+            25
+        """
+        if self.Data:
+            return self.Data.fl
+        else:
+            return None
+ 
+    cpdef int characteristic(self):
+        """
+        Return the characteristic of the field over which the matrix is defined
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(25,'z'),1,4)([4,3,2,1])
+            sage: M.characteristic()
+            5
+        """
+        if self.Data:
+            return factor(self.Data.fl)[0][0]
+        else:
+            return None
+
+    def matrix_from_rows(self, rows):
+        """
+        Matrix formed by a given list of rows from self.
+
+        INPUT:
+
+        ``rows`` -- list of integers, providing the rows that shall be put together
+
+        OUTPUT:
+
+        A matrix that is formed by the indicated rows of self, in the given order.
+
+        EXAMPLES::
+
+            sage: MS = MatrixSpace(GF(5^3,'y'),3,2)
+            sage: A = MS(range(6))
+            sage: A
+            [0 1]
+            [2 3]
+            [4 0]
+            sage: A.matrix_from_rows([2,1])
+            [4 0]
+            [2 3]
+            sage: A.matrix_from_rows([1,1])
+            [2 3]
+            [2 3]
+
+        """
+        if not (PY_TYPE_CHECK(rows, list) or PY_TYPE_CHECK(rows, tuple)):
+            raise TypeError, "rows must be a list of integers"
+        cdef Matrix_modpn_dense A
+        cdef int k
+        cdef Py_ssize_t nrows,c
+
+        nrows = PyList_GET_SIZE(rows)
+        A = self.new_matrix(nrows = nrows)
+
+        cdef PTR src,dest
+        dest = A.Data.d
+        for k in rows:
+            src = self.Data.d
+            zadvance(&(src),k)
+            zcpyrow(dest,src,zrowsize)
+            zsteprow(&(dest))
+        return A
+
+    def __getitem__(self,key):
+        """
+        Return element, row, or slice of self.
+
+        NOTE:
+
+        We could simply inherit the __getitem__ method form
+        the base class of dense matrices. Hover, this would
+        be a lot slower.
+
+        INPUT:
+
+        - a tuple (i,j) with i, j non-negative integers, or
+        - a non-negative integer, or
+        - a slice object, created via [i:j]
+
+        USAGE:
+
+        - A[i, j] -- the i,j of A, or
+        - A[i]    -- the i-th row of A, or
+        - A[i:j]  -- the i-th through (j-1)-st rows of A.
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(9,'z'),3)(list(GF(9,'z')))
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M[1,2]    # indirect doctest
+            z
+            sage: M[1]
+            (z + 2, 2, z)
+            sage: M[1:3]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+        """
+        if not self.Data:
+            raise IndexError, "Matrix is empty"
+        zsetfield(self.Data.fl)
+        cdef matrix_t *mtxrow
+        cdef PTR p
+        cdef Matrix_modpn_dense A
+        cdef int i,j
+        cdef tuple key_tuple
+        cdef int nr=self.Data.nor
+        cdef int nc=self.Data.noc
+        zsetlen(nc)
+        if PyTuple_CheckExact(key):
+            # key is a tuple, so we get i and j efficiently
+            # TODO: Support pairs of slices
+            key_tuple = <tuple>key
+            if len(key_tuple) != 2:
+                raise IndexError("index must be an integer or pair of integers")
+            try:
+                i = <object>PyTuple_GET_ITEM(key_tuple, 0)
+                j = <object>PyTuple_GET_ITEM(key_tuple, 1)
+            except TypeError:
+                raise IndexError("index must be an integer or pair of integers")
+            if (i<0):
+                i += nr
+            if (i<0) or (i>=nr):
+                raise IndexError, "matrix index out of range"
+            if j<0:
+                j+=nc
+            if j<0 or j>=nc:
+                raise IndexError, "matrix index out of range"
+            p = self.Data.d
+            zadvance(&(p),i)
+            return self._converter.int_to_field(zftoi(zextract(p,j+1)))
+        if PySlice_Check(key):
+            # Slice interpretation is passed to the sequence
+            # constructed by range
+            if (key.start>=0) and (key.stop<=nr) and \
+               (key.start<=nr) and (key.stop>=0):
+                A=Matrix_modpn_dense(None)
+                A.Data = _matextract(self.Data, long(key.start+1), long(key.stop-key.start))
+                if A.Data:
+                    A._mutability = copy(self._mutability)
+                    from sage.all import MatrixSpace
+                    A._parent = MatrixSpace(self._base_ring, key.stop-key.start, self._ncols)
+                    A._converter = self._converter
+                    A._base_ring = self._base_ring
+                    A._ncols = self._ncols
+                    A._nrows = key.stop-key.start
+                    A._cache = {}
+                    return A
+                else:
+                    raise IndexError, "matrix index out of range"
+            else:
+                raise IndexError, "matrix index out of range"
+        try:
+            if (key>=0) and (key<nr):
+                return self.row(key)  # this is inherited
+        except:
+            raise IndexError, "matrix index out of range"
+        raise IndexError, "matrix index out of range"
+
+    def __setitem__(self, key, el):
+        """
+        x.__setitem__(key, el) <==> x[key]=el
+
+        USAGE:
+
+        - key= pair of integers, el= integer
+        - key= integer, el= list/tuple of integers 
+          row number <key> is filled with the entries of <el>,
+          until the end of the row or the end of <el> is reached
+        - key= integer, el= Matrix_modpn_dense with the same number of columns as self
+          The rows of self starting with number <key> are filled with the 
+          contents of el.
+
+        EXAMPLES::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F)); M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M[0] = [1,2,3]   # indirect doctest
+            sage: M
+            [      1       2       0]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M[1,2] = 0
+            sage: M
+            [      1       2       0]
+            [  z + 2       2       0]
+            [2*z + 2 2*z + 1       1]
+            
+            sage: N = MatrixSpace(F,2,3)([z,z^2,z^3,z^3,z^4,z^5])
+            sage: N
+            [      z   z + 1 2*z + 1]
+            [2*z + 1       2     2*z]
+            sage: M[1] = N
+            sage: M
+            [      1       2       0]
+            [      z   z + 1 2*z + 1]
+            [2*z + 1       2     2*z]
+
+        """
+        self.check_mutability()
+        cdef PTR x
+        cdef int j,mini,idx
+        cdef Matrix_modpn_dense tmp
+        cdef list value_list
+        if self.Data:
+            zsetfield(self.Data.fl)
+        else:
+            raise IndexError, "Matrix is empty"
+        if isinstance(key,tuple):
+            try:
+                j,idx = key
+            except (TypeError, ValueError):
+                raise TypeError, "pair of integers expected"
+            if (0<=j<self.Data.nor) and (0<=idx<self.Data.noc):
+                x = self.Data.d
+                zsetlen(self.Data.noc)
+                zadvance(&(x),j)
+                zinsert(x, idx+1, zitof(self._converter.field_to_int(self._coerce_element(el))))
+            else:
+                raise IndexError, "index out of range"
+        elif (isinstance(key,Integer) or isinstance(key,int)): # insert a list or another MTX
+            j = key
+            if (j>=0) and (j<self.Data.nor):
+                x = self.Data.d
+                zsetlen(self.Data.noc)
+                zadvance(&(x),j)
+                if isinstance(el,Matrix_modpn_dense):
+                    tmp = el
+                    if tmp.Data.fl!=self.Data.fl:
+                        raise ArithmeticError, "Fields are incompatible"
+                    if tmp.Data.nor+j>self.Data.nor:
+                        raise IndexError, "Matrix to be inserted has too many rows"
+                    if tmp.Data.noc!=self.Data.noc:
+                        raise IndexError, "Matrix to be inserted must have the same number of columns as self"
+                    zsetlen(tmp.Data.noc)
+                    zcpyrow(x,tmp.Data.d,zsize(tmp.Data.nor))
+                    zsetlen(self.Data.noc)
+                    return
+                value_list = el  # type error, if no list.
+                if len(value_list)!=self.Data.noc:
+                    raise ValueError, "List length differs from number of columns"
+                idx = 0
+                memset(x, 0, zrowsize)
+                for v in value_list:
+                    x[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(v)))]
+                    idx += 1
+                    if idx == MPB:
+                        x+=1
+                        idx=0
+            else:
+                raise IndexError, "index out of range"
+        else:
+            raise IndexError, "index must be an integer or a pair of integers"
+
+    def _rowlist_(self, i, j=-1):
+        """
+        M._rowlist_(i): Return row <i> as a list of integers
+        M._rowlist_(i,j): Return rows <i> up to and including <j> in a single list
+        
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: M._rowlist_(0)
+            [0, 6, 4]
+
+        We verify that these integers really correspond
+        to the actual elements of the top row::
+
+            sage: M[0]
+            (0, 2*z, z + 1)
+            sage: (2*z).integer_representation()
+            6
+            sage: (z+1).integer_representation()
+            4
+            sage: M._rowlist_(1,2)
+            [5, 2, 3, 8, 7, 1]
+        """
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+        else:
+            return []
+        if not(isinstance(i,int) or isinstance(i,Integer)):
+            raise TypeError, "Index must be an integer"
+        if (i<0) or (i>=self.Data.nor):
+            raise IndexError, "Index out of range"
+        cdef PTR p,q
+        cdef int k,l,ii,jj, idx
+        p = self.Data.d
+        zadvance(&(p),i)
+        idx = 0
+        q = p
+        cdef list L = [zftoi(zextract_step(&(q),&(idx))) for k in range(self.Data.noc)]
+        if j!=-1:
+            if not(isinstance(j,int) or isinstance(j,Integer)):
+                raise TypeError, "Second index must be an integer"
+            if j >= self.Data.nor:
+                raise IndexError, "Index out of range"
+            ii=i
+            jj=j
+            for k from ii < k <= jj:
+                zsteprow(&(p))
+                idx = 0
+                q = p
+                for l from 0<=l<self.Data.noc:
+                    L.append(zftoi(zextract_step(&(q),&(idx))))
+        return L
+
+    cpdef list _matlist_(self):
+        """
+        M._matlist_(): Return the entries of M as a list of lists of integers
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: M._matlist_()
+            [[0, 6, 4], [5, 2, 3], [8, 7, 1]]
+        """
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+        else:
+            raise IndexError, "Matrix is empty"
+        cdef PTR p,q
+        cdef int i,j
+        p = self.Data.d
+        cdef int idx
+        cdef list l_out=[]
+        for i from 1<=i<self.Data.nor:
+            q = p
+            idx = 0
+            l_out.append([zftoi(zextract_step(&(q),&(idx))) for j in range(self.Data.noc)])
+            zsteprow(&(p))
+        idx = 0
+        l_out.append([zftoi(zextract_step(&(p),&(idx))) for j in range(self.Data.noc)])
+        return l_out
+
+#########################
+## Arithmetics
+    cpdef ModuleElement _add_(self, ModuleElement other):
+        """
+        Add two MeatAxe matrices of equal size
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(25)
+            sage: MS = MatrixSpace(F, 2,4)
+            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
+            sage: N = MS([1,z,z^2,z^3, z^4,z^5,z^5,z^7])
+            sage: M + N    # indirect doctest
+            [      2   z + 2   z + 1 4*z + 2]
+            [3*z + 2   z + 1 2*z + 1       z]
+            
+        We verify that the answer is correct::
+
+            sage: (M+N).list() == [a+b for a,b in zip(M.list(), N.list())]
+            True
+
+        """
+        cdef Matrix_modpn_dense right = other
+        if ((self.Data.noc==0) or (right.Data.noc==0)):
+            raise NotImplementedError, "The matrices must not be empty"
+        cdef Matrix_modpn_dense left = self.__copy__()
+        if not right._mutability.is_mutable():
+            left._mutability.set_immutable()
+        if matadd(left.Data,right.Data):
+            return left
+        else:
+            raise ArithmeticError, "Matrix sizes or fields not compatible"
+
+    cpdef ModuleElement _sub_(self, ModuleElement other):
+        """
+        Subtract two MeatAxe matrices of equal size
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(25)
+            sage: MS = MatrixSpace(F, 2,4)
+            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
+            sage: N = MS([1,z,z^2,z^3, z^4,z^5,z^5,z^7])
+            sage: M - N    # indirect doctest
+            [      0 4*z + 2     4*z   z + 1]
+            [4*z + 3 3*z + 4 4*z + 4     2*z]
+
+        We verify that the answer is correct::
+
+            sage: (M-N).list() == [a-b for a,b in zip(M.list(), N.list())]
+            True
+
+        """
+        cdef Matrix_modpn_dense right = other
+        if (self.Data.fl==2):
+            return self.__add__(right)
+        if ((self.Data.noc==0) or (right.Data.noc==0)):
+            raise NotImplementedError, "The matrices must not be empty"
+        cdef Matrix_modpn_dense left = self.__copy__()
+        if right.is_immutable():
+            left.set_immutable()
+        if matsub(left.Data,right.Data):
+            return left
+        else:
+            raise ArithmeticError, "Matrix sizes or fields not compatible"
+
+    cpdef ModuleElement _neg_(self):
+        """
+        Negation of a MeatAxe matrix
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(25)
+            sage: MS = MatrixSpace(F, 2,4)
+            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
+            sage: -M    # indirect doctest
+            [  4   3   2   1]
+            [4*z 3*z 2*z   z]
+            sage: (-M) + M
+            [0 0 0 0]
+            [0 0 0 0]
+
+        """
+        if (self.Data.fl==2):
+            return self.__copy__()
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        cdef Matrix_modpn_dense result = self.__copy__()
+        if matneg(result.Data):
+            return result
+        else:
+            raise ArithmeticError, "Something went wrong"
+
+    # The default __mul__ relies on _multiply_classical, _multiply_strassen,
+    #_strassen_default_cutoff, _will_use_strassen, and _lmul_
+    def _multiply_classical(left, Matrix_modpn_dense other):
+        """
+        Multiply two MeatAxe matrices.
+
+        NOTE:
+
+        This method is implicitly called when multiplying matrices.
+        Large matrices will in instead use the asymptotically fast
+        Strassen-Winograd multiplication algorithm.
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: N = MatrixSpace(F,3)(list(reversed(list(F))))
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: N
+            [      1 2*z + 1 2*z + 2]
+            [      z       2   z + 2]
+            [  z + 1     2*z       0]
+            sage: M*N                   # indirect doctest
+            [2*z + 1 2*z + 2       2]
+            [    2*z       1   z + 1]
+            [      2   z + 2       z]
+
+        The result is mutable if and only if at least one of the factors
+        is mutable::
+
+            sage: M.is_mutable()
+            True
+            sage: (M*N).is_mutable()
+            True
+            sage: M.set_immutable()
+            sage: (M*N).is_mutable()
+            True
+            sage: N.set_immutable()
+            sage: (M*N).is_mutable()
+            False
+
+        """
+        cdef Matrix_modpn_dense Left = left
+        if (Left.ncols()==0) or (other.ncols()==0):
+            raise ValueError, "The matrices must not be empty"
+        if Left.Data.fl != other.Data.fl:
+            raise ValueError, "Different base fields"
+        from sage.all import MatrixSpace
+        parent = MatrixSpace(Left._base_ring, Left.Data.nor, other.Data.noc)
+        cdef Matrix_modpn_dense out = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        out.Data = matmalloc(Left.Data.fl, Left.Data.nor, other.Data.noc)
+        if matmul_result(Left.Data, other.Data, out.Data) != 0:
+            raise ValueError, "Matrix sizes are not compatible"
+        out._mutability = Mutability(Left._mutability._is_immutable and other._mutability._is_immutable)
+        out._parent = parent
+        out._base_ring = Left._base_ring
+        out._converter = Left._converter
+        out._ncols = other._ncols
+        out._nrows = Left._nrows
+        out._cache = {}
+        return out
+
+    cdef bint _will_use_strassen(self, Matrix right) except -2:
+        cutoff = self._strassen_default_cutoff(right)
+        if cutoff<0:
+            return False
+        if self._nrows < 2*cutoff:
+            return False
+        zsetfield(self.base_ring().order())
+        zsetlen(self._ncols)
+        if zrowsize_io<2*cutoff:
+            return False
+        return True
+
+    def _multiply_strassen(self, Matrix_modpn_dense right, long cutoff=0):
+        """
+        Multiplication using the Strassen-Winograd algorithm.
+
+        INPUT:
+
+        - Another matrix -- the second factor of the product
+        - ``cutoff`` -- optional integer, giving the size
+          below which the algorithm aborts the recursion
+          and reverts to classical (school book) multiplication.
+
+        NOTE:
+
+        By default, Strassen-Winograd multiplication is used
+        if the number of rows and the byte length of a row
+        exceeds 220.
+
+        EXAMPLE::
+
+            sage: MS = MatrixSpace(GF(9,'z'),1500,1500)
+            sage: A = MS.random_element()
+            sage: B = MS.random_element()
+            sage: type(A)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: C = A*B   # Here, Strassen-Winograd multiplication is the default
+            sage: C == A._multiply_classical(B)
+            True
+
+        In the next example, we force the use of MeatAxe matrices,
+        in order to demonstrate that Strassen-Winograd also works
+        in a case where Sage would use M4RIE matrices by default::
+
+            sage: MS = MatrixSpace(GF(16,'a'),200,200)
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: MS._MatrixSpace_generic__matrix_class = Matrix_modpn_dense
+            sage: A = MS.random_element()
+            sage: B = MS.random_element()
+            sage: C = A*B   # Here, school book multiplication is the default.
+            sage: type(C)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: C == A._multiply_strassen(B, cutoff=10)
+            True
+
+        """
+        if self._ncols != right._nrows:
+            raise ArithmeticError, "Number of columns of self must equal number of rows of right."
+        if not self._base_ring is right.base_ring():
+            raise TypeError, "Base rings must be the same."
+        if cutoff == 0:
+            cutoff = self._strassen_default_cutoff(right)
+        if cutoff <= 0:
+            raise ValueError, "cutoff must be at least 1"
+        set_cutoff(cutoff,cutoff)
+        from sage.all import MatrixSpace
+        MS = MatrixSpace(self.base_ring(), self._nrows, right._ncols)
+        result = Matrix_modpn_dense(MS,None)
+        if multiply_strassen(self.Data, right.Data, result.Data)!=0:
+            return RuntimeError, "Strassen multiplication did not work"
+        return result
+
+    cdef int _strassen_default_cutoff(self, Matrix right) except -2:
+        return 110
+
+    cpdef ModuleElement _lmul_(self, RingElement right):
+        """
+        Multiply an MeatAxe matrix with a field element.
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M*z                  # indirect doctest
+            [      0 2*z + 2 2*z + 1]
+            [      1     2*z   z + 1]
+            [  z + 2       2       z]
+            sage: 2*M
+            [      0       z 2*z + 2]
+            [2*z + 1       1     2*z]
+            [  z + 1   z + 2       2]
+            
+        """
+        cdef Matrix_modpn_dense left
+        if right==0:
+            return Matrix_modpn_dense(self._parent, entries=None, mutable=self._mutability._is_immutable) # this is a zero matrix
+        elif right==1:
+            return self.__copy__()
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        left = self.__copy__()
+        if matmulF(left.Data,zitof(self._converter.field_to_int(right))):
+            return left
+        else:
+            raise ArithmeticError, "Matrix sizes or fields not compatible"
+            
+    cpdef Matrix_modpn_dense _mulInt_(self, long right):
+        """
+        Multiply a MeatAxe matrix with a field element represented by an integer
+
+        NOTE:
+
+        This method is of internal use only
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M._mulInt_(4)
+            [      0   z + 2       2]
+            [      z 2*z + 2 2*z + 1]
+            [      1     2*z   z + 1]
+            sage: M*F.fetch_int(4)
+            [      0   z + 2       2]
+            [      z 2*z + 2 2*z + 1]
+            [      1     2*z   z + 1]
+
+        """
+        cdef Matrix_modpn_dense left
+        if right==0:
+            return Matrix_modpn_dense(self._parent, entries=None, mutable=self._mutability._is_immutable) # this is a zero matrix
+        elif right==1:
+            return self.__copy__()
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        left = self.__copy__()
+        if matmulF(left.Data,zitof(right)):
+            return left
+        else:
+            raise ArithmeticError, "Matrix sizes or fields not compatible"
+
+    def __div__(Matrix_modpn_dense self, long p):
+        """
+        Divide this matrix by a field element represented by an integer
+
+        NOTE:
+
+        If the field is not of prime order then the integer -1 does not
+        represent the additive inverse of 1.
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M/3    # indirect doctest
+            [      0       2       z]
+            [2*z + 2 2*z + 1       1]
+            [    2*z   z + 1   z + 2]
+            sage: M/-1
+            [      0   z + 2       2]
+            [      z 2*z + 2 2*z + 1]
+            [      1     2*z   z + 1]
+
+        We verify the results::
+
+            sage: M/3 == M*~F.fetch_int(3)
+            True
+            sage: M/-1 == M*~F.fetch_int(8)
+            True
+
+        """
+        if ((p%self.Data.fl)==0):
+            raise ZeroDivisionError, "Cannot divide by zero"
+        elif p==1:
+            return self.__copy__()
+        elif not self.Data:
+            return self
+        return self._mulInt_(tmultinv[zitof(p)])
+
+    def __pow__(Matrix_modpn_dense self,n,m):
+        """
+        M.__pow__(n): return M^n
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: M[2,2] = 0
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       0]
+            sage: M^2
+            [      0       0 2*z + 2]
+            [      0       2       0]
+            [  z + 1       0       0]
+            sage: M^4
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+            sage: M^-3
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       0]
+            sage: M^-3*M^3
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        if self._ncols != self._nrows:
+            raise TypeError, "Matrix is no square matrix"
+        cdef Matrix_modpn_dense OUT
+        cdef Matrix_modpn_dense SELFINV
+        OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        if n>=0:
+            _sig_on
+            OUT.Data = matpower(self.Data,n)
+            _sig_off
+        else:
+            SELFINV = self.__invert__()
+            _sig_on
+            OUT.Data = matpower(SELFINV.Data,-n)
+            _sig_off
+        if OUT.Data:
+            OUT._mutability = copy(self._mutability)
+            OUT._parent = self._parent
+            OUT._converter = self._converter
+            OUT._base_ring = self._base_ring
+            OUT._ncols = self._ncols
+            OUT._nrows = self._nrows
+            OUT._cache = {}
+            return OUT
+        else:
+            raise RuntimeError, "Something went wrong"
+
+    def __invert__(Matrix_modpn_dense self):
+        """
+        M__invert__(): return M^(-1), if defined
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: ~M
+            Traceback (most recent call last):
+            ...
+            ArithmeticError: This MeatAxe matrix is not invertible
+            sage: M[2,2] = 0
+            sage: ~M    # indirect doctest
+            [      2       z       0]
+            [2*z + 1       1     2*z]
+            [      0   z + 2       2]
+            sage: M*~M
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+        
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        _sig_on
+        OUT.Data = matinv(self.Data)
+        _sig_off
+        if OUT.Data:
+            OUT._mutability = copy(self._mutability)
+            OUT._parent = self._parent
+            OUT._converter = self._converter
+            OUT._base_ring = self._base_ring
+            OUT._ncols = self._ncols
+            OUT._nrows = self._nrows
+            OUT._cache = {}
+            return OUT
+        else:
+            raise ArithmeticError, "This MeatAxe matrix is not invertible"
+
+    cpdef Matrix_modpn_dense transpose(Matrix_modpn_dense self):
+        """
+        Return the transposed matrix
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(121,'z'),3)([2,3,4,5,6,7,8,9,0])
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: M
+            [2 3 4]
+            [5 6 7]
+            [8 9 0]
+            sage: M.transpose()
+            [2 5 8]
+            [3 6 9]
+            [4 7 0]
+            sage: transpose(M)   # indirect doctest
+            [2 5 8]
+            [3 6 9]
+            [4 7 0]
+
+        TESTS::
+
+            sage: M = MatrixSpace(GF(121,'z'),2,3)(range(6))
+            sage: M
+            [0 1 2]
+            [3 4 5]
+            sage: M.parent()
+            Full MatrixSpace of 2 by 3 dense matrices over Finite Field in z of size 11^2
+            sage: M.transpose()
+            [0 3]
+            [1 4]
+            [2 5]
+            sage: M.transpose().parent()
+            Full MatrixSpace of 3 by 2 dense matrices over Finite Field in z of size 11^2
+
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        OUT.Data = mattr(self.Data)
+        if OUT.Data:
+            OUT._mutability = copy(self._mutability)
+            from sage.all import MatrixSpace
+            OUT._parent = MatrixSpace(self._base_ring, self._ncols, self._nrows)
+            OUT._base_ring = self._base_ring
+            OUT._converter = self._converter
+            OUT._ncols = self._nrows
+            OUT._nrows = self._ncols
+            OUT._cache = {}
+            return OUT
+        else:
+            raise ArithmeticError, "Error allocating memory"
+
+    cpdef int order(Matrix_modpn_dense self) except -1:
+        """
+        Return multiplicative order of self 
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)([F.fetch_int(i) for i in [2,3,4,5,6,7,8,9,0]])
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: order(M)   # indirect doctest
+            728
+            sage: M.order()
+            728
+            sage: M^728
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+
+        We verify the result::
+
+            sage: E = MatrixSpace(F,3)(1)
+            sage: any([M^i==E for i in range(1,728)])
+            False
+            sage: M^728
+            [1 0 0]
+            [0 1 0]
+            [0 0 1]
+
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        if (self.Data.nor <> self.Data.noc):
+            raise NotImplementedError, "only defined for square matrices"
+        cdef int o = matorder(self.Data)
+        if o==-1:
+            raise ArithmeticError, "order too large"
+        else:
+            return o
+
+    def stack(self, *L):
+        """
+        Stack one matrix over one or several others, provided the number of columns coincides
+
+        EXAMPLES::
+
+            sage: MS = MatrixSpace(GF(25,'z'),2,4)
+            sage: M = MS([[4,3,2,1],[3,2,1,0]])
+            sage: N = MS([[1,1,2,2],[3,3,4,4]])
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: print M.stack(N)
+            [4 3 2 1]
+            [3 2 1 0]
+            [1 1 2 2]
+            [3 3 4 4]
+            sage: print M.stack(N,M)
+            [4 3 2 1]
+            [3 2 1 0]
+            [1 1 2 2]
+            [3 3 4 4]
+            [4 3 2 1]
+            [3 2 1 0]
+
+        """
+        cdef Matrix_modpn_dense N
+        cdef Matrix_modpn_dense OUT
+        cdef PTR cur
+        cdef int totalrows = self.Data.nor
+        for N in L:
+            if self.Data.fl!=N.Data.fl:
+                raise TypeError, "Matrices must be defined over the same field"
+            if self.Data.noc!=N.Data.noc:
+                raise RuntimeError, "Matrices must have the same number of columns"
+            totalrows += N.Data.nor
+        from sage.all import MatrixSpace
+        OUT = Matrix_modpn_dense(MatrixSpace(self._base_ring, totalrows, self.Data.noc), entries=None)
+        zsetfield(self.Data.fl)
+        zsetlen(self.Data.noc)
+        cur = OUT.Data.d
+        if not copyrows(cur,self.Data.d,self.Data.nor):
+            raise MemoryError, "Error copying first matrix"
+        zadvance(&(cur),self.Data.nor)
+        for i in range(len(L)):
+            N = L[i]
+            copyrows(cur,N.Data.d,N.Data.nor)
+            zadvance(&(cur),N.Data.nor)
+        return OUT
+
+
+###################
+## Gauss algorithm
+
+    cpdef Matrix_modpn_dense normalized(Matrix_modpn_dense self):
+        """
+        Return a normalized copy of self,  i.e., all leading coefficients are 1
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)([F.fetch_int(i) for i in [2,3,4,5,6,7,8,9,0]])
+            sage: M.normalized()
+            [      1     2*z 2*z + 2]
+            [      1 2*z + 2       2]
+            [      1       0       0]
+            sage: M
+            [      2       z   z + 1]
+            [  z + 2     2*z 2*z + 1]
+            [2*z + 2       0       0]
+
+        """
+        if self.Data.fl==2:
+            return self.__copy__()
+        cdef Matrix_modpn_dense OUT = self.__copy__()
+        cdef int i
+        cdef PTR p = OUT.Data.d
+        cdef piv_t P
+        cdef FEL f
+        for i from 0<=i<self.Data.nor:
+            P = _zfindpiv_(p)
+            if P.m:
+                if P.m!=1:
+                    f = tmultinv[zitof(P.m)]
+                    zmulrow(p,f)
+            zsteprow(&(p))
+        return OUT
+
+    cpdef Matrix_modpn_dense semi_echelon(Matrix_modpn_dense self):
+        """
+        Return a normalized semi echelon form of self with zero rows removed
+
+        EXAMPLES::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: M
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+            sage: M.semi_echelon()
+            [  0   1 2*z]
+            [  1   0   0]
+
+        TESTS::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: F.<z> = GF(29)
+            sage: M = Matrix_modpn_dense(MatrixSpace(F,5), range(1,26))
+            sage: M.semi_echelon()
+            [1 2 3 4 5]
+            [0 1 2 3 4]
+            sage: M.echelon_form()
+            [ 1  0 28 27 26]
+            [ 0  1  2  3  4]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        cdef Matrix_modpn_dense OUT
+        cdef int i
+        OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
+        _sig_on
+        OUT.Data = echelon(self.Data)
+        _sig_off
+        cdef PTR p = OUT.Data.d
+        cdef piv_t P
+        cdef FEL f
+        if self.Data.fl != 2:
+            for i from 0<=i<OUT.Data.nor:
+                P = _zfindpiv_(p)
+                if P.m:
+                    if P.m!=1:
+                        f = tmultinv[zitof(P.m)]
+                        zmulrow(p,f) # OUT[i] = OUT[i]/f
+                zsteprow(&(p))
+        OUT._mutability = copy(self._mutability)
+        from sage.all import MatrixSpace
+        OUT._parent = MatrixSpace(self._base_ring, OUT.Data.nor, OUT.Data.noc)
+        OUT._base_ring = self._base_ring
+        OUT._converter = self._converter
+        OUT._ncols = OUT.Data.noc
+        OUT._nrows = OUT.Data.nor
+        OUT._cache = {}
+        return OUT
+
+    def _echelon_in_place_classical(self, reduce=True):
+        """
+        Change self into row echelon form.
+
+        INPUT:
+
+        - ``reduce`` (optional, default True): Compute the *reduced* row echelon form.
+
+        EXAMPLES:
+
+        Note that in our example we use MeatAxe matrices in a case where
+        one would usually have Linbox as a backend::
+
+            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+            sage: F.<z> = GF(29)
+            sage: M = Matrix_modpn_dense(MatrixSpace(F,5), range(1,26))
+            sage: M
+            [ 1  2  3  4  5]
+            [ 6  7  8  9 10]
+            [11 12 13 14 15]
+            [16 17 18 19 20]
+            [21 22 23 24 25]
+            sage: M.echelon_form()
+            [ 1  0 28 27 26]
+            [ 0  1  2  3  4]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+            sage: M._echelon_in_place_classical(reduce=False)
+            sage: M
+            [ 1  2  3  4  5]
+            [ 0 24 19 14  9]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+            [ 0  0  0  0  0]
+
+        TESTS:
+
+        The following gave rise to a segfault::
+
+            sage: M = matrix(GF(125,'z'), 0,2000)
+            sage: M.echelon_form()
+            0 x 2000 dense matrix over Finite Field in z of size 5^3
+            sage: M = matrix(GF(125,'z'), 2000, 0)
+            sage: M.echelon_form()
+            2000 x 0 dense matrix over Finite Field in z of size 5^3
+
+        """
+        if (self.Data.nor == 0) or (self.Data.noc == 0):
+            self.cache('pivots', ())
+            return
+        zsetfield(self.Data.fl)
+        zsetlen(self.Data.noc)
+        cdef int nr = self.Data.nor
+        cdef long* piv = <long*>calloc(nr+1,sizeof(long))
+        zmkechelon(self.Data.d, nr, piv)
+        cdef int i
+        cdef int rk = piv[0]
+        self.cache('pivots', tuple(piv[i]-1 for i in xrange(1,rk+1)))
+        if rk == 0:
+            free(piv)
+            return
+        # Need to empty the remaining rows - zmkechelon was moving them
+        # up, but did not clear them.
+        cdef PTR p = self.Data.d
+        zadvance(&(p), rk)
+        for i from rk<=i<nr:
+            memset(p, 0, zrowsize)
+            zsteprow(&(p))
+        # Sage wants that the reduced row echelon form is created.
+        # Hence, we need to normalize and to clear the parts on top
+        # of the pivots.
+        cdef PTR q
+        cdef FEL f
+        cdef int j
+        if reduce:
+            p = self.Data.d
+            f = tmultinv[zitof(zextract(p,piv[1]))]
+            zmulrow(p,f)
+            q = self.Data.d
+            for i from 2<=i<=rk:
+                zsteprow(&(q))
+                # normalize the row
+                f = tmultinv[zitof(zextract(q,piv[i]))]
+                zmulrow(q,f)
+                # clear all rows on top of it
+                p = self.Data.d
+                for j from 1<=j<i:
+                    f = taddinv[zitof(zextract(p,piv[i]))]
+                    zaddmulrow(p, q, f)
+                    zsteprow(&(p))
+        free(piv)
+
+    cpdef int nullity(Matrix_modpn_dense self):
+        """
+        Return the nullity of a matrix
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(25,'z'),2,4)([[4,3,2,1],[3,2,1,0]])
+            sage: M.nullity()
+            2
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        return nullity(self.Data)
+
+    cpdef tuple lead(self):
+        """
+        ``(f,i) = M.lead()`` `\\iff` ``f=M[0,i]`` is the first non-zero coefficient in the first row of M.
+        
+        If the first row of M has no non-zero entry then ``f==0`` and ``i`` is the number of columns
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(GF(25,'z'),1,4)([0,2,4,1])
+            sage: M.lead()
+            (2, 1)
+            sage: N = MatrixSpace(GF(25,'z'),3,15)(0)
+            sage: N.lead()
+            (0, 15)
+            
+        """
+        if (self.ncols()==0):
+            raise NotImplementedError, "The matrix must not be empty"
+        zsetfield(self.Data.fl)
+        zsetlen(self.Data.noc)
+        cdef piv_t P = _zfindpiv_(self.Data.d)
+        return (P.m, P.i)
+
+    # set/get "fast and dirty"
+
+    cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
+        """
+        Get an element without checking.
+
+        TEST::
+
+            sage: F.<z> = GF(9)
+            sage: M = MatrixSpace(F,3)(list(F))
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: M    # indirect doctest
+            [      0     2*z   z + 1]
+            [  z + 2       2       z]
+            [2*z + 2 2*z + 1       1]
+
+        """
+        cdef PTR p = self.Data.d
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+        else:
+            raise IndexError, "Matrix is empty"
+        zadvance(&(p),i)
+        return self._converter.int_to_field(zftoi(zextract(p,j+1)))
+
+    cdef int get_unsafe_int(self, Py_ssize_t i, Py_ssize_t j):
+        # NOTE:
+        # It is essential that you call zsetfield and zsetlen YOURSELF
+        # and that you assert that the matrix is not empty!
+        # This method is here for speed!
+        cdef PTR p = self.Data.d
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+        else:
+            raise IndexError, "Matrix is empty"
+        zadvance(&(p),i)
+        return zftoi(zextract(p,j+1))
+
+    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
+        # ASSUMPTION: value's parent is the base ring
+        if self.Data:
+            zsetfield(self.Data.fl)
+            zsetlen(self.Data.noc)
+        else:
+            raise IndexError, "Matrix is empty"
+        cdef PTR x = self.Data.d
+        zadvance(&(x),i)
+        zinsert(x, j+1, zitof(self._converter.field_to_int(value)))
+
+    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value):
+        # NOTE:
+        # It is essential that you call zsetfield and zsetlen YOURSELF
+        # and that you assert that the matrix is not empty!
+        # This method is here for speed!
+        cdef PTR x = self.Data.d
+        zadvance(&(x),i)
+        zinsert(x, j+1, zitof(value))
+
+    def randomize(self, density=None, nonzero=False, *args, **kwds):
+        """
+        Fill the matrix with random values.
+
+        INPUT:
+
+        - ``density`` (optional real number between zero and one) --
+          the expected density of the resulting matrix
+        - ``nonzero`` (optional bool, default ``False``) -- 
+          If true, all inserted marks are non-zero.
+
+        EXAMPLE::
+
+            sage: MS = MatrixSpace(GF(27,'z'),6,6)
+            sage: M = MS.random_element(); M    # indirect doctest
+            [              1           z + 1     z^2 + z + 1             z^2       2*z^2 + z           z + 1]
+            [2*z^2 + 2*z + 2   2*z^2 + z + 2         z^2 + 1 2*z^2 + 2*z + 2         z^2 + z   2*z^2 + z + 1]
+            [        2*z + 2     z^2 + z + 2           z + 2 2*z^2 + 2*z + 2           2*z^2           2*z^2]
+            [  2*z^2 + z + 2             z^2           z + 2         z^2 + z       2*z^2 + 2         z^2 + 2]
+            [      2*z^2 + z             2*z 2*z^2 + 2*z + 1       2*z^2 + 1 2*z^2 + 2*z + 1       2*z^2 + z]
+            [        2*z + 1         z^2 + z             z^2             z^2     2*z^2 + 2*z           z + 1]
+            sage: type(M)
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+            sage: MS.random_element(nonzero=True)
+            [            2*z               1   z^2 + 2*z + 1   2*z^2 + z + 1             z^2     z^2 + z + 1]
+            [    2*z^2 + 2*z   2*z^2 + z + 2         2*z + 1       z^2 + 2*z     2*z^2 + 2*z             z^2]
+            [        z^2 + z     z^2 + z + 2 2*z^2 + 2*z + 1         z^2 + 2               1           2*z^2]
+            [              z     2*z^2 + 2*z           2*z^2         2*z + 1           z + 2           z + 2]
+            [        z^2 + z             z^2           z + 2     2*z^2 + 2*z         2*z + 1         z^2 + z]
+            [    z^2 + z + 2       2*z^2 + z             z^2           z + 1     2*z^2 + 2*z   z^2 + 2*z + 1]
+            sage: MS.random_element(density=0.5)
+            [        z^2 + 2               0   z^2 + 2*z + 2       2*z^2 + z               0     z^2 + z + 2]
+            [              0               1               0               0               0               0]
+            [  2*z^2 + z + 1   2*z^2 + z + 2               0     z^2 + z + 2               0     z^2 + z + 1]
+            [              0               0               0               0               0               0]
+            [2*z^2 + 2*z + 2               0               0   2*z^2 + z + 2               0         2*z + 1]
+            [              0       2*z^2 + z               0               1               0   2*z^2 + z + 1]
+
+        """
+        self.check_mutability()
+        cdef int fl = self.Data.fl
+        density = float(density)
+        if density <= 0:
+            return
+        if density > 1:
+            density = float(1)
+
+        self.clear_cache()
+
+        cdef unsigned char *x,*y
+        x = <unsigned char*>(self.Data.d)
+        cdef int nr = self.Data.nor
+        cdef int nc = self.Data.noc
+        cdef int i, j, k
+
+        zsetfield(fl)
+        zsetlen(nc)
+        cdef int O
+        randint = current_randstate().c_random
+        randdouble = current_randstate().c_rand_double
+
+        if not nonzero:
+            if density == 1:
+                O = (self.base_ring().order()**MPB)
+                sig_on()
+                for i from 0 <= i < nr:
+                    y = x
+                    for j from 0 <= j < zrowsize_io:
+                        y[j] = randint()%O
+                    zsteprow(&(x))
+                sig_off()
+            else:
+                sig_on()
+                for i from 0 <= i < nr:
+                    y = x
+                    for j from 0 < j <= nc:
+                        if randdouble() < density:
+                            zinsert(y, j, zitof( (randint()%fl) ))
+                    zsteprow(&(x))
+                sig_off()
+        else:
+            if density == 1:
+                fl -= 1
+                sig_on()
+                for i from 0 <= i < nr:
+                    y = x
+                    for j from 0 < j <= nc:
+                        zinsert(y, j, zitof( (randint()%fl)+1 ))
+                    zsteprow(&(x))
+                sig_off()
+            else:
+                fl -= 1
+                sig_on()
+                for i from 0 <= i < nr:
+                    y = x
+                    for j from 0 < j <= nc:
+                        if randdouble() < density:
+                            zinsert(y, j, zitof( (randint()%fl)+1 ))
+                    zsteprow(&(x))
+                sig_off()

sage/matrix/matrix_space.py

@@ -952 +952 @@
             <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
             sage: type(matrix(GF(16007), 2, range(4)))
             <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
+            sage: type(matrix(GF(2), 2, range(4)))
+            <type 'sage.matrix.matrix_mod2_dense.Matrix_mod2_dense'>
+            sage: type(matrix(GF(64,'z'), 2, range(4)))
+            <type 'sage.matrix.matrix_mod2e_dense.Matrix_mod2e_dense'>
+            sage: type(matrix(GF(125,'z'), 2, range(4)))
+            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
+
         """
         R = self.base_ring()
         if self.is_dense():
@@ -981 +988 @@
                 return matrix_generic_dense.Matrix_generic_dense
             elif sage.rings.finite_rings.all.is_FiniteField(R) and R.characteristic() == 2 and R.order() <= 1024:
                 return matrix_mod2e_dense.Matrix_mod2e_dense
+            elif sage.rings.finite_rings.all.is_FiniteField(R) and R.characteristic() != 2 and R.order() <= 255:
+                from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
+                return Matrix_modpn_dense
             elif sage.rings.polynomial.multi_polynomial_ring_generic.is_MPolynomialRing(R) and R.base_ring() in _Fields:
                 return matrix_mpolynomial_dense.Matrix_mpolynomial_dense
             #elif isinstance(R, sage.rings.padics.padic_ring_capped_relative.pAdicRingCappedRelative):

sage/modules/quotient_module.py

@@ -25 +25 @@
 
     To obtain V or W use \code{self.V()} and \code{self.W()}.
 
-    EXAMPLES:
+    EXAMPLES::
+
         sage: k.<i> = QuadraticField(-1)
         sage: A = k^3; V = A.span([[1,0,i], [2,i,0]])
         sage: W = A.span([[3,i,i]])
@@ -72 +73 @@
         """
         Create this quotient space, from the given domain, sub-module, and quotient_matrix.
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: A = QQ^5; V = A.span_of_basis([[1,0,-1,1,1], [1,-1,0,2/3,3/4]]); V
             Vector space of degree 5 and dimension 2 over Rational Field
             User basis matrix:
@@ -83 +85 @@
             User basis matrix:
             [1/3 2/3  -1 5/9 1/2]
 
-        This creates a quotient vector space, which calls the init method:
+        This creates a quotient vector space, which calls the init method::
+
             sage: Q = V / W
 
-        Behold the type of Q:
+        Behold the type of Q::
+
             sage: type(Q)
             <class 'sage.modules.quotient_module.FreeModule_ambient_field_quotient_with_category'>
             
         We do some consistency checks on the extra quotient and
         lifting structure of Q.
+        ::
+
             sage: Q(V.0)
             (1)
             sage: Q( V.0 - 2/3*V.1 )
@@ -116 +122 @@
         Return the rather verbose string representation of this quotient space V/W.
 
         EXAMPLES:
-        We create a quotient vector space over a finite field:
+
+        We create a quotient vector space over a finite field::
+
             sage: k.<a> = GF(9); A = k^3; V = A.span_of_basis([[1,0,a], [a,a,1]]); W = V.span([V.1])
             sage: Q = V/W
 
-        Note the type:
+        Note the type::
+
             sage: type(Q)
             <class 'sage.modules.quotient_module.FreeModule_ambient_field_quotient_with_category'>
 
         The string representation mentions that this is a quotient V/W, that the quotient
-        has dimension 1 and is over a finite field, and also describes V and W:
+        has dimension 1 and is over a finite field, and also describes V and W::
+
             sage: Q._repr_()
             'Vector space quotient V/W of dimension 1 over Finite Field in a of size 3^2 where\nV: Vector space of degree 3 and dimension 2 over Finite Field in a of size 3^2\nUser basis matrix:\n[1 0 a]\n[a a 1]\nW: Vector space of degree 3 and dimension 1 over Finite Field in a of size 3^2\nBasis matrix:\n[    1     1 a + 2]'
         """
@@ -140 +150 @@
         the tuple (V,W).
 
         EXAMPLES:
-        We compute the hash of a certain 0-dimension quotient vector space:
+
+        We compute the hash of a certain 0-dimension quotient vector space::
+
             sage: A = QQ^2; V = A.span_of_basis([[1,0], [1,1]]); W = V.span([V.1, V.0])
             sage: Q = V/W; Q.dimension()
             0
@@ -149 +161 @@
             -5856620741060301410    # 64-bit
 
         The hash is just got by hashing both V and W.
+        ::
+
             sage: hash((V, W))
             954887582             # 32-bit
             -5856620741060301410  # 64-bit
@@ -165 +179 @@
         self and other.
 
         EXAMPLES:
-        We create three quotient spaces and compare them:
+
+        We create three quotient spaces and compare them::
+
             sage: A = QQ^2; V = A.span_of_basis([[1,0], [1,1]]);
             sage: W0 = V.span([V.1, V.0]); W1 = V.span([V.1]); W2 = V.span([V.1])
             sage: Q0 = V/W0; Q1 = V/W1; Q2 = V/W2
@@ -191 +207 @@
         if it can be made sense of as a list of length the dimension of self.
                 
         EXAMPLES:
+
         We create a 2-dimensional quotient of a 3-dimension ambient vector space.
+        ::
+
             sage: M = QQ^3 / [[1,2,3]]
 
         A list of length 3 coerces into QQ^3, so it coerces into M.
+        ::
+
             sage: M([1,2,4])
             (-1/3, -2/3)
             sage: M([1,2,3])
@@ -202 +223 @@
 
         A list of length 2 at least coerces into M, where here it just gives
         the corresponding linear combination of the basis for M.
+        ::
+            
             sage: M([1,2])
             (1, 2)
             sage: M.0 + 2*M.1
@@ -225 +248 @@
         Elements canonically coerce into self if they canonically
         coerce into V.
         
-        EXAMPLES:
+        EXAMPLES::
+
             sage: V = QQ^3; W = V.span([[1,0,0]]); Q = V/W
             sage: Q._coerce_impl(V.0)
             (0, 0)
@@ -239 +263 @@
             (2, 0)
 
         Here we coerce in something that is over ZZ, so it canonically coerce into V hence self.
+        ::
+
             sage: Q._coerce_impl((ZZ^3)([1,2,3]))
             (2, 3)
 
-        Here there is no canonical coercion:
+        Here there is no canonical coercion::
+
             sage: Q._coerce_impl((GF(17)^3)([1,2,3]))
             Traceback (most recent call last):
             ...
@@ -259 +286 @@
         """
         Given this quotient space $Q = V/W$, return the natural quotient map from V to Q.
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^3 / [[1,2,3]]
             sage: M.quotient_map()
             Free module morphism defined by the matrix
@@ -279 +307 @@
         Given this quotient space $Q = V/W$, return a fixed choice of linear homomorphism
         (a section) from Q to V.
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^3 / [[1,2,3]]
             sage: M.lift_map()
             Free module morphism defined by the matrix
@@ -296 +325 @@
 
         The lift is a fixed homomorphism.  
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^3 / [[1,2,3]]
             sage: M.lift(M.0)
             (1, 0, 0)
@@ -311 +341 @@
         """
         Given this quotient space $Q = V/W$, return W.
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^10 / [range(10), range(2,12)]
             sage: M.W()
             Vector space of degree 10 and dimension 2 over Rational Field
@@ -325 +356 @@
         """
         Given this quotient space $Q = V/W$, return $V$.
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^10 / [range(10), range(2,12)]
             sage: M.V()
             Vector space of dimension 10 over Rational Field
@@ -336 +368 @@
         """
         Given this quotient space $Q = V/W$, return $V$.  This is the same as self.V().
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^10 / [range(10), range(2,12)]
             sage: M.cover()
             Vector space of dimension 10 over Rational Field
@@ -347 +380 @@
         """
         Given this quotient space $Q = V/W$, return $W$.  This is the same as self.W().
 
-        EXAMPLES:
+        EXAMPLES::
+
             sage: M = QQ^10 / [range(10), range(2,12)]
             sage: M.relations()
             Vector space of degree 10 and dimension 2 over Rational Field