# HG changeset patch
# User Martin Albrecht <martinralbrecht@googlemail.com>
# Date 1306771873 3600
# Node ID 632260cc66e03df6c20be829b28bee185cae9c66
# Parent a48af224839e4eb2555175d15d703ecd81d0b4c1
#11574: updated interface for M4RI version 20110601 and up.
diff git a/module_list.py b/module_list.py
a

b


886  886  # TODO  change to use BLAS at some point. 
887  887  Extension('sage.matrix.matrix_integer_dense', 
888  888  sources = ['sage/matrix/matrix_integer_dense.pyx'], 
 889  extra_compile_args = ['std=c99'], 
889  890  # order matters for cygwin!! 
890  891  libraries = ['iml', 'pari', 'm', 'gmp', BLAS, BLAS2]), 
891  892  
… 
… 

896  897  Extension('sage.matrix.matrix_mod2_dense', 
897  898  sources = ['sage/matrix/matrix_mod2_dense.pyx'], 
898  899  libraries = ['gmp','m4ri', 'gd', 'png12', 'z'], 
 900  extra_compile_args = ['std=c99'], 
899  901  depends = [SAGE_INC + "png.h", SAGE_INC + "m4ri/m4ri.h"]), 
900  902  
901  903  Extension('sage.matrix.matrix_modn_dense', 
… 
… 

1097  1099  Extension('sage.modules.vector_mod2_dense', 
1098  1100  sources = ['sage/modules/vector_mod2_dense.pyx'], 
1099  1101  libraries = ['gmp','m4ri', 'png12', 'gd'], 
 1102  extra_compile_args = ['std=c99'], 
1100  1103  depends = [SAGE_INC + "png.h", SAGE_INC + "m4ri/m4ri.h"]), 
1101  1104  
1102  1105  Extension('sage.modules.vector_rational_dense', 
diff git a/sage/libs/m4ri.pxd b/sage/libs/m4ri.pxd
a

b


1  1  cdef extern from "m4ri/m4ri.h": 
 2  ctypedef int rci_t 
 3  ctypedef int wi_t 
2  4  ctypedef unsigned long long m4ri_word "word" 
 5  ctypedef int BIT 
3  6  
4  7  ctypedef struct mzd_t: 
5   int nrows 
6   int ncols 
7   int width 
 8  rci_t nrows 
 9  rci_t ncols 
 10  wi_t width 
 11  int offset 
8  12  m4ri_word **rows 
9  13  
10  14  ctypedef struct mzp_t: 
11   int *values 
12   int size 
 15  rci_t *values 
 16  rci_t size 
13  17  
14   ctypedef int BIT 
15   
16   cdef int RADIX 
 18  cdef int m4ri_radix 
17  19  
18  20  ############## 
19  21  # Maintainance 
… 
… 

28  30  ############## 
29  31  
30  32  # create empty matrix 
31   cdef mzd_t *mzd_init(int,int) 
 33  cdef mzd_t *mzd_init(rci_t , rci_t) 
32  34  
33  35  # create the identity permutation 
34   cdef mzp_t *mzp_init(long) 
 36  cdef mzp_t *mzp_init(rci_t) 
35  37  
36  38  # free memory for the matrix 
37  39  cdef void mzd_free(mzd_t *) 
… 
… 

43  45  cdef void mzd_randomize(mzd_t *) 
44  46  
45  47  # identity matrix if i%2 
46   cdef void mzd_set_ui(mzd_t *, unsigned int i) 
 48  cdef void mzd_set_ui(mzd_t *, unsigned int ) 
47  49  
48  50  # [A],[B] > [AB] 
49   cdef mzd_t *mzd_concat(mzd_t *C, mzd_t *A, mzd_t *B) 
 51  cdef mzd_t *mzd_concat(mzd_t *, mzd_t *, mzd_t *) 
50  52  
51  53  # [A],[B] >  A  
52  54  #  B  
53   cdef mzd_t *mzd_stack(mzd_t *C, mzd_t *A, mzd_t *B) 
 55  cdef mzd_t *mzd_stack(mzd_t *, mzd_t *, mzd_t *) 
54  56  
55  57  # returns a submatrix from a 
56   cdef mzd_t *mzd_submatrix(mzd_t *S, mzd_t *A, int lowr, int lowc, int highr, int highc) 
 58  cdef mzd_t *mzd_submatrix(mzd_t *, mzd_t *, rci_t lowr, rci_t lowc, rci_t highr, rci_t highc) 
57  59  
58  60  # return a matrix window to A 
59   cdef mzd_t *mzd_init_window( mzd_t *A, int lowr, int lowc, int highr, int highc) 
 61  cdef mzd_t *mzd_init_window(mzd_t *, rci_t lowr, rci_t lowc, rci_t highr, rci_t highc) 
60  62  
61  63  cdef void mzd_free_window(mzd_t *) 
62  64  
63  65  # deep copy 
64  66  cdef mzd_t *mzd_copy(mzd_t *, mzd_t *) 
65  67  
 68  # printing 
 69  cdef void mzd_print(mzd_t *) 
 70  
66  71  ############## 
67  72  # Bit Level IO 
68  73  ############## 
69  74  
70  75  # set BIT 
71   cdef void mzd_write_bit( mzd_t *m, int row, int col, BIT value) 
 76  cdef void mzd_write_bit( mzd_t *m, rci_t row, rci_t col, BIT value) 
72  77  
73  78  # get BIT 
74   cdef BIT mzd_read_bit( mzd_t *m, int row, int col ) 
 79  cdef BIT mzd_read_bit( mzd_t *m, rci_t row, rci_t col ) 
75  80  
76  81  # get BITs (n<=64) 
77   cdef m4ri_word mzd_read_bits( mzd_t *m, int row, int col, int n) 
 82  cdef m4ri_word mzd_read_bits( mzd_t *m, rci_t row, rci_t col, int n) 
78  83  
79  84  ##################### 
80  85  # Row/Column Based IO 
81  86  ##################### 
82  87  
83   cdef void mzd_row_swap(mzd_t *, int, int) 
 88  cdef void mzd_row_swap(mzd_t *, rci_t, rci_t) 
84  89  
85   cdef void mzd_col_swap(mzd_t *, int, int) 
 90  cdef void mzd_col_swap(mzd_t *, rci_t, rci_t) 
86  91  
87   cdef void mzd_row_clear_offset(mzd_t *m, int, int) 
 92  cdef void mzd_row_clear_offset(mzd_t *m, rci_t, rci_t) 
88  93  
89   cdef void mzd_row_add_offset(mzd_t *m, int, int, int) 
 94  cdef void mzd_row_add_offset(mzd_t *m, rci_t, rci_t, rci_t) 
90  95  
91  96  ############ 
92  97  # Arithmetic 
… 
… 

171  176  # reduced row echelon form using PLUQ factorization 
172  177  cdef mzd_t *mzd_kernel_left_pluq(mzd_t *A, int cutoff) 
173  178  
 179  ######################## 
 180  # Bit operations 
 181  ######################## 
 182  
 183  cdef m4ri_word __M4RI_LEFT_BITMASK(int) 
 184  
 185  cdef m4ri_word m4ri_swap_bits(m4ri_word) 
 186  
174  187  ################################## 
175  188  # Internal Functions for Debugging 
176  189  ################################## 
177  190  
178   cdef void mzd_write_zeroed_bits(mzd_t *m, int x, int y, int n, m4ri_word values) 
179   
180  191  cdef void mzd_clear_bits(mzd_t *m, int x, int y, int n) 
181  192  
diff git a/sage/matrix/matrix_mod2_dense.pyx b/sage/matrix/matrix_mod2_dense.pyx
a

b


350  350  cdef unsigned long i, j, truerow 
351  351  cdef unsigned long start, shift 
352  352  cdef m4ri_word row_xor 
353   cdef m4ri_word end_mask = ~(((<m4ri_word>1)<<(RADIX  self._ncols%RADIX))1) 
 353  cdef m4ri_word end_mask = __M4RI_LEFT_BITMASK(self._ncols%m4ri_radix) 
354  354  cdef m4ri_word top_mask, bot_mask 
355  355  cdef m4ri_word cur 
356  356  cdef m4ri_word* row 
… 
… 

370  370  row = self._entries.rows[i] 
371  371  start = (i*self._entries.ncols) >> 6 
372  372  shift = (i*self._entries.ncols) & 0x3F 
373   bot_mask = (~(<m4ri_word>0)) << shift 
 373  bot_mask = __M4RI_LEFT_BITMASK(m4ri_radix  shift) 
374  374  top_mask = ~bot_mask 
375  375  
376  376  if self._entries.width > 1: 
… 
… 

379  379  
380  380  for j from 1 <= j < self._entries.width  1: 
381  381  row_xor ^= row[j] 
382   cur = ((row[j1] << (63shift)) << 1) ^ (row[j] >> shift) 
 382  cur = ((row[j1] >> (63shift)) >> 1) ^ (row[j] << shift) 
383  383  running_parity ^= (start+j) & parity_mask(cur) 
384  384  
385  385  running_parity ^= (start+j) & parity_mask(row[j1] & top_mask) 
… 
… 

393  393  running_parity ^= (start+j) & parity_mask(cur & bot_mask) 
394  394  running_parity ^= (start+j+1) & parity_mask(cur & top_mask) 
395  395  
396   start = (i*self._entries.ncols) & (RADIX1) 
397   running_xor ^= (row_xor >> start) ^ ((row_xor << (RADIX1start)) << 1) 
 396  running_xor ^= (row_xor << shift) ^ ((row_xor >> (63shift)) >> 1) 
398  397  
399  398  cdef unsigned long bit_is_set 
400  399  cdef unsigned long h 
401  400  
402  401  # Now we assemble the running_parity and running_xor to get the hash. 
403  402  # Viewing the flattened matrix as a list of a_i, the hash is the xor 
404   # of the i for which a_i is nonzero. We split i into the lower RADIX 
405   # bits and the rest, so i = i1 << RADIX + i0. Now two matching i0 
 403  # of the i for which a_i is nonzero. We split i into the lower m4ri_radix 
 404  # bits and the rest, so i = i1 << m4ri_radix + i0. Now two matching i0 
406  405  # would cancel, so we only need the parity of how many of each 
407  406  # possible i0 occur. This is stored in the bits of running_xor. 
408  407  # 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 
411  410  # word (but because two cancel, we only need keep track of the 
412  411  # parity. 
413   h = RADIX * running_parity 
414   for i from 0 <= i < RADIX: 
415   bit_is_set = (running_xor >> (RADIX1i)) & 1 
416   h ^= (RADIX1) & ~(bit_is_set1) & i 
 412  
 413  h = m4ri_radix * running_parity 
 414  for i from 0 <= i < m4ri_radix: 
 415  bit_is_set = (running_xor >> i) & 1 
 416  h ^= (m4ri_radix1) & ~(bit_is_set1) & i 
417  417  
418  418  if h == 1: 
419  419  h = 2 
… 
… 

665  665  cdef Vector_mod2_dense c = PY_NEW(Vector_mod2_dense) 
666  666  c._init(self._nrows, VS) 
667  667  c._entries = mzd_init(1, self._nrows) 
668   tmp = mzd_init(self._nrows, 1) 
669   _mzd_mul_naive(tmp, self._entries, (<Vector_mod2_dense>v)._entries, 0) 
670   mzd_transpose(c._entries, tmp) 
671   mzd_free(tmp) 
 668  if c._entries.nrows and c._entries.ncols: 
 669  tmp = mzd_init(self._nrows, 1) 
 670  _mzd_mul_naive(tmp, self._entries, (<Vector_mod2_dense>v)._entries, 0) 
 671  mzd_transpose(c._entries, tmp) 
 672  mzd_free(tmp) 
672  673  return c 
673  674  
674  675  cdef Matrix _matrix_times_matrix_(self, Matrix right): 
… 
… 

1266  1267  cdef int i, j, k 
1267  1268  cdef int nc 
1268  1269  cdef unsigned int low, high 
1269   cdef m4ri_word mask = 1 
 1270  cdef m4ri_word mask = 0 
1270  1271  
1271  1272  # Original code, before adding the ``nonzero`` option. 
1272  1273  if not nonzero: 
1273  1274  if density == 1: 
1274  1275  assert(sizeof(m4ri_word) == 8) 
1275   mask = ~((mask<<(RADIX  self._entries.ncols%RADIX))1) 
 1276  mask = __M4RI_LEFT_BITMASK(self._entries.ncols % m4ri_radix) 
1276  1277  for i from 0 <= i < self._nrows: 
1277  1278  for j from 0 <= j < self._entries.width: 
1278  1279  # for portability we get 32bit twice rather than 64bit once 
1279  1280  low = gmp_urandomb_ui(rstate.gmp_state, 32) 
1280  1281  high = gmp_urandomb_ui(rstate.gmp_state, 32) 
1281   self._entries.rows[i][j] = ((<unsigned long long>high)<<32)  (<unsigned long long>low) 
 1282  self._entries.rows[i][j] = m4ri_swap_bits( ((<unsigned long long>high)<<32)  (<unsigned long long>low) ) 
1282  1283  self._entries.rows[i][self._entries.width  1] &= mask 
1283  1284  else: 
1284  1285  nc = self._ncols 