Ticket #9562: m4rie_for_sage.patch

File m4rie_for_sage.patch, 94.3 KB (added by malb, 8 years ago)

rebased to 4.7.3.alpha0

  • module_list.py 

    # HG changeset patch
# User Martin Albrecht <martinralbrecht@googlemail.com>
# Date 1279820685 -3600
# Node ID 533b016f33a66cff78c4a7e6025fe8ba65b06a93
# Parent  8040f2122c9bb4a68c7a75c0f4b53e2d51d46633
#9562 Matrices over GF(2^n) for 2 <= n <= 10 using the M4RIE library

diff --git a/module_list.py b/module_list.py
    a b  
    917917              extra_compile_args = ['-std=c99'] + m4ri_extra_compile_args,
    918918              depends = [SAGE_INC + "png.h", SAGE_INC + "m4ri/m4ri.h"]),
    919919
     920    Extension('sage.matrix.matrix_mod2e_dense',
     921              sources = ['sage/matrix/matrix_mod2e_dense.pyx'],
     922              libraries = ['m4rie', 'm4ri', 'givaro', 'ntl', 'gmpxx', 'gmp', 'm', 'stdc++'],
     923              depends = [SAGE_INC + "m4rie/m4rie.h"],
     924              include_dirs = [SAGE_INC + 'm4rie'],
     925              extra_compile_args = m4ri_extra_compile_args,
     926              language="c++"),
     927
    920928    Extension('sage.matrix.matrix_modn_dense',
    921929              sources = ['sage/matrix/matrix_modn_dense.pyx'],
    922930              libraries = ['gmp']),

  • sage/crypto/block_cipher/miniaes.py 

    diff --git a/sage/crypto/block_cipher/miniaes.py b/sage/crypto/block_cipher/miniaes.py
    a b  
    2626# http://www.gnu.org/licenses/
    2727###########################################################################
    2828
    29 from sage.matrix.matrix_generic_dense import Matrix_generic_dense
     29from sage.matrix.matrix_dense import Matrix_dense
    3030from sage.matrix.matrix_space import MatrixSpace, MatrixSpace_generic
    3131from sage.monoids.string_monoid import BinaryStrings
    3232from sage.monoids.string_monoid_element import StringMonoidElement
    … …  
    536536            ...
    537537            TypeError: round key must be a 2 x 2 matrix over GF(16)
    538538        """
    539         if not isinstance(block, (MatrixSpace_generic, Matrix_generic_dense)):
     539        if not isinstance(block, Matrix_dense) or \
     540                not (block.base_ring().order() == 16 and block.base_ring().is_field()):
    540541            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    541542        if not (block.nrows() == block.ncols() == 2):
    542543            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    543         if not isinstance(rkey, (MatrixSpace_generic, Matrix_generic_dense)):
     544
     545        if not isinstance(rkey, Matrix_dense) or \
     546                not (rkey.base_ring().order() == 16 and rkey.base_ring().is_field()):
    544547            raise TypeError, "round key must be a 2 x 2 matrix over GF(16)"
    545548        if not (rkey.nrows() == rkey.ncols() == 2):
    546549            raise TypeError, "round key must be a 2 x 2 matrix over GF(16)"
    … …  
    705708            ...
    706709            TypeError: secret key must be a 2 x 2 matrix over GF(16)
    707710        """
    708         if not isinstance(C, (MatrixSpace_generic, Matrix_generic_dense)):
     711        if not isinstance(C, Matrix_dense) or \
     712                not (C.base_ring().order() == 16 and C.base_ring().is_field()):
    709713            raise TypeError, "ciphertext block must be a 2 x 2 matrix over GF(16)"
    710714        if not (C.nrows() == C.ncols() == 2):
    711715            raise TypeError, "ciphertext block must be a 2 x 2 matrix over GF(16)"
    712         if not isinstance(key, (MatrixSpace_generic, Matrix_generic_dense)):
     716        if not isinstance(key, Matrix_dense) or \
     717                not (key.base_ring().order() == 16 and key.base_ring().is_field()):
    713718            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    714719        if not (key.nrows() == key.ncols() == 2):
    715720            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    … …  
    864869            ...
    865870            TypeError: secret key must be a 2 x 2 matrix over GF(16)
    866871        """
    867         if not isinstance(P, (MatrixSpace_generic, Matrix_generic_dense)):
     872        if not isinstance(P, Matrix_dense) or \
     873                not (P.base_ring().order() == 16 and P.base_ring().is_field()):
    868874            raise TypeError, "plaintext block must be a 2 x 2 matrix over GF(16)"
    869875        if not (P.nrows() == P.ncols() == 2):
    870876            raise TypeError, "plaintext block must be a 2 x 2 matrix over GF(16)"
    871         if not isinstance(key, (MatrixSpace_generic, Matrix_generic_dense)):
     877        if not isinstance(key, Matrix_dense) or \
     878                not (key.base_ring().order() == 16 and key.base_ring().is_field()):
    872879            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    873880        if not (key.nrows() == key.ncols() == 2):
    874881            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    … …  
    10161023            ...
    10171024            TypeError: input block must be a 2 x 2 matrix over GF(16)
    10181025        """
    1019         if not isinstance(block, (MatrixSpace_generic, Matrix_generic_dense)):
     1026        if not isinstance(block, Matrix_dense) or \
     1027                not (block.base_ring().order() == 16 and block.base_ring().is_field()):
    10201028            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    10211029        if not (block.nrows() == block.ncols() == 2):
    10221030            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    … …  
    12061214            ...
    12071215            ValueError: the algorithm for nibble-sub must be either 'encrypt' or 'decrypt'
    12081216        """
    1209         if not isinstance(block, (MatrixSpace_generic, Matrix_generic_dense)):
     1217        if not isinstance(block, Matrix_dense) or \
     1218                not (block.base_ring().order() == 16 and block.base_ring().is_field()):
    12101219            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    12111220        if not (block.nrows() == block.ncols() == 2):
    12121221            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    … …  
    13491358            ...
    13501359            TypeError: secret key must be a 2 x 2 matrix over GF(16)
    13511360        """
    1352         if not isinstance(key, (MatrixSpace_generic, Matrix_generic_dense)):
     1361        if not isinstance(key, Matrix_dense) or \
     1362                not (key.base_ring().order() == 16 and key.base_ring().is_field()):
    13531363            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    13541364        if not (key.nrows() == key.ncols() == 2):
    13551365            raise TypeError, "secret key must be a 2 x 2 matrix over GF(16)"
    … …  
    14981508            ...
    14991509            TypeError: input block must be a 2 x 2 matrix over GF(16)
    15001510        """
    1501         if not isinstance(block, (MatrixSpace_generic, Matrix_generic_dense)):
     1511        if not isinstance(block, Matrix_dense) or \
     1512                not (block.base_ring().order() == 16 and block.base_ring().is_field()):
    15021513            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    15031514        if not (block.nrows() == block.ncols() == 2):
    15041515            raise TypeError, "input block must be a 2 x 2 matrix over GF(16)"
    … …  
    16521663            S = "".join([str(self._GF_to_bin[g]) for g in G])
    16531664            return B(S)
    16541665        # G is a matrix over GF(16)
    1655         elif isinstance(G, (MatrixSpace_generic, Matrix_generic_dense)):
     1666        elif isinstance(G, Matrix_dense):
    16561667            if not (G.base_ring() is K):
    16571668                raise TypeError, "input G must be an element of GF(16), a list of elements of GF(16), or a matrix over GF(16)"
    16581669            S = "".join([str(self._GF_to_bin[G[i][j]]) for i in xrange(G.nrows()) for j in xrange(G.ncols())])
    … …  
    17831794                raise ValueError, "input G must be an element of GF(16), a list of elements of GF(16), or a matrix over GF(16)"
    17841795            return [self._GF_to_int[g] for g in G]
    17851796        # G is a matrix over GF(16)
    1786         elif isinstance(G, (MatrixSpace_generic, Matrix_generic_dense)):
     1797        elif isinstance(G, Matrix_dense):
    17871798            if not (G.base_ring() is K):
    17881799                raise TypeError, "input G must be an element of GF(16), a list of elements of GF(16), or a matrix over GF(16)"
    17891800            return [self._GF_to_int[G[i][j]] for i in xrange(G.nrows()) for j in xrange(G.ncols())]

  • sage/crypto/mq/mpolynomialsystem.py 

    diff --git a/sage/crypto/mq/mpolynomialsystem.py b/sage/crypto/mq/mpolynomialsystem.py
    a b  
    3030We can construct a polynomial system for a random plaintext-ciphertext
    3131pair and study it::
    3232
     33   sage: set_random_seed(1)
    3334   sage: F,s = sr.polynomial_system()
    3435   sage: F
    3536   Polynomial Sequence with 112 Polynomials in 64 Variables
    … …  
    102103    sage: A.rank()
    103104    4056
    104105    sage: A[4055]*v
    105     (k002*k003)
    106 
    107 The last output tells us that ``k002`` and ``k003`` can't both be 1.
     106    (k001*k002 + k001*k003)
    108107
    109108TEST::
    110109
    … …  
    296295            sage: F,s = sr.polynomial_system()
    297296            sage: r = F.part(0)
    298297            sage: copy(r)
    299             (w100 + k000 + (a^3 + a + 1), w101 + k001 + (a^3 + 1),
    300             w102 + k002 + (a^3 + a^2 + 1), w103 + k003 + (a^3 + a^2 +
    301             a))
     298            (w100 + k000 + (a^2 + 1), w101 + k001 + (a), w102 + k002 + (a^2), w103 + k003 + (a + 1))
    302299            sage: type(r) == type(copy(r))
    303300            True
    304301        """
    … …  
    826823            sage: P = F.ring()
    827824            sage: I = F.ideal()
    828825            sage: I.elimination_ideal(P('s000*s001*s002*s003*w100*w101*w102*w103*x100*x101*x102*x103'))
    829             Ideal (k002 + (a^2)*k003 + 1,
    830                    k001 + (a^3)*k003 + (a + 1),
    831                    k000 + (a^3 + a^2 + a)*k003 + (a^3 + a^2 + a),
    832                    k103 + (a^3)*k003 + (a^2 + 1),
    833                    k102 + (a^3 + a^2 + a)*k003 + (a^3 + 1),
    834                    k101 + k003 + (a^2 + a),
    835                    k100 + (a^2)*k003 + (a^2 + a),
    836                    k003^2 + (a^3 + a^2 + a)*k003 + (a^3 + a^2 + a))
    837             of Multivariate Polynomial Ring in k100, k101, k102, k103,
    838             x100, x101, x102, x103, w100, w101, w102, w103, s000,
    839             s001, s002, s003, k000, k001, k002, k003 over Finite Field
    840             in a of size 2^4
     826            Ideal (k002 + (a^3 + a^2 + 1)*k003 + (a^3),
     827                   k001 + (a^3 + a^2 + a + 1)*k003 + (a^3 + a^2 + a),
     828                   k000 + (a + 1)*k003 + (a^2 + a + 1),
     829                   k103 + (a^3 + a)*k003 + (a^3 + a),
     830                   k102 + (a^2 + a + 1)*k003 + (a^3 + a^2 + a),
     831                   k101 + (a^3)*k003 + (a^3),
     832                   k100 + (a^3 + a + 1)*k003 + (a^2 + 1),
     833                   k003^2 + (a + 1)*k003 + (a^2 + a + 1))
     834            of Multivariate Polynomial Ring in
     835            k100, k101, k102, k103, x100, x101, x102, x103, w100, w101, w102, w103,
     836            s000, s001, s002, s003, k000, k001, k002, k003 over Finite Field in a of size 2^4
    841837        """
    842838        return self._ring.ideal(self.gens())
    843839

  • sage/crypto/mq/sr.py 

    diff --git a/sage/crypto/mq/sr.py b/sage/crypto/mq/sr.py
    a b  
    10581058       
    10591059            sage: sr = mq.SR(2, 2, 2, 4)
    10601060            sage: sr.random_state_array()
    1061             [a^3 + a + 1       a + 1]
    1062             [      a + 1         a^2]
     1061            [      a^2 + 1 a^3 + a^2 + a]
     1062            [  a^3 + a + 1         a + 1]
    10631063        """
    10641064        return random_matrix(self.base_ring(), self._r, self._c, *args, **kwds)
    10651065
    … …  
    10721072       
    10731073            sage: sr = mq.SR(2, 2, 2, 4)
    10741074            sage: sr.random_vector()
     1075            [      a^2 + 1]
     1076            [            a]
     1077            [          a^2]
     1078            [        a + 1]
    10751079            [  a^3 + a + 1]
    10761080            [      a^3 + 1]
    10771081            [a^3 + a^2 + 1]
    10781082            [a^3 + a^2 + a]
     1083            [a^3 + a^2 + a]
     1084            [  a^3 + a + 1]
     1085            [      a^3 + 1]
     1086            [a^3 + a^2 + 1]
    10791087            [        a + 1]
    10801088            [      a^2 + 1]
    10811089            [            a]
    10821090            [          a^2]
    1083             [        a + 1]
    1084             [      a^2 + 1]
    1085             [            a]
    1086             [          a^2]
    1087             [          a^2]
    1088             [        a + 1]
    1089             [      a^2 + 1]
    1090             [            a]
    10911091       
    10921092        .. note::
    10931093
    … …  
    11111111       
    11121112            sage: sr = mq.SR()
    11131113            sage: sr.random_element()
    1114             [  a^3 + a + 1]
    1115             [      a^3 + 1]
    1116             [a^3 + a^2 + 1]
    1117             [a^3 + a^2 + a]
     1114            [a^2 + 1]
     1115            [      a]
     1116            [    a^2]
     1117            [  a + 1]
    11181118            sage: sr.random_element('state_array')
    1119             [a + 1]
     1119            [a^3 + 1]
    11201120
    11211121        Passes extra positional or keyword arguments through::
    11221122
    … …  
    20352035            sage: P = sr.vars("P",0)
    20362036            sage: F,s = sr.polynomial_system(P=P,C=C)
    20372037            sage: [(k,v) for k,v in sorted(s.iteritems())] # this can be ignored
    2038             [(k003, 1), (k002, 1), (k001, 0), (k000, 0)]
     2038            [(k003, 1), (k002, 0), (k001, 0), (k000, 1)]
    20392039            sage: F
    20402040            Polynomial Sequence with 36 Polynomials in 28 Variables
    20412041            sage: F.part(0)
    … …  
    22222222            sage: sr = mq.SR()
    22232223            sage: A = sr.random_state_array()
    22242224            sage: A
    2225             [a^3 + a + 1]
     2225            [a^2 + 1]
    22262226            sage: sr.antiphi(sr.phi(A)) == A
    22272227            True
    22282228        """
    … …  
    26272627            sage: sr = mq.SR(gf2=True)
    26282628            sage: A = sr.random_state_array()
    26292629            sage: A
    2630             [a^3 + a + 1]
     2630            [a^2 + 1]
    26312631            sage: sr.antiphi(sr.phi(A)) == A
    26322632            True
    26332633        """

  • sage/libs/m4ri.pxd 

    diff --git a/sage/libs/m4ri.pxd b/sage/libs/m4ri.pxd
    a b  
     1
    12cdef extern from "m4ri/m4ri.h":
    23    ctypedef int rci_t
    34    ctypedef int wi_t

  • new file sage/libs/m4rie.pxd 

    diff --git a/sage/libs/m4rie.pxd b/sage/libs/m4rie.pxd
