Ticket #11574: m4ri_20110601.patch

module_list.py

@@ -886 +886 @@
     # TODO -- change to use BLAS at some point.
     Extension('sage.matrix.matrix_integer_dense',
               sources = ['sage/matrix/matrix_integer_dense.pyx'],
+              extra_compile_args = ['-std=c99'],
               # order matters for cygwin!!
               libraries = ['iml', 'pari', 'm', 'gmp', BLAS, BLAS2]),
 
@@ -896 +897 @@
     Extension('sage.matrix.matrix_mod2_dense',
               sources = ['sage/matrix/matrix_mod2_dense.pyx'],
               libraries = ['gmp','m4ri', 'gd', 'png12', 'z'],
+              extra_compile_args = ['-std=c99'],
               depends = [SAGE_INC + "png.h", SAGE_INC + "m4ri/m4ri.h"]),
 
     Extension('sage.matrix.matrix_modn_dense',
@@ -1097 +1099 @@
     Extension('sage.modules.vector_mod2_dense',
               sources = ['sage/modules/vector_mod2_dense.pyx'],
               libraries = ['gmp','m4ri', 'png12', 'gd'],
+              extra_compile_args = ['-std=c99'],
               depends = [SAGE_INC + "png.h", SAGE_INC + "m4ri/m4ri.h"]),
     
     Extension('sage.modules.vector_rational_dense',

sage/libs/m4ri.pxd

@@ -1 +1 @@
 cdef extern from "m4ri/m4ri.h":
+    ctypedef int rci_t
+    ctypedef int wi_t
     ctypedef unsigned long long m4ri_word "word"
+    ctypedef int BIT
 
     ctypedef struct mzd_t:
-        int nrows
-        int ncols
-        int width
+        rci_t nrows
+        rci_t ncols
+        wi_t width
+        int offset
         m4ri_word **rows
 
     ctypedef struct mzp_t:
-        int *values
-        int size
+        rci_t *values
+        rci_t size
 
-    ctypedef int BIT
-
-    cdef int RADIX
+    cdef int m4ri_radix
 
     ##############
     # Maintainance
@@ -28 +30 @@
     ##############
   
     # create empty matrix
-    cdef mzd_t *mzd_init(int,int)
+    cdef mzd_t *mzd_init(rci_t , rci_t)
 
     # create the identity permutation
-    cdef mzp_t *mzp_init(long)
+    cdef mzp_t *mzp_init(rci_t)
     
     # free memory for the matrix
     cdef void mzd_free(mzd_t *)
@@ -43 +45 @@
     cdef void mzd_randomize(mzd_t *)
 
     # identity matrix if i%2
-    cdef void mzd_set_ui(mzd_t *, unsigned int i)
+    cdef void mzd_set_ui(mzd_t *, unsigned int )
 
     # [A],[B] -> [AB]
-    cdef mzd_t *mzd_concat(mzd_t *C, mzd_t *A, mzd_t *B)
+    cdef mzd_t *mzd_concat(mzd_t *, mzd_t *, mzd_t *)
 
     # [A],[B] -> | A |
     #            | B |
-    cdef mzd_t *mzd_stack(mzd_t *C, mzd_t *A, mzd_t *B)
+    cdef mzd_t *mzd_stack(mzd_t *, mzd_t *, mzd_t *)
     
     # returns a submatrix from a
-    cdef mzd_t *mzd_submatrix(mzd_t *S, mzd_t *A, int lowr, int lowc, int highr, int highc)
+    cdef mzd_t *mzd_submatrix(mzd_t *, mzd_t *, rci_t lowr, rci_t lowc, rci_t highr, rci_t highc)
 
     # return a matrix window to A
-    cdef mzd_t *mzd_init_window( mzd_t *A, int lowr, int lowc, int highr, int highc)
+    cdef mzd_t *mzd_init_window(mzd_t *, rci_t lowr, rci_t lowc, rci_t highr, rci_t highc)
 
     cdef void mzd_free_window(mzd_t *)
 
     # deep copy
     cdef mzd_t *mzd_copy(mzd_t *, mzd_t *)
 
+    # printing
+    cdef void mzd_print(mzd_t *)
+
     ##############
     # Bit Level IO
     ##############
 
     # set BIT
-    cdef void mzd_write_bit( mzd_t *m, int row, int col, BIT value)
+    cdef void mzd_write_bit( mzd_t *m, rci_t row, rci_t col, BIT value)
 
     # get BIT
-    cdef BIT mzd_read_bit( mzd_t *m, int row, int col )
+    cdef BIT mzd_read_bit( mzd_t *m, rci_t row, rci_t col )
 
     # get BITs (n<=64)
-    cdef m4ri_word mzd_read_bits( mzd_t *m, int row, int col, int n)
+    cdef m4ri_word mzd_read_bits( mzd_t *m, rci_t row, rci_t col, int n)
     
     #####################
     # Row/Column Based IO
     #####################
     
-    cdef void mzd_row_swap(mzd_t *, int, int)
+    cdef void mzd_row_swap(mzd_t *, rci_t, rci_t)
 
-    cdef void mzd_col_swap(mzd_t *, int, int)
+    cdef void mzd_col_swap(mzd_t *, rci_t, rci_t)
 
-    cdef void mzd_row_clear_offset(mzd_t *m, int, int)
+    cdef void mzd_row_clear_offset(mzd_t *m, rci_t, rci_t)
 
-    cdef void mzd_row_add_offset(mzd_t *m, int, int, int)
+    cdef void mzd_row_add_offset(mzd_t *m, rci_t, rci_t, rci_t)
 
     ############
     # Arithmetic
@@ -171 +176 @@
     # reduced row echelon form using PLUQ factorization
     cdef mzd_t *mzd_kernel_left_pluq(mzd_t *A, int cutoff)
 
+    ########################
+    # Bit operations
+    ########################
+
+    cdef m4ri_word __M4RI_LEFT_BITMASK(int)
+
+    cdef m4ri_word m4ri_swap_bits(m4ri_word)
+
     ##################################
     # Internal Functions for Debugging
     ##################################
 
-    cdef void mzd_write_zeroed_bits(mzd_t *m, int x, int y, int n, m4ri_word values)
-     
     cdef void mzd_clear_bits(mzd_t *m, int x, int y, int n)
 

sage/matrix/matrix_mod2_dense.pyx

@@ -350 +350 @@
         cdef unsigned long i, j, truerow
         cdef unsigned long start, shift
         cdef m4ri_word row_xor
-        cdef m4ri_word end_mask = ~(((<m4ri_word>1)<<(RADIX - self._ncols%RADIX))-1)
+        cdef m4ri_word end_mask = __M4RI_LEFT_BITMASK(self._ncols%m4ri_radix)
         cdef m4ri_word top_mask, bot_mask
         cdef m4ri_word cur
         cdef m4ri_word* row
@@ -370 +370 @@
             row = self._entries.rows[i]
             start = (i*self._entries.ncols) >> 6
             shift = (i*self._entries.ncols) & 0x3F
-            bot_mask = (~(<m4ri_word>0)) << shift
+            bot_mask = __M4RI_LEFT_BITMASK(m4ri_radix - shift)
             top_mask = ~bot_mask
 
             if self._entries.width > 1:
@@ -379 +379 @@
 
                 for j from 1 <= j < self._entries.width - 1:
                     row_xor ^= row[j]
-                    cur = ((row[j-1] << (63-shift)) << 1) ^ (row[j] >> shift)
+                    cur = ((row[j-1] >> (63-shift)) >> 1) ^ (row[j] << shift)
                     running_parity ^= (start+j) & parity_mask(cur)
 
                 running_parity ^= (start+j) & parity_mask(row[j-1] & top_mask)
@@ -393 +393 @@
             running_parity ^= (start+j) & parity_mask(cur & bot_mask)
             running_parity ^= (start+j+1) & parity_mask(cur & top_mask)
 
-            start = (i*self._entries.ncols) & (RADIX-1)
-            running_xor ^= (row_xor >> start) ^ ((row_xor << (RADIX-1-start)) << 1)
+            running_xor ^= (row_xor << shift) ^ ((row_xor >> (63-shift)) >> 1)
 
         cdef unsigned long bit_is_set
         cdef unsigned long h
 
         # Now we assemble the running_parity and running_xor to get the hash. 
         # Viewing the flattened matrix as a list of a_i, the hash is the xor
-        # of the i for which a_i is non-zero. We split i into the lower RADIX 
-        # bits and the rest, so i = i1 << RADIX + i0. Now two matching i0 
+        # of the i for which a_i is non-zero. We split i into the lower m4ri_radix 
+        # bits and the rest, so i = i1 << m4ri_radix + i0. Now two matching i0 
         # would cancel, so we only need the parity of how many of each 
         # possible i0 occur. This is stored in the bits of running_xor.
         # Similarly, running_parity is the xor of the i1 needed. It's called
@@ -410 +409 @@
         # the number of i1 to add is equal to the number of set bits in that 
         # word (but because two cancel, we only need keep track of the 
         # parity. 
-        h = RADIX * running_parity
-        for i from 0 <= i < RADIX:
-            bit_is_set = (running_xor >> (RADIX-1-i)) & 1
-            h ^= (RADIX-1) & ~(bit_is_set-1) & i
+
+        h = m4ri_radix * running_parity
+        for i from 0 <= i < m4ri_radix:
+            bit_is_set = (running_xor >> i) & 1
+            h ^= (m4ri_radix-1) & ~(bit_is_set-1) & i
 
         if h == -1:
             h = -2
@@ -665 +665 @@
         cdef Vector_mod2_dense c = PY_NEW(Vector_mod2_dense)
         c._init(self._nrows, VS)
         c._entries = mzd_init(1, self._nrows)
-        tmp = mzd_init(self._nrows, 1)
-        _mzd_mul_naive(tmp, self._entries, (<Vector_mod2_dense>v)._entries, 0)
-        mzd_transpose(c._entries, tmp)
-        mzd_free(tmp)
+        if c._entries.nrows and c._entries.ncols:
+            tmp = mzd_init(self._nrows, 1)
+            _mzd_mul_naive(tmp, self._entries, (<Vector_mod2_dense>v)._entries, 0)
+            mzd_transpose(c._entries, tmp)
+            mzd_free(tmp)
         return c
         
     cdef Matrix _matrix_times_matrix_(self, Matrix right):
@@ -1266 +1267 @@
         cdef int i, j, k
         cdef int nc
         cdef unsigned int low, high
-        cdef m4ri_word mask = 1
+        cdef m4ri_word mask = 0
         
         # Original code, before adding the ``nonzero`` option.
         if not nonzero:
             if density == 1:
                 assert(sizeof(m4ri_word) == 8)
-                mask = ~((mask<<(RADIX - self._entries.ncols%RADIX))-1)
+                mask = __M4RI_LEFT_BITMASK(self._entries.ncols % m4ri_radix)
                 for i from 0 <= i < self._nrows:
                     for j from 0 <= j < self._entries.width:
                         # for portability we get 32-bit twice rather than 64-bit once
                         low = gmp_urandomb_ui(rstate.gmp_state, 32)
                         high = gmp_urandomb_ui(rstate.gmp_state, 32)
-                        self._entries.rows[i][j] = ((<unsigned long long>high)<<32) | (<unsigned long long>low)
+                        self._entries.rows[i][j] = m4ri_swap_bits( ((<unsigned long long>high)<<32) | (<unsigned long long>low) )
                     self._entries.rows[i][self._entries.width - 1] &= mask
             else:
                 nc = self._ncols