trac #17665: Merged with 6.5.beta6

src/sage/data_structures/binary_matrix.pxd

@@ -0,0 +1,8 @@
+from sage.data_structures.bitset cimport bitset_t
+
+cdef struct binary_matrix_s:
+    long n_cols
+    long n_rows
+    bitset_t* rows
+
+ctypedef binary_matrix_s[1] binary_matrix_t

src/sage/misc/binary_matrix.pxi → src/sage/data_structures/binary_matrix.pxi

@@ -13,16 +13,14 @@ a ``binary_matrix_t`` structure contains :
 
 - ``long width`` -- number of ``unsigned long`` per row
 
-- ``unsigned long ** rows`` -- ``rows[i]`` points toward the ``width`` blocks of
-  type ``unsigned long`` containing the bits of row `i`.
+- ``bitset_t * rows`` -- ``rows[i]`` points toward a block of type ``bitset_t``
+  containing the bits of row `i`.
 
-.. NOTE::
-
-    The rows are stored contiguously in memory, i.e. ``row[i] = row[i-1] +
-    width``.
 """
-include "binary_matrix_pxd.pxi"
-include "sage/ext/stdsage.pxi"
+
+from sage.data_structures.binary_matrix cimport *
+include 'sage/data_structures/bitset.pxi'
+
 
 cdef inline binary_matrix_init(binary_matrix_t m, long n_rows, long n_cols):
     r"""
@@ -32,82 +30,78 @@ cdef inline binary_matrix_init(binary_matrix_t m, long n_rows, long n_cols):
 
     m.n_cols = n_cols
     m.n_rows = n_rows
-    m.width = (n_cols - 1)/(8*sizeof(unsigned long)) + 1
-    m.rows = <unsigned long**>sage_malloc(n_rows * sizeof(unsigned long *))
+    m.rows = <bitset_t *>sage_malloc(n_rows * sizeof(bitset_t))
     if m.rows == NULL:
         raise MemoryError
 
-    m.rows[0] = <unsigned long*>sage_malloc(n_rows * m.width * sizeof(unsigned long))
-    if m.rows[0] == NULL:
-        sage_free(m.rows)
-        raise MemoryError
-
-    for i in range(1,n_rows):
-        m.rows[i] = m.rows[i-1] + m.width
-        m.rows[i][m.width-1] = 0
+    for i from 0 <= i < n_rows:
+        bitset_init(m.rows[i], n_cols)
 
 cdef inline binary_matrix_free(binary_matrix_t m):
     r"""
     Frees the memory allocated by the matrix
     """
-    sage_free(m.rows[0])
+    cdef int i
+
+    for i from 0 <= i < m.n_rows:
+        bitset_free(m.rows[i])
     sage_free(m.rows)
 
 cdef inline binary_matrix_fill(binary_matrix_t m, bint bit):
     r"""
     Fill the whole matrix with a bit
     """
-    memset(m.rows[0],-(<unsigned char> bit),m.width * m.n_rows * sizeof(unsigned long))
+    cdef int i
+
+    if bit: # set the matrix to 1
+        for i from 0 <= i < m.n_rows:
+            bitset_set_first_n(m.rows[i], m.n_cols)
+    else:
+        for i from 0 <= i < m.n_rows:
+            bitset_clear(m.rows[i])
 
 cdef inline binary_matrix_complement(binary_matrix_t m):
     r"""
     Complements all of the matrix' bits.
     """
-    cdef int i,j
-    for i in range(m.n_rows):
-        for j in range(m.width):
-            m.rows[i][j] = ~m.rows[i][j]
-
-        # Set the "useless" bits to 0 too.
-        if ((m.n_cols)  & offset_mask):
-            m.rows[i][m.width-1] &= (<unsigned long> 1<<((m.n_cols)  & offset_mask))-1
+    cdef int i
+    for i from 0 <= i < m.n_rows:
+        bitset_complement(m.rows[i], m.rows[i])
 
 cdef inline binary_matrix_set1(binary_matrix_t m, long row, long col):
     r"""
     Sets an entry to 1
     """
-    m.rows[row][col >> index_shift] |= (<unsigned long>1) << (col & offset_mask)
+    bitset_add(m.rows[row], col)
 
 cdef inline binary_matrix_set0(binary_matrix_t m, long row, long col):
     r"""
     Sets an entry to 0
     """
-    m.rows[row][col >> index_shift] &= ~((<unsigned long>1) << (col & offset_mask))
+    bitset_discard(m.rows[row], col)
 
 cdef inline binary_matrix_set(binary_matrix_t m, long row, long col, bint value):
     r"""
     Sets an entry
     """
-    binary_matrix_set0(m,row,col)
-    m.rows[row][col >> index_shift] |= (<unsigned long>1) << (value & offset_mask)
+    if value:
+        binary_matrix_set1(m,row,col)
+    else:
+        binary_matrix_set0(m,row,col)
 
 cdef inline bint binary_matrix_get(binary_matrix_t m, long row, long col):
     r"""
     Returns the value of a given entry
     """
