Ticket #12103: trac12103_meataxe_rel11900.patch

File trac12103_meataxe_rel11900.patch, 103.6 KB (added by SimonKing, 9 years ago)

Use MeatAxe as back end for linear algebra over some finite fields. Relative to #11900

  • module_list.py

    # HG changeset patch
    # User Simon King <simon.king@uni-jena.de>
    # Date 1322738482 -3600
    # Node ID c4a1b1e7b3ad7b23fd7a4563ecf120cc63875cc5
    # Parent  2d7c2948fc20d042dfd43a4f92a4a44561f88496
    #12103: MeatAxe wrapper, fully compatible with other matrices in Sage
    The spkg includes Strassen-Winograd multiplication
    So far missing: Strassen algorithm for computation of echelon forms.
    
    diff --git a/module_list.py b/module_list.py
    a b  
    946946              extra_compile_args = m4ri_extra_compile_args,
    947947              language="c++"),
    948948
     949    Extension('sage.matrix.matrix_modpn_dense',
     950              sources = ['sage/matrix/matrix_modpn_dense.pyx'],
     951              libraries = ['meataxe', 'csage'],
     952              depends = [SAGE_INC + "meataxe.h"]),
     953
    949954    Extension('sage.matrix.matrix_modn_dense',
    950955              sources = ['sage/matrix/matrix_modn_dense.pyx'],
    951956              libraries = ['gmp']),
  • sage/matrix/constructor.py

    diff --git a/sage/matrix/constructor.py b/sage/matrix/constructor.py
    a b  
    10621062
    10631063        sage: K.<a>=FiniteField(3^2)
    10641064        sage: random_matrix(K, 2, 5)
    1065         [      1     2*a       1   a + 2       2]
    1066         [      a 2*a + 1       0       2       1]
    1067        
     1065        [      a       0       2       2   a + 1]
     1066        [      a   a + 2   a + 1 2*a + 1     2*a]
     1067
    10681068        sage: random_matrix(RR, 3, 4, density=0.66)
    1069         [ 0.000000000000000  0.566500636438206 0.0870635178173962  0.000000000000000]
    1070         [-0.662290145671671  0.000000000000000  0.475667133865666  0.000000000000000]
    1071         [-0.276405104068647  0.000000000000000  0.000000000000000 -0.636689607643642]
    1072        
     1069        [ 0.000000000000000 -0.615392177288427 -0.636689607643642  0.000000000000000]
     1070        [ 0.575979457284225  0.000000000000000 0.0870635178173962  0.000000000000000]
     1071        [-0.662290145671671  0.000000000000000  0.000000000000000 0.0362592516947566]
     1072
    10731073        sage: A = random_matrix(ComplexField(32), 3, density=0.8, sparse=True); A
    1074         [                             0   -0.193242896 + 0.460619713*I                              0]
    1075         [   0.909948633 + 0.611092515*I 0.00128001347 + 0.0659103142*I   0.199796580 + 0.0955426861*I]
    1076         [                             0                              0    0.354729790 - 0.184624095*I]
     1074        [                             0    0.909948633 + 0.611092515*I                              0]
     1075        [  -0.728572036 - 0.455087671*I    0.354729790 - 0.184624095*I 0.00128001347 + 0.0659103142*I]
     1076        [                             0                              0   -0.360952012 - 0.698699617*I]
    10771077        sage: A.is_sparse()
    10781078        True
    10791079
  • sage/matrix/matrix0.pyx

    diff --git a/sage/matrix/matrix0.pyx b/sage/matrix/matrix0.pyx
    a b  
    33353335        x = self.fetch('rank')
    33363336        if not x is None: return x
    33373337        if self._nrows == 0 or self._ncols == 0:
    3338             return 0       
     3338            return 0
    33393339        r = len(self.pivots())
    33403340        self.cache('rank', r)
    33413341        return r
  • sage/matrix/matrix2.pyx

    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    26282628        We test that the generic code is called for matrices over fields,
    26292629        lacking any more specific routine. ::
    26302630
    2631             sage: F.<a> = FiniteField(5^2)
     2631            sage: F.<a> = FiniteField(5^4)
    26322632            sage: A = matrix(F, 3, 4, [[  1,   a,     1+a,  a^3+a^5],
    26332633            ...                        [  a, a^4,   a+a^4,  a^4+a^8],
    26342634            ...                        [a^2, a^6, a^2+a^6, a^5+a^10]])
    26352635            sage: set_verbose(1)
    26362636            sage: A.right_kernel(algorithm='default')
    2637             verbose ...
     2637            verbose 1 ...
    26382638            verbose 1 (<module>) computing right kernel matrix over an arbitrary field for 3x4 matrix
    2639             ...
    2640             verbose 1 (<module>) done computing right kernel matrix over an arbitrary field for 3x4 matrix
    2641             ...
    2642             Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^2
     2639            verbose 1 (<module>) generic in-place Gauss elimination on 3 x 4 matrix
     2640            verbose 1 (<module>) done with gauss echelon form (time = ...)
     2641            verbose 1 (<module>) done computing right kernel matrix over an arbitrary field for 3x4 matrix (time = ...)
     2642            verbose 1 (<module>) generic in-place Gauss elimination on 2 x 4 matrix
     2643            verbose 1 (<module>) done with gauss echelon form (time = ...)
     2644            verbose 1 ...
     2645            Vector space of degree 4 and dimension 2 over Finite Field in a of size 5^4
    26432646            Basis matrix:
    2644             [      1       0 3*a + 4 2*a + 2]
    2645             [      0       1     2*a 3*a + 3]
     2647            [                1                 0     a^3 + a^2 + 3   2*a^3 + a^2 + 1]
     2648            [                0                 1 4*a^3 + 4*a^2 + 1 3*a^3 + 4*a^2 + 4]
    26462649            sage: set_verbose(0)
    26472650
    26482651        Over the Integers:
  • new file sage/matrix/matrix_modpn_dense.pxd

    diff --git a/sage/matrix/matrix_modpn_dense.pxd b/sage/matrix/matrix_modpn_dense.pxd
    new file mode 100644
    - +  
     1#*****************************************************************************
     2#
     3#    Dense matrices over GF(p^n), p^n<255, using LibMeataxe as backend
     4#
     5#    Copyright (C) 2011 Simon A. King <simon.king@uni-jena.de>
     6#
     7#  Distributed under the terms of the GNU General Public License (GPL),
     8#  version 2 or later (at your choice)
     9#
     10#    This code is distributed in the hope that it will be useful,
     11#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     12#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     13#    General Public License for more details.
     14#
     15#  The full text of the GPL is available at:
     16#
     17#                  http://www.gnu.org/licenses/
     18#*****************************************************************************
     19
     20include "../ext/interrupt.pxi"
     21
     22#
     23# First, something from Python.h and string.h
     24#
     25
     26cdef extern from "Python.h":
     27    object PyString_FromStringAndSize(char *s, Py_ssize_t len)
     28    char* PyString_AsString(object string)
     29
     30#
     31# Next, import SOME types from meataxe.h
     32# (most types are not needed, but listed here
     33# in the comments, for completeness)
     34#
     35
     36cdef extern from "stdlib.h":
     37    void* calloc ( size_t num, size_t size )
     38    void free ( void * ptr )
     39    void * memset ( void * ptr, int value, size_t num )
     40
     41cdef extern from "meataxe.h":
     42    # general ctype emulations
     43    ctypedef int size_t
     44    ctypedef unsigned long Ulong
     45    ctypedef unsigned short Ushort
     46    ctypedef unsigned char Uchar
     47    ctypedef struct FILE
     48   
     49    # assume small fields (default)
     50    ctypedef unsigned char FEL
     51    ctypedef FEL *PTR
     52    cdef extern FEL tmult[256][256]
     53    cdef extern FEL tadd[256][256]
     54    cdef extern FEL taddinv[256], tmultinv[256]
     55    cdef extern FEL tinsert[8][256], textract[8][256]
     56
     57    # pivots
     58    ctypedef struct piv_t:
     59        long i,m
     60
     61    #
     62    # cdef extern long zfl, zchar       #/* Field order and characteristic */
     63    # cdef extern FEL zgen              #/* Generator */
     64    # cdef extern long znoc             #/* Number of columns for row ops */
     65    cdef extern size_t zrowsize         #/* Row size in memory */
     66    cdef extern size_t zrowsize_io      #/* Row size in files */
     67    cdef extern size_t MPB              # number of FEL stored in one byte
     68    # cdef extern char zzzversion[]
     69    # cdef extern char zzz_cc[]
     70       
     71#########################################
     72# function prototypes
     73    # FEL zadd(FEL a,FEL b)
     74    # FEL zneg(FEL a)
     75    # FEL zmul(FEL a,FEL b)
     76    # FEL zinv(FEL a)
     77    #
     78    # FEL zembed(FEL a, long subfield)  #/* Embed from subfield */
     79    # FEL zrestrict(FEL a, long subfield)       #/* Restrict to subfield */
     80    void *copyrows(PTR dest, PTR src, size_t nor)
     81
     82    FEL zitof(long l)
     83    long zftoi(FEL f)
     84
     85    int zsetlen(long ncols)             #/* Set row size */
     86    int zsetfield(long field)
     87    void zinsert(PTR row, long col, FEL mark)
     88    FEL zextract(PTR row, long col)
     89    void zinsert_step(PTR *pos, int *idx, FEL mark)
     90    FEL zextract_step(PTR *pos, int *idx)
     91
     92    piv_t _zfindpiv_(PTR row)
     93    # PTR zextractcol(PTR mat, long nor, long col, PTR result)
     94    # PTR zaddrow(PTR dest, PTR src)
     95    # PTR zsubrow(PTR dest, PTR src)
     96    # PTR znegrow(PTR row)
     97    PTR zmulrow(PTR row, FEL mark)
     98    PTR zaddmulrow(PTR dest, PTR src, FEL f)
     99    # PTR zmaprow(PTR row, PTR matrix, long nor, PTR result)
     100    # PTR zpermrow(PTR row, long *perm, PTR result)
     101
     102    void zmoverow(PTR dest, PTR src)
     103    void zcpyrow(PTR dest, PTR src, long sz)
     104    # void zswaprow(PTR dest, PTR src)
     105    # int zcmprow(PTR dest, PTR src)
     106    size_t zsize(long nrows)
     107    # PTR zalloc(long NROWS)
     108    void zadvance(PTR *ptr, long nrows)
     109    void zsteprow(PTR *ptr)
     110    void zbackrow(PTR *ptr)
     111    PTR ptrPlus(PTR base, long offset)
     112    # int mtxinit()
     113
     114    # skip gap-conversion
     115    # bitstring operations
     116    # ctypedef struct bitstring_t:
     117    #     Ushort id     #/* always = T_BITS */
     118    #     char dummy[sizeof(long)-sizeof(Ushort)]
     119    #     unsigned char buf[1]
     120    # cdef extern size_t bs_size
     121    # void bs_print(bitstring_t *b)
     122    # char *bs_desc(bitstring_t *b)
     123    # void bs_setlen(int l)
     124    # void bs_reset(bitstring_t *x)
     125    # bitstring_t *bs_alloc()
     126    # bitstring_t *bs_dup(bitstring_t *s)
     127    # void bs_free(bitstring_t *b)
     128    # void bs_set(bitstring_t *b, int i)
     129    # void bs_clear(bitstring_t *b, int i)
     130    # int bs_test(bitstring_t *b, int i)
     131    # bitstring_t *bs_read(FILE *f)
     132    # int bs_write(FILE *f, bitstring_t *b)
     133    # void bs_and(bitstring_t *dest, bitstring_t *src)
     134    # void bs_minus(bitstring_t *dest, bitstring_t *src)
     135    # void bs_or(bitstring_t *dest, bitstring_t *src)
     136    # int bs_match(bitstring_t *x, bitstring_t *y)
     137    # int bs_issub(bitstring_t *a, bitstring_t *b)
     138    # int bs_cmp(bitstring_t *a, bitstring_t *b)
     139    # void bs_cpy(bitstring_t *a, bitstring_t *b)
     140
     141    # file i/o
     142    # size_t zreadlong(FILE *f, long *buf, size_t n)
     143    # size_t zwritelong(FILE *f, long *buf, size_t n)
     144    # size_t zreadvec(FILE *f, PTR buf, size_t n)
     145    # size_t zwritevec(FILE *f, PTR buf, size_t n)
     146    # int zseek(FILE *f, long pos)
     147    # FILE *zreadhdr(char *name, long *fld, long *nor, long *noc)
     148    # FILE *zwritehdr(char *name, long fld, long nor, long noc)
     149
     150    # Gauss elimination and related functions
     151    # int zmkpivot(PTR matrix, long nor, long *piv)
     152    long zmkechelon(PTR matrix, long nor, long *piv)
     153    # int zcleanrow(PTR row, PTR matrix, long nor, long *piv)
     154    # int zcleanrow2(PTR row,PTR matrix,long nor,long *piv,PTR row2)
     155    # long znullsp(PTR matrix, long nor, long *piv, PTR nsp)
     156    # int zmatinv(PTR mat, PTR result)
     157
     158    #  Projection on the quotient
     159    # int zquotinit(PTR subspace, long dim, long *piv)
     160    # int zquot(PTR space, long dim, PTR quot)
     161    # int zquotop(PTR matrix, PTR quot)
     162
     163    # Seed vector generation
     164    # cdef extern PTR zpseed_vec                #/* The seed vector */
     165    # void zpseed_free()
     166    # int zpseed_init(long sdim, PTR sbasis)
     167    # long zpseed_make(long num)
     168    #long zpseed_next()
     169
     170    # Command line parsing
     171    # ctypedef struct proginfo_t:
     172    #     char *name
     173    #     char *shortdesc
     174    #     char *rcsrev
     175    #     char **helptext
     176    # cdef extern char opt_char         #/* Current option */
     177    # cdef extern char opt_text[50]     #/* Option text */
     178    # cdef extern char *opt_text_ptr    #/* Current position in option text */
     179    # cdef extern int opt_ind           #/* Index in argv[] list */
     180    # cdef extern char MeatAxeBinDir[250]       #/* MeatAxe program directory */
     181
     182    # void initargs(int argc, char **argv, proginfo_t *pi)
     183    # int zgetopt(char *pattern)
     184    # long getint()
     185
     186    # Messages
     187    # cdef extern int mtxerrno
     188    # cdef extern int mtxerraction
     189    # void mtxerror(char *text)
     190    # void errexit(int code, char *text)
     191    # int errhandler(int errno, char *errfile, int errline)
     192    # int Message(FILE *f, char *fmt, ...)
     193    # int ErrorExit(char *fmt, ...)
     194    # cdef extern int msg_level
     195
     196    # CPU time
     197    # void prtimes()
     198
     199    # Random numbers
     200    # void RandInit(unsigned seed)
     201    # long int Random()
     202
     203    # Miscellaneous
     204    # long gcd(long a, long b)
     205    # long lcm(long a, long b)
     206
     207    # Matrices
     208    ctypedef struct matrix_t:
     209        Ushort id                       #/* Always = T_MATRIX */
     210        long fl, nor, noc               #/* Field, #rows, #columns */
     211        PTR d                           #/* Pointer to data area */
     212
     213    # void matprint(char *name, matrix_t *p)
     214    # char *matdesc(matrix_t *s)
     215    matrix_t *matalloc(long fl, long nor, long noc)
     216    matrix_t *matmalloc(long fl, long nor, long noc)
     217    matrix_t *MTXmatid(long fl, long nor)
     218    void matfree(matrix_t *m)
     219
     220    long matcmp(matrix_t *m1, matrix_t *m2)
     221    matrix_t *matdup(matrix_t *src)
     222    # int  matmove(matrix_t *dest, matrix_t *src)
     223    matrix_t *_matextract(matrix_t *src, long first, long n)
     224    matrix_t *matread(FILE *f)
     225    int matwrite(FILE *f, matrix_t *mat)
     226    matrix_t *matload(char *fn)
     227    int matsave(matrix_t *mat, char *fn)
     228
     229    matrix_t *matadd(matrix_t *dest, matrix_t *src)
     230    matrix_t *matmul(matrix_t *dest, matrix_t *src)
     231    int matmul_result(matrix_t *dest, matrix_t *src, matrix_t *res)
     232    matrix_t *matmulF(matrix_t *src, FEL mark)
     233    matrix_t *matsub(matrix_t *dest, matrix_t *src)
     234    matrix_t *matneg(matrix_t *src)
     235    matrix_t *mattr(matrix_t *src)
     236    matrix_t *matinv(matrix_t *src)
     237    # int chkechelon(matrix_t *mat)
     238    matrix_t *matpower(matrix_t *mat, long n)
     239    long matorder(matrix_t *mat)
     240   
     241    # Higher-level gauss functions
     242    matrix_t *echelon(matrix_t *mat)
     243    long nullity(matrix_t *mat)
     244    matrix_t *nullspace(matrix_t *mat)
     245
     246    # ycomp.c
     247    # int spccomp(matrix_t *m1, matrix_t *m2, long n)
     248    # int spcequal(matrix_t *m1, matrix_t *m2)
     249    # int spccontains(matrix_t *m1, matrix_t *m2)
     250
     251    # yspin.c
     252    # matrix_t *matspin(matrix_t *seed, int ngen, matrix_t *gen[])
     253    # matrix_t *matspin_f(matrix_t *seed, int ngen, matrix_t *gen[], long *sdim)
     254
     255    # Spin-up, split, and standard basis
     256    # long zspinup(PTR space, long nseed, long *piv, int ngen, PTR gen[],       int gentype)
     257    # int zsbasis(PTR seed, long nseed, int ngen, PTR gen[], PTR space, long *piv, PTR basis)
     258    # cdef extern long *split_pivot
     259
     260    # ... various options in meataxe.h
     261    # int spinup(matrix_t *seed,int ngen,matrix_t *gen[],int options, matrix_t **subspace)
     262    # int split(matrix_t *subspace, int ngen, matrix_t *gen[], matrix_t **sub, matrix_t **quot)
     263    # matrix_t *quotproj(matrix_t *subspace, matrix_t *vectors)
     264    # matrix_t *sbasis(matrix_t *seed, int ngen, matrix_t **gen)
     265    # int chbasis(matrix_t *basis, int ngen, matrix_t **gen,matrix_t **newgen)
     266
     267    # Polynomials
     268    # ctypedef struct poly_t:
     269    #     Ushort id                     #/* Always = T_POLY */
     270    #     long fl
     271    #     long deg
     272    #     size_t size  #/* deg <= size-1 */
     273    #     FEL *buf
     274
     275    # char *poldesc(poly_t *s)
     276    # poly_t *polalloc(long fl, long degree)
     277    # void polfree(poly_t *x)
     278    # poly_t *poldup(poly_t *x)
     279    # poly_t *poladd(poly_t *dest, poly_t *src)
     280    # poly_t *polmul(poly_t *dest, poly_t *src)
     281    # poly_t *poldivmod(poly_t *a, poly_t *b)
     282    # poly_t *polmod(poly_t *a, poly_t *b)
     283    # void polprint(char *name, poly_t *p)
     284    # poly_t *polread(FILE *f)
     285    # int polwrite(FILE *f, poly_t *p)
     286    # int polcmp(poly_t *a, poly_t *b)
     287    # poly_t *polgcd(poly_t *a, poly_t *b)
     288    # poly_t *polderive(poly_t *p)
     289    # int vec2pol(PTR vec, poly_t *pol)
     290    # int pol2vec(poly_t *pol, PTR vec)
     291    # int polpack(poly_t *pol, PTR vec)
     292    # poly_t *polshiftmod(poly_t *p, long n, poly_t *q)
     293
     294    # The word generator
     295    # ctypedef struct wgdata_t:
     296    #     int ngen
     297    #     size_t len, max
     298    #     matrix_t **b
     299    # wgdata_t *WGInit(int ngen, matrix_t *gen[])
     300    # int WGFree(wgdata_t *b)
     301    # matrix_t *MakeWord(wgdata_t *b, long n)
     302    # void makefp(wgdata_t *b, long fp[])
     303    # char *SymbolicName(wgdata_t *b, long n)
     304
     305    # Permutations
     306    # ctypedef struct perm_t:
     307    #     Ushort id                     #/* Always = T_PERM */
     308    #     long deg
     309    #     PTR d
     310    # char *permdesc(perm_t *s)
     311    # void permprint(char *name, perm_t *x)
     312    # perm_t *permalloc(long deg)
     313    # void permfree(perm_t *m)
     314    # perm_t *permdup(perm_t *src)
     315    # perm_t *permmove(perm_t *dest, perm_t *src)
     316    # perm_t *permread(FILE *f)
     317    # int permwrite(FILE *f, perm_t *perm)
     318    # perm_t *permload(char *fn)
     319    # int permsave(perm_t *perm, char *fn)
     320    # perm_t *permmul(perm_t *dest, perm_t *src)
     321    # long permorder(perm_t *perm)
     322    # perm_t *permpower(perm_t *p, long n)
     323
     324    # Insertion
     325    # matrix_t *matinsert(matrix_t *mat, poly_t *pol)
     326    # matrix_t *matinsert_(matrix_t *mat, poly_t *pol)
     327
     328    # Factored polynomials
     329    # ctypedef struct fpoly_t:
     330    #     Ushort id                     #/* Always T_FPOLY */
     331    #     size_t len, max
     332    #     poly_t **p                    #/* Irreducible factors */
     333    #     long *e                       #/* Multiplicities */
     334    # char *fpoldesc(fpoly_t *p)
     335    # void fpolprint(char *name, fpoly_t *p)
     336    # fpoly_t *fpolalloc()
     337    # fpoly_t *fpoldup(fpoly_t *s)
     338    # void fpolfree(fpoly_t *p)
     339    # fpoly_t *fpolread(FILE *f)
     340    # int fpolwrite(FILE *f, fpoly_t *p)
     341    # fpoly_t *fpolmulp(fpoly_t *dest, poly_t *src, long pwr)
     342    # fpoly_t *fpolmul(fpoly_t *dest, fpoly_t *src)
     343
     344    # Characteristic and minimal polynomials
     345    # cdef extern long CharPolSeed
     346    # poly_t *charpolfactor(matrix_t *mat)
     347    # fpoly_t *charpol(matrix_t *mat)
     348    # poly_t *minpolfactor(matrix_t *mat)
     349    # fpoly_t *minpol(matrix_t *mat)
     350
     351    # Berlekamp factorization
     352    # fpoly_t *factorization(poly_t *pol)
     353
     354    # Sets
     355    # ctypedef struct set_t:
     356    #     Ushort id                     #/* Always = T_SET */
     357    #     size_t len, max
     358    #     long *buf
     359    # int set_allocstrategy(size_t first,size_t blocksize)
     360    # char *set_desc(set_t *s)
     361    # void set_print(char *name, set_t *x)
     362    # set_t *set_alloc()
     363    # set_t *set_dup(set_t *s)
     364    # void set_free(set_t *x)
     365    # int set_insert(set_t *set, long elem)
     366    # int set_contains(set_t *set, long elem)
     367    # set_t *set_read(FILE *f)
     368    # int set_write(FILE *f, set_t *set)
     369
     370    # Sequences
     371    # ctypedef struct sequence_t:
     372    #     Ushort id                     #/* Always = T_SEQUENCE */
     373    #     size_t len, max
     374    #     void **buf
     375    # char *seq_desc(sequence_t *s)
     376    # void seq_print(char *name, sequence_t *s)
     377    # sequence_t *seq_alloc(size_t size)
     378    # sequence_t *seq_dup(sequence_t *s)
     379    # void seq_free(sequence_t *x)
     380    # int seq_insert(sequence_t *s, int pos, void *x)
     381    # int seq_remove(sequence_t *s, int pos)
     382    # sequence_t *seq_read(FILE *f)
     383    # int seq_write(FILE *f, sequence_t *x)
     384
     385    # Integers and strings
     386    # ctypedef struct integer_t:
     387    #     Ushort id                     #/* Always = T_INT */
     388    #     long l
     389    # integer_t *intalloc(long l)
     390    # integer_t *intdup(integer_t *s)
     391    # char *intdesc(integer_t *s)
     392    # void intprint(char *name, integer_t *x)
     393    # void intfree(integer_t *x)
     394    # integer_t *intread(FILE *f)
     395    # int intwrite(FILE *f, integer_t *x)
     396    # integer_t *intadd(integer_t *d, integer_t *s)
     397    # integer_t *intsub(integer_t *d, integer_t *s)
     398    # integer_t *intmul(integer_t *d, integer_t *s)
     399    # integer_t *intdiv(integer_t *d, integer_t *s)
     400    # integer_t *intneg(integer_t *x)
     401
     402    # ctypedef struct string_t:
     403    #     Ushort id                     #/* Always = T_STR */
     404    #     char *s
     405    # string_t *stringalloc(char *c)
     406    # string_t *stringdup(string_t *s)
     407    # char *stringdesc(string_t *s)
     408    # void stringprint(char *name, string_t *x)
     409    # void stringfree(string_t *x)
     410    # string_t *stringread(FILE *f)
     411    # int stringwrite(FILE *f, string_t *x)
     412
     413    # MeatAxe object methods
     414    # void mtxfree(void *x)
     415    # void *mtxdup(void *x)
     416    # char *mtxgetname(void *x)
     417    # char *mtxgetdesc(void *x)
     418    # void *mtxread(FILE *f)
     419    # void mtxprint(void *x)
     420    # int mtxwrite(FILE *f, void *x)
     421    # void *mtxadd(void *l, void *r)
     422    # void *mtxsub(void *l, void *r)
     423    # void *mtxmul(void *l, void *r)
     424    # void *mtxdiv(void *l, void *r)
     425    # void *mtxpwr(void *l, void *r)
     426    # void *mtxneg(void *x)
     427    # void *mtxorder(void *x)
     428
     429    # Profiling
     430    # cdef extern long ProfFileIO               #/* File I/O */
     431    # cdef extern long ProfNullSpace    #/* Matrix null space */
     432    # cdef extern long ProfSpinUp               #/* Spin-up */
     433    # cdef extern long ProfMakeWord     #/* Make word */
     434    # cdef extern long ProfMatInsert    #/* Matrix insertion */
     435    # cdef extern long ProfPolFactor    #/* Polynomial factorization */
     436    # cdef extern long ProfCharPol      #/* Char. and minimal polynomial */
     437
     438    # Added by Simon King:
     439    cdef char* mat2str(matrix_t* m)           # Convert m->d into a char*
     440    cdef void str2mat(matrix_t* m, long l, char* x) # *COPY* the data from x to m->d
     441    cdef char* zrow2str(PTR d)       # Convert d into a char*
     442    cdef void zstr2row(PTR d, char* x) # *COPY* the data from x to d, assuming
     443                                       # that d has been assigned to at least
     444                                       # zrowsize_io byte.
     445
     446    cdef int multiply_strassen(matrix_t* A, matrix_t* B, matrix_t *C)
     447    cdef void set_cutoff(long rows, long size)
     448
     449# thats's all of meataxe.h !
     450
     451##############################################################
     452
     453#
     454# In the .pxd-declaration file, we merely need to define the
     455# cdef-parts of the extension class. Everything else is
     456# in the .pyx-definition file.
     457
     458from sage.matrix.matrix_dense cimport Matrix_dense
     459from sage.matrix.matrix0 cimport Matrix
     460from sage.structure.element cimport Element, ModuleElement, RingElement
     461cimport sage.structure.element
     462
     463cdef class FieldConverter_class:
     464    cdef object field
     465    cpdef inline object int_to_field(self, int x)
     466    cpdef inline int field_to_int(self, x)
     467
     468cdef class Matrix_modpn_dense(Matrix_dense):
     469    cdef matrix_t *Data
     470    cdef FieldConverter_class _converter
     471    cpdef Matrix_modpn_dense _mulInt_(self, long right)
     472    cpdef int base(self)
     473    cpdef int characteristic(self)
     474    cpdef list _matlist_(self)
     475    cpdef Matrix_modpn_dense transpose(Matrix_modpn_dense self)
     476    cpdef int order(Matrix_modpn_dense self) except -1
     477    cpdef Matrix_modpn_dense normalized(Matrix_modpn_dense self)
     478    cpdef Matrix_modpn_dense semi_echelon(Matrix_modpn_dense self)
     479    cpdef int nullity(Matrix_modpn_dense self)
     480    cpdef tuple lead(self)
     481    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value)
     482    cdef int get_unsafe_int(self, Py_ssize_t i, Py_ssize_t j)
  • new file sage/matrix/matrix_modpn_dense.pyx

    diff --git a/sage/matrix/matrix_modpn_dense.pyx b/sage/matrix/matrix_modpn_dense.pyx
    new file mode 100644
    - +  
     1#*****************************************************************************
     2#
     3#    Dense matrices over GF(p^n), p^n<255, using LibMeataxe as backend
     4#
     5#    Copyright (C) 2011 Simon A. King <simon.king@uni-jena.de>
     6#
     7#  Distributed under the terms of the GNU General Public License (GPL),
     8#  version 2 or later (at your choice)
     9#
     10#    This code is distributed in the hope that it will be useful,
     11#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     12#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     13#    General Public License for more details.
     14#
     15#  The full text of the GPL is available at:
     16#
     17#                  http://www.gnu.org/licenses/
     18#*****************************************************************************
     19r"""
     20Dense Matrices over `\mathbb F_q`, with `q<255` odd and not prime
     21
     22This module is a wrapper for a major modification of version 2.2.4 of the
     23Aachen `C-MeatAxe <http://www.math.rwth-aachen.de/homes/MTX/download.html>`_
     24and provides matrices over the finite field `\mathbb F_q`, where
     25`q\le 255` is odd and not prime.
     26
     27The modification strips down MeatAxe, which is a collection of executables,
     28to a C-library for linear algebra over fields of order at most 255. In
     29addition to some smaller improvements, the modifications include an
     30implementation of Winograd-Strassen multiplication.
     31
     32AUTHORS:
     33
     34- Simon King <simon.king@uni-jena.de>
     35
     36"""
     37
     38####################
     39#
     40# import sage types
     41#
     42####################
     43include "../ext/python.pxi"
     44include "../ext/python_object.pxi"
     45include "../ext/interrupt.pxi"
     46include "../ext/stdsage.pxi"
     47include "../ext/python_slice.pxi"
     48
     49import sage
     50import sage.all
     51from warnings import showwarning
     52import sys
     53from sage.structure.mutability import Mutability
     54from sage.all import Integer
     55from sage.all import cputime
     56from sage.all import GF
     57from sage.all import is_prime_power
     58from sage.all import factor
     59from sage.all import save, copy
     60from sage.all import cached_method, cached_function
     61from sage.all import SAGE_ROOT
     62from sage.all import current_randstate
     63
     64## The tables for MeatAxe computations will be found in the package data folder
     65import os
     66os.environ['MTXLIB'] = SAGE_ROOT + '/local/share/meataxe/'
     67
     68#########################################################################
     69#########################################################################
     70#
     71# Definition of the Matrix_modpn_dense extension type
     72#
     73#########################################################################
     74#########################################################################
     75from sage.rings.finite_rings.integer_mod import IntegerMod_int
     76
     77# Fast conversion from field to int and int to field
     78cdef class FieldConverter_class:
     79    """
     80    An auxiliary class, used to convert between <int> and finite field element
     81
     82    This class is for non-prime fields only. The method
     83    :meth:`int_to_field` exists for speed. The method
     84    :meth:`field_to_int` exists in order to have a common interface
     85    for elements of prime and non-prime fields; see
     86    :class:`PrimeFieldConverter_class`.
     87
     88    EXAMPLE::
     89
     90        sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
     91        sage: F.<y> = GF(125)
     92        sage: C = FieldConverter_class(F)
     93        sage: C.int_to_field(15)
     94        3*y
     95        sage: F.fetch_int(15)
     96        3*y
     97        sage: %timeit C.int_to_field(15)    #not tested
     98        625 loops, best of 3: 1.04 µs per loop
     99        sage: %timeit F.fetch_int(15)       #not tested
     100        625 loops, best of 3: 3.97 µs per loop
     101        sage: C.field_to_int(y)
     102        5
     103        sage: y.integer_representation()
     104        5
     105
     106    """
     107    def __init__(self, field):
     108        """
     109        INPUT:
     110
     111        A finite *non-prime* field. This assumption is not tested.
     112
     113        EXAMPLE::
     114
     115            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
     116            sage: F.<y> = GF(125)
     117            sage: C = FieldConverter_class(F)
     118            sage: C.int_to_field(15)
     119            3*y
     120            sage: F.fetch_int(15)
     121            3*y
     122            sage: C.field_to_int(y)
     123            5
     124            sage: y.integer_representation()
     125            5
     126
     127        """
     128        self.field = field._cache.fetch_int
     129    cpdef inline object int_to_field(self, int x):
     130        """
     131        Fetch a python int into the field.
     132
     133        EXAMPLE::
     134
     135            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
     136            sage: F.<y> = GF(125)
     137            sage: C = FieldConverter_class(F)
     138            sage: C.int_to_field(15)
     139            3*y
     140            sage: F.fetch_int(15)
     141            3*y
     142
     143        """
     144        return self.field(x)
     145    cpdef inline int field_to_int(self, x):
     146        """
     147        Represent a field element by a python int.
     148
     149        EXAMPLE::
     150
     151            sage: from sage.matrix.matrix_modpn_dense import FieldConverter_class
     152            sage: F.<y> = GF(125)
     153            sage: C = FieldConverter_class(F)
     154            sage: C.field_to_int(y)
     155            5
     156            sage: y.integer_representation()
     157            5
     158
     159        """
     160        return x.integer_representation()
     161
     162cdef class PrimeFieldConverter_class(FieldConverter_class):
     163    """
     164    An auxiliary class, used to convert between <int> and finite field element
     165
     166    This class is for prime fields only. The methods
     167    :meth:`int_to_field` and :meth:`field_to_int` exist in order to
     168    have a common interface for elements of prime and non-prime fields;
     169    see :class:`FieldConverter_class`.
     170
     171    EXAMPLE::
     172
     173        sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
     174        sage: F = GF(5)
     175        sage: C = PrimeFieldConverter_class(F)
     176        sage: C.int_to_field(int(2))
     177        2
     178        sage: F(2)
     179        2
     180        sage: C.field_to_int(F(2))
     181        2
     182        sage: int(F(2))
     183        2
     184
     185    """
     186    def __init__(self, field):
     187        """
     188        INPUT:
     189
     190        A finite *prime* field. This assumption is not tested.
     191
     192        EXAMPLE::
     193
     194            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
     195            sage: F = GF(5)
     196            sage: C = PrimeFieldConverter_class(F)
     197            sage: C.int_to_field(int(2))
     198            2
     199            sage: F(2)
     200            2
     201            sage: C.field_to_int(F(2))
     202            2
     203            sage: int(F(2))
     204            2
     205
     206        """
     207        self.field = field
     208    cpdef inline object int_to_field(self, int x):
     209        """
     210        Fetch a python int into the field.
     211
     212        EXAMPLE::
     213
     214            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
     215            sage: F = GF(5)
     216            sage: C = PrimeFieldConverter_class(F)
     217            sage: C.int_to_field(int(2))
     218            2
     219            sage: F(2)
     220            2
     221
     222        """
     223        return IntegerMod_int(self.field, x)
     224    cpdef inline int field_to_int(self, x):
     225        """
     226        Represent a field element by a python int.
     227
     228        EXAMPLE::
     229
     230            sage: from sage.matrix.matrix_modpn_dense import PrimeFieldConverter_class
     231            sage: F = GF(5)
     232            sage: C = PrimeFieldConverter_class(F)
     233            sage: C.field_to_int(F(2))
     234            2
     235            sage: int(F(2))
     236            2
     237
     238        """
     239        return int(x)
     240
     241cdef public dict _converter_cache = {}
     242cdef FieldConverter_class FieldConverter(field):
     243    """
     244    Return a :class:`FieldConverter_class` or :class:`PrimeFieldConverter_class` instance,
     245    depending whether the field is prime or not.
     246
     247    EXAMPLE::
     248
     249        sage: MS = MatrixSpace(GF(5^3,'y'),2)
     250        sage: A = MS.random_element()
     251        sage: A*2 == A+A    # indirect doctest
     252        True
     253        sage: MS = MatrixSpace(GF(5),2)
     254        sage: A = MS.random_element()
     255        sage: A*2 == A+A
     256        True
     257
     258    """
     259    try:
     260        return _converter_cache[field]
     261    except KeyError:
     262        if field.is_prime_field():
     263            return _converter_cache.setdefault(field, PrimeFieldConverter_class(field))
     264        return _converter_cache.setdefault(field, FieldConverter_class(field))
     265
     266class MTX_unpickle_class:
     267    """
     268    Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
     269
     270    EXAMPLES::
     271
     272        sage: M = MatrixSpace(GF(25,'z'),2)([1,2,3,4])
     273        sage: M == loads(dumps(M))   # indirect doctest
     274        True
     275
     276    """
     277    def __init__(self):
     278        """
     279        Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
     280
     281        TESTS::
     282
     283            sage: from sage.matrix.matrix_modpn_dense import MTX_unpickle_class
     284            sage: unpickle = MTX_unpickle_class()   # indirect doctest
     285            sage: M = MatrixSpace(GF(25,'z'),2)([1,2,3,4])
     286            sage: L = M.__reduce__()[1]
     287            sage: N = unpickle(*L)
     288            sage: print N
     289            [1 2]
     290            [3 4]
     291
     292        """
     293        self.__safe_for_unpickling__=True
     294
     295    def __call__(self,f,nr,nc,object Data,m):
     296        """
     297        Unpickling :class:`~sage.matrix.matrix_modpn_dense.Matrix_modpn_dense` matrices
     298
     299        TESTS::
     300
     301            sage: from sage.matrix.matrix_modpn_dense import MTX_unpickle_class
     302            sage: unpickle = MTX_unpickle_class()
     303            sage: M = MatrixSpace(GF(125,'z'),1,4)([1,2,3,4])
     304            sage: L = M.__reduce__()[1]
     305            sage: N = unpickle(*L)   # indirect doctest
     306            sage: print N
     307            [1 2 3 4]
     308
     309        """
     310        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     311        cdef char* x
     312        #cdef list L
     313        cdef int i
     314        cdef long ff
     315        cdef int NR
     316        cdef PTR p
     317        if isinstance(f, (int,long)):
     318            ff = f
     319            B = GF(f,'z')
     320            # TODO: MatrixSpace with "implementation"
     321            from sage.all import MatrixSpace
     322            parent = MatrixSpace(B, nr, nc)
     323        else:
     324            parent = f
     325            B = parent.base_ring()
     326            ff = B.order()
     327        OUT._parent = parent
     328        OUT._base_ring = B
     329        OUT._converter = FieldConverter(B)
     330        OUT._ncols = nc
     331        OUT._nrows = nr
     332        OUT._cache = {}
     333        if Data:
     334            OUT.Data = matalloc(ff, nr, nc)
     335            if not OUT.Data:
     336                raise MemoryError, "Not enough memory to allocate %dx%d MTY matrix"%(nr,nc)
     337            x=PyString_AsString(Data)
     338            str2mat(OUT.Data,len(Data),x)
     339        OUT._mutability = Mutability(m)
     340        return OUT
     341
     342mtx_unpickle = MTX_unpickle_class()
     343
     344
     345cdef class Matrix_modpn_dense(Matrix_dense):
     346    r"""
     347    Dense matrices over `\mathbb F_q`, `q<255` odd and not prime.
     348
     349    NOTE:
     350
     351    This class uses a major modification of the Aachen C-MeatAxe
     352    as backend. In principle, it would also work for prime fields
     353    and in characteristic two. However, other matrices in Sage,
     354    relying on linbox, m4ri or m4rie, are more efficient in these
     355    cases.
     356
     357    EXAMPLES::
     358
     359        sage: M = MatrixSpace(GF(25,'z'),2,3)([1,2,3,4,5,6])
     360        sage: print M
     361        [1 2 3]
     362        [4 0 1]
     363        sage: type(M)
     364        <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     365
     366    The documentation of the ``__init__`` methods shows further
     367    ways of creating a :class:`Matrix_modpn_dense` instance.
     368    However, these should only be of internal use.
     369
     370    """
     371
     372###################
     373## Init, Dealloc, Copy
     374    def __cinit__(self, parent=None, entries=None, **kwds):
     375        """
     376        TESTS::
     377
     378            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     379            sage: Matrix_modpn_dense.__new__(Matrix_modpn_dense)   # indirect doctest
     380            []
     381            sage: Matrix_modpn_dense(MatrixSpace(GF(64,'z'),4), None)
     382            [0 0 0 0]
     383            [0 0 0 0]
     384            [0 0 0 0]
     385            [0 0 0 0]
     386
     387        """
     388        if parent is None:  # this makes Matrix_modpn_dense.__new__(Matrix_modpn_dense) work,
     389                            # returning a non-initialised matrix
     390            return
     391        if isinstance(parent,basestring): # this allows to provide a file when initialising a matrix
     392            return
     393        cdef long f = parent.base_ring().order()
     394        cdef long nrows = parent.nrows()
     395        cdef long ncols = parent.ncols()
     396        if entries is None or entries==0:
     397            self.Data = matalloc(f, nrows, ncols)
     398            return
     399        if isinstance(entries, list):
     400            self.Data = matalloc(f, nrows, ncols)
     401            return
     402        if nrows!=ncols:
     403            raise TypeError, "Matrix ought to be a square matrix"
     404        self.Data = MTXmatid(f, <long>nrows)
     405
     406    def __init__(self, parent, entries=None, mutable=True, copy=False, coerce=False):
     407        """
     408        Matrix extension class using libmeataxe as backend
     409
     410        INPUT:
     411
     412        Instances of this class can be created by providing one of
     413        the following input data, where ``q<255`` is a prime power,
     414        ``m,n`` are non-negative integers, and `a_{11},...,a_{mn}`
     415        can be coerced into ``GF(q)``. Note that a user should
     416        create these instances via the matrix constructors; what
     417        we explain here is for internal use only!
     418
     419        - None => empty matrix over an unspecified field (used for unpickling)
     420        - a string ``f`` ==> load matrix from the file named ``f``
     421        - A matrix space of `m\\times n` matrices over GF(q) and either
     422
     423          - a list `[a_{11},a_{12},...,a_{1n},a_{21},...,a_{m1},...,a_{mn}]`,
     424            which results in a matrix with the given marks
     425          - ``None``, which is the fastest way to creata a zero matrix.
     426          - an element of GF(q), which results in a diagonal matrix with the
     427            given element on the diagonal.
     428
     429        If the optional parameter ``mutable`` is ``False`` (by default,
     430        it is ``True``), the resulting matrix can not be changed, and
     431        it can be used as key in a Python dictionary.
     432
     433        NOTE:
     434
     435        If ``q`` is not a prime power, a ``MemoryError`` will be raised.
     436
     437        EXAMPLES::
     438
     439            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     440
     441        1. Creating an empty matrix::
     442
     443            sage: Matrix_modpn_dense(None)
     444            []
     445
     446        2. Creating a zero (3x2)-matrix::
     447
     448            sage: Matrix_modpn_dense(MatrixSpace(GF(4,'z'),3,2))
     449            [0 0]
     450            [0 0]
     451            [0 0]
     452
     453        3. Creating a matrix from a list or list of lists::
     454
     455            sage: Matrix_modpn_dense(MatrixSpace(GF(5),2,3),[1,2,3,4,5,6])
     456            [1 2 3]
     457            [4 0 1]
     458            sage: Matrix_modpn_dense(MatrixSpace(GF(5),2,3),[[1,2,3],[4,5,6]])  # indirect doctest
     459            [1 2 3]
     460            [4 0 1]
     461
     462        4. Creating a diagonal matrix::
     463
     464            sage: M = Matrix_modpn_dense(MatrixSpace(GF(7),5),2); M
     465            [2 0 0 0 0]
     466            [0 2 0 0 0]
     467            [0 0 2 0 0]
     468            [0 0 0 2 0]
     469            [0 0 0 0 2]
     470
     471        5. Creating a matrix from a file in MeatAxe format.
     472
     473           First, we have to create that file; we use a temporary file,
     474           that will be removed when leaving Sage. Note that the method
     475           :meth:`msave` must be used, which does not use Python pickling
     476           but relies on the intrinsic C--MeatAxe way of saving.
     477           ::
     478
     479            sage: f = tmp_filename()
     480            sage: M.msave(f)
     481            sage: Matrix_modpn_dense(f)
     482            [2 0 0 0 0]
     483            [0 2 0 0 0]
     484            [0 0 2 0 0]
     485            [0 0 0 2 0]
     486            [0 0 0 0 2]
     487
     488        TESTS::
     489
     490            sage: MS = MatrixSpace(GF(125,'y'),2)
     491            sage: A = MS(0)
     492            sage: A.left_kernel()
     493            Vector space of degree 2 and dimension 2 over Finite Field in y of size 5^3
     494            Basis matrix:
     495            [1 0]
     496            [0 1]
     497            sage: A.right_kernel()
     498            Vector space of degree 2 and dimension 2 over Finite Field in y of size 5^3
     499            Basis matrix:
     500            [1 0]
     501            [0 1]
     502
     503        """
     504        if parent is None:
     505            self._mutability = Mutability(False)
     506            self._ncols = 0
     507            self._nrows = 0
     508            self._cache = {}
     509            return
     510        if isinstance(parent, basestring): # load from file
     511            FILE = os.path.realpath(parent)
     512            try:
     513                fsock = open(FILE,"rb",0)
     514                fsock.close()
     515            except (OSError,IOError):
     516                return
     517            self.Data = matload(FILE)
     518            if zsetfield(self.Data.fl):
     519                raise ValueError, "Invalid saved date"
     520            B = GF(self.Data.fl, 'z')
     521            # TODO: MatrixSpace with "implementation"
     522            from sage.all import MatrixSpace
     523            parent = MatrixSpace(B, self.Data.nor, self.Data.noc)
     524            self._mutability = Mutability(False)
     525            self._parent = parent
     526            self._base_ring = B
     527            self._converter = FieldConverter(B)
     528            self._ncols = self.Data.noc
     529            self._nrows = self.Data.nor
     530            self._cache = {}
     531            return
     532
     533        if not self.Data: # should be initialised by __cinit__
     534            raise MemoryError, "Error allocating memory for MeatAxe matrix"
     535        Matrix_dense.__init__(self, parent)
     536        self._mutability=Mutability(not mutable)
     537        B = self._base_ring
     538        self._converter = FieldConverter(B)
     539        if entries is None:
     540            return
     541
     542        cdef int i,j
     543        cdef int f
     544        if not isinstance(entries,list): # __cinit__ initialized the unit matrix
     545            f = self._converter.field_to_int(self._coerce_element(entries))
     546            if matmulF(self.Data,zitof(f)):
     547                return
     548            else:
     549                raise ArithmeticError, "Matrix sizes or fields not compatible"
     550        cdef PTR x
     551        cdef PTR y
     552        x = self.Data.d
     553        cdef int nr = self.Data.nor
     554        cdef int nc = self.Data.noc
     555        if nr==0 or nc==0:
     556            return
     557        if len(entries)<nr:
     558            raise ValueError, "Expected a list of size at least the number of rows"
     559        cdef int idx
     560        cdef list dt, dt_i
     561        if isinstance(entries[0],list):
     562            dt = entries
     563            for i from 0 <= i < nr:
     564                y = x
     565                idx = 0
     566                dt_i = dt[i]
     567                #dtnext = dt[i].__iter__().next
     568                for j from 0 <= j < nc:
     569                    y[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(dt_i[j])))]
     570                    idx += 1
     571                    if idx == MPB:
     572                        y+=1
     573                        idx=0
     574                    # zinsert_step is not the fastest, since tnull[*idx] and tinsert[*idx] are looked up over and over again.
     575                zsteprow(&(x))
     576        else:
     577            dtnext = entries.__iter__().next
     578            for i from 0 <= i < nr:
     579                y = x
     580                idx = 0
     581                for j from 0 <= j < nc:
     582                    bla = dtnext()
     583
     584                    y[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(bla)))] #dtnext())))]
     585                    idx += 1
     586                    if idx == MPB:
     587                        y+=1
     588                        idx=0
     589                zsteprow(&(x))
     590                           
     591    def __dealloc__(self):
     592        """
     593        Deallocate this matrix
     594
     595        TESTS::
     596
     597            sage: M = MatrixSpace(GF(25,'z'),2,4)([1,2,3,4,4,3,2,1])
     598            sage: del M    # indirect doctest
     599            sage: M
     600            Traceback (most recent call last):
     601            ...
     602            NameError: name 'M' is not defined
     603        """
     604        if self.Data:
     605            matfree(self.Data)
     606        #cdef matrix_t *p
     607        self.Data = NULL
     608
     609    def __copy__(self):
     610        """
     611        Return a mutable copy of this matrix
     612
     613        EXAMPLES:
     614
     615        The following example caused a problem in a very
     616        early stage of wrapping C-MeatAxe::
     617
     618            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     619            sage: M = Matrix_modpn_dense(MatrixSpace(GF(2),3,20), [20*[0],20*[0],[1]+19*[0]], mutable=False)
     620            sage: M.is_mutable()
     621            False
     622            sage: N = copy(M)    # indirect doctest
     623            sage: N.is_mutable()
     624            True
     625            sage: print N
     626            [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
     627            [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
     628            [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
     629            sage: N==M
     630            True
     631            sage: N is M
     632            False
     633
     634        """       
     635        cdef Matrix_modpn_dense retval = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     636        # Do the initialisation "manually"
     637        retval._mutability = Mutability(False)
     638        retval._parent = self._parent
     639        retval._base_ring = self._base_ring
     640        retval._converter = self._converter
     641        retval._ncols = self._ncols
     642        retval._nrows = self._nrows
     643        retval._cache = {}
     644        if self.Data:
     645            retval.Data = matdup(self.Data)
     646            if not retval.Data:
     647                raise MemoryError, "Error copying a %s instance"%repr(type(self))
     648        return(retval)
     649   
     650    def __reduce__(self):
     651        r"""
     652        Pickle/save this matrix.
     653
     654        TESTS::
     655
     656            sage: M = MatrixSpace(GF(25,'z'),1,4)([1,2,3,4])
     657            sage: M.__reduce__()[0]
     658            <sage.matrix.matrix_modpn_dense.MTX_unpickle_class instance at ...>
     659            sage: M.__reduce__()[1]
     660            (Full MatrixSpace of 1 by 4 dense matrices over Finite Field in z of size 5^2, 1, 4, '\x01\x02\x03\x04', False)
     661            sage: loads(dumps(M))
     662            [1 2 3 4]
     663            sage: loads(dumps(M)) == M  # indirect doctest
     664            True
     665
     666        """
     667        cdef char* d
     668        cdef L = []
     669        Lnext = L.append
     670        cdef int i,NR
     671        cdef PTR p
     672        if self.Data:
     673            zsetfield(self.Data.fl)
     674            zsetlen(self.Data.noc)
     675
     676            NR = self.Data.nor
     677            p = self.Data.d
     678            for i from 0 <= i < NR:
     679                d = zrow2str(p)
     680                Lnext(PyString_FromStringAndSize(d,zrowsize_io))
     681                zsteprow(&(p))
     682            return mtx_unpickle, (self._parent,self.Data.nor,self.Data.noc,''.join(L),self._mutability._is_immutable) # or ...,L,self._mutable)
     683        else:
     684            return mtx_unpickle, (0,0,0,'',self._mutability._is_immutable)
     685
     686
     687#######################
     688## String representation is taken care of by implementing get_unsafe
     689
     690
     691##################
     692## ==,<,>, dictionary
     693    def __richcmp__(left, right, int op):
     694        return (<Element>left)._richcmp(right, op)
     695   
     696    def __cmp__(left, right):
     697        return (<Element>left)._cmp(right)
     698
     699    cdef int _cmp_c_impl(self, Element right) except -2:
     700        """
     701        Compare two MeatAxe matrices
     702       
     703        Of course, '<' and '>' doesn't make much sense for matrices.
     704
     705        EXAMPLES::
     706
     707            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     708            sage: M = Matrix_modpn_dense(MatrixSpace(GF(2), 3,20), [20*[0],20*[0],[1]+19*[0]])
     709            sage: N = copy(M)
     710            sage: M == N   # indirect doctest
     711            True
     712            sage: M != N
     713            False
     714            sage: M<N
     715            False
     716            sage: N[2,19] = 1
     717            sage: M == N
     718            False
     719            sage: M != N
     720            True
     721        """
     722        return matcmp((<Matrix_modpn_dense>self).Data, (<Matrix_modpn_dense>right).Data)
     723
     724    def __hash__(self):
     725        """
     726        Return a hash value for M, provided M is immutable
     727
     728        EXAMPLES::
     729
     730            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     731            sage: M = Matrix_modpn_dense(MatrixSpace(GF(5),1,4),[[4,3,2,1]],mutable=False)
     732            sage: if sys.byteorder == 'little':
     733            ...     print hash(M) == 7606091044269354279  # indirect doctest
     734            ... else:
     735            ...     print hash(M) == 1060097699   # indirect doctest
     736            True
     737
     738        """
     739        cdef char* d
     740        if not self.Data:
     741            return 0
     742        else:
     743            zsetfield(self.Data.fl)
     744            zsetlen(self.Data.noc)
     745            d = mat2str(self.Data)
     746            return hash((self.Data.fl,self.Data.nor,self.Data.noc,PyString_FromStringAndSize(d,zsize(self.Data.nor))))
     747
     748##########################
     749## Saving should be done via pickling
     750## However, we kept a method that relies on MeatAxe matsave:
     751    def msave(self,f):
     752        """
     753        Save self into file, using the original C--MeatAxe format
     754     
     755        INPUT:
     756
     757        ``f`` - a file name (string)
     758
     759        NOTE:
     760
     761        A matrix ``M`` should usually be saved by ``save(M,'filename')``
     762        and reloaded with ``load('filename')``. Saving with
     763        ``M.msave('filename')`` uses a totally different format, so that
     764        the data can not be reloaded with ``load`` but with
     765        ``Matrix_modpn_dense('filename')``.
     766
     767        EXAMPLES:
     768
     769        First, we create a temporary file name, that will be removed
     770        when leaving Sage. Of course, in an application a permanent
     771        file would usually be chosen.
     772        ::
     773       
     774            sage: f = tmp_filename()+'.sobj'
     775
     776        Next, we create a matrix and use msave::
     777       
     778            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     779            sage: M = MatrixSpace(GF(25,'z'),2,4)([1,2,3,4,4,3,2,1])
     780            sage: M.msave(f)
     781
     782        It can not be reloaded with 'load'::
     783       
     784            sage: N = load(f)
     785            Traceback (most recent call last):
     786            ...
     787            UnpicklingError: invalid load key, ...
     788
     789        But it can be reloaded with the class constructor::
     790       
     791            sage: N = Matrix_modpn_dense(f)
     792            sage: M == N
     793            True
     794        """
     795        matsave(self.Data,f)
     796
     797##########################
     798## Structural parts
     799
     800    cpdef int base(self):
     801        """
     802        Return the order of the field over which the matrix is defined
     803
     804        EXAMPLES::
     805
     806            sage: M = MatrixSpace(GF(25,'z'),1,4)([4,3,2,1])
     807            sage: M.base()
     808            25
     809        """
     810        if self.Data:
     811            return self.Data.fl
     812        else:
     813            return None
     814 
     815    cpdef int characteristic(self):
     816        """
     817        Return the characteristic of the field over which the matrix is defined
     818
     819        EXAMPLES::
     820
     821            sage: M = MatrixSpace(GF(25,'z'),1,4)([4,3,2,1])
     822            sage: M.characteristic()
     823            5
     824        """
     825        if self.Data:
     826            return factor(self.Data.fl)[0][0]
     827        else:
     828            return None
     829
     830    def matrix_from_rows(self, rows):
     831        """
     832        Matrix formed by a given list of rows from self.
     833
     834        INPUT:
     835
     836        ``rows`` -- list of integers, providing the rows that shall be put together
     837
     838        OUTPUT:
     839
     840        A matrix that is formed by the indicated rows of self, in the given order.
     841
     842        EXAMPLES::
     843
     844            sage: MS = MatrixSpace(GF(5^3,'y'),3,2)
     845            sage: A = MS(range(6))
     846            sage: A
     847            [0 1]
     848            [2 3]
     849            [4 0]
     850            sage: A.matrix_from_rows([2,1])
     851            [4 0]
     852            [2 3]
     853            sage: A.matrix_from_rows([1,1])
     854            [2 3]
     855            [2 3]
     856
     857        """
     858        if not (PY_TYPE_CHECK(rows, list) or PY_TYPE_CHECK(rows, tuple)):
     859            raise TypeError, "rows must be a list of integers"
     860        cdef Matrix_modpn_dense A
     861        cdef int k
     862        cdef Py_ssize_t nrows,c
     863
     864        nrows = PyList_GET_SIZE(rows)
     865        A = self.new_matrix(nrows = nrows)
     866
     867        cdef PTR src,dest
     868        dest = A.Data.d
     869        for k in rows:
     870            src = self.Data.d
     871            zadvance(&(src),k)
     872            zcpyrow(dest,src,zrowsize)
     873            zsteprow(&(dest))
     874        return A
     875
     876    def __getitem__(self,key):
     877        """
     878        Return element, row, or slice of self.
     879
     880        NOTE:
     881
     882        We could simply inherit the __getitem__ method form
     883        the base class of dense matrices. Hover, this would
     884        be a lot slower.
     885
     886        INPUT:
     887
     888        - a tuple (i,j) with i, j non-negative integers, or
     889        - a non-negative integer, or
     890        - a slice object, created via [i:j]
     891
     892        USAGE:
     893
     894        - A[i, j] -- the i,j of A, or
     895        - A[i]    -- the i-th row of A, or
     896        - A[i:j]  -- the i-th through (j-1)-st rows of A.
     897
     898        EXAMPLES::
     899
     900            sage: M = MatrixSpace(GF(9,'z'),3)(list(GF(9,'z')))
     901            sage: M
     902            [      0     2*z   z + 1]
     903            [  z + 2       2       z]
     904            [2*z + 2 2*z + 1       1]
     905            sage: M[1,2]    # indirect doctest
     906            z
     907            sage: M[1]
     908            (z + 2, 2, z)
     909            sage: M[1:3]
     910            [  z + 2       2       z]
     911            [2*z + 2 2*z + 1       1]
     912        """
     913        if not self.Data:
     914            raise IndexError, "Matrix is empty"
     915        zsetfield(self.Data.fl)
     916        cdef matrix_t *mtxrow
     917        cdef PTR p
     918        cdef Matrix_modpn_dense A
     919        cdef int i,j
     920        cdef tuple key_tuple
     921        cdef int nr=self.Data.nor
     922        cdef int nc=self.Data.noc
     923        zsetlen(nc)
     924        if PyTuple_CheckExact(key):
     925            # key is a tuple, so we get i and j efficiently
     926            # TODO: Support pairs of slices
     927            key_tuple = <tuple>key
     928            if len(key_tuple) != 2:
     929                raise IndexError("index must be an integer or pair of integers")
     930            try:
     931                i = <object>PyTuple_GET_ITEM(key_tuple, 0)
     932                j = <object>PyTuple_GET_ITEM(key_tuple, 1)
     933            except TypeError:
     934                raise IndexError("index must be an integer or pair of integers")
     935            if (i<0):
     936                i += nr
     937            if (i<0) or (i>=nr):
     938                raise IndexError, "matrix index out of range"
     939            if j<0:
     940                j+=nc
     941            if j<0 or j>=nc:
     942                raise IndexError, "matrix index out of range"
     943            p = self.Data.d
     944            zadvance(&(p),i)
     945            return self._converter.int_to_field(zftoi(zextract(p,j+1)))
     946        if PySlice_Check(key):
     947            # Slice interpretation is passed to the sequence
     948            # constructed by range
     949            if (key.start>=0) and (key.stop<=nr) and \
     950               (key.start<=nr) and (key.stop>=0):
     951                A=Matrix_modpn_dense(None)
     952                A.Data = _matextract(self.Data, long(key.start+1), long(key.stop-key.start))
     953                if A.Data:
     954                    A._mutability = copy(self._mutability)
     955                    from sage.all import MatrixSpace
     956                    A._parent = MatrixSpace(self._base_ring, key.stop-key.start, self._ncols)
     957                    A._converter = self._converter
     958                    A._base_ring = self._base_ring
     959                    A._ncols = self._ncols
     960                    A._nrows = key.stop-key.start
     961                    A._cache = {}
     962                    return A
     963                else:
     964                    raise IndexError, "matrix index out of range"
     965            else:
     966                raise IndexError, "matrix index out of range"
     967        try:
     968            if (key>=0) and (key<nr):
     969                return self.row(key)  # this is inherited
     970        except:
     971            raise IndexError, "matrix index out of range"
     972        raise IndexError, "matrix index out of range"
     973
     974    def __setitem__(self, key, el):
     975        """
     976        x.__setitem__(key, el) <==> x[key]=el
     977
     978        USAGE:
     979
     980        - key= pair of integers, el= integer
     981        - key= integer, el= list/tuple of integers
     982          row number <key> is filled with the entries of <el>,
     983          until the end of the row or the end of <el> is reached
     984        - key= integer, el= Matrix_modpn_dense with the same number of columns as self
     985          The rows of self starting with number <key> are filled with the
     986          contents of el.
     987
     988        EXAMPLES::
     989
     990            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     991            sage: F.<z> = GF(9)
     992            sage: M = MatrixSpace(F,3)(list(F)); M
     993            [      0     2*z   z + 1]
     994            [  z + 2       2       z]
     995            [2*z + 2 2*z + 1       1]
     996            sage: M[0] = [1,2,3]   # indirect doctest
     997            sage: M
     998            [      1       2       0]
     999            [  z + 2       2       z]
     1000            [2*z + 2 2*z + 1       1]
     1001            sage: M[1,2] = 0
     1002            sage: M
     1003            [      1       2       0]
     1004            [  z + 2       2       0]
     1005            [2*z + 2 2*z + 1       1]
     1006           
     1007            sage: N = MatrixSpace(F,2,3)([z,z^2,z^3,z^3,z^4,z^5])
     1008            sage: N
     1009            [      z   z + 1 2*z + 1]
     1010            [2*z + 1       2     2*z]
     1011            sage: M[1] = N
     1012            sage: M
     1013            [      1       2       0]
     1014            [      z   z + 1 2*z + 1]
     1015            [2*z + 1       2     2*z]
     1016
     1017        """
     1018        self.check_mutability()
     1019        cdef PTR x
     1020        cdef int j,mini,idx
     1021        cdef Matrix_modpn_dense tmp
     1022        cdef list value_list
     1023        if self.Data:
     1024            zsetfield(self.Data.fl)
     1025        else:
     1026            raise IndexError, "Matrix is empty"
     1027        if isinstance(key,tuple):
     1028            try:
     1029                j,idx = key
     1030            except (TypeError, ValueError):
     1031                raise TypeError, "pair of integers expected"
     1032            if (0<=j<self.Data.nor) and (0<=idx<self.Data.noc):
     1033                x = self.Data.d
     1034                zsetlen(self.Data.noc)
     1035                zadvance(&(x),j)
     1036                zinsert(x, idx+1, zitof(self._converter.field_to_int(self._coerce_element(el))))
     1037            else:
     1038                raise IndexError, "index out of range"
     1039        elif (isinstance(key,Integer) or isinstance(key,int)): # insert a list or another MTX
     1040            j = key
     1041            if (j>=0) and (j<self.Data.nor):
     1042                x = self.Data.d
     1043                zsetlen(self.Data.noc)
     1044                zadvance(&(x),j)
     1045                if isinstance(el,Matrix_modpn_dense):
     1046                    tmp = el
     1047                    if tmp.Data.fl!=self.Data.fl:
     1048                        raise ArithmeticError, "Fields are incompatible"
     1049                    if tmp.Data.nor+j>self.Data.nor:
     1050                        raise IndexError, "Matrix to be inserted has too many rows"
     1051                    if tmp.Data.noc!=self.Data.noc:
     1052                        raise IndexError, "Matrix to be inserted must have the same number of columns as self"
     1053                    zsetlen(tmp.Data.noc)
     1054                    zcpyrow(x,tmp.Data.d,zsize(tmp.Data.nor))
     1055                    zsetlen(self.Data.noc)
     1056                    return
     1057                value_list = el  # type error, if no list.
     1058                if len(value_list)!=self.Data.noc:
     1059                    raise ValueError, "List length differs from number of columns"
     1060                idx = 0
     1061                memset(x, 0, zrowsize)
     1062                for v in value_list:
     1063                    x[0] += tinsert[idx][zitof(self._converter.field_to_int(self._coerce_element(v)))]
     1064                    idx += 1
     1065                    if idx == MPB:
     1066                        x+=1
     1067                        idx=0
     1068            else:
     1069                raise IndexError, "index out of range"
     1070        else:
     1071            raise IndexError, "index must be an integer or a pair of integers"
     1072
     1073    def _rowlist_(self, i, j=-1):
     1074        """
     1075        M._rowlist_(i): Return row <i> as a list of integers
     1076        M._rowlist_(i,j): Return rows <i> up to and including <j> in a single list
     1077       
     1078        EXAMPLES::
     1079
     1080            sage: F.<z> = GF(9)
     1081            sage: M = MatrixSpace(F,3)(list(F))
     1082            sage: M._rowlist_(0)
     1083            [0, 6, 4]
     1084
     1085        We verify that these integers really correspond
     1086        to the actual elements of the top row::
     1087
     1088            sage: M[0]
     1089            (0, 2*z, z + 1)
     1090            sage: (2*z).integer_representation()
     1091            6
     1092            sage: (z+1).integer_representation()
     1093            4
     1094            sage: M._rowlist_(1,2)
     1095            [5, 2, 3, 8, 7, 1]
     1096        """
     1097        if self.Data:
     1098            zsetfield(self.Data.fl)
     1099            zsetlen(self.Data.noc)
     1100        else:
     1101            return []
     1102        if not(isinstance(i,int) or isinstance(i,Integer)):
     1103            raise TypeError, "Index must be an integer"
     1104        if (i<0) or (i>=self.Data.nor):
     1105            raise IndexError, "Index out of range"
     1106        cdef PTR p,q
     1107        cdef int k,l,ii,jj, idx
     1108        p = self.Data.d
     1109        zadvance(&(p),i)
     1110        idx = 0
     1111        q = p
     1112        cdef list L = [zftoi(zextract_step(&(q),&(idx))) for k in range(self.Data.noc)]
     1113        if j!=-1:
     1114            if not(isinstance(j,int) or isinstance(j,Integer)):
     1115                raise TypeError, "Second index must be an integer"
     1116            if j >= self.Data.nor:
     1117                raise IndexError, "Index out of range"
     1118            ii=i
     1119            jj=j
     1120            for k from ii < k <= jj:
     1121                zsteprow(&(p))
     1122                idx = 0
     1123                q = p
     1124                for l from 0<=l<self.Data.noc:
     1125                    L.append(zftoi(zextract_step(&(q),&(idx))))
     1126        return L
     1127
     1128    cpdef list _matlist_(self):
     1129        """
     1130        M._matlist_(): Return the entries of M as a list of lists of integers
     1131
     1132        EXAMPLES::
     1133
     1134            sage: F.<z> = GF(9)
     1135            sage: M = MatrixSpace(F,3)(list(F))
     1136            sage: M._matlist_()
     1137            [[0, 6, 4], [5, 2, 3], [8, 7, 1]]
     1138        """
     1139        if self.Data:
     1140            zsetfield(self.Data.fl)
     1141            zsetlen(self.Data.noc)
     1142        else:
     1143            raise IndexError, "Matrix is empty"
     1144        cdef PTR p,q
     1145        cdef int i,j
     1146        p = self.Data.d
     1147        cdef int idx
     1148        cdef list l_out=[]
     1149        for i from 1<=i<self.Data.nor:
     1150            q = p
     1151            idx = 0
     1152            l_out.append([zftoi(zextract_step(&(q),&(idx))) for j in range(self.Data.noc)])
     1153            zsteprow(&(p))
     1154        idx = 0
     1155        l_out.append([zftoi(zextract_step(&(p),&(idx))) for j in range(self.Data.noc)])
     1156        return l_out
     1157
     1158#########################
     1159## Arithmetics
     1160    cpdef ModuleElement _add_(self, ModuleElement other):
     1161        """
     1162        Add two MeatAxe matrices of equal size
     1163
     1164        EXAMPLES::
     1165
     1166            sage: F.<z> = GF(25)
     1167            sage: MS = MatrixSpace(F, 2,4)
     1168            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
     1169            sage: N = MS([1,z,z^2,z^3, z^4,z^5,z^5,z^7])
     1170            sage: M + N    # indirect doctest
     1171            [      2   z + 2   z + 1 4*z + 2]
     1172            [3*z + 2   z + 1 2*z + 1       z]
     1173           
     1174        We verify that the answer is correct::
     1175
     1176            sage: (M+N).list() == [a+b for a,b in zip(M.list(), N.list())]
     1177            True
     1178
     1179        """
     1180        cdef Matrix_modpn_dense right = other
     1181        if ((self.Data.noc==0) or (right.Data.noc==0)):
     1182            raise NotImplementedError, "The matrices must not be empty"
     1183        cdef Matrix_modpn_dense left = self.__copy__()
     1184        if not right._mutability.is_mutable():
     1185            left._mutability.set_immutable()
     1186        if matadd(left.Data,right.Data):
     1187            return left
     1188        else:
     1189            raise ArithmeticError, "Matrix sizes or fields not compatible"
     1190
     1191    cpdef ModuleElement _sub_(self, ModuleElement other):
     1192        """
     1193        Subtract two MeatAxe matrices of equal size
     1194
     1195        EXAMPLES::
     1196
     1197            sage: F.<z> = GF(25)
     1198            sage: MS = MatrixSpace(F, 2,4)
     1199            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
     1200            sage: N = MS([1,z,z^2,z^3, z^4,z^5,z^5,z^7])
     1201            sage: M - N    # indirect doctest
     1202            [      0 4*z + 2     4*z   z + 1]
     1203            [4*z + 3 3*z + 4 4*z + 4     2*z]
     1204
     1205        We verify that the answer is correct::
     1206
     1207            sage: (M-N).list() == [a-b for a,b in zip(M.list(), N.list())]
     1208            True
     1209
     1210        """
     1211        cdef Matrix_modpn_dense right = other
     1212        if (self.Data.fl==2):
     1213            return self.__add__(right)
     1214        if ((self.Data.noc==0) or (right.Data.noc==0)):
     1215            raise NotImplementedError, "The matrices must not be empty"
     1216        cdef Matrix_modpn_dense left = self.__copy__()
     1217        if right.is_immutable():
     1218            left.set_immutable()
     1219        if matsub(left.Data,right.Data):
     1220            return left
     1221        else:
     1222            raise ArithmeticError, "Matrix sizes or fields not compatible"
     1223
     1224    cpdef ModuleElement _neg_(self):
     1225        """
     1226        Negation of a MeatAxe matrix
     1227
     1228        EXAMPLES::
     1229
     1230            sage: F.<z> = GF(25)
     1231            sage: MS = MatrixSpace(F, 2,4)
     1232            sage: M = MS([1,2,3,4, z,2*z,3*z,4*z])
     1233            sage: -M    # indirect doctest
     1234            [  4   3   2   1]
     1235            [4*z 3*z 2*z   z]
     1236            sage: (-M) + M
     1237            [0 0 0 0]
     1238            [0 0 0 0]
     1239
     1240        """
     1241        if (self.Data.fl==2):
     1242            return self.__copy__()
     1243        if (self.ncols()==0):
     1244            raise NotImplementedError, "The matrix must not be empty"
     1245        cdef Matrix_modpn_dense result = self.__copy__()
     1246        if matneg(result.Data):
     1247            return result
     1248        else:
     1249            raise ArithmeticError, "Something went wrong"
     1250
     1251    # The default __mul__ relies on _multiply_classical, _multiply_strassen,
     1252    #_strassen_default_cutoff, _will_use_strassen, and _lmul_
     1253    def _multiply_classical(left, Matrix_modpn_dense other):
     1254        """
     1255        Multiply two MeatAxe matrices.
     1256
     1257        NOTE:
     1258
     1259        This method is implicitly called when multiplying matrices.
     1260        Large matrices will in instead use the asymptotically fast
     1261        Strassen-Winograd multiplication algorithm.
     1262
     1263        EXAMPLES::
     1264
     1265            sage: F.<z> = GF(9)
     1266            sage: M = MatrixSpace(F,3)(list(F))
     1267            sage: N = MatrixSpace(F,3)(list(reversed(list(F))))
     1268            sage: type(M)
     1269            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1270            sage: M
     1271            [      0     2*z   z + 1]
     1272            [  z + 2       2       z]
     1273            [2*z + 2 2*z + 1       1]
     1274            sage: N
     1275            [      1 2*z + 1 2*z + 2]
     1276            [      z       2   z + 2]
     1277            [  z + 1     2*z       0]
     1278            sage: M*N                   # indirect doctest
     1279            [2*z + 1 2*z + 2       2]
     1280            [    2*z       1   z + 1]
     1281            [      2   z + 2       z]
     1282
     1283        The result is mutable if and only if at least one of the factors
     1284        is mutable::
     1285
     1286            sage: M.is_mutable()
     1287            True
     1288            sage: (M*N).is_mutable()
     1289            True
     1290            sage: M.set_immutable()
     1291            sage: (M*N).is_mutable()
     1292            True
     1293            sage: N.set_immutable()
     1294            sage: (M*N).is_mutable()
     1295            False
     1296
     1297        """
     1298        cdef Matrix_modpn_dense Left = left
     1299        if (Left.ncols()==0) or (other.ncols()==0):
     1300            raise ValueError, "The matrices must not be empty"
     1301        if Left.Data.fl != other.Data.fl:
     1302            raise ValueError, "Different base fields"
     1303        from sage.all import MatrixSpace
     1304        parent = MatrixSpace(Left._base_ring, Left.Data.nor, other.Data.noc)
     1305        cdef Matrix_modpn_dense out = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     1306        out.Data = matmalloc(Left.Data.fl, Left.Data.nor, other.Data.noc)
     1307        if matmul_result(Left.Data, other.Data, out.Data) != 0:
     1308            raise ValueError, "Matrix sizes are not compatible"
     1309        out._mutability = Mutability(Left._mutability._is_immutable and other._mutability._is_immutable)
     1310        out._parent = parent
     1311        out._base_ring = Left._base_ring
     1312        out._converter = Left._converter
     1313        out._ncols = other._ncols
     1314        out._nrows = Left._nrows
     1315        out._cache = {}
     1316        return out
     1317
     1318    cdef bint _will_use_strassen(self, Matrix right) except -2:
     1319        cutoff = self._strassen_default_cutoff(right)
     1320        if cutoff<0:
     1321            return False
     1322        if self._nrows < 2*cutoff:
     1323            return False
     1324        zsetfield(self.base_ring().order())
     1325        zsetlen(self._ncols)
     1326        if zrowsize_io<2*cutoff:
     1327            return False
     1328        return True
     1329
     1330    def _multiply_strassen(self, Matrix_modpn_dense right, long cutoff=0):
     1331        """
     1332        Multiplication using the Strassen-Winograd algorithm.
     1333
     1334        INPUT:
     1335
     1336        - Another matrix -- the second factor of the product
     1337        - ``cutoff`` -- optional integer, giving the size
     1338          below which the algorithm aborts the recursion
     1339          and reverts to classical (school book) multiplication.
     1340
     1341        NOTE:
     1342
     1343        By default, Strassen-Winograd multiplication is used
     1344        if the number of rows and the byte length of a row
     1345        exceeds 220.
     1346
     1347        EXAMPLE::
     1348
     1349            sage: MS = MatrixSpace(GF(9,'z'),1500,1500)
     1350            sage: A = MS.random_element()
     1351            sage: B = MS.random_element()
     1352            sage: type(A)
     1353            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1354            sage: C = A*B   # Here, Strassen-Winograd multiplication is the default
     1355            sage: C == A._multiply_classical(B)
     1356            True
     1357
     1358        In the next example, we force the use of MeatAxe matrices,
     1359        in order to demonstrate that Strassen-Winograd also works
     1360        in a case where Sage would use M4RIE matrices by default::
     1361
     1362            sage: MS = MatrixSpace(GF(16,'a'),200,200)
     1363            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     1364            sage: MS._MatrixSpace_generic__matrix_class = Matrix_modpn_dense
     1365            sage: A = MS.random_element()
     1366            sage: B = MS.random_element()
     1367            sage: C = A*B   # Here, school book multiplication is the default.
     1368            sage: type(C)
     1369            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1370            sage: C == A._multiply_strassen(B, cutoff=10)
     1371            True
     1372
     1373        """
     1374        if self._ncols != right._nrows:
     1375            raise ArithmeticError, "Number of columns of self must equal number of rows of right."
     1376        if not self._base_ring is right.base_ring():
     1377            raise TypeError, "Base rings must be the same."
     1378        if cutoff == 0:
     1379            cutoff = self._strassen_default_cutoff(right)
     1380        if cutoff <= 0:
     1381            raise ValueError, "cutoff must be at least 1"
     1382        set_cutoff(cutoff,cutoff)
     1383        from sage.all import MatrixSpace
     1384        MS = MatrixSpace(self.base_ring(), self._nrows, right._ncols)
     1385        result = Matrix_modpn_dense(MS,None)
     1386        if multiply_strassen(self.Data, right.Data, result.Data)!=0:
     1387            return RuntimeError, "Strassen multiplication did not work"
     1388        return result
     1389
     1390    cdef int _strassen_default_cutoff(self, Matrix right) except -2:
     1391        return 110
     1392
     1393    cpdef ModuleElement _lmul_(self, RingElement right):
     1394        """
     1395        Multiply an MeatAxe matrix with a field element.
     1396
     1397        EXAMPLES::
     1398
     1399            sage: F.<z> = GF(9)
     1400            sage: M = MatrixSpace(F,3)(list(F))
     1401            sage: type(M)
     1402            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1403            sage: M
     1404            [      0     2*z   z + 1]
     1405            [  z + 2       2       z]
     1406            [2*z + 2 2*z + 1       1]
     1407            sage: M*z                  # indirect doctest
     1408            [      0 2*z + 2 2*z + 1]
     1409            [      1     2*z   z + 1]
     1410            [  z + 2       2       z]
     1411            sage: 2*M
     1412            [      0       z 2*z + 2]
     1413            [2*z + 1       1     2*z]
     1414            [  z + 1   z + 2       2]
     1415           
     1416        """
     1417        cdef Matrix_modpn_dense left
     1418        if right==0:
     1419            return Matrix_modpn_dense(self._parent, entries=None, mutable=self._mutability._is_immutable) # this is a zero matrix
     1420        elif right==1:
     1421            return self.__copy__()
     1422        if (self.ncols()==0):
     1423            raise NotImplementedError, "The matrix must not be empty"
     1424        left = self.__copy__()
     1425        if matmulF(left.Data,zitof(self._converter.field_to_int(right))):
     1426            return left
     1427        else:
     1428            raise ArithmeticError, "Matrix sizes or fields not compatible"
     1429           
     1430    cpdef Matrix_modpn_dense _mulInt_(self, long right):
     1431        """
     1432        Multiply a MeatAxe matrix with a field element represented by an integer
     1433
     1434        NOTE:
     1435
     1436        This method is of internal use only
     1437
     1438        EXAMPLES::
     1439
     1440            sage: F.<z> = GF(9)
     1441            sage: M = MatrixSpace(F,3)(list(F))
     1442            sage: type(M)
     1443            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1444            sage: M
     1445            [      0     2*z   z + 1]
     1446            [  z + 2       2       z]
     1447            [2*z + 2 2*z + 1       1]
     1448            sage: M._mulInt_(4)
     1449            [      0   z + 2       2]
     1450            [      z 2*z + 2 2*z + 1]
     1451            [      1     2*z   z + 1]
     1452            sage: M*F.fetch_int(4)
     1453            [      0   z + 2       2]
     1454            [      z 2*z + 2 2*z + 1]
     1455            [      1     2*z   z + 1]
     1456
     1457        """
     1458        cdef Matrix_modpn_dense left
     1459        if right==0:
     1460            return Matrix_modpn_dense(self._parent, entries=None, mutable=self._mutability._is_immutable) # this is a zero matrix
     1461        elif right==1:
     1462            return self.__copy__()
     1463        if (self.ncols()==0):
     1464            raise NotImplementedError, "The matrix must not be empty"
     1465        left = self.__copy__()
     1466        if matmulF(left.Data,zitof(right)):
     1467            return left
     1468        else:
     1469            raise ArithmeticError, "Matrix sizes or fields not compatible"
     1470
     1471    def __div__(Matrix_modpn_dense self, long p):
     1472        """
     1473        Divide this matrix by a field element represented by an integer
     1474
     1475        NOTE:
     1476
     1477        If the field is not of prime order then the integer -1 does not
     1478        represent the additive inverse of 1.
     1479
     1480        EXAMPLES::
     1481
     1482            sage: F.<z> = GF(9)
     1483            sage: M = MatrixSpace(F,3)(list(F))
     1484            sage: M
     1485            [      0     2*z   z + 1]
     1486            [  z + 2       2       z]
     1487            [2*z + 2 2*z + 1       1]
     1488            sage: M/3    # indirect doctest
     1489            [      0       2       z]
     1490            [2*z + 2 2*z + 1       1]
     1491            [    2*z   z + 1   z + 2]
     1492            sage: M/-1
     1493            [      0   z + 2       2]
     1494            [      z 2*z + 2 2*z + 1]
     1495            [      1     2*z   z + 1]
     1496
     1497        We verify the results::
     1498
     1499            sage: M/3 == M*~F.fetch_int(3)
     1500            True
     1501            sage: M/-1 == M*~F.fetch_int(8)
     1502            True
     1503
     1504        """
     1505        if ((p%self.Data.fl)==0):
     1506            raise ZeroDivisionError, "Cannot divide by zero"
     1507        elif p==1:
     1508            return self.__copy__()
     1509        elif not self.Data:
     1510            return self
     1511        return self._mulInt_(tmultinv[zitof(p)])
     1512
     1513    def __pow__(Matrix_modpn_dense self,n,m):
     1514        """
     1515        M.__pow__(n): return M^n
     1516
     1517        EXAMPLES::
     1518
     1519            sage: F.<z> = GF(9)
     1520            sage: M = MatrixSpace(F,3)(list(F))
     1521            sage: M[2,2] = 0
     1522            sage: M
     1523            [      0     2*z   z + 1]
     1524            [  z + 2       2       z]
     1525            [2*z + 2 2*z + 1       0]
     1526            sage: M^2
     1527            [      0       0 2*z + 2]
     1528            [      0       2       0]
     1529            [  z + 1       0       0]
     1530            sage: M^4
     1531            [1 0 0]
     1532            [0 1 0]
     1533            [0 0 1]
     1534            sage: M^-3
     1535            [      0     2*z   z + 1]
     1536            [  z + 2       2       z]
     1537            [2*z + 2 2*z + 1       0]
     1538            sage: M^-3*M^3
     1539            [1 0 0]
     1540            [0 1 0]
     1541            [0 0 1]
     1542
     1543        """
     1544        if (self.ncols()==0):
     1545            raise NotImplementedError, "The matrix must not be empty"
     1546        if self._ncols != self._nrows:
     1547            raise TypeError, "Matrix is no square matrix"
     1548        cdef Matrix_modpn_dense OUT
     1549        cdef Matrix_modpn_dense SELFINV
     1550        OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     1551        if n>=0:
     1552            _sig_on
     1553            OUT.Data = matpower(self.Data,n)
     1554            _sig_off
     1555        else:
     1556            SELFINV = self.__invert__()
     1557            _sig_on
     1558            OUT.Data = matpower(SELFINV.Data,-n)
     1559            _sig_off
     1560        if OUT.Data:
     1561            OUT._mutability = copy(self._mutability)
     1562            OUT._parent = self._parent
     1563            OUT._converter = self._converter
     1564            OUT._base_ring = self._base_ring
     1565            OUT._ncols = self._ncols
     1566            OUT._nrows = self._nrows
     1567            OUT._cache = {}
     1568            return OUT
     1569        else:
     1570            raise RuntimeError, "Something went wrong"
     1571
     1572    def __invert__(Matrix_modpn_dense self):
     1573        """
     1574        M__invert__(): return M^(-1), if defined
     1575
     1576        EXAMPLES::
     1577
     1578            sage: F.<z> = GF(9)
     1579            sage: M = MatrixSpace(F,3)(list(F))
     1580            sage: ~M
     1581            Traceback (most recent call last):
     1582            ...
     1583            ArithmeticError: This MeatAxe matrix is not invertible
     1584            sage: M[2,2] = 0
     1585            sage: ~M    # indirect doctest
     1586            [      2       z       0]
     1587            [2*z + 1       1     2*z]
     1588            [      0   z + 2       2]
     1589            sage: M*~M
     1590            [1 0 0]
     1591            [0 1 0]
     1592            [0 0 1]
     1593       
     1594        """
     1595        if (self.ncols()==0):
     1596            raise NotImplementedError, "The matrix must not be empty"
     1597        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     1598        _sig_on
     1599        OUT.Data = matinv(self.Data)
     1600        _sig_off
     1601        if OUT.Data:
     1602            OUT._mutability = copy(self._mutability)
     1603            OUT._parent = self._parent
     1604            OUT._converter = self._converter
     1605            OUT._base_ring = self._base_ring
     1606            OUT._ncols = self._ncols
     1607            OUT._nrows = self._nrows
     1608            OUT._cache = {}
     1609            return OUT
     1610        else:
     1611            raise ArithmeticError, "This MeatAxe matrix is not invertible"
     1612
     1613    cpdef Matrix_modpn_dense transpose(Matrix_modpn_dense self):
     1614        """
     1615        Return the transposed matrix
     1616
     1617        EXAMPLES::
     1618
     1619            sage: M = MatrixSpace(GF(121,'z'),3)([2,3,4,5,6,7,8,9,0])
     1620            sage: type(M)
     1621            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1622            sage: M
     1623            [2 3 4]
     1624            [5 6 7]
     1625            [8 9 0]
     1626            sage: M.transpose()
     1627            [2 5 8]
     1628            [3 6 9]
     1629            [4 7 0]
     1630            sage: transpose(M)   # indirect doctest
     1631            [2 5 8]
     1632            [3 6 9]
     1633            [4 7 0]
     1634
     1635        TESTS::
     1636
     1637            sage: M = MatrixSpace(GF(121,'z'),2,3)(range(6))
     1638            sage: M
     1639            [0 1 2]
     1640            [3 4 5]
     1641            sage: M.parent()
     1642            Full MatrixSpace of 2 by 3 dense matrices over Finite Field in z of size 11^2
     1643            sage: M.transpose()
     1644            [0 3]
     1645            [1 4]
     1646            [2 5]
     1647            sage: M.transpose().parent()
     1648            Full MatrixSpace of 3 by 2 dense matrices over Finite Field in z of size 11^2
     1649
     1650        """
     1651        if (self.ncols()==0):
     1652            raise NotImplementedError, "The matrix must not be empty"
     1653        cdef Matrix_modpn_dense OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     1654        OUT.Data = mattr(self.Data)
     1655        if OUT.Data:
     1656            OUT._mutability = copy(self._mutability)
     1657            from sage.all import MatrixSpace
     1658            OUT._parent = MatrixSpace(self._base_ring, self._ncols, self._nrows)
     1659            OUT._base_ring = self._base_ring
     1660            OUT._converter = self._converter
     1661            OUT._ncols = self._nrows
     1662            OUT._nrows = self._ncols
     1663            OUT._cache = {}
     1664            return OUT
     1665        else:
     1666            raise ArithmeticError, "Error allocating memory"
     1667
     1668    cpdef int order(Matrix_modpn_dense self) except -1:
     1669        """
     1670        Return multiplicative order of self
     1671
     1672        EXAMPLES::
     1673
     1674            sage: F.<z> = GF(9)
     1675            sage: M = MatrixSpace(F,3)([F.fetch_int(i) for i in [2,3,4,5,6,7,8,9,0]])
     1676            sage: type(M)
     1677            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1678            sage: order(M)   # indirect doctest
     1679            728
     1680            sage: M.order()
     1681            728
     1682            sage: M^728
     1683            [1 0 0]
     1684            [0 1 0]
     1685            [0 0 1]
     1686
     1687        We verify the result::
     1688
     1689            sage: E = MatrixSpace(F,3)(1)
     1690            sage: any([M^i==E for i in range(1,728)])
     1691            False
     1692            sage: M^728
     1693            [1 0 0]
     1694            [0 1 0]
     1695            [0 0 1]
     1696
     1697        """
     1698        if (self.ncols()==0):
     1699            raise NotImplementedError, "The matrix must not be empty"
     1700        if (self.Data.nor <> self.Data.noc):
     1701            raise NotImplementedError, "only defined for square matrices"
     1702        cdef int o = matorder(self.Data)
     1703        if o==-1:
     1704            raise ArithmeticError, "order too large"
     1705        else:
     1706            return o
     1707
     1708    def stack(self, *L):
     1709        """
     1710        Stack one matrix over one or several others, provided the number of columns coincides
     1711
     1712        EXAMPLES::
     1713
     1714            sage: MS = MatrixSpace(GF(25,'z'),2,4)
     1715            sage: M = MS([[4,3,2,1],[3,2,1,0]])
     1716            sage: N = MS([[1,1,2,2],[3,3,4,4]])
     1717            sage: type(M)
     1718            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1719            sage: print M.stack(N)
     1720            [4 3 2 1]
     1721            [3 2 1 0]
     1722            [1 1 2 2]
     1723            [3 3 4 4]
     1724            sage: print M.stack(N,M)
     1725            [4 3 2 1]
     1726            [3 2 1 0]
     1727            [1 1 2 2]
     1728            [3 3 4 4]
     1729            [4 3 2 1]
     1730            [3 2 1 0]
     1731
     1732        """
     1733        cdef Matrix_modpn_dense N
     1734        cdef Matrix_modpn_dense OUT
     1735        cdef PTR cur
     1736        cdef int totalrows = self.Data.nor
     1737        for N in L:
     1738            if self.Data.fl!=N.Data.fl:
     1739                raise TypeError, "Matrices must be defined over the same field"
     1740            if self.Data.noc!=N.Data.noc:
     1741                raise RuntimeError, "Matrices must have the same number of columns"
     1742            totalrows += N.Data.nor
     1743        from sage.all import MatrixSpace
     1744        OUT = Matrix_modpn_dense(MatrixSpace(self._base_ring, totalrows, self.Data.noc), entries=None)
     1745        zsetfield(self.Data.fl)
     1746        zsetlen(self.Data.noc)
     1747        cur = OUT.Data.d
     1748        if not copyrows(cur,self.Data.d,self.Data.nor):
     1749            raise MemoryError, "Error copying first matrix"
     1750        zadvance(&(cur),self.Data.nor)
     1751        for i in range(len(L)):
     1752            N = L[i]
     1753            copyrows(cur,N.Data.d,N.Data.nor)
     1754            zadvance(&(cur),N.Data.nor)
     1755        return OUT
     1756
     1757
     1758###################
     1759## Gauss algorithm
     1760
     1761    cpdef Matrix_modpn_dense normalized(Matrix_modpn_dense self):
     1762        """
     1763        Return a normalized copy of self,  i.e., all leading coefficients are 1
     1764
     1765        EXAMPLES::
     1766
     1767            sage: F.<z> = GF(9)
     1768            sage: M = MatrixSpace(F,3)([F.fetch_int(i) for i in [2,3,4,5,6,7,8,9,0]])
     1769            sage: M.normalized()
     1770            [      1     2*z 2*z + 2]
     1771            [      1 2*z + 2       2]
     1772            [      1       0       0]
     1773            sage: M
     1774            [      2       z   z + 1]
     1775            [  z + 2     2*z 2*z + 1]
     1776            [2*z + 2       0       0]
     1777
     1778        """
     1779        if self.Data.fl==2:
     1780            return self.__copy__()
     1781        cdef Matrix_modpn_dense OUT = self.__copy__()
     1782        cdef int i
     1783        cdef PTR p = OUT.Data.d
     1784        cdef piv_t P
     1785        cdef FEL f
     1786        for i from 0<=i<self.Data.nor:
     1787            P = _zfindpiv_(p)
     1788            if P.m:
     1789                if P.m!=1:
     1790                    f = tmultinv[zitof(P.m)]
     1791                    zmulrow(p,f)
     1792            zsteprow(&(p))
     1793        return OUT
     1794
     1795    cpdef Matrix_modpn_dense semi_echelon(Matrix_modpn_dense self):
     1796        """
     1797        Return a normalized semi echelon form of self with zero rows removed
     1798
     1799        EXAMPLES::
     1800
     1801            sage: F.<z> = GF(9)
     1802            sage: M = MatrixSpace(F,3)(list(F))
     1803            sage: M
     1804            [      0     2*z   z + 1]
     1805            [  z + 2       2       z]
     1806            [2*z + 2 2*z + 1       1]
     1807            sage: M.semi_echelon()
     1808            [  0   1 2*z]
     1809            [  1   0   0]
     1810
     1811        TESTS::
     1812
     1813            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     1814            sage: F.<z> = GF(29)
     1815            sage: M = Matrix_modpn_dense(MatrixSpace(F,5), range(1,26))
     1816            sage: M.semi_echelon()
     1817            [1 2 3 4 5]
     1818            [0 1 2 3 4]
     1819            sage: M.echelon_form()
     1820            [ 1  0 28 27 26]
     1821            [ 0  1  2  3  4]
     1822            [ 0  0  0  0  0]
     1823            [ 0  0  0  0  0]
     1824            [ 0  0  0  0  0]
     1825
     1826        """
     1827        if (self.ncols()==0):
     1828            raise NotImplementedError, "The matrix must not be empty"
     1829        cdef Matrix_modpn_dense OUT
     1830        cdef int i
     1831        OUT = Matrix_modpn_dense.__new__(Matrix_modpn_dense)
     1832        _sig_on
     1833        OUT.Data = echelon(self.Data)
     1834        _sig_off
     1835        cdef PTR p = OUT.Data.d
     1836        cdef piv_t P
     1837        cdef FEL f
     1838        if self.Data.fl != 2:
     1839            for i from 0<=i<OUT.Data.nor:
     1840                P = _zfindpiv_(p)
     1841                if P.m:
     1842                    if P.m!=1:
     1843                        f = tmultinv[zitof(P.m)]
     1844                        zmulrow(p,f) # OUT[i] = OUT[i]/f
     1845                zsteprow(&(p))
     1846        OUT._mutability = copy(self._mutability)
     1847        from sage.all import MatrixSpace
     1848        OUT._parent = MatrixSpace(self._base_ring, OUT.Data.nor, OUT.Data.noc)
     1849        OUT._base_ring = self._base_ring
     1850        OUT._converter = self._converter
     1851        OUT._ncols = OUT.Data.noc
     1852        OUT._nrows = OUT.Data.nor
     1853        OUT._cache = {}
     1854        return OUT
     1855
     1856    def _echelon_in_place_classical(self, reduce=True):
     1857        """
     1858        Change self into row echelon form.
     1859
     1860        INPUT:
     1861
     1862        - ``reduce`` (optional, default True): Compute the *reduced* row echelon form.
     1863
     1864        EXAMPLES:
     1865
     1866        Note that in our example we use MeatAxe matrices in a case where
     1867        one would usually have Linbox as a backend::
     1868
     1869            sage: from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     1870            sage: F.<z> = GF(29)
     1871            sage: M = Matrix_modpn_dense(MatrixSpace(F,5), range(1,26))
     1872            sage: M
     1873            [ 1  2  3  4  5]
     1874            [ 6  7  8  9 10]
     1875            [11 12 13 14 15]
     1876            [16 17 18 19 20]
     1877            [21 22 23 24 25]
     1878            sage: M.echelon_form()
     1879            [ 1  0 28 27 26]
     1880            [ 0  1  2  3  4]
     1881            [ 0  0  0  0  0]
     1882            [ 0  0  0  0  0]
     1883            [ 0  0  0  0  0]
     1884            sage: M._echelon_in_place_classical(reduce=False)
     1885            sage: M
     1886            [ 1  2  3  4  5]
     1887            [ 0 24 19 14  9]
     1888            [ 0  0  0  0  0]
     1889            [ 0  0  0  0  0]
     1890            [ 0  0  0  0  0]
     1891
     1892        TESTS:
     1893
     1894        The following gave rise to a segfault::
     1895
     1896            sage: M = matrix(GF(125,'z'), 0,2000)
     1897            sage: M.echelon_form()
     1898            0 x 2000 dense matrix over Finite Field in z of size 5^3
     1899            sage: M = matrix(GF(125,'z'), 2000, 0)
     1900            sage: M.echelon_form()
     1901            2000 x 0 dense matrix over Finite Field in z of size 5^3
     1902
     1903        """
     1904        if (self.Data.nor == 0) or (self.Data.noc == 0):
     1905            self.cache('pivots', ())
     1906            return
     1907        zsetfield(self.Data.fl)
     1908        zsetlen(self.Data.noc)
     1909        cdef int nr = self.Data.nor
     1910        cdef long* piv = <long*>calloc(nr+1,sizeof(long))
     1911        zmkechelon(self.Data.d, nr, piv)
     1912        cdef int i
     1913        cdef int rk = piv[0]
     1914        self.cache('pivots', tuple(piv[i]-1 for i in xrange(1,rk+1)))
     1915        if rk == 0:
     1916            free(piv)
     1917            return
     1918        # Need to empty the remaining rows - zmkechelon was moving them
     1919        # up, but did not clear them.
     1920        cdef PTR p = self.Data.d
     1921        zadvance(&(p), rk)
     1922        for i from rk<=i<nr:
     1923            memset(p, 0, zrowsize)
     1924            zsteprow(&(p))
     1925        # Sage wants that the reduced row echelon form is created.
     1926        # Hence, we need to normalize and to clear the parts on top
     1927        # of the pivots.
     1928        cdef PTR q
     1929        cdef FEL f
     1930        cdef int j
     1931        if reduce:
     1932            p = self.Data.d
     1933            f = tmultinv[zitof(zextract(p,piv[1]))]
     1934            zmulrow(p,f)
     1935            q = self.Data.d
     1936            for i from 2<=i<=rk:
     1937                zsteprow(&(q))
     1938                # normalize the row
     1939                f = tmultinv[zitof(zextract(q,piv[i]))]
     1940                zmulrow(q,f)
     1941                # clear all rows on top of it
     1942                p = self.Data.d
     1943                for j from 1<=j<i:
     1944                    f = taddinv[zitof(zextract(p,piv[i]))]
     1945                    zaddmulrow(p, q, f)
     1946                    zsteprow(&(p))
     1947        free(piv)
     1948
     1949    cpdef int nullity(Matrix_modpn_dense self):
     1950        """
     1951        Return the nullity of a matrix
     1952
     1953        EXAMPLES::
     1954
     1955            sage: M = MatrixSpace(GF(25,'z'),2,4)([[4,3,2,1],[3,2,1,0]])
     1956            sage: M.nullity()
     1957            2
     1958        """
     1959        if (self.ncols()==0):
     1960            raise NotImplementedError, "The matrix must not be empty"
     1961        return nullity(self.Data)
     1962
     1963    cpdef tuple lead(self):
     1964        """
     1965        ``(f,i) = M.lead()`` `\\iff` ``f=M[0,i]`` is the first non-zero coefficient in the first row of M.
     1966       
     1967        If the first row of M has no non-zero entry then ``f==0`` and ``i`` is the number of columns
     1968
     1969        EXAMPLES::
     1970
     1971            sage: M = MatrixSpace(GF(25,'z'),1,4)([0,2,4,1])
     1972            sage: M.lead()
     1973            (2, 1)
     1974            sage: N = MatrixSpace(GF(25,'z'),3,15)(0)
     1975            sage: N.lead()
     1976            (0, 15)
     1977           
     1978        """
     1979        if (self.ncols()==0):
     1980            raise NotImplementedError, "The matrix must not be empty"
     1981        zsetfield(self.Data.fl)
     1982        zsetlen(self.Data.noc)
     1983        cdef piv_t P = _zfindpiv_(self.Data.d)
     1984        return (P.m, P.i)
     1985
     1986    # set/get "fast and dirty"
     1987
     1988    cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
     1989        """
     1990        Get an element without checking.
     1991
     1992        TEST::
     1993
     1994            sage: F.<z> = GF(9)
     1995            sage: M = MatrixSpace(F,3)(list(F))
     1996            sage: type(M)
     1997            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     1998            sage: M    # indirect doctest
     1999            [      0     2*z   z + 1]
     2000            [  z + 2       2       z]
     2001            [2*z + 2 2*z + 1       1]
     2002
     2003        """
     2004        cdef PTR p = self.Data.d
     2005        if self.Data:
     2006            zsetfield(self.Data.fl)
     2007            zsetlen(self.Data.noc)
     2008        else:
     2009            raise IndexError, "Matrix is empty"
     2010        zadvance(&(p),i)
     2011        return self._converter.int_to_field(zftoi(zextract(p,j+1)))
     2012
     2013    cdef int get_unsafe_int(self, Py_ssize_t i, Py_ssize_t j):
     2014        # NOTE:
     2015        # It is essential that you call zsetfield and zsetlen YOURSELF
     2016        # and that you assert that the matrix is not empty!
     2017        # This method is here for speed!
     2018        cdef PTR p = self.Data.d
     2019        if self.Data:
     2020            zsetfield(self.Data.fl)
     2021            zsetlen(self.Data.noc)
     2022        else:
     2023            raise IndexError, "Matrix is empty"
     2024        zadvance(&(p),i)
     2025        return zftoi(zextract(p,j+1))
     2026
     2027    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
     2028        # ASSUMPTION: value's parent is the base ring
     2029        if self.Data:
     2030            zsetfield(self.Data.fl)
     2031            zsetlen(self.Data.noc)
     2032        else:
     2033            raise IndexError, "Matrix is empty"
     2034        cdef PTR x = self.Data.d
     2035        zadvance(&(x),i)
     2036        zinsert(x, j+1, zitof(self._converter.field_to_int(value)))
     2037
     2038    cdef set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value):
     2039        # NOTE:
     2040        # It is essential that you call zsetfield and zsetlen YOURSELF
     2041        # and that you assert that the matrix is not empty!
     2042        # This method is here for speed!
     2043        cdef PTR x = self.Data.d
     2044        zadvance(&(x),i)
     2045        zinsert(x, j+1, zitof(value))
     2046
     2047    def randomize(self, density=None, nonzero=False, *args, **kwds):
     2048        """
     2049        Fill the matrix with random values.
     2050
     2051        INPUT:
     2052
     2053        - ``density`` (optional real number between zero and one) --
     2054          the expected density of the resulting matrix
     2055        - ``nonzero`` (optional bool, default ``False``) --
     2056          If true, all inserted marks are non-zero.
     2057
     2058        EXAMPLE::
     2059
     2060            sage: MS = MatrixSpace(GF(27,'z'),6,6)
     2061            sage: M = MS.random_element(); M    # indirect doctest
     2062            [              1           z + 1     z^2 + z + 1             z^2       2*z^2 + z           z + 1]
     2063            [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]
     2064            [        2*z + 2     z^2 + z + 2           z + 2 2*z^2 + 2*z + 2           2*z^2           2*z^2]
     2065            [  2*z^2 + z + 2             z^2           z + 2         z^2 + z       2*z^2 + 2         z^2 + 2]
     2066            [      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]
     2067            [        2*z + 1         z^2 + z             z^2             z^2     2*z^2 + 2*z           z + 1]
     2068            sage: type(M)
     2069            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     2070            sage: MS.random_element(nonzero=True)
     2071            [            2*z               1   z^2 + 2*z + 1   2*z^2 + z + 1             z^2     z^2 + z + 1]
     2072            [    2*z^2 + 2*z   2*z^2 + z + 2         2*z + 1       z^2 + 2*z     2*z^2 + 2*z             z^2]
     2073            [        z^2 + z     z^2 + z + 2 2*z^2 + 2*z + 1         z^2 + 2               1           2*z^2]
     2074            [              z     2*z^2 + 2*z           2*z^2         2*z + 1           z + 2           z + 2]
     2075            [        z^2 + z             z^2           z + 2     2*z^2 + 2*z         2*z + 1         z^2 + z]
     2076            [    z^2 + z + 2       2*z^2 + z             z^2           z + 1     2*z^2 + 2*z   z^2 + 2*z + 1]
     2077            sage: MS.random_element(density=0.5)
     2078            [        z^2 + 2               0   z^2 + 2*z + 2       2*z^2 + z               0     z^2 + z + 2]
     2079            [              0               1               0               0               0               0]
     2080            [  2*z^2 + z + 1   2*z^2 + z + 2               0     z^2 + z + 2               0     z^2 + z + 1]
     2081            [              0               0               0               0               0               0]
     2082            [2*z^2 + 2*z + 2               0               0   2*z^2 + z + 2               0         2*z + 1]
     2083            [              0       2*z^2 + z               0               1               0   2*z^2 + z + 1]
     2084
     2085        """
     2086        self.check_mutability()
     2087        cdef int fl = self.Data.fl
     2088        density = float(density)
     2089        if density <= 0:
     2090            return
     2091        if density > 1:
     2092            density = float(1)
     2093
     2094        self.clear_cache()
     2095
     2096        cdef unsigned char *x,*y
     2097        x = <unsigned char*>(self.Data.d)
     2098        cdef int nr = self.Data.nor
     2099        cdef int nc = self.Data.noc
     2100        cdef int i, j, k
     2101
     2102        zsetfield(fl)
     2103        zsetlen(nc)
     2104        cdef int O
     2105        randint = current_randstate().c_random
     2106        randdouble = current_randstate().c_rand_double
     2107
     2108        if not nonzero:
     2109            if density == 1:
     2110                O = (self.base_ring().order()**MPB)
     2111                sig_on()
     2112                for i from 0 <= i < nr:
     2113                    y = x
     2114                    for j from 0 <= j < zrowsize_io:
     2115                        y[j] = randint()%O
     2116                    zsteprow(&(x))
     2117                sig_off()
     2118            else:
     2119                sig_on()
     2120                for i from 0 <= i < nr:
     2121                    y = x
     2122                    for j from 0 < j <= nc:
     2123                        if randdouble() < density:
     2124                            zinsert(y, j, zitof( (randint()%fl) ))
     2125                    zsteprow(&(x))
     2126                sig_off()
     2127        else:
     2128            if density == 1:
     2129                fl -= 1
     2130                sig_on()
     2131                for i from 0 <= i < nr:
     2132                    y = x
     2133                    for j from 0 < j <= nc:
     2134                        zinsert(y, j, zitof( (randint()%fl)+1 ))
     2135                    zsteprow(&(x))
     2136                sig_off()
     2137            else:
     2138                fl -= 1
     2139                sig_on()
     2140                for i from 0 <= i < nr:
     2141                    y = x
     2142                    for j from 0 < j <= nc:
     2143                        if randdouble() < density:
     2144                            zinsert(y, j, zitof( (randint()%fl)+1 ))
     2145                    zsteprow(&(x))
     2146                sig_off()
  • sage/matrix/matrix_space.py

    diff --git a/sage/matrix/matrix_space.py b/sage/matrix/matrix_space.py
    a b  
    952952            <type 'sage.matrix.matrix_modn_dense_float.Matrix_modn_dense_float'>
    953953            sage: type(matrix(GF(16007), 2, range(4)))
    954954            <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
     955            sage: type(matrix(GF(2), 2, range(4)))
     956            <type 'sage.matrix.matrix_mod2_dense.Matrix_mod2_dense'>
     957            sage: type(matrix(GF(64,'z'), 2, range(4)))
     958            <type 'sage.matrix.matrix_mod2e_dense.Matrix_mod2e_dense'>
     959            sage: type(matrix(GF(125,'z'), 2, range(4)))
     960            <type 'sage.matrix.matrix_modpn_dense.Matrix_modpn_dense'>
     961
    955962        """
    956963        R = self.base_ring()
    957964        if self.is_dense():
     
    981988                return matrix_generic_dense.Matrix_generic_dense
    982989            elif sage.rings.finite_rings.all.is_FiniteField(R) and R.characteristic() == 2 and R.order() <= 1024:
    983990                return matrix_mod2e_dense.Matrix_mod2e_dense
     991            elif sage.rings.finite_rings.all.is_FiniteField(R) and R.characteristic() != 2 and R.order() <= 255:
     992                from sage.matrix.matrix_modpn_dense import Matrix_modpn_dense
     993                return Matrix_modpn_dense
    984994            elif sage.rings.polynomial.multi_polynomial_ring_generic.is_MPolynomialRing(R) and R.base_ring() in _Fields:
    985995                return matrix_mpolynomial_dense.Matrix_mpolynomial_dense
    986996            #elif isinstance(R, sage.rings.padics.padic_ring_capped_relative.pAdicRingCappedRelative):
  • sage/modules/quotient_module.py

    diff --git a/sage/modules/quotient_module.py b/sage/modules/quotient_module.py
    a b  
    2525
    2626    To obtain V or W use \code{self.V()} and \code{self.W()}.
    2727
    28     EXAMPLES:
     28    EXAMPLES::
     29
    2930        sage: k.<i> = QuadraticField(-1)
    3031        sage: A = k^3; V = A.span([[1,0,i], [2,i,0]])
    3132        sage: W = A.span([[3,i,i]])
     
    7273        """
    7374        Create this quotient space, from the given domain, sub-module, and quotient_matrix.
    7475
    75         EXAMPLES:
     76        EXAMPLES::
     77
    7678            sage: A = QQ^5; V = A.span_of_basis([[1,0,-1,1,1], [1,-1,0,2/3,3/4]]); V
    7779            Vector space of degree 5 and dimension 2 over Rational Field
    7880            User basis matrix:
     
    8385            User basis matrix:
    8486            [1/3 2/3  -1 5/9 1/2]
    8587
    86         This creates a quotient vector space, which calls the init method:
     88        This creates a quotient vector space, which calls the init method::
     89
    8790            sage: Q = V / W
    8891
    89         Behold the type of Q:
     92        Behold the type of Q::
     93
    9094            sage: type(Q)
    9195            <class 'sage.modules.quotient_module.FreeModule_ambient_field_quotient_with_category'>
    9296           
    9397        We do some consistency checks on the extra quotient and
    9498        lifting structure of Q.
     99        ::
     100
    95101            sage: Q(V.0)
    96102            (1)
    97103            sage: Q( V.0 - 2/3*V.1 )
     
    116122        Return the rather verbose string representation of this quotient space V/W.
    117123
    118124        EXAMPLES:
    119         We create a quotient vector space over a finite field:
     125
     126        We create a quotient vector space over a finite field::
     127
    120128            sage: k.<a> = GF(9); A = k^3; V = A.span_of_basis([[1,0,a], [a,a,1]]); W = V.span([V.1])
    121129            sage: Q = V/W
    122130
    123         Note the type:
     131        Note the type::
     132
    124133            sage: type(Q)
    125134            <class 'sage.modules.quotient_module.FreeModule_ambient_field_quotient_with_category'>
    126135
    127136        The string representation mentions that this is a quotient V/W, that the quotient
    128         has dimension 1 and is over a finite field, and also describes V and W:
     137        has dimension 1 and is over a finite field, and also describes V and W::
     138
    129139            sage: Q._repr_()
    130140            '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]'
    131141        """
     
    140150        the tuple (V,W).
    141151
    142152        EXAMPLES:
    143         We compute the hash of a certain 0-dimension quotient vector space:
     153
     154        We compute the hash of a certain 0-dimension quotient vector space::
     155
    144156            sage: A = QQ^2; V = A.span_of_basis([[1,0], [1,1]]); W = V.span([V.1, V.0])
    145157            sage: Q = V/W; Q.dimension()
    146158            0
     
    149161            -5856620741060301410    # 64-bit
    150162
    151163        The hash is just got by hashing both V and W.
     164        ::
     165
    152166            sage: hash((V, W))
    153167            954887582             # 32-bit
    154168            -5856620741060301410  # 64-bit
     
    165179        self and other.
    166180
    167181        EXAMPLES:
    168         We create three quotient spaces and compare them:
     182
     183        We create three quotient spaces and compare them::
     184
    169185            sage: A = QQ^2; V = A.span_of_basis([[1,0], [1,1]]);
    170186            sage: W0 = V.span([V.1, V.0]); W1 = V.span([V.1]); W2 = V.span([V.1])
    171187            sage: Q0 = V/W0; Q1 = V/W1; Q2 = V/W2
     
    191207        if it can be made sense of as a list of length the dimension of self.
    192208               
    193209        EXAMPLES:
     210
    194211        We create a 2-dimensional quotient of a 3-dimension ambient vector space.
     212        ::
     213
    195214            sage: M = QQ^3 / [[1,2,3]]
    196215
    197216        A list of length 3 coerces into QQ^3, so it coerces into M.
     217        ::
     218
    198219            sage: M([1,2,4])
    199220            (-1/3, -2/3)
    200221            sage: M([1,2,3])
     
    202223
    203224        A list of length 2 at least coerces into M, where here it just gives
    204225        the corresponding linear combination of the basis for M.
     226        ::
     227           
    205228            sage: M([1,2])
    206229            (1, 2)
    207230            sage: M.0 + 2*M.1
     
    225248        Elements canonically coerce into self if they canonically
    226249        coerce into V.
    227250       
    228         EXAMPLES:
     251        EXAMPLES::
     252
    229253            sage: V = QQ^3; W = V.span([[1,0,0]]); Q = V/W
    230254            sage: Q._coerce_impl(V.0)
    231255            (0, 0)
     
    239263            (2, 0)
    240264
    241265        Here we coerce in something that is over ZZ, so it canonically coerce into V hence self.
     266        ::
     267
    242268            sage: Q._coerce_impl((ZZ^3)([1,2,3]))
    243269            (2, 3)
    244270
    245         Here there is no canonical coercion:
     271        Here there is no canonical coercion::
     272
    246273            sage: Q._coerce_impl((GF(17)^3)([1,2,3]))
    247274            Traceback (most recent call last):
    248275            ...
     
    259286        """
    260287        Given this quotient space $Q = V/W$, return the natural quotient map from V to Q.
    261288
    262         EXAMPLES:
     289        EXAMPLES::
     290
    263291            sage: M = QQ^3 / [[1,2,3]]
    264292            sage: M.quotient_map()
    265293            Free module morphism defined by the matrix
     
    279307        Given this quotient space $Q = V/W$, return a fixed choice of linear homomorphism
    280308        (a section) from Q to V.
    281309
    282         EXAMPLES:
     310        EXAMPLES::
     311
    283312            sage: M = QQ^3 / [[1,2,3]]
    284313            sage: M.lift_map()
    285314            Free module morphism defined by the matrix
     
    296325
    297326        The lift is a fixed homomorphism. 
    298327
    299         EXAMPLES:
     328        EXAMPLES::
     329
    300330            sage: M = QQ^3 / [[1,2,3]]
    301331            sage: M.lift(M.0)
    302332            (1, 0, 0)
     
    311341        """
    312342        Given this quotient space $Q = V/W$, return W.
    313343
    314         EXAMPLES:
     344        EXAMPLES::
     345
    315346            sage: M = QQ^10 / [range(10), range(2,12)]
    316347            sage: M.W()
    317348            Vector space of degree 10 and dimension 2 over Rational Field
     
    325356        """
    326357        Given this quotient space $Q = V/W$, return $V$.
    327358
    328         EXAMPLES:
     359        EXAMPLES::
     360
    329361            sage: M = QQ^10 / [range(10), range(2,12)]
    330362            sage: M.V()
    331363            Vector space of dimension 10 over Rational Field
     
    336368        """
    337369        Given this quotient space $Q = V/W$, return $V$.  This is the same as self.V().
    338370
    339         EXAMPLES:
     371        EXAMPLES::
     372
    340373            sage: M = QQ^10 / [range(10), range(2,12)]
    341374            sage: M.cover()
    342375            Vector space of dimension 10 over Rational Field
     
    347380        """
    348381        Given this quotient space $Q = V/W$, return $W$.  This is the same as self.W().
    349382
    350         EXAMPLES:
     383        EXAMPLES::
     384
    351385            sage: M = QQ^10 / [range(10), range(2,12)]
    352386            sage: M.relations()
    353387            Vector space of degree 10 and dimension 2 over Rational Field