new file mode 100644
    - +  
     1##############################################################################
     2#       Copyright (C) 2010 Martin Albrecht <martinralbrecht@googlemail.com>
     3#  Distributed under the terms of the GNU General Public License (GPL)
     4#  The full text of the GPL is available at:
     5#                  http://www.gnu.org/licenses/
     6##############################################################################
     7
     8from sage.libs.m4ri cimport mzd_t, m4ri_word
     9
     10
     11
     12cdef extern from "m4rie/m4rie.h":
     13    ctypedef struct gf2e:
     14        m4ri_word **mul
     15        m4ri_word *inv
     16        size_t degree
     17        m4ri_word minpoly
     18
     19    void gf2e_free(gf2e *ff)
     20
     21#cdef extern from "m4rie/gf2e_matrix.h":
     22    ctypedef struct mzed_t:
     23        mzd_t *x
     24        gf2e *finite_field
     25        int nrows
     26        int ncols
     27        int w
     28
     29    ctypedef int const_int "const int"
     30    ctypedef size_t const_size_t "const size_t"
     31    ctypedef mzed_t const_mzed_t "const mzed_t"
     32
     33    mzed_t *mzed_init(gf2e *, size_t m, size_t n)
     34
     35    void mzed_free(mzed_t *)
     36
     37    int mzed_read_elem(const_mzed_t *M, const_size_t row, const_size_t col)
     38
     39    void mzed_write_elem(mzed_t *, const_size_t row, const_size_t col, const_int elem)
     40
     41    void mzed_set_ui(mzed_t *A, m4ri_word value)
     42
     43    mzed_t *mzed_copy(mzed_t *o, const_mzed_t *i)
     44
     45    int mzed_cmp(mzed_t *l, mzed_t *r)
     46
     47    mzed_t *mzed_randomize(mzed_t *)
     48
     49    mzed_t *mzed_add(mzed_t *, mzed_t *, mzed_t *)
     50
     51    size_t mzed_echelonize_naive(mzed_t *, size_t)
     52
     53    void mzed_add_elem(mzed_t *a, const_size_t row, const_size_t col, const_int elem)
     54
     55    void mzed_add_multiple_of_row(mzed_t *A, size_t ar, mzed_t *B, size_t br, m4ri_word *X, size_t start_col)
     56
     57    void mzed_rescale_row(mzed_t *A, size_t r, size_t c, m4ri_word *X)
     58
     59    void mzed_row_swap(mzed_t *M, const_size_t rowa, const_size_t rowb)
     60
     61    void mzed_copy_row(mzed_t* B, size_t i, const_mzed_t* A, size_t j)
     62
     63    void mzed_col_swap(mzed_t *M, const_size_t cola, const_size_t colb)
     64
     65    void mzed_row_add(mzed_t *M, const_size_t sourcerow, const_size_t destrow)
     66
     67    size_t mzed_first_zero_row(mzed_t *A)
     68
     69    int mzed_is_zero(mzed_t *A)
     70
     71    void mzed_row_clear_offset(mzed_t *M, const_size_t row, const_size_t coloffset)
     72
     73    mzed_t *mzed_concat(mzed_t *C, const_mzed_t *A, const_mzed_t *B)
     74
     75    mzed_t *mzed_stack(mzed_t *C, const_mzed_t *A, const_mzed_t *B)
     76
     77    mzed_t *mzed_submatrix(mzed_t *S, const_mzed_t *M, size_t lowr, size_t lowc, size_t highr, size_t highc)
     78
     79    mzed_t *mzed_mul(mzed_t *C, const_mzed_t *A, const_mzed_t *B)
     80
     81    mzed_t *mzed_mul_naive(mzed_t *C, const_mzed_t *A, const_mzed_t *B)
     82
     83    mzed_t *mzed_mul_scalar(mzed_t *C, m4ri_word a, const_mzed_t *B)
     84
     85    # TODO: not implemented yet in m4rie
     86    mzed_t *mzed_transpose(mzed_t *DST, const_mzed_t *A )
     87
     88    void mzed_print(const_mzed_t *M)
     89 
     90    mzed_t *mzed_invert_travolta(mzed_t *A, mzed_t *B)
     91
     92    # TODO: not implemented yet in m4rie
     93    double mzed_density(mzed_t *A, int res)
     94
     95    # TODO: not implemented yet in m4rie
     96    double _mzed_density(mzed_t *A, int res, size_t r, size_t c)
     97
     98#cdef extern from "m4rie/travolta.h":
     99    size_t mzed_echelonize_travolta(mzed_t *, size_t)
     100
     101    mzed_t *mzed_mul_travolta(mzed_t *, mzed_t *, mzed_t *)
     102
     103#cdef extern from "m4rie/echelonform.h":
     104    size_t mzed_echelonize(mzed_t *, size_t)
     105
     106    size_t mzed_echelonize_ple(mzed_t *, size_t)
     107
     108#cdef extern from "m4rie/strassen.h":
     109    mzed_t *mzed_mul_strassen(mzed_t *, mzed_t *, mzed_t *, size_t cutoff)
     110
     111    size_t _mzed_strassen_cutoff(mzed_t *C, mzed_t *A, mzed_t *B)
     112
     113#cdef extern from "m4rie/bitslice.h":
     114
     115    int __M4RIE_MAX_KARATSUBA_DEGREE
     116
     117    ctypedef struct mzd_slice_t:
     118        mzd_t *x[16]
     119        gf2e *finite_field
     120        int nrows
     121        int ncols
     122        int depth
     123
     124    mzd_slice_t *mzd_slice_init(gf2e *ff, size_t m, size_t n)
     125
     126    void mzd_slice_free(mzd_slice_t *A)
     127
     128    mzed_t *mzed_cling(mzed_t *A, mzd_slice_t *Z)
     129
     130    mzd_slice_t *mzed_slice(mzd_slice_t *A, mzed_t *Z)
     131
     132    mzd_slice_t *mzd_slice_concat(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     133
     134    mzd_slice_t *mzd_slice_stack(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     135
     136    mzd_slice_t *mzd_slice_submatrix(mzd_slice_t *S, mzd_slice_t *A, size_t lowr, size_t lowc, size_t highr, size_t highc)
     137
     138    mzd_slice_t *mzd_slice_init_window(mzd_slice_t *A, size_t lowr, size_t lowc, size_t highr, size_t highc)
     139
     140    void mzd_slice_free_window(mzd_slice_t *A)
     141
     142    mzd_slice_t *_mzd_slice_add(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     143
     144    mzd_slice_t *mzd_slice_add(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     145
     146    void mzd_slice_randomize(mzd_slice_t *A)
     147
     148    mzd_slice_t *mzd_slice_copy(mzd_slice_t *B, mzd_slice_t *A)
     149
     150    void mzd_slice_set_ui(mzd_slice_t *A, m4ri_word value)
     151
     152    m4ri_word mzd_slice_read_elem(mzd_slice_t *A, size_t row, size_t col)
     153
     154    void mzd_slice_add_elem(mzd_slice_t *A, size_t row, size_t col, m4ri_word elem)
     155
     156    void mzd_slice_write_elem(mzd_slice_t *A, size_t row, size_t col, m4ri_word elem)
     157
     158    int mzd_slice_cmp(mzd_slice_t *A, mzd_slice_t *B)
     159
     160    int mzd_slice_is_zero(mzd_slice_t *A)
     161
     162    void mzd_slice_rescale_row(mzd_slice_t *A, size_t r, size_t c, m4ri_word *X)
     163
     164    void mzd_slice_row_swap(mzd_slice_t *A, size_t rowa, size_t rowb)
     165
     166    void mzd_slice_copy_row(mzd_slice_t* B, size_t i, mzd_slice_t* A, size_t j)
     167 
     168    void mzd_slice_col_swap(mzd_slice_t *A, size_t cola, size_t colb)
     169
     170    void mzd_slice_row_add(mzd_slice_t *A, size_t sourcerow, size_t destrow)
     171
     172
     173    void mzd_slice_row_clear_offset(mzd_slice_t *A, size_t row, size_t coloffset)
     174
     175    void mzd_slice_print(mzd_slice_t *A)
     176
     177    mzd_slice_t *_mzed_slice2(mzd_slice_t *A, mzed_t *Z)
     178
     179    mzd_slice_t *mzed_slice2(mzd_slice_t *A, mzed_t *Z)
     180
     181    mzd_slice_t *_mzed_slice4(mzd_slice_t *A, mzed_t *Z)
     182
     183    mzed_t *_mzed_cling2(mzed_t *A, mzd_slice_t *Z)
     184
     185    mzed_t* mzed_cling2(mzed_t *A, mzd_slice_t *Z)
     186
     187    mzed_t *_mzed_cling4(mzed_t *A, mzd_slice_t *Z)
     188
     189    mzed_t* mzed_cling4(mzed_t *A, mzd_slice_t *Z)
     190
     191    mzed_t *mzed_mul_karatsuba(mzed_t *C, mzed_t *A, mzed_t *B)
     192
     193    mzed_t *mzed_addmul_karatsuba(mzed_t *C, mzed_t *A, mzed_t *B)
     194
     195    mzed_t *_mzed_mul_karatsuba(mzed_t *C, mzed_t *A, mzed_t *B)
     196
     197    mzd_slice_t *_mzd_slice_mul_karatsuba2(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     198
     199    mzd_slice_t *_mzd_slice_mul_karatsuba3(mzd_slice_t *C, mzd_slice_t *A, mzd_slice_t *B)
     200
     201cdef extern from "m4rie/finite_field_givaro.h":
     202    ctypedef struct M4RIE__FiniteField "M4RIE::FiniteField":
     203        int (* pol2log)(int r)
     204        int (* log2pol)(int r)
     205       
     206    gf2e *gf2e_init_givgfq(M4RIE__FiniteField *givgfq)
     207
     208    int mzed_read_elem_log(const_mzed_t *a, const_size_t row, const_size_t col, M4RIE__FiniteField *ff)
     209    void mzed_write_elem_log(mzed_t *a, const_size_t row, const_size_t col, const_int elem, M4RIE__FiniteField *ff)
     210    void mzed_add_elem_log(mzed_t *a, const_size_t row, const_size_t col, const_int elem, M4RIE__FiniteField *ff)
     211       

  • sage/libs/ntl/ntl_mat_GF2E.pyx 

    diff --git a/sage/libs/ntl/ntl_mat_GF2E.pyx b/sage/libs/ntl/ntl_mat_GF2E.pyx
    a b  
    2525
    2626include "../../ext/interrupt.pxi"
    2727include "../../ext/stdsage.pxi"
     28include '../../ext/random.pxi'
    2829include 'misc.pxi'
    2930include 'decl.pxi'
    3031
    … …  
    4950            ncols -- number of columns
    5051            v     -- either a list or a matrix over GF(2^x)
    5152
    52         EXAMPLES:
     53        EXAMPLES::
     54
    5355            sage: k.<a> = GF(2^4)
    5456            sage: ctx = ntl.GF2EContext(k)
    5557            sage: ntl.GF2XHexOutput(1)
    … …  
    146148        """
    147149        Returns the structure that holds the underlying NTL GF2E modulus.
    148150
    149         EXAMPLES:
     151        EXAMPLES::
     152
    150153            sage: ntl.GF2XHexOutput(0)
    151154            sage: ctx = ntl.GF2EContext( ntl.GF2X([1,1,0,1,1,0,0,0,1]) )
    152155            sage: a = ntl.GF2E(ntl.ZZ_pX([1,1,3],2), ctx)
    … …  
    163166
    164167    def __reduce__(self):
    165168        """
     169        EXAMPLE::
     170
    166171            sage: k.<a> = GF(2^4)
    167172            sage: ctx = ntl.GF2EContext(k)
    168173            sage: A = ntl.mat_GF2E(ctx, 5,5, [0..24])
    … …  
    175180        """
    176181        Return the string representation of self.
    177182       
    178         EXAMPLES:
     183        EXAMPLES::
     184
    179185            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    180186            sage: ntl.GF2XHexOutput(1)
    181187            sage: ntl.mat_GF2E(ctx, 2,2,range(4)).__repr__()
    … …  
    189195
    190196    def __mul__(ntl_mat_GF2E self, other):
    191197        """
    192         EXAMPLES:
     198        EXAMPLES::
     199
    193200            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    194201            sage: ntl.GF2XHexOutput(1)
    195202            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    … …  
    214221
    215222    def __sub__(ntl_mat_GF2E self, other):
    216223        """
    217         EXAMPLES:
     224        EXAMPLES::
     225
    218226            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    219227            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    220228            sage: n = ntl.mat_GF2E(ctx, 5,5,[3..27])
    … …  
    239247
    240248    def __add__(ntl_mat_GF2E self, other):
    241249        """
    242         EXAMPLES:
     250        EXAMPLES::
     251
    243252            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    244253            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    245254            sage: n = ntl.mat_GF2E(ctx, 5,5,[3..27])
    … …  
    263272
    264273    def __neg__(ntl_mat_GF2E self):
    265274        """
    266         EXAMPLES:
     275        EXAMPLES::
     276
    267277            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    268278            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    269279            sage: -m == m ## indirect doctest
    … …  
    277287
    278288    def __pow__(ntl_mat_GF2E self, long e, ignored):
    279289        """
    280         EXAMPLES:
     290        EXAMPLES::
     291
    281292            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    282293            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    283294            sage: m**2 == m*m ## indirect doctest
    … …  
    291302       
    292303    def __richcmp__(ntl_mat_GF2E self, other, op):
    293304        """
    294         EXAMPLES:
     305        EXAMPLES::
     306
    295307            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    296308            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    297309            sage: n = ntl.mat_GF2E(ctx, 5,5,[3..27])
    … …  
    319331        """
    320332        Return the number of rows in self.
    321333
    322         EXAMPLES:
     334        EXAMPLES::
     335
    323336            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    324337            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24]) ; m.NumRows()
    325338            5
    … …  
    330343        """
    331344        Return the number of columns in self.
    332345
    333         EXAMPLES:
     346        EXAMPLES::
     347
    334348            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    335349            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24]) ; m.NumCols()
    336350            5
    … …  
    339353
    340354    def __setitem__(self, ij, x):
    341355        """
     356        EXAMPLES::
     357
    342358            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    343359            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    344360            sage: ntl.GF2XHexOutput(0)
    … …  
    375391
    376392    def __getitem__(self, ij):
    377393        """
     394        EXAMPLES::
     395
    378396            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    379397            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    380398            sage: m[0,1]
    … …  
    406424        """
    407425        Returns the determinant.
    408426       
    409         EXAMPLES:
     427        EXAMPLES::
     428
    410429            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    411430            sage: ntl.GF2XHexOutput(0)
    412431            sage: ntl.mat_GF2E(ctx, 5,5,[0..24]).determinant()
    … …  
    429448        ncols columns).
    430449
    431450        INPUT:
    432            ncols -- number of columns to process (default: all)
    433451
    434         EXAMPLES:
     452        - ``ncols`` - number of columns to process (default: all)
     453
     454        EXAMPLES::
     455
    435456            sage: m = ctx = ntl.GF2EContext(ntl.GF2X([1,1,0,1,1,0,0,0,1]))
    436457            sage: ntl.mat_GF2E(ctx, 5,5,[3..27]).gauss()
    437458            5
    … …  
    448469        """
    449470        Returns a list of the entries in this matrix
    450471
    451         EXAMPLES:
     472        EXAMPLES::
     473
    452474            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    453475            sage: m = ntl.mat_GF2E(ctx, 2,2,[ntl.GF2E_random(ctx) for x in xrange(2*2)])
    454476            sage: ntl.GF2XHexOutput(0)
    … …  
    461483        """
    462484        Return True if self is zero, and false otherwise.
    463485       
    464         EXAMPLES:
     486        EXAMPLES::
     487
    465488            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    466489            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    467490            sage: n = ntl.mat_GF2E(ctx, 5,5)
    … …  
    478501
    479502    def _sage_(ntl_mat_GF2E self, k=None):
    480503        """
    481         Returns a \class{Matrix} over a FiniteField representation
     504        Returns a ``Matrix`` over a ``FiniteField`` representation
    482505        of this element.
    483506
    484507        INPUT:
    485             self  -- \class{mat_GF2E} element
    486             k     -- optional GF(2**deg)
     508
     509        - ``k`` - optional GF(2**deg)
    487510
    488511        OUTPUT:
    489512            Matrix over k
    490513
    491         EXAMPLES:
     514        EXAMPLES::
     515
    492516            sage: ctx = ntl.GF2EContext([1,1,0,1])
    493517            sage: m = ntl.mat_GF2E(ctx, 2,2,[3..6])
    494518            sage: ntl.GF2XHexOutput(0)
    … …  
    518542        OUTPUT:
    519543            transposed Matrix
    520544
    521         EXAMPLES:
     545        EXAMPLES::
     546
    522547            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    523548            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    524549            sage: n = m.transpose()
    … …  
    537562        """
    538563        Return $X = A^{-1}$; an error is raised if A is singular.
    539564
    540         EXAMPLES:
     565        EXAMPLES::
     566
    541567            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    542568            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    543569            sage: n = ~m
    … …  
    555581        """
    556582        test if A is the n x n identity matrix
    557583
    558         EXAMPLES:
     584        EXAMPLES::
     585
    559586            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    560587            sage: m = ntl.mat_GF2E(ctx, 5,5,[0..24])
    561588            sage: n = ~m
    … …  
    569596
    570597    def IsDiag(self, long n, ntl_GF2E d):
    571598        """
    572         test if X is an  n x n diagonal matrix with d on diagonal
     599        Test if X is an  n x n diagonal matrix with d on diagonal.
    573600
    574         EXAMPLES:
     601        EXAMPLES::
     602
    575603            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    576604            sage: m = ntl.mat_GF2E(ctx, 3,3,[[0,1],0,0, 0,[0,1],0, 0,0,[0,1]])
    577605            sage: m.IsDiag(2, ntl.GF2E([0,1],ctx))
    … …  
    586614        The rows of X are computed as basis of A's row space.  X is is
    587615        row echelon form.
    588616
    589         EXAMPLE:
     617        EXAMPLE::
     618
    590619            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    591620            sage: m = ntl.mat_GF2E(ctx, 3,3,[0..24])
    592621            sage: ntl.GF2XHexOutput(1)
    … …  
    604633       
    605634    def kernel(self):
    606635        """
    607         Computes a basis for the kernel of the map x -> x*A. where x is a
    608         row vector.
     636        Computes a basis for the kernel of the map ``x -> x*A``, where
     637        ``x`` is a row vector.
    609638
    610         EXAMPLE:
     639        EXAMPLE::
     640
    611641            sage: ctx = ntl.GF2EContext([1,1,0,1,1,0,0,0,1])
    612642            sage: m = ntl.mat_GF2E(ctx, 3,3,[0..24])
    613643            sage: ntl.GF2XHexOutput(1)
    … …  
    619649        mat_GF2E_kernel(X.x, self.x)
    620650        sig_off()
    621651        return X
     652
     653    def randomize(self, density=1, nonzero=False):
     654        """
     655        Randomize ``density`` proportion of the entries of this matrix,
     656        leaving the rest unchanged.
     657
     658        INPUT:
     659       
     660        -  ``density`` - float; proportion (roughly) to be considered for
     661           changes
     662        -  ``nonzero`` - Bool (default: ``False``); whether the new entries
     663           are forced to be non-zero
     664
     665        EXAMPLES::
     666
     667            sage: k.<a> = GF(2^4)
     668            sage: ctx = ntl.GF2EContext(k)
     669            sage: ntl.GF2XHexOutput(1)
     670            sage: A = ntl.mat_GF2E(ctx, 100,100)
     671            sage: A.randomize()
     672            sage: len([e for e in A.list() if e!=0])
     673            9389
     674
     675            sage: A = ntl.mat_GF2E(ctx, 100,100)
     676            sage: A.randomize(nonzero=True)
     677            sage: len([e for e in A.list() if e!=0])
     678            10000
     679
     680            sage: A = ntl.mat_GF2E(ctx, 100,100)
     681            sage: A.randomize(nonzero=True, density=0.1)
     682            sage: len([e for e in A.list() if e!=0])
     683            994
     684
     685        """
     686        cdef long i,j
     687        cdef GF2E_c tmp
     688
     689        cdef float _density = density
     690        cdef randstate rstate = current_randstate()
     691
     692        if _density <= 0:
     693            return
     694        if _density > 1:
     695            _density = 1.0
     696
     697        if not nonzero:
     698            if _density == 1.0:
     699                for i in xrange(self.x.NumRows()):
     700                    for j in xrange(self.x.NumCols()):
     701                        tmp = GF2E_random()
     702                        mat_GF2E_setitem(&self.x, i, j, &tmp)
     703            else:
     704                for i in xrange(self.x.NumRows()):
     705                    for j in xrange(self.x.NumCols()):
     706                        if rstate.c_rand_double() <= _density:
     707                            tmp = GF2E_random()
     708                            mat_GF2E_setitem(&self.x, i, j, &tmp)
     709        else:
     710            if _density == 1.0:
     711                for i in xrange(self.x.NumRows()):
     712                    for j in xrange(self.x.NumCols()):
     713                        tmp = GF2E_random()
     714                        while GF2E_IsZero(tmp):
     715                            tmp = GF2E_random()                           
     716                        mat_GF2E_setitem(&self.x, i, j, &tmp)
     717            else:
     718                for i in xrange(self.x.NumRows()):
     719                    for j in xrange(self.x.NumCols()):
     720                        if rstate.c_rand_double() <= _density:
     721                            tmp = GF2E_random()
     722                            while GF2E_IsZero(tmp):
     723                                tmp = GF2E_random()                           
     724                            mat_GF2E_setitem(&self.x, i, j, &tmp)

  • sage/matrix/matrix_mod2_dense.pxd 

    diff --git a/sage/matrix/matrix_mod2_dense.pxd b/sage/matrix/matrix_mod2_dense.pxd
    a b  
    11# choose: dense or sparse
    22
    3 cimport matrix_dense 
     3cimport matrix_dense
    44
    55from sage.libs.m4ri cimport *
    66

  • new file sage/matrix/matrix_mod2e_dense.pxd 

    diff --git a/sage/matrix/matrix_mod2e_dense.pxd b/sage/matrix/matrix_mod2e_dense.pxd
new file mode 100644
    - +  
     1# choose: dense or sparse
     2
     3from sage.rings.finite_rings.element_givaro cimport GivaroGfq, Cache_givaro
     4
     5from sage.libs.m4rie cimport mzed_t
     6
     7cimport matrix_dense
     8
     9cdef class Matrix_mod2e_dense(matrix_dense.Matrix_dense):
     10    cdef mzed_t *_entries
     11    cdef Cache_givaro cc
     12    cdef object _one
     13    cdef object _zero
     14
     15    cpdef Matrix_mod2e_dense _multiply_travolta(Matrix_mod2e_dense self, Matrix_mod2e_dense right)
     16    cpdef Matrix_mod2e_dense _multiply_karatsuba(Matrix_mod2e_dense self, Matrix_mod2e_dense right)
     17    cpdef Matrix_mod2e_dense _multiply_strassen(Matrix_mod2e_dense self, Matrix_mod2e_dense right, cutoff=*)

  • new file sage/matrix/matrix_mod2e_dense.pyx 

    diff --git a/sage/matrix/matrix_mod2e_dense.pyx b/sage/matrix/matrix_mod2e_dense.pyx
new file mode 100644
    - +  
     1"""
     2Dense matrices over `\GF{2^e}` for `2 <= e <= 10` using the M4RIE library.
     3
     4The M4RIE library offers two matrix representations:
     5
     61) ``mzed_t``
     7
     8  m x n matrices over `\GF{2^e}` are internally represented roughly as
     9  m x (en) matrices over `\GF{2}`. Several elements are packed into
     10  words such that each element is filled with zeroes until the next
     11  power of two. Thus, for example, elements of `\GF{2^3}` are
     12  represented as ``[0xxx|0xxx|0xxx|0xxx|...]``. This representation is
     13  wrapped as :class:`Matrix_mod2e_dense` in Sage.
     14
     15  Multiplication and elimination both use "Travolta" tables. These
     16  tables are simply all possible multiples of a given row in a matrix
     17  such that a scale+add operation is reduced to a table lookup +
     18  add. On top of Travolta multiplication M4RIE implements
     19  asymptotically fast Strassen-Winograd multiplication. Elimination
     20  uses simple Gaussian elimination which requires `O(n^3)` additions
     21  but only `O(n^2 * 2^e)` multiplications.
     22
     232) ``mzd_slice_t``
     24
     25  m x n matrices over `\GF{2^e}` are internally represented as slices
     26  of m x n matrices over `\GF{2}`. This representation allows for very
     27  fast matrix times matrix products using Karatsuba's polynomial
     28  multiplication for polynomials over matrices. However, it is not
     29  feature complete yet and hence not wrapped in Sage for now.
     30
     31See http://m4ri.sagemath.org for more details on the M4RIE library.
     32
     33EXAMPLE::
     34
     35    sage: K.<a> = GF(2^8)
     36    sage: A = random_matrix(K, 3,4)
     37    sage: A
     38    [                  a^3 + a^2 + 1         a^7 + a^5 + a^4 + a + 1                       a^7 + a^3                   a^6 + a^4 + a]
     39    [      a^5 + a^4 + a^3 + a^2 + 1 a^7 + a^6 + a^5 + a^4 + a^3 + a             a^5 + a^4 + a^3 + a             a^6 + a^4 + a^2 + 1]
     40    [  a^6 + a^5 + a^4 + a^2 + a + 1   a^7 + a^5 + a^3 + a^2 + a + 1       a^7 + a^6 + a^3 + a^2 + a             a^7 + a^5 + a^3 + a]
     41
     42    sage: A.echelon_form()
     43    [                        1                         0                         0       a^7 + a^6 + a^2 + a]
     44    [                        0                         1                         0     a^6 + a^4 + a^3 + a^2]
     45    [                        0                         0                         1 a^6 + a^4 + a^3 + a^2 + a]
     46
     47AUTHOR:
     48
     49 * Martin Albrecht <martinralbrecht@googlemail.com>
     50
     51TESTS::
     52
     53    sage: TestSuite(sage.matrix.matrix_mod2e_dense.Matrix_mod2e_dense).run(verbose=True)
     54    running ._test_pickling() . . . pass
     55
     56TODO:
     57
     58 - wrap ``mzd_slice_t``
     59
     60           
     61REFERENCES:
     62
     63.. [BB09] Tomas J. Boothby and Robert W. Bradshaw. *Bitslicing
     64and the Method of Four Russians Over Larger Finite
     65Fields*. arXiv:0901.1413v1,
     662009. http://arxiv.org/abs/0901.1413
     67"""
     68
     69include "../ext/interrupt.pxi"
     70include "../ext/cdefs.pxi"
     71include '../ext/stdsage.pxi'
     72include '../ext/random.pxi'
     73
     74cimport matrix_dense
     75from sage.structure.element cimport Matrix, Vector
     76from sage.structure.element cimport ModuleElement, Element, RingElement
     77
     78from sage.rings.all import FiniteField as GF
     79from sage.rings.finite_rings.element_givaro cimport FiniteField_givaroElement, GivRandom, GivRandomSeeded
     80from sage.misc.randstate cimport randstate, current_randstate
     81
     82from sage.matrix.matrix_mod2_dense cimport Matrix_mod2_dense
     83
     84from sage.libs.m4ri cimport m4ri_word, mzd_copy, mzd_init
     85from sage.libs.m4rie cimport *
     86from sage.libs.m4rie cimport mzed_t, M4RIE__FiniteField
     87
     88
     89# we must keep a copy of the internal finite field representation
     90# around to avoid re-creating it over and over again. Furthermore,
     91# M4RIE assumes pointer equivalence of identical fields.
     92
     93_givaro_cache = {}
     94
     95cdef class M4RIE_finite_field:
     96    """
     97    A thin wrapper around the M4RIE finite field class such that we
     98    can put it in a hash table. This class is not meant for public
     99    consumption.
     100    """
     101    cdef gf2e *ff
     102
     103    def __cinit__(self):
     104        """
     105        EXAMPLE::
     106
     107            sage: from sage.matrix.matrix_mod2e_dense import M4RIE_finite_field
     108            sage: K = M4RIE_finite_field(); K
     109            <sage.matrix.matrix_mod2e_dense.M4RIE_finite_field object at 0x...>
     110        """
     111        pass
     112
     113    def __dealloc__(self):
     114        """
     115        EXAMPLE::
     116
     117            sage: from sage.matrix.matrix_mod2e_dense import M4RIE_finite_field
     118            sage: K = M4RIE_finite_field()
     119            sage: del K
     120        """
     121        if self.ff:
     122            gf2e_free(self.ff)
     123
     124cdef class Matrix_mod2e_dense(matrix_dense.Matrix_dense):
     125    ########################################################################
     126    # LEVEL 1 functionality
     127    ########################################################################
     128    def __cinit__(self, parent, entries, copy, coerce, alloc=True):
     129        """
     130        Create new matrix over `GF(2^e)` for 2<=e<=10.
     131
     132        INPUT:
     133
     134        - ``parent`` - a :class:`MatrixSpace`.
     135        - ``entries`` - may be list or a finite field element.
     136        - ``copy`` - ignored, elements are always copied
     137        - ``coerce`` - ignored, elements are always coerced
     138        - ``alloc`` - if ``True`` the matrix is allocated first (default: ``True``)
     139
     140        EXAMPLES::
     141
     142            sage: K.<a> = GF(2^4)
     143            sage: A = Matrix(K, 3, 4); A
     144            [0 0 0 0]
     145            [0 0 0 0]
     146            [0 0 0 0]
     147
     148            sage: A.randomize(); A
     149            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     150            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     151            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     152
     153            sage: K.<a> = GF(2^3)
     154            sage: A = Matrix(K,3,4); A
     155            [0 0 0 0]
     156            [0 0 0 0]
     157            [0 0 0 0]
     158
     159            sage: A.randomize(); A
     160            [a^2 + 1       0   a + 1       0]
     161            [      a       a   a + 1       a]
     162            [      1     a^2 a^2 + a       0]
     163        """
     164        matrix_dense.Matrix_dense.__init__(self, parent)
     165
     166        cdef M4RIE_finite_field FF
     167
     168        R = parent.base_ring()
     169
     170        self.cc = <Cache_givaro>R._cache
     171       
     172        if alloc and self._nrows and self._ncols:
     173            if self.cc in _givaro_cache:
     174                self._entries = mzed_init((<M4RIE_finite_field>_givaro_cache[self.cc]).ff, self._nrows, self._ncols)
     175            else:
     176                FF = PY_NEW(M4RIE_finite_field)
     177                FF.ff = gf2e_init_givgfq(<M4RIE__FiniteField*>self.cc.objectptr)
     178                self._entries = mzed_init(FF.ff, self._nrows, self._ncols)
     179                _givaro_cache[self.cc] = FF
     180
     181        # cache elements
     182        self._zero = self._base_ring(0)
     183        self._one = self._base_ring(1)
     184
     185    def __dealloc__(self):
     186        """
     187        TESTS::
     188       
     189            sage: K.<a> = GF(2^4)
     190            sage: A = Matrix(K, 1000, 1000)
     191            sage: del A
     192            sage: A = Matrix(K, 1000, 1000)
     193            sage: del A
     194        """
     195        if self._entries:
     196            mzed_free(self._entries)
     197            self._entries = NULL
     198
     199    def __init__(self, parent, entries, copy, coerce):
     200        """
     201        Create new matrix over `GF(2^e)` for 2<=e<=10.
     202
     203        INPUT:
     204
     205        - ``parent`` - a :class:`MatrixSpace`.
     206        - ``entries`` - may be list or a finite field element.
     207        - ``copy`` - ignored, elements are always copied
     208        - ``coerce`` - ignored, elements are always coerced
     209
     210        EXAMPLE::
     211       
     212            sage: K.<a> = GF(2^4)
     213            sage: l = [K.random_element() for _ in range(3*4)]; l
     214            [a^2 + 1, a^3 + 1, 0, 0, a, a^3 + a + 1, a + 1, a + 1, a^2, a^3 + a + 1, a^3 + a, a^3 + a]
     215
     216            sage: A = Matrix(K, 3, 4, l); A
     217            [    a^2 + 1     a^3 + 1           0           0]
     218            [          a a^3 + a + 1       a + 1       a + 1]
     219            [        a^2 a^3 + a + 1     a^3 + a     a^3 + a]
     220
     221            sage: A.list()
     222            [a^2 + 1, a^3 + 1, 0, 0, a, a^3 + a + 1, a + 1, a + 1, a^2, a^3 + a + 1, a^3 + a, a^3 + a]
     223
     224            sage: l[0], A[0,0]
     225            (a^2 + 1, a^2 + 1)
     226
     227            sage: A = Matrix(K, 3, 3, a); A
     228            [a 0 0]
     229            [0 a 0]
     230            [0 0 a]
     231        """
     232        cdef int i,j
     233        cdef FiniteField_givaroElement e
     234
     235        if entries is None:
     236            return
     237
     238        R = self.base_ring()
     239       
     240        # scalar ? 
     241        if not isinstance(entries, list):
     242            if entries != 0:
     243                if self.nrows() != self.ncols():
     244                    raise TypeError("self must be a  square matrices for scalar assignment")
     245                for i in range(self.nrows()):
     246                    self.set_unsafe(i,i, R(entries))
     247            return
     248
     249        # all entries are given as a long list
     250        if len(entries) != self._nrows * self._ncols:
     251            raise IndexError("The vector of entries has the wrong length.")
     252       
     253        k = 0
     254
     255        for i from 0 <= i < self._nrows:
     256            if PyErr_CheckSignals(): raise KeyboardInterrupt
     257            for j from 0 <= j < self._ncols:
     258                e = R(entries[k])
     259               
     260                mzed_write_elem_log(self._entries,i,j, e.element, <M4RIE__FiniteField*>self.cc.objectptr)
     261                k = k + 1
     262           
     263    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
     264        """
     265        A[i,j] = value without bound checks
     266
     267        INPUT:
     268        - ``i`` - row index
     269        - ``j`` - column index
     270        - ``value`` - a finite field element (not checked but assumed)
     271
     272        EXAMPLE::
     273       
     274            sage: K.<a> = GF(2^4)
     275            sage: A = Matrix(K,3,4,[K.random_element() for _ in range(3*4)]); A
     276            [    a^2 + 1     a^3 + 1           0           0]
     277            [          a a^3 + a + 1       a + 1       a + 1]
     278            [        a^2 a^3 + a + 1     a^3 + a     a^3 + a]
     279
     280            sage: A[0,0] = a                     # indirect doctest
     281            sage: A
     282            [          a     a^3 + 1           0           0]
     283            [          a a^3 + a + 1       a + 1       a + 1]
     284            [        a^2 a^3 + a + 1     a^3 + a     a^3 + a]
     285        """
     286        mzed_write_elem_log(self._entries, i, j, (<FiniteField_givaroElement>value).element, <M4RIE__FiniteField*>self.cc.objectptr)
     287
     288    cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
     289        """
     290        Get A[i,j] without bound checks.
     291
     292        INPUT:
     293        - ``i`` - row index
     294        - ``j`` - column index
     295
     296        EXAMPLE::
     297       
     298            sage: K.<a> = GF(2^4)
     299            sage: A = random_matrix(K,3,4)
     300            sage: A[2,3]                         # indirect doctest
     301            a^3 + a + 1
     302            sage: K.<a> = GF(2^3)
     303            sage: m,n  = 3, 4
     304            sage: A = random_matrix(K,3,4); A
     305            [a^2 + 1       0   a + 1       0]
     306            [      a       a   a + 1       a]
     307            [      1     a^2 a^2 + a       0]
     308        """
     309        cdef int r = mzed_read_elem_log(self._entries, i, j, <M4RIE__FiniteField*>self.cc.objectptr)
     310        return self.cc._new_c(r)
     311
     312    cpdef ModuleElement _add_(self, ModuleElement right):
     313        """
     314        Return A+B
     315
     316        INPUT:
     317
     318        - ``right`` - a matrix
     319
     320        EXAMPLE::
     321       
     322            sage: K.<a> = GF(2^4)
     323            sage: A = random_matrix(K,3,4); A
     324            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     325            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     326            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     327           
     328            sage: B = random_matrix(K,3,4); B
     329            [          a^3 + 1                 0               a^3                 0]
     330            [          a^3 + a                 a     a^3 + a^2 + a           a^3 + a]
     331            [      a^3 + a + 1               a^2 a^3 + a^2 + a + 1           a^2 + 1]
     332
     333            sage: C = A + B; C # indirect doctest
     334            [    a^3 + a^2 a^3 + a^2 + a         a + 1         a + 1]
     335            [      a^3 + 1             1           a^2       a^2 + a]
     336            [      a^2 + 1       a^2 + a a^3 + a^2 + a a^3 + a^2 + a]
     337        """
     338        cdef Matrix_mod2e_dense A
     339        A = Matrix_mod2e_dense.__new__(Matrix_mod2e_dense, self._parent, 0, 0, 0, alloc=False)
     340        if self._nrows == 0 or self._ncols == 0:
     341            return A
     342        A._entries = mzed_add(NULL, self._entries, (<Matrix_mod2e_dense>right)._entries)
     343
     344        return A
     345
     346    cpdef ModuleElement _sub_(self, ModuleElement right):
     347        """
     348        EXAMPLE::
     349       
     350            sage: from sage.matrix.matrix_mod2e_dense import Matrix_mod2e_dense
     351            sage: K.<a> = GF(2^4)
     352            sage: m,n  = 3, 4
     353            sage: MS = MatrixSpace(K,m,n)
     354            sage: A = random_matrix(K,3,4); A
     355            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     356            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     357            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     358
     359            sage: B = random_matrix(K,3,4); B
     360            [          a^3 + 1                 0               a^3                 0]
     361            [          a^3 + a                 a     a^3 + a^2 + a           a^3 + a]
     362            [      a^3 + a + 1               a^2 a^3 + a^2 + a + 1           a^2 + 1]
     363
     364            sage: C = A - B; C  # indirect doctest
     365            [    a^3 + a^2 a^3 + a^2 + a         a + 1         a + 1]
     366            [      a^3 + 1             1           a^2       a^2 + a]
     367            [      a^2 + 1       a^2 + a a^3 + a^2 + a a^3 + a^2 + a]
     368        """
     369        return self._add_(right)
     370
     371    def _multiply_classical(self, Matrix right):
     372        """
     373        Classical cubic matrix multiplication.
     374
     375        EXAMPLES::
     376
     377            sage: K.<a> = GF(2^2)
     378            sage: A = random_matrix(K, 50, 50)
     379            sage: B = random_matrix(K, 50, 50)
     380            sage: A*B == A._multiply_classical(B)
     381            True
     382
     383            sage: K.<a> = GF(2^4)
     384            sage: A = random_matrix(K, 50, 50)
     385            sage: B = random_matrix(K, 50, 50)
     386            sage: A*B == A._multiply_classical(B)
     387            True
     388
     389            sage: K.<a> = GF(2^8)
     390            sage: A = random_matrix(K, 50, 50)
     391            sage: B = random_matrix(K, 50, 50)
     392            sage: A*B == A._multiply_classical(B)
     393            True
     394
     395            sage: K.<a> = GF(2^10)
     396            sage: A = random_matrix(K, 50, 50)
     397            sage: B = random_matrix(K, 50, 50)
     398            sage: A*B == A._multiply_classical(B)
     399            True
     400
     401        .. note::
     402
     403            This function is very slow. Use ``*`` operator instead.
     404        """
     405        if self._ncols != right._nrows:
     406            raise ArithmeticError("left ncols must match right nrows")
     407
     408        cdef Matrix_mod2e_dense ans
     409       
     410        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     411        if self._nrows == 0 or self._ncols == 0 or right._ncols == 0:
     412            return ans
     413        _sig_on
     414        ans._entries = mzed_mul_naive(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries)
     415        _sig_off
     416        return ans
     417
     418    cdef Matrix _matrix_times_matrix_(self, Matrix right):
     419        """
     420        Return A*B
     421
     422        Uses the M4RIE machinery to decide which function to call.
     423
     424        INPUT:
     425
     426        - ``right`` - a matrix
     427
     428        EXAMPLES::
     429
     430            sage: K.<a> = GF(2^3)
     431            sage: A = random_matrix(K, 51, 50)
     432            sage: B = random_matrix(K, 50, 52)
     433            sage: A*B == A._multiply_travolta(B) # indirect doctest
     434            True
     435
     436            sage: K.<a> = GF(2^5)
     437            sage: A = random_matrix(K, 10, 50)
     438            sage: B = random_matrix(K, 50, 12)
     439            sage: A*B == A._multiply_travolta(B)
     440            True
     441
     442            sage: K.<a> = GF(2^7)
     443            sage: A = random_matrix(K,100, 50)
     444            sage: B = random_matrix(K, 50, 17)
     445            sage: A*B == A._multiply_classical(B)
     446            True
     447        """
     448        if self._ncols != right._nrows:
     449            raise ArithmeticError("left ncols must match right nrows")
     450
     451        cdef Matrix_mod2e_dense ans
     452       
     453        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     454        if self._nrows == 0 or self._ncols == 0 or right._ncols == 0:
     455            return ans
     456        _sig_on
     457        ans._entries = mzed_mul(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries)
     458        _sig_off
     459        return ans
     460
     461    cpdef Matrix_mod2e_dense _multiply_travolta(Matrix_mod2e_dense self, Matrix_mod2e_dense right):
     462        """
     463        Return A*B using Travolta tables.
     464
     465        We can write classical cubic multiplication (``C=A*B``) as::
     466
     467        for i in range(A.ncols()):
     468           for j in range(A.nrows()):
     469             C[j] += A[j,i] * B[j]
     470
     471        Hence, in the inner-most loop we compute multiples of ``B[j]``
     472        by the values ``A[j,i]``. If the matrix ``A`` is big and the
     473        finite field is small, there is a very high chance that
     474        ``e * B[j]`` is computed more than once for any ``e`` in the finite
     475        field. Instead, we compute all possible
     476        multiples of ``B[j]`` and re-use this data in the inner loop.
     477        This is what is called a "Travolta" table in M4RIE.
     478
     479        INPUT:
     480
     481        - ``right`` - a matrix
     482
     483        EXAMPLES::
     484
     485            sage: K.<a> = GF(2^2)
     486            sage: A = random_matrix(K, 50, 50)
     487            sage: B = random_matrix(K, 50, 50)
     488            sage: A._multiply_travolta(B) == A._multiply_classical(B) # indirect doctest
     489            True
     490
     491            sage: K.<a> = GF(2^4)
     492            sage: A = random_matrix(K, 50, 50)
     493            sage: B = random_matrix(K, 50, 50)
     494            sage: A._multiply_travolta(B) == A._multiply_classical(B)
     495            True
     496
     497            sage: K.<a> = GF(2^8)
     498            sage: A = random_matrix(K, 50, 50)
     499            sage: B = random_matrix(K, 50, 50)
     500            sage: A._multiply_travolta(B) == A._multiply_classical(B)
     501            True
     502
     503            sage: K.<a> = GF(2^10)
     504            sage: A = random_matrix(K, 50, 50)
     505            sage: B = random_matrix(K, 50, 50)
     506            sage: A._multiply_travolta(B) == A._multiply_classical(B)
     507            True
     508        """
     509        if self._ncols != right._nrows:
     510            raise ArithmeticError("left ncols must match right nrows")
     511
     512        cdef Matrix_mod2e_dense ans
     513       
     514        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     515        if self._nrows == 0 or self._ncols == 0 or right._ncols == 0:
     516            return ans
     517
     518        _sig_on
     519        ans._entries = mzed_mul_travolta(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries)
     520        _sig_off
     521        return ans
     522
     523    cpdef Matrix_mod2e_dense _multiply_karatsuba(Matrix_mod2e_dense self, Matrix_mod2e_dense right):
     524        """
     525        Matrix multiplication using Karatsuba over polynomials with
     526        matrix coefficients over GF(2).
     527
     528        The idea behind Karatsuba multiplication for matrices over
     529        `\GF{p^n}` is to treat these matrices as polynomials with
     530        coefficients of matrices over `\GF{p}`. Then, Karatsuba-style
     531        formulas can be used to perform multiplication, cf. [BB09]_.
     532
     533        INPUT:
     534
     535        - ``right`` - a matrix
     536       
     537        EXAMPLES::
     538
     539            sage: K.<a> = GF(2^2)
     540            sage: A = random_matrix(K, 50, 50)
     541            sage: B = random_matrix(K, 50, 50)
     542            sage: A._multiply_karatsuba(B) == A._multiply_classical(B) # indirect doctest
     543            True
     544
     545            sage: K.<a> = GF(2^2)
     546            sage: A = random_matrix(K, 137, 11)
     547            sage: B = random_matrix(K, 11, 23)
     548            sage: A._multiply_karatsuba(B) == A._multiply_classical(B)
     549            True
     550
     551        TESTS::
     552
     553            sage: K.<a> = GF(2^5)
     554            sage: A = random_matrix(K, 50, 50)
     555            sage: B = random_matrix(K, 50, 50)
     556            sage: A._multiply_karatsuba(B) == A._multiply_classical(B)
     557            Traceback (most recent call last):
     558            ...
     559            NotImplementedError: Karatsuba multiplication is only implemented for matrices over GF(2^e) with e <= 4.
     560        """
     561        if self._ncols != right._nrows:
     562            raise ArithmeticError("left ncols must match right nrows")
     563
     564        if self._entries.finite_field.degree > __M4RIE_MAX_KARATSUBA_DEGREE:
     565            raise NotImplementedError("Karatsuba multiplication is only implemented for matrices over GF(2^e) with e <= %d."%__M4RIE_MAX_KARATSUBA_DEGREE)
     566
     567        cdef Matrix_mod2e_dense ans
     568       
     569        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     570        if self._nrows == 0 or self._ncols == 0 or right._ncols == 0:
     571            return ans
     572
     573        _sig_on
     574        ans._entries = mzed_mul_karatsuba(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries)
     575        _sig_off
     576        return ans
     577
     578    cpdef Matrix_mod2e_dense _multiply_strassen(Matrix_mod2e_dense self, Matrix_mod2e_dense right, cutoff=0):
     579        """
     580        Winograd-Strassen matrix multiplication with Travolta
     581        multiplication as base case.
     582
     583        INPUT:
     584
     585        - ``right`` - a matrix
     586        - ``cutoff`` - row or column dimension to switch over to
     587          Travolta multiplication (default: 64)
     588
     589        EXAMPLES::
     590
     591            sage: K.<a> = GF(2^2)
     592            sage: A = random_matrix(K, 50, 50)
     593            sage: B = random_matrix(K, 50, 50)
     594            sage: A._multiply_strassen(B) == A._multiply_classical(B) # indirect doctest
     595            True
     596
     597            sage: K.<a> = GF(2^4)
     598            sage: A = random_matrix(K, 50, 50)
     599            sage: B = random_matrix(K, 50, 50)
     600            sage: A._multiply_strassen(B) == A._multiply_classical(B)
     601            True
     602
     603            sage: K.<a> = GF(2^8)
     604            sage: A = random_matrix(K, 50, 50)
     605            sage: B = random_matrix(K, 50, 50)
     606            sage: A._multiply_strassen(B) == A._multiply_classical(B)
     607            True
     608
     609            sage: K.<a> = GF(2^10)
     610            sage: A = random_matrix(K, 50, 50)
     611            sage: B = random_matrix(K, 50, 50)
     612            sage: A._multiply_strassen(B) == A._multiply_classical(B)
     613            True
     614        """
     615        if self._ncols != right._nrows:
     616            raise ArithmeticError("left ncols must match right nrows")
     617
     618        cdef Matrix_mod2e_dense ans
     619       
     620        ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
     621        if self._nrows == 0 or self._ncols == 0 or right._ncols == 0:
     622            return ans
     623
     624        if cutoff == 0:
     625            cutoff = _mzed_strassen_cutoff(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries)
     626
     627        _sig_on
     628        ans._entries = mzed_mul_strassen(ans._entries, self._entries, (<Matrix_mod2e_dense>right)._entries, cutoff)
     629        _sig_off
     630        return ans
     631
     632    cpdef ModuleElement _lmul_(self, RingElement right):
     633        """
     634        Return ``a*B`` for ``a`` an element of the base field.
     635
     636        INPUT:
     637
     638        - ``right`` - an element of the base field
     639
     640        EXAMPLES::
     641
     642             sage: K.<a> = GF(4)
     643             sage: A = random_matrix(K,10,10)
     644             sage: A
     645             [    0     1     1     0     0     0     a a + 1     1     a]
     646             [    1     1 a + 1     a a + 1     0     1     1     1 a + 1]
     647             [    1     0     a     1     1 a + 1     a     1 a + 1 a + 1]
     648             [    0     a     a     a     0     1     1     1     0     1]
     649             [    1     1     a     a     a a + 1 a + 1 a + 1     1     0]
     650             [a + 1 a + 1     a     a     0 a + 1     0     a     1     1]
     651             [a + 1 a + 1     0 a + 1     0 a + 1 a + 1     a     a     a]
     652             [    0     1     a     0     1     a a + 1 a + 1     a a + 1]
     653             [    0     1     a     1     1 a + 1     1 a + 1     a     a]
     654             [a + 1     a     0     a     a     a     0     a     0 a + 1]
     655             
     656             sage: a*A # indirect doctest
     657             [    0     a     a     0     0     0 a + 1     1     a a + 1]
     658             [    a     a     1 a + 1     1     0     a     a     a     1]
     659             [    a     0 a + 1     a     a     1 a + 1     a     1     1]
     660             [    0 a + 1 a + 1 a + 1     0     a     a     a     0     a]
     661             [    a     a a + 1 a + 1 a + 1     1     1     1     a     0]
     662             [    1     1 a + 1 a + 1     0     1     0 a + 1     a     a]
     663             [    1     1     0     1     0     1     1 a + 1 a + 1 a + 1]
     664             [    0     a a + 1     0     a a + 1     1     1 a + 1     1]
     665             [    0     a a + 1     a     a     1     a     1 a + 1 a + 1]
     666             [    1 a + 1     0 a + 1 a + 1 a + 1     0 a + 1     0     1]
     667        """
     668        if not isinstance(right, FiniteField_givaroElement):
     669            raise TypeError("right must be element of type 'FiniteField_givaroElement'")
     670
     671        cdef Matrix_mod2e_dense C = Matrix_mod2e_dense.__new__(Matrix_mod2e_dense, self._parent, 0, 0, 0)
     672        cdef m4ri_word a = <m4ri_word>(<M4RIE__FiniteField*>self.cc.objectptr).log2pol((<FiniteField_givaroElement>right).element)
     673        mzed_mul_scalar(C._entries, a, self._entries)
     674        return C
     675
     676    def __neg__(self):
     677        """
     678        EXAMPLE::
     679       
     680            sage: K.<a> = GF(2^4)
     681            sage: A = random_matrix(K, 3, 4); A
     682            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     683            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     684            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     685
     686            sage: -A
     687            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     688            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     689            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     690        """
     691        return self.__copy__()
     692
     693    def __richcmp__(Matrix self, right, int op):  # always need for mysterious reasons.
     694        """
     695        EXAMPLE::
     696       
     697            sage: K.<a> = GF(2^4)
     698            sage: A = random_matrix(K,3,4)
     699            sage: B = copy(A)
     700            sage: A == B
     701            True
     702            sage: A[0,0] = a
     703            sage: A == B
     704            False
     705        """
     706        return self._richcmp(right, op)
     707
     708    cdef int _cmp_c_impl(self, Element right) except -2:
     709        if self._nrows == 0 or self._ncols == 0:
     710            return 0
     711        return mzed_cmp(self._entries, (<Matrix_mod2e_dense>right)._entries)
     712
     713    def __copy__(self):
     714        """
     715        EXAMPLE::
     716       
     717            sage: K.<a> = GF(2^4)
     718            sage: m,n  = 3, 4
     719            sage: A = random_matrix(K,3,4); A
     720            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     721            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     722            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     723
     724            sage: A2 = copy(A); A2
     725            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     726            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     727            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     728
     729            sage: A[0,0] = 1
     730            sage: A2[0,0]
     731            a^2 + 1
     732        """
     733        cdef Matrix_mod2e_dense A
     734        A = Matrix_mod2e_dense.__new__(Matrix_mod2e_dense, self._parent, 0, 0, 0)
     735       
     736        if self._nrows and self._ncols:
     737            mzed_copy(A._entries, <const_mzed_t *>self._entries)
     738
     739        return A
     740       
     741    def _list(self):
     742        """
     743        EXAMPLE::
     744       
     745            sage: K.<a> = GF(2^4)
     746            sage: m,n  = 3, 4
     747            sage: A = random_matrix(K,3,4); A
     748            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1]
     749            [        a + 1         a + 1       a^3 + a     a^3 + a^2]
     750            [a^3 + a^2 + a             a             1   a^3 + a + 1]
     751
     752            sage: A.list() # indirect doctest
     753            [a^2 + 1, a^3 + a^2 + a, a^3 + a + 1, a + 1, a + 1, a + 1, a^3 + a, a^3 + a^2, a^3 + a^2 + a, a, 1, a^3 + a + 1]
     754        """
     755        cdef int i,j
     756        l = []
     757        for i from 0 <= i < self._nrows:
     758            for j from 0 <= j < self._ncols:
     759                l.append(self.get_unsafe(i,j))
     760        return l
     761
     762    def randomize(self, density=1, nonzero=False, *args, **kwds):
     763        """
     764        Randomize ``density`` proportion of the entries of this matrix,
     765        leaving the rest unchanged.
     766       
     767        INPUT:
     768
     769        -  ``density`` - float; proportion (roughly) to be considered for
     770           changes
     771        -  ``nonzero`` - Bool (default: ``False``); whether the new entries
     772           are forced to be non-zero
     773
     774        OUTPUT:
     775       
     776        -  None, the matrix is modified in-place
     777       
     778        EXAMPLE::
     779       
     780            sage: K.<a> = GF(2^4)
     781            sage: A = Matrix(K,3,3)
     782
     783            sage: A.randomize(); A
     784            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1]
     785            [        a + 1         a + 1         a + 1]
     786            [      a^3 + a     a^3 + a^2 a^3 + a^2 + a]
     787
     788            sage: K.<a> = GF(2^4)
     789            sage: A = random_matrix(K,1000,1000,density=0.1)
     790            sage: float(A.density())
     791            0.099738...
     792
     793            sage: A = random_matrix(K,1000,1000,density=1.0)
     794            sage: float(A.density())
     795            1.0
     796
     797            sage: A = random_matrix(K,1000,1000,density=0.5)
     798            sage: float(A.density())
     799            0.499759...
     800
     801        Note, that the matrix is updated and not zero-ed out before
     802        being randomized::
     803
     804            sage: A = matrix(K,1000,1000)
     805            sage: A.randomize(nonzero=False, density=0.1)
     806            sage: float(A.density())
     807            0.0937679...
     808           
     809            sage: A.randomize(nonzero=False, density=0.05)
     810            sage: float(A.density())
     811            0.13626...
     812        """
     813        if self._ncols == 0 or self._nrows == 0:
     814            return
     815       
     816        cdef Py_ssize_t i,j
     817        cdef int seed = current_randstate().c_random()
     818        cdef GivRandom generator = GivRandomSeeded(seed)
     819        cdef int tmp
     820
     821        self.check_mutability()
     822        self.clear_cache()
     823
     824        cdef randstate rstate = current_randstate()
     825
     826        if self._ncols == 0 or self._nrows == 0:
     827            return
     828
     829        cdef float _density = density
     830        if _density <= 0:
     831            return
     832        if _density > 1:
     833            _density = 1.0
     834
     835        if _density == 1:
     836            if nonzero == False:
     837                _sig_on
     838                for i in range(self._nrows):
     839                    for j in range(self._ncols):
     840                        tmp = self.cc.objectptr.random(generator, tmp)
     841                        mzed_write_elem_log(self._entries, i, j, tmp, <M4RIE__FiniteField*>self.cc.objectptr)
     842                _sig_off
     843            else:
     844                _sig_on
     845                for i in range(self._nrows):
     846                    for j in range(self._ncols):
     847                        tmp = self.cc.objectptr.random(generator,tmp)
     848                        while self.cc.objectptr.isZero(tmp):
     849                            tmp = self.cc.objectptr.random(generator,tmp)
     850                        mzed_write_elem_log(self._entries, i, j, tmp, <M4RIE__FiniteField*>self.cc.objectptr)
     851                _sig_off
     852        else:
     853            if nonzero == False:
     854                _sig_on
     855                for i in range(self._nrows):
     856                    for j in range(self._ncols):
     857                        if rstate.c_rand_double() <= _density:
     858                            tmp = self.cc.objectptr.random(generator, tmp)
     859                            mzed_write_elem_log(self._entries, i, j, tmp, <M4RIE__FiniteField*>self.cc.objectptr)
     860                _sig_off
     861            else:
     862                _sig_on
     863                for i in range(self._nrows):
     864                    for j in range(self._ncols):
     865                        if rstate.c_rand_double() <= _density:
     866                            tmp = self.cc.objectptr.random(generator,tmp)
     867                            while self.cc.objectptr.isZero(tmp):
     868                                tmp = self.cc.objectptr.random(generator,tmp)
     869                            mzed_write_elem_log(self._entries, i, j, tmp, <M4RIE__FiniteField*>self.cc.objectptr)
     870                _sig_off
     871
     872    def echelonize(self, algorithm='heuristic', reduced=True, **kwds):
     873        """
     874        Compute the row echelon form of ``self`` in place.
     875
     876        INPUT:
     877
     878        - ``algorithm`` - one of the following
     879          - ``heuristic`` - let M4RIE decide (default)
     880          - ``travolta`` - use travolta table based algorithm
     881          - ``ple`` - use PLE decomposition
     882          - ``naive`` - use naive cubic Gaussian elimination (M4RIE implementation)
     883          - ``builtin`` - use naive cubic Gaussian elimination (Sage implementation)
     884        - ``reduced`` - if ``True`` return reduced echelon form. No
     885          guarantee is given that the matrix is *not* reduced if
     886          ``False`` (default: ``True``)
     887
     888        EXAMPLE::
     889       
     890            sage: K.<a> = GF(2^4)
     891            sage: m,n  = 3, 5
     892            sage: A = random_matrix(K, 3, 5); A
     893            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1         a + 1         a + 1]
     894            [        a + 1       a^3 + a     a^3 + a^2 a^3 + a^2 + a             a]
     895            [            1   a^3 + a + 1     a^3 + a^2 a^3 + a^2 + 1           a^2]
     896
     897            sage: A.echelonize(); A
     898            [            1             0             0             a             a]
     899            [            0             1             0   a^2 + a + 1             a]
     900            [            0             0             1             a a^3 + a^2 + 1]
     901
     902            sage: K.<a> = GF(2^3)
     903            sage: m,n  = 3, 5
     904            sage: MS = MatrixSpace(K,m,n)
     905            sage: A = random_matrix(K, 3, 5)
     906
     907            sage: copy(A).echelon_form('travolta')
     908            [          1           0           0         a^2     a^2 + 1]
     909            [          0           1           0         a^2           1]
     910            [          0           0           1 a^2 + a + 1       a + 1]
     911           
     912            sage: copy(A).echelon_form('naive');
     913            [          1           0           0         a^2     a^2 + 1]
     914            [          0           1           0         a^2           1]
     915            [          0           0           1 a^2 + a + 1       a + 1]
     916
     917            sage: copy(A).echelon_form('builtin');
     918            [          1           0           0         a^2     a^2 + 1]
     919            [          0           1           0         a^2           1]
     920            [          0           0           1 a^2 + a + 1       a + 1]
     921        """
     922        if self._nrows == 0 or self._ncols == 0:
     923            self.cache('in_echelon_form',True)
     924            self.cache('rank', 0)
     925            self.cache('pivots', [])
     926            return self
     927
     928        cdef int k, n, full
     929
     930        full = int(reduced)
     931
     932        x = self.fetch('in_echelon_form')
     933        if not x is None: return  # already known to be in echelon form
     934
     935        self.check_mutability()
     936        self.clear_cache()       
     937
     938        if algorithm == 'naive':
     939            _sig_on
     940            r =  mzed_echelonize_naive(self._entries, full)
     941            _sig_off
     942           
     943        elif algorithm == 'travolta':
     944            _sig_on
     945            r =  mzed_echelonize_travolta(self._entries, full)
     946            _sig_off
     947           
     948        elif algorithm == 'ple':
     949            _sig_on
     950            r =  mzed_echelonize_ple(self._entries, full)
     951            _sig_off
     952           
     953        elif algorithm == 'heuristic':
     954            _sig_on
     955            r =  mzed_echelonize(self._entries, full)
     956            _sig_off
     957           
     958        elif algorithm == 'builtin':
     959            self._echelon_in_place_classical()
     960
     961        else:
     962            raise ValueError("No algorithm '%s'."%algorithm)
     963
     964        self.cache('in_echelon_form',True)
     965        self.cache('rank', r)
     966        self.cache('pivots', self._pivots())
     967
     968    def _pivots(self):
     969        """
     970        EXAMPLE::
     971
     972            sage: K.<a> = GF(2^8)
     973            sage: A = random_matrix(K, 15, 15)
     974            sage: A.pivots() # indirect doctest
     975            (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)
     976        """
     977        if not self.fetch('in_echelon_form'):
     978            raise ValueError("self must be in reduced row echelon form first.")
     979        pivots = []
     980        cdef Py_ssize_t i, j, nc
     981        nc = self._ncols
     982        i = 0
     983        while i < self._nrows:
     984            for j from i <= j < nc:
     985                if self.get_unsafe(i,j):
     986                    pivots.append(j)
     987                    i += 1
     988                    break
     989            if j == nc:
     990                break
     991        return pivots
     992
     993    def __invert__(self):
     994        """
     995        EXAMPLE::
     996
     997            sage: K.<a> = GF(2^3)
     998            sage: A = random_matrix(K,3,3); A
     999            [      0   a + 1       1]
     1000            [a^2 + a a^2 + a a^2 + a]
     1001            [      a a^2 + 1   a + 1]
     1002
     1003            sage: B = ~A; B
     1004            [        a^2           0     a^2 + 1]
     1005            [a^2 + a + 1         a^2 a^2 + a + 1]
     1006            [      a + 1 a^2 + a + 1           a]
     1007           
     1008            sage: A*B
     1009            [1 0 0]
     1010            [0 1 0]
     1011            [0 0 1]
     1012        """
     1013        cdef Matrix_mod2e_dense A
     1014        A = Matrix_mod2e_dense.__new__(Matrix_mod2e_dense, self._parent, 0, 0, 0)
     1015       
     1016        if self._nrows and self._nrows == self._ncols:
     1017            mzed_invert_travolta(A._entries, self._entries)
     1018
     1019        return A
     1020       
     1021    cdef rescale_row_c(self, Py_ssize_t row, multiple, Py_ssize_t start_col):
     1022        """
     1023        Return ``multiple * self[row][start_col:]``
     1024
     1025        INPUT:
     1026
     1027        - ``row`` - row index for row to rescale
     1028        - ``multiple`` - finite field element to scale by
     1029        - ``start_col`` - only start at this column index.
     1030
     1031        EXAMPLE::
     1032
     1033            sage: K.<a> = GF(2^3)
     1034            sage: A = random_matrix(K,3,3); A
     1035            [      0   a + 1       1]
     1036            [a^2 + a a^2 + a a^2 + a]
     1037            [      a a^2 + 1   a + 1]
     1038
     1039            sage: A.rescale_row(0, a , 0); A
     1040            [      0 a^2 + a       a]
     1041            [a^2 + a a^2 + a a^2 + a]
     1042            [      a a^2 + 1   a + 1]
     1043
     1044            sage: A.rescale_row(0,0,0); A
     1045            [      0       0       0]
     1046            [a^2 + a a^2 + a a^2 + a]
     1047            [      a a^2 + 1   a + 1]
     1048        """       
     1049        cdef m4ri_word x = <m4ri_word>(<M4RIE__FiniteField*>self.cc.objectptr).log2pol((<FiniteField_givaroElement>multiple).element)
     1050        cdef m4ri_word *X = self._entries.finite_field.mul[x]
     1051        mzed_rescale_row(self._entries, row, start_col, X)
     1052
     1053
     1054    cdef add_multiple_of_row_c(self,  Py_ssize_t row_to, Py_ssize_t row_from, multiple, Py_ssize_t start_col):
     1055        """
     1056        Compute ``self[row_to][start_col:] += multiple * self[row_from][start_col:]``.
     1057       
     1058        INPUT:
     1059
     1060        - ``row_to`` - row index of source
     1061        - ``row_from`` - row index of destination
     1062        - ``multiple`` -  finite field element
     1063        - ``start_col`` - only start at this column index
     1064
     1065        EXAMPLE::
     1066
     1067            sage: K.<a> = GF(2^3)
     1068            sage: A = random_matrix(K,3,3); A
     1069            [      0   a + 1       1]
     1070            [a^2 + a a^2 + a a^2 + a]
     1071            [      a a^2 + 1   a + 1]
     1072
     1073            sage: A.add_multiple_of_row(0,1,a,0); A
     1074            [a^2 + a + 1         a^2     a^2 + a]
     1075            [    a^2 + a     a^2 + a     a^2 + a]
     1076            [          a     a^2 + 1       a + 1]
     1077        """
     1078       
     1079        cdef m4ri_word x = <m4ri_word>(<M4RIE__FiniteField*>self.cc.objectptr).log2pol((<FiniteField_givaroElement>multiple).element)
     1080        cdef m4ri_word *X = self._entries.finite_field.mul[x]
     1081        mzed_add_multiple_of_row(self._entries, row_to, self._entries, row_from, X, start_col)
     1082
     1083
     1084    cdef swap_rows_c(self, Py_ssize_t row1, Py_ssize_t row2):
     1085        """
     1086        Swap rows ``row1`` and ``row2``.
     1087
     1088        INPUT:
     1089
     1090        - ``row1`` - row index
     1091        - ``row2`` - row index
     1092
     1093        EXAMPLE::
     1094
     1095            sage: K.<a> = GF(2^3)
     1096            sage: A = random_matrix(K,3,3)
     1097            sage: A
     1098            [      0   a + 1       1]
     1099            [a^2 + a a^2 + a a^2 + a]
     1100            [      a a^2 + 1   a + 1]
     1101
     1102            sage: A.swap_rows(0,1); A
     1103            [a^2 + a a^2 + a a^2 + a]
     1104            [      0   a + 1       1]
     1105            [      a a^2 + 1   a + 1]
     1106
     1107        """
     1108        mzed_row_swap(self._entries, row1, row2)
     1109
     1110    cdef swap_columns_c(self, Py_ssize_t col1, Py_ssize_t col2):
     1111        """
     1112        Swap columns ``col1`` and ``col2``.
     1113
     1114        INPUT:
     1115
     1116        - ``col1`` - column index
     1117        - ``col2`` - column index
     1118
     1119        EXAMPLE::
     1120
     1121            sage: K.<a> = GF(2^3)
     1122            sage: A = random_matrix(K,3,3)
     1123            sage: A
     1124            [      0   a + 1       1]
     1125            [a^2 + a a^2 + a a^2 + a]
     1126            [      a a^2 + 1   a + 1]
     1127
     1128            sage: A.swap_columns(0,1); A
     1129            [  a + 1       0       1]
     1130            [a^2 + a a^2 + a a^2 + a]
     1131            [a^2 + 1       a   a + 1]
     1132
     1133            sage: A = random_matrix(K,4,16)
     1134
     1135            sage: B = copy(A)
     1136            sage: B.swap_columns(0,1)
     1137            sage: B.swap_columns(0,1)
     1138            sage: A == B
     1139            True
     1140
     1141            sage: A.swap_columns(0,15)
     1142            sage: A.column(0) == B.column(15)
     1143            True
     1144            sage: A.swap_columns(14,15)
     1145            sage: A.column(14) == B.column(0)
     1146            True
     1147        """
     1148        mzed_col_swap(self._entries, col1, col2)
     1149
     1150    def augment(self, Matrix_mod2e_dense right):
     1151        """
     1152        Augments ``self`` with ``right``.
     1153
     1154        INPUT:
     1155
     1156        - ``right`` - a matrix
     1157
     1158        EXAMPLE::
     1159
     1160            sage: K.<a> = GF(2^4)
     1161            sage: MS = MatrixSpace(K,3,3)
     1162            sage: A = random_matrix(K,3,3)
     1163            sage: B = A.augment(MS(1)); B
     1164            [      a^2 + 1 a^3 + a^2 + a   a^3 + a + 1             1             0             0]
     1165            [        a + 1         a + 1         a + 1             0             1             0]
     1166            [      a^3 + a     a^3 + a^2 a^3 + a^2 + a             0             0             1]
     1167
     1168            sage: B.echelonize(); B
     1169            [            1             0             0 a^3 + a^2 + 1 a^3 + a^2 + 1       a^2 + a]
     1170            [            0             1             0       a^3 + 1       a^2 + 1       a^2 + 1]
     1171            [            0             0             1           a^2       a^2 + a         a + 1]
     1172
     1173            sage: C = B.matrix_from_columns([3,4,5]); C
     1174            [a^3 + a^2 + 1 a^3 + a^2 + 1       a^2 + a]
     1175            [      a^3 + 1       a^2 + 1       a^2 + 1]
     1176            [          a^2       a^2 + a         a + 1]
     1177
     1178            sage: C == ~A
     1179            True
     1180
     1181            sage: C*A == MS(1)
     1182            True
     1183
     1184        TESTS::
     1185
     1186            sage: K.<a> = GF(2^4)
     1187            sage: A = random_matrix(K,2,3)
     1188            sage: B = random_matrix(K,2,0)
     1189            sage: A.augment(B)
     1190            [a^3 + 1       0     a^3]
     1191            [      0 a^3 + a       a]
     1192
     1193            sage: B.augment(A)
     1194            [a^3 + 1       0     a^3]
     1195            [      0 a^3 + a       a]
     1196
     1197            sage: M = Matrix(K, 0, 0, 0)
     1198            sage: N = Matrix(K, 0, 19, 0)
     1199            sage: W = M.augment(N)
     1200            sage: W.ncols()
     1201            19
     1202
     1203            sage: M = Matrix(K, 0, 1, 0)
     1204            sage: N = Matrix(K, 0, 1, 0)
     1205            sage: M.augment(N)
     1206            []
     1207        """
     1208        cdef Matrix_mod2e_dense A
     1209
     1210        if self._nrows != right._nrows:
     1211            raise TypeError, "Both numbers of rows must match."
     1212
     1213        if self._ncols == 0:
     1214            return right.__copy__()
     1215        if right._ncols == 0:
     1216            return self.__copy__()
     1217
     1218        A = self.new_matrix(ncols = self._ncols + right._ncols)
     1219        if self._nrows == 0:
     1220            return A
     1221        A._entries = mzed_concat(A._entries, self._entries, right._entries)
     1222        return A
     1223
     1224    def stack(self, Matrix_mod2e_dense other):
     1225        """
     1226        Stack ``self`` on top of ``other``.
     1227
     1228        INPUT:
     1229
     1230        - ``other`` - a matrix
     1231
     1232        EXAMPLE::
     1233
     1234            sage: K.<a> = GF(2^4)
     1235            sage: A = random_matrix(K,2,2); A
     1236            [      a^2 + 1 a^3 + a^2 + a]
     1237            [  a^3 + a + 1         a + 1]
     1238           
     1239            sage: B = random_matrix(K,2,2); B
     1240            [a^3 + 1       0]
     1241            [    a^3       0]
     1242
     1243            sage: A.stack(B)
     1244            [      a^2 + 1 a^3 + a^2 + a]
     1245            [  a^3 + a + 1         a + 1]
     1246            [      a^3 + 1             0]
     1247            [          a^3             0]
     1248
     1249            sage: B.stack(A)
     1250            [      a^3 + 1             0]
     1251            [          a^3             0]
     1252            [      a^2 + 1 a^3 + a^2 + a]
     1253            [  a^3 + a + 1         a + 1]
     1254
     1255        TESTS::
     1256
     1257            sage: A = random_matrix(K,0,3)
     1258            sage: B = random_matrix(K,3,3)
     1259            sage: A.stack(B)
     1260            [                0     a^3 + a^2 + a     a^3 + a^2 + a]
     1261            [    a^3 + a^2 + a         a^3 + a^2       a^3 + a + 1]
     1262            [a^3 + a^2 + a + 1               a^3               a^2]
     1263
     1264            sage: B.stack(A)
     1265            [                0     a^3 + a^2 + a     a^3 + a^2 + a]
     1266            [    a^3 + a^2 + a         a^3 + a^2       a^3 + a + 1]
     1267            [a^3 + a^2 + a + 1               a^3               a^2]
     1268
     1269            sage: M = Matrix(K, 0, 0, 0)
     1270            sage: N = Matrix(K, 19, 0, 0)
     1271            sage: W = M.stack(N)
     1272            sage: W.nrows()
     1273            19
     1274            sage: M = Matrix(K, 1, 0, 0)
     1275            sage: N = Matrix(K, 1, 0, 0)
     1276            sage: M.stack(N)
     1277            []
     1278        """
     1279        if self._ncols != other._ncols:
     1280            raise TypeError, "Both numbers of columns must match."
     1281
     1282        if self._nrows == 0:
     1283            return other.__copy__()
     1284        if other._nrows == 0:
     1285            return self.__copy__()
     1286
     1287        cdef Matrix_mod2e_dense A
     1288        A = self.new_matrix(nrows = self._nrows + other._nrows)
     1289        if self._ncols == 0:
     1290            return A
     1291        A._entries = mzed_stack(A._entries, self._entries, other._entries)
     1292        return A
     1293
     1294    def submatrix(self, lowr, lowc, nrows , ncols):
     1295        """
     1296        Return submatrix from the index ``lowr,lowc`` (inclusive) with
     1297        dimension ``nrows x ncols``.
     1298
     1299        INPUT:
     1300 
     1301        - ``lowr`` -- index of start row
     1302        - ``lowc`` -- index of start column
     1303        - ``nrows`` -- number of rows of submatrix
     1304        - ``ncols`` -- number of columns of submatrix
     1305       
     1306        EXAMPLES::
     1307
     1308             sage: K.<a> = GF(2^10)
     1309             sage: A = random_matrix(K,200,200)
     1310             sage: A[0:2,0:2] == A.submatrix(0,0,2,2)
     1311             True
     1312             sage: A[0:100,0:100] == A.submatrix(0,0,100,100)
     1313             True
     1314             sage: A == A.submatrix(0,0,200,200)
     1315             True
     1316
     1317             sage: A[1:3,1:3] == A.submatrix(1,1,2,2)
     1318             True
     1319             sage: A[1:100,1:100] == A.submatrix(1,1,99,99)
     1320             True
     1321             sage: A[1:200,1:200] == A.submatrix(1,1,199,199)
     1322             True
     1323        """
     1324        cdef int highr = lowr + nrows
     1325        cdef int highc = lowc + ncols
     1326
     1327        if nrows <= 0 or ncols <= 0:
     1328            raise TypeError("Expected nrows, ncols to be > 0, but got %d,%d instead."%(nrows, ncols))
     1329 
     1330        if highc > self._entries.ncols:
     1331            raise TypeError("Expected highc <= self.ncols(), but got %d > %d instead."%(highc, self._entries.ncols))
     1332
     1333        if highr > self._entries.nrows:
     1334            raise TypeError("Expected highr <= self.nrows(), but got %d > %d instead."%(highr, self._entries.nrows))
     1335
     1336        if lowr < 0:
     1337            raise TypeError("Expected lowr >= 0, but got %d instead."%lowr)
     1338
     1339        if lowc < 0:
     1340            raise TypeError("Expected lowc >= 0, but got %d instead."%lowc)
     1341
     1342        cdef Matrix_mod2e_dense A = self.new_matrix(nrows = nrows, ncols = ncols)
     1343        if self._ncols == 0 or self._nrows == 0:
     1344            return A
     1345        A._entries = mzed_submatrix(A._entries, self._entries, lowr, lowc, highr, highc)
     1346        return A
     1347
     1348    def rank(self):
     1349        """
     1350        Return the rank of this matrix (cached).
     1351
     1352        EXAMPLE::
     1353
     1354            sage: K.<a> = GF(2^4)
     1355            sage: A = random_matrix(K, 1000, 1000)
     1356            sage: A.rank()
     1357            1000
     1358
     1359            sage: A = matrix(K, 10, 0)
     1360            sage: A.rank()
     1361            0
     1362        """
     1363        x = self.fetch('rank')
     1364        if not x is None:
     1365            return x
     1366        if self._nrows == 0 or self._ncols == 0:
     1367            return 0
     1368        cdef mzed_t *A = mzed_copy(NULL, self._entries)
     1369
     1370        cdef size_t r = mzed_echelonize(A, 0)
     1371        mzed_free(A)
     1372        self.cache('rank', r)
     1373        return r
     1374
     1375    def __hash__(self):
     1376        """
     1377        EXAMPLE::
     1378       
     1379            sage: K.<a> = GF(2^4)
     1380            sage: A = random_matrix(K, 1000, 1000)
     1381            sage: A.set_immutable()
     1382            sage: {A:1} #indirect doctest
     1383            {1000 x 1000 dense matrix over Finite Field in a of size 2^4 (type 'print A.str()' to see all of the entries): 1}
     1384           
     1385        """
     1386        return self._hash()
     1387
     1388    def __reduce__(self):
     1389        """
     1390        EXAMPLE::
     1391
     1392            sage: K.<a> = GF(2^8)
     1393            sage: A = random_matrix(K,70,70)
     1394            sage: f, s= A.__reduce__()
     1395            sage: from sage.matrix.matrix_mod2e_dense import unpickle_matrix_mod2e_dense_v0
     1396            sage: f == unpickle_matrix_mod2e_dense_v0
     1397            True
     1398            sage: f(*s) == A
     1399            True
     1400        """
     1401        from sage.matrix.matrix_space import MatrixSpace
     1402
     1403        cdef Matrix_mod2_dense A
     1404        MS = MatrixSpace(GF(2), self._entries.x.nrows, self._entries.x.ncols)
     1405        A = Matrix_mod2_dense.__new__(Matrix_mod2_dense, MS, 0, 0, 0, alloc = False)
     1406        A._entries = mzd_copy( NULL, self._entries.x)
     1407        return unpickle_matrix_mod2e_dense_v0, (A, self.base_ring(), self.nrows(), self.ncols())
     1408
     1409    def slice(self):
     1410        r"""
     1411        Unpack this matrix into matrices over `\GF{2}`.
     1412
     1413        Elements in `\GF{2^e}` can be represented as `\sum c_i a^i`
     1414        where `a` is a root the minimal polynomial. This function
     1415        returns a tuple of matrices `C` whose entry `C_i[x,y]` is the
     1416        coefficient of `c_i` in `A[x,y]` if this matrix is `A`.
     1417
     1418        EXAMPLE::
     1419       
     1420            sage: K.<a> = GF(2^2)
     1421            sage: A = random_matrix(K, 5, 5); A
     1422            [    0     1     1     0     0]
     1423            [    0     a a + 1     1     a]
     1424            [    1     1 a + 1     a a + 1]
     1425            [    0     1     1     1 a + 1]
     1426            [    1     0     a     1     1]
     1427
     1428            sage: A1,A0 = A.slice()
     1429            sage: A0
     1430            [0 0 0 0 0]
     1431            [0 1 1 0 1]
     1432            [0 0 1 1 1]
     1433            [0 0 0 0 1]
     1434            [0 0 1 0 0]
     1435
     1436            sage: A1
     1437            [0 1 1 0 0]
     1438            [0 0 1 1 0]
     1439            [1 1 1 0 1]
     1440            [0 1 1 1 1]
     1441            [1 0 0 1 1]
     1442
     1443            sage: A0[2,4]*a + A1[2,4], A[2,4]
     1444            (a + 1, a + 1)
     1445
     1446            sage: K.<a> = GF(2^3)
     1447            sage: A = random_matrix(K, 5, 5); A
     1448            [    a^2 + 1           0       a + 1           0           a]
     1449            [          a       a + 1           a           1         a^2]
     1450            [    a^2 + a           0           0     a^2 + 1       a + 1]
     1451            [    a^2 + 1           0     a^2 + 1       a + 1           1]
     1452            [a^2 + a + 1           1 a^2 + a + 1           0     a^2 + 1]
     1453
     1454            sage: A0,A1,A2 = A.slice()
     1455            sage: A0
     1456            [1 0 1 0 0]
     1457            [0 1 0 1 0]
     1458            [0 0 0 1 1]
     1459            [1 0 1 1 1]
     1460            [1 1 1 0 1]
     1461
     1462        Slicing and clinging are inverse operations::
     1463
     1464            sage: B = matrix(K, 5, 5)
     1465            sage: B.cling(A0,A1,A2)
     1466            sage: B == A
     1467            True
     1468        """
     1469        if self._entries.finite_field.degree > 4:
     1470            raise NotImplementedError("Slicing is only implemented for degree <= 4.")
     1471
     1472        from sage.matrix.matrix_space import MatrixSpace
     1473
     1474        MS = MatrixSpace(GF(2), self._nrows, self._ncols)
     1475        cdef mzd_slice_t *a = mzed_slice(NULL, self._entries)
     1476
     1477        cdef Matrix_mod2_dense A0, A1, A2, A3
     1478        A0 = Matrix_mod2_dense.__new__(Matrix_mod2_dense, MS, 0, 0, 0, alloc = True)
     1479        A1 = Matrix_mod2_dense.__new__(Matrix_mod2_dense, MS, 0, 0, 0, alloc = True)
     1480        mzd_copy(A0._entries, a.x[0])
     1481        mzd_copy(A1._entries, a.x[1])
     1482        if self._entries.finite_field.degree > 2:
     1483            A2 = Matrix_mod2_dense.__new__(Matrix_mod2_dense, MS, 0, 0, 0, alloc = True)
     1484            mzd_copy(A2._entries, a.x[2])
     1485        if self._entries.finite_field.degree > 3:
     1486            A3 = Matrix_mod2_dense.__new__(Matrix_mod2_dense, MS, 0, 0, 0, alloc = True)
     1487            mzd_copy(A3._entries, a.x[3])
     1488
     1489        mzd_slice_free(a)
     1490        if self._entries.finite_field.degree == 2:
     1491            return A0,A1
     1492        elif self._entries.finite_field.degree == 3:
     1493            return A0,A1,A2
     1494        elif self._entries.finite_field.degree == 4:
     1495            return A0,A1,A2,A3
     1496
     1497    def cling(self, *C):
     1498        r"""
     1499        Pack the matrices over `\GF{2}` into this matrix over `\GF{2^e}`.
     1500
     1501        Elements in `\GF{2^e}` can be represented as `\sum c_i a^i` where
     1502        `a` is a root the minimal polynomial. If this matrix is `A`
     1503        then this function writes `c_i a^i` to the entry `A[x,y]`
     1504        where `c_i` is the entry `C_i[x,y]`.
     1505
     1506        INPUT:
     1507
     1508        - ``C`` - a list of matrices over GF(2)
     1509
     1510        EXAMPLE::
     1511       
     1512            sage: K.<a> = GF(2^2)
     1513            sage: A = matrix(K, 5, 5)
     1514            sage: A0 = random_matrix(GF(2), 5, 5); A0
     1515            [0 1 0 1 1]
     1516            [0 1 1 1 0]
     1517            [0 0 0 1 0]
     1518            [0 1 1 0 0]
     1519            [0 0 0 1 1]
     1520
     1521            sage: A1 = random_matrix(GF(2), 5, 5); A1
     1522            [0 0 1 1 1]
     1523            [1 1 1 1 0]
     1524            [0 0 0 1 1]
     1525            [1 0 0 0 1]
     1526            [1 0 0 1 1]
     1527
     1528            sage: A.cling(A1, A0); A
     1529            [    0     a     1 a + 1 a + 1]
     1530            [    1 a + 1 a + 1 a + 1     0]
     1531            [    0     0     0 a + 1     1]
     1532            [    1     a     a     0     1]
     1533            [    1     0     0 a + 1 a + 1]
     1534
     1535            sage: A0[0,3]*a + A1[0,3], A[0,3]
     1536            (a + 1, a + 1)
     1537
     1538        Slicing and clinging are inverse operations::
     1539
     1540            sage: B1, B0 = A.slice()
     1541            sage: B0 == A0 and B1 == A1
     1542            True
     1543
     1544        TESTS::
     1545
     1546            sage: K.<a> = GF(2^2)
     1547            sage: A = matrix(K, 5, 5)
     1548            sage: A0 = random_matrix(GF(2), 5, 5)
     1549            sage: A1 = random_matrix(GF(2), 5, 5)
     1550            sage: A.cling(A0, A1)
     1551            sage: B = copy(A)
     1552            sage: A.cling(A0, A1)
     1553            sage: A == B
     1554            True
     1555
     1556            sage: A.cling(A0)
     1557            Traceback (most recent call last):
     1558            ...
     1559            ValueError: The number of input matrices must be equal to the degree of the base field.
     1560       
     1561            sage: K.<a> = GF(2^5)
     1562            sage: A = matrix(K, 5, 5)
     1563            sage: A0 = random_matrix(GF(2), 5, 5)
     1564            sage: A1 = random_matrix(GF(2), 5, 5)
     1565            sage: A2 = random_matrix(GF(2), 5, 5)
     1566            sage: A3 = random_matrix(GF(2), 5, 5)
     1567            sage: A4 = random_matrix(GF(2), 5, 5)
     1568            sage: A.cling(A0, A1, A2, A3, A4)
     1569            Traceback (most recent call last):
     1570            ...
     1571            NotImplementedError: Cling is only implemented for degree <= 4.
     1572        """
     1573        cdef Py_ssize_t i
     1574
     1575        if self._entries.finite_field.degree > 4:
     1576            raise NotImplementedError("Cling is only implemented for degree <= 4.")
     1577
     1578        if self._mutability._is_immutable:
     1579            raise TypeError("Mutable matrices cannot be modified.")
     1580
     1581        if len(C) != self._entries.finite_field.degree:
     1582            raise ValueError("The number of input matrices must be equal to the degree of the base field.")
     1583
     1584        cdef mzd_slice_t *v = mzd_slice_init(self._entries.finite_field, self._nrows, self._ncols)
     1585        for i in range(self._entries.finite_field.degree):
     1586            if not PY_TYPE_CHECK(C[i], Matrix_mod2_dense):
     1587                mzd_slice_free(v)
     1588                raise TypeError("All input matrices must be over GF(2).")
     1589            mzd_copy(v.x[i], (<Matrix_mod2_dense>C[i])._entries)
     1590        mzed_set_ui(self._entries, 0)
     1591        mzed_cling(self._entries, v)
     1592        mzd_slice_free(v)
     1593
     1594def unpickle_matrix_mod2e_dense_v0(Matrix_mod2_dense a, base_ring, nrows, ncols):
     1595    """
     1596    EXAMPLE::
     1597
     1598        sage: K.<a> = GF(2^2)
     1599        sage: A = random_matrix(K,10,10)
     1600        sage: f, s= A.__reduce__()
     1601        sage: from sage.matrix.matrix_mod2e_dense import unpickle_matrix_mod2e_dense_v0
     1602        sage: f == unpickle_matrix_mod2e_dense_v0
     1603        True
     1604        sage: f(*s) == A
     1605        True
     1606    """
     1607    from sage.matrix.matrix_space import MatrixSpace
     1608
     1609    MS = MatrixSpace(base_ring, nrows, ncols)
     1610    cdef Matrix_mod2e_dense A  = Matrix_mod2e_dense.__new__(Matrix_mod2e_dense, MS, 0, 0, 0)
     1611    mzd_copy(A._entries.x, a._entries)
     1612    return A

  • sage/matrix/matrix_space.py 

    diff --git a/sage/matrix/matrix_space.py b/sage/matrix/matrix_space.py
    a b  
    3838import matrix_modn_sparse
    3939
    4040import matrix_mod2_dense
    41 #import matrix_mod2_sparse
     41import matrix_mod2e_dense
    4242
    4343import matrix_integer_dense
    4444import matrix_integer_sparse
    … …  
    906906                # elif R.order() < matrix_modn_dense.MAX_MODULUS:
    907907                #     return matrix_modn_dense.Matrix_modn_dense
    908908                return matrix_generic_dense.Matrix_generic_dense
     909            elif sage.rings.finite_rings.all.is_FiniteField(R) and R.characteristic() == 2 and R.order() <= 1024:
     910                return matrix_mod2e_dense.Matrix_mod2e_dense
    909911            elif sage.rings.polynomial.multi_polynomial_ring_generic.is_MPolynomialRing(R) and R.base_ring().is_field():
    910912                return matrix_mpolynomial_dense.Matrix_mpolynomial_dense
    911913            #elif isinstance(R, sage.rings.padics.padic_ring_capped_relative.pAdicRingCappedRelative):

  • sage/rings/finite_rings/element_givaro.pxd 

    diff --git a/sage/rings/finite_rings/element_givaro.pxd b/sage/rings/finite_rings/element_givaro.pxd
    a b  
    7878    cpdef int characteristic(self)
    7979    cpdef FiniteField_givaroElement gen(self)
    8080    cpdef FiniteField_givaroElement element_from_data(self, e)
     81    cdef FiniteField_givaroElement _new_c(self, int value)
    8182
    8283cdef class FiniteField_givaro_iterator:
    8384    cdef int iterator

  • sage/rings/finite_rings/element_givaro.pyx 

    diff --git a/sage/rings/finite_rings/element_givaro.pyx b/sage/rings/finite_rings/element_givaro.pyx
    a b  
    766766            rep = 'int'
    767767        return unpickle_Cache_givaro, (self.parent, p, k, self.parent.polynomial(), rep, self._has_array)
    768768
     769    cdef FiniteField_givaroElement _new_c(self, int value):
     770        return make_FiniteField_givaroElement(self, value)
     771
    769772
    770773def unpickle_Cache_givaro(parent, p, k, modulus, rep, cache):
    771774    """

  • sage/rings/polynomial/multi_polynomial_sequence.py 

    diff --git a/sage/rings/polynomial/multi_polynomial_sequence.py b/sage/rings/polynomial/multi_polynomial_sequence.py
    a b  
    130130    sage: A.rank()
    131131    4056
    132132    sage: A[4055]*v
    133     (k001*k003)
    134 
    135 The last output tells us that ``k001`` and ``k003`` can't both be 1.
     133    (k001*k002 + k001*k003)
    136134
    137135TEST::
    138136
    … …  
    439437            sage: P = F.ring()
    440438            sage: I = F.ideal()
    441439            sage: I.elimination_ideal(P('s000*s001*s002*s003*w100*w101*w102*w103*x100*x101*x102*x103'))
    442             Ideal (k002 + (a^2)*k003 + 1,
    443                    k001 + (a^3)*k003 + (a + 1),
    444                    k000 + (a^3 + a^2 + a)*k003 + (a^3 + a^2 + a),
    445                    k103 + (a^3)*k003 + (a^2 + 1),
    446                    k102 + (a^3 + a^2 + a)*k003 + (a^3 + 1),
    447                    k101 + k003 + (a^2 + a),
    448                    k100 + (a^2)*k003 + (a^2 + a),
    449                    k003^2 + (a^3 + a^2 + a)*k003 + (a^3 + a^2 + a))
     440            Ideal (k002 + (a^3 + a^2 + 1)*k003 + (a^3),
     441                   k001 + (a^3 + a^2 + a + 1)*k003 + (a^3 + a^2 + a),
     442                   k000 + (a + 1)*k003 + (a^2 + a + 1),
     443                   k103 + (a^3 + a)*k003 + (a^3 + a),
     444                   k102 + (a^2 + a + 1)*k003 + (a^3 + a^2 + a),
     445                   k101 + (a^3)*k003 + (a^3),
     446                   k100 + (a^3 + a + 1)*k003 + (a^2 + 1),
     447                   k003^2 + (a + 1)*k003 + (a^2 + a + 1))
    450448            of Multivariate Polynomial Ring in
    451449            k100, k101, k102, k103, x100, x101, x102, x103, w100, w101, w102, w103,
    452450            s000, s001, s002, s003, k000, k001, k002, k003 over Finite Field in a of size 2^4

  • sage/rings/polynomial/pbori.pyx 

    diff --git a/sage/rings/polynomial/pbori.pyx b/sage/rings/polynomial/pbori.pyx
    a b  
    44904490            sage: F,s = sr.polynomial_system()
    44914491            sage: I = F.ideal()
    44924492            sage: I.interreduced_basis()
    4493             [k100 + k003 + 1,
    4494              k101 + k003,
    4495              k102,
    4496              k103 + k003,
    4497              x100 + 1,
    4498              x101 + 1,
    4499              x102 + 1,
    4500              x103 + k003,
    4501              w100 + k003,
    4502              w101,
    4503              w102 + k003 + 1,
    4504              w103 + k003 + 1,
    4505              s000 + 1,
    4506              s001 + 1,
    4507              s002 + 1,
    4508              s003 + k003 + 1,
    4509              k000 + k003 + 1,
    4510              k001,
    4511              k002 + k003]
     4493            [k100 + k003,
     4494            k101 + k003,
     4495            k102 + k003 + 1,
     4496            k103, x100 + k003,
     4497            x101, x102 + 1,
     4498            x103, w100 + 1,
     4499            w101 + k003, w102,
     4500            w103 + k003 + 1,
     4501            s000 + k003 + 1,
     4502            s001, s002 + 1,
     4503            s003,
     4504            k000 + 1,
     4505            k001 + k003 + 1,
     4506            k002]
    45124507        """
    45134508        R = self.ring()
    45144509        set_cring(R)