-    return (m.rows[row][col >> index_shift] >> (col & offset_mask)) & 1
+    return bitset_in(m.rows[row], col)
 
 cdef inline binary_matrix_print(binary_matrix_t m):
     r"""
     Prints the binary matrix
     """
     cdef int i,j
     import sys
-    for i in range(m.n_rows):
-        # If you want to print the *whole* matrix, including the useless bits,
-        # use the following line instead
-        #
-        # for j in (m.width*8*sizeof(unsigned long)):
-        for j in range(m.n_cols):
+    for i from 0 <= i < m.n_rows:
+        for j from 0 <= j < m.n_cols:
             sys.stdout.write("1" if binary_matrix_get(m, i, j) else ".",)
         print ""

src/sage/data_structures/bitset.pxi

@@ -289,6 +289,20 @@ cdef inline bint bitset_issuperset(bitset_t a, bitset_t b):
     """
     return bitset_issubset(b, a)
 
+
+cdef inline bint bitset_are_disjoint(bitset_t a, bitset_t b):
+    """
+    Tests whether ``a`` and ``b`` have an empty intersection.
+
+    We assume ``a.limbs == b.limbs``.
+    """
+    cdef mp_size_t i
+    for i from 0 <= i < a.limbs:
+        if (a.bits[i]&b.bits[i]) != 0:
+            return False
+    return True
+
+
 #############################################################################
 # Bitset Bit Manipulation
 #############################################################################

src/sage/graphs/base/static_dense_graph.pxd

@@ -1,3 +1,4 @@
-include "sage/misc/binary_matrix.pxi"
+from sage.data_structures.binary_matrix cimport binary_matrix_t
+
 
 cdef dict dense_graph_init(binary_matrix_t m, g, translation = ?)

src/sage/graphs/base/static_dense_graph.pyx

@@ -40,6 +40,7 @@ Index
 Functions
 ---------
 """
+include "sage/data_structures/binary_matrix.pxi"
 
 cdef dict dense_graph_init(binary_matrix_t m, g, translation=False):
     r"""
@@ -59,25 +60,25 @@ cdef dict dense_graph_init(binary_matrix_t m, g, translation=False):
     cdef int is_undirected = isinstance(g, Graph)
     cdef int n = g.order()
 
-    binary_matrix_init(m,n,n)
-    binary_matrix_fill(m,0)
+    binary_matrix_init(m, n, n)
+    binary_matrix_fill(m, 0)
 
     # If the vertices are 0...n-1, let's avoid an unnecessary dictionary
     if g.vertices() == range(n):
         if translation:
             d_translation = {i:i for i in range(n)}
 
         for i,j in g.edge_iterator(labels = False):
-            binary_matrix_set1(m,i,j)
+            binary_matrix_set1(m, i, j)
             if is_undirected:
-                binary_matrix_set1(m,j,i)
+                binary_matrix_set1(m, j, i)
     else:
         d_translation = {v:i for i,v in enumerate(g.vertices())}
 
         for u,v in g.edge_iterator(labels = False):
-            binary_matrix_set1(m,d_translation[u],d_translation[v])
+            binary_matrix_set1(m, d_translation[u], d_translation[v])
             if is_undirected:
-                binary_matrix_set1(m,d_translation[v],d_translation[u])
+                binary_matrix_set1(m, d_translation[v], d_translation[u])
 
     if translation:
         return d_translation
@@ -171,6 +172,7 @@ def is_strongly_regular(g, parameters = False):
     """
     g._scream_if_not_simple()
     cdef binary_matrix_t m
+    cdef bitset_t b_tmp
     cdef int n = g.order()
     cdef int inter
     cdef int i,j,l, k
@@ -186,6 +188,8 @@ def is_strongly_regular(g, parameters = False):
     if not all(d == k for d in degree):
         return False
 
+    bitset_init(b_tmp, n)
+
     # m i now our copy of the graph
     dense_graph_init(m, g)
 
@@ -197,25 +201,27 @@ def is_strongly_regular(g, parameters = False):
 
             # The intersection of the common neighbors of i and j is a AND of
             # their respective rows. A popcount then returns its cardinality.
-            inter = 0
-            for l in range(m.width):
-                inter += __builtin_popcountl(m.rows[i][l] & m.rows[j][l])
+            bitset_and(b_tmp, m.rows[i], m.rows[j])
+            inter = bitset_len(b_tmp)
 
             # Check that this cardinality is correct according to the values of lambda and mu
             if binary_matrix_get(m,i,j):
                 if llambda == -1:
                     llambda = inter
                 elif llambda != inter:
                     binary_matrix_free(m)
+                    bitset_free(b_tmp)
                     return False
             else:
                 if mu == -1:
                     mu = inter
                 elif mu != inter:
                     binary_matrix_free(m)
+                    bitset_free(b_tmp)
                     return False
 
     binary_matrix_free(m)
+    bitset_free(b_tmp)
 
     if parameters:
         return (n,k,llambda,mu)

src/sage/graphs/distances_all_pairs.pyx

@@ -143,8 +143,7 @@ Functions
 #                  http://www.gnu.org/licenses/
 ##############################################################################
 
-include "sage/data_structures/bitset.pxi"
-include "sage/misc/binary_matrix.pxi"
+include "sage/data_structures/binary_matrix.pxi"
 from libc.stdint cimport uint64_t, uint32_t, INT32_MAX, UINT32_MAX
 from sage.graphs.base.c_graph cimport CGraph
 from sage.graphs.base.c_graph cimport vertex_label
@@ -538,22 +537,26 @@ def is_distance_regular(G, parameters = False):
         if distance_matrix[i] > diameter and distance_matrix[i] != infinity:
             diameter = distance_matrix[i]
 
+    cdef bitset_t b_tmp
+    bitset_init(b_tmp, n)
+
     # b_distance_matrix[d*n+v] is the set of vertices at distance d from v.
     cdef binary_matrix_t b_distance_matrix
     try:
         binary_matrix_init(b_distance_matrix,n*(diameter+2),n)
         binary_matrix_fill(b_distance_matrix, 0)
     except MemoryError:
         sage_free(distance_matrix)
+        bitset_free(b_tmp)
         raise
 
     # Fills b_distance_matrix
     for u in range(n):
         for v in range(u,n):
             d = distance_matrix[u*n+v]
             if d != infinity:
-                binary_matrix_set1(b_distance_matrix, d*n+u,v)
-                binary_matrix_set1(b_distance_matrix, d*n+v,u)
+                binary_matrix_set1(b_distance_matrix, d*n+u, v)
+                binary_matrix_set1(b_distance_matrix, d*n+v, u)
 
     cdef list bi = [-1 for i in range(diameter +1)]
     cdef list ci = [-1 for i in range(diameter +1)]
@@ -570,11 +573,10 @@ def is_distance_regular(G, parameters = False):
 
             # Computations of b_d and c_d for u,v. We intersect sets stored in
             # b_distance_matrix.
-            b = 0
-            c = 0
-            for l in range(b_distance_matrix.width):
-                b += __builtin_popcountl(b_distance_matrix.rows[(d+1)*n+u][l] & b_distance_matrix.rows[1*n+v][l])
-                c += __builtin_popcountl(b_distance_matrix.rows[(d-1)*n+u][l] & b_distance_matrix.rows[1*n+v][l])
+            bitset_and(b_tmp, b_distance_matrix.rows[(d+1)*n+u], b_distance_matrix.rows[n+v])
+            b = bitset_len(b_tmp)
+            bitset_and(b_tmp, b_distance_matrix.rows[(d-1)*n+u], b_distance_matrix.rows[n+v])
+            c = bitset_len(b_tmp)
 
             # Consistency of b_d and c_d
             if bi[d] == -1:
@@ -584,10 +586,12 @@ def is_distance_regular(G, parameters = False):
             elif bi[d] != b or ci[d] != c:
                 sage_free(distance_matrix)
                 binary_matrix_free(b_distance_matrix)
+                bitset_free(b_tmp)
                 return False
 
     sage_free(distance_matrix)
     binary_matrix_free(b_distance_matrix)
+    bitset_free(b_tmp)
 
     if parameters:
         bi[0] = k

src/sage/graphs/generic_graph_pyx.pyx

@@ -18,7 +18,7 @@ AUTHORS:
 include "sage/ext/interrupt.pxi"
 include 'sage/ext/cdefs.pxi'
 include 'sage/ext/stdsage.pxi'
-include "sage/misc/binary_matrix.pxi"
+include "sage/data_structures/binary_matrix.pxi"
 
 # import from Python standard library
 from sage.misc.prandom import random
@@ -1321,8 +1321,7 @@ def transitive_reduction_acyclic(G):
         for vv in G.neighbors_out(uu):
             v = v_to_int[vv]
             binary_matrix_set1(closure, u, v)
-            for i in range(closure.width):
-                closure.rows[u][i] |= closure.rows[v][i]
+            bitset_or(closure.rows[u], closure.rows[u], closure.rows[v])
 
     # Build the transitive reduction of G
     #
@@ -1334,8 +1333,7 @@ def transitive_reduction_acyclic(G):
         u = v_to_int[uu]
         for vv in G.neighbors_out(uu):
             v = v_to_int[vv]
-            for i in range(closure.width):
-                closure.rows[u][i] &= ~closure.rows[v][i]
+            bitset_difference(closure.rows[u], closure.rows[u], closure.rows[v])
         for vv in G.neighbors_out(uu):
             v = v_to_int[vv]
             if binary_matrix_get(closure, u, v):

src/sage/graphs/independent_sets.pxd

@@ -1,5 +1,6 @@
 from libc.stdint cimport uint32_t, uint64_t
-include "sage/misc/binary_matrix.pxi"
+from sage.data_structures.binary_matrix cimport *
+
 
 cdef class IndependentSets:
     cdef binary_matrix_t g