Use mpn functions for bitsets

src/sage/crypto/boolean_function.pyx

@@ -105,7 +105,7 @@ cdef long yellow_code(unsigned long a):
         m ^= (m<<s)
     return r
 
-cdef reed_muller(unsigned long *f, int ldn):
+cdef reed_muller(mp_limb_t* f, int ldn):
     r"""
     The Reed Muller transform (also known as binary Moebius transform)
     is an orthogonal transform. For a function `f` defined by

src/sage/graphs/independent_sets.pyx

@@ -19,12 +19,12 @@ include "sage/misc/bitset.pxi"
 from sage.misc.cachefunc import cached_method
 from sage.graphs.base.static_dense_graph cimport dense_graph_init
 
-cdef inline int bitset_are_disjoint(unsigned long * b1, unsigned long * b2, int width):
+cdef inline bint bitset_are_disjoint(mp_limb_t* b1, mp_limb_t* b2, mp_size_t width):
     r"""
     Tests whether two bitsets of length width*sizeof(unsigned int) have an empty
     intersection.
     """
-    cdef int i
+    cdef mp_size_t i
     for i in range(width):
         if b1[i]&b2[i]:
             return False
@@ -33,7 +33,7 @@ cdef inline int bitset_are_disjoint(unsigned long * b1, unsigned long * b2, int
 cdef inline int ismaximal(binary_matrix_t g, int n, bitset_t s):
     cdef int i
     for i in range(n):
-        if (not bitset_in(s,i)) and bitset_are_disjoint(g.rows[i],s.bits,g.width):
+        if (not bitset_in(s,i)) and bitset_are_disjoint(<mp_limb_t*>(g.rows[i]), s.bits, g.width):
             return False
 
     return True
@@ -245,7 +245,7 @@ cdef class IndependentSets:
             if bitset_in(current_set,i):
 
                 # We have found an independent set !
-                if bitset_are_disjoint(self.g.rows[i],current_set.bits,self.g.width):
+                if bitset_are_disjoint(<mp_limb_t*>(self.g.rows[i]), current_set.bits, self.g.width):
 
                     # Saving that set
                     bitset_copy(tmp, current_set)
@@ -391,7 +391,7 @@ cdef class IndependentSets:
             bitset_add(s, i)
 
             # Checking that the set s is independent
-            if not bitset_are_disjoint(self.g.rows[i],s.bits,self.g.width):
+            if not bitset_are_disjoint(<mp_limb_t*>(self.g.rows[i]), s.bits, self.g.width):
                 return False
 
         if self.maximal and not ismaximal(self.g, self.n,s):

src/sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi

@@ -783,8 +783,8 @@ cdef StabilizerChain *SC_new(int n, bint init_gens=True):
     SC.gen_is_id.size  = default_num_bits
     SC.gen_used.limbs  = limbs
     SC.gen_is_id.limbs = limbs
-    SC.gen_used.bits   = <unsigned long*>sage_malloc(limbs * sizeof(unsigned long))
-    SC.gen_is_id.bits  = <unsigned long*>sage_malloc(limbs * sizeof(unsigned long))
+    SC.gen_used.bits   = <mp_limb_t*>sage_malloc(limbs * sizeof(mp_limb_t))
+    SC.gen_is_id.bits  = <mp_limb_t*>sage_malloc(limbs * sizeof(mp_limb_t))
 
     # check for allocation failures
     if int_array        is NULL or int_ptrs          is NULL or \
@@ -945,12 +945,12 @@ cdef int SC_realloc_bitsets(StabilizerChain *SC, long size):
         new_size *= 2
     cdef unsigned long limbs_old = SC.gen_used.limbs
     cdef long limbs = (new_size - 1)/(8*sizeof(unsigned long)) + 1
-    cdef unsigned long *tmp = <unsigned long *> sage_realloc(SC.gen_used.bits, limbs * sizeof(unsigned long))
+    cdef mp_limb_t *tmp = <mp_limb_t*> sage_realloc(SC.gen_used.bits, limbs * sizeof(mp_limb_t))
     if tmp is not NULL:
         SC.gen_used.bits = tmp
     else:
         return 1
-    tmp = <unsigned long *> sage_realloc(SC.gen_is_id.bits, limbs * sizeof(unsigned long))
+    tmp = <mp_limb_t*> sage_realloc(SC.gen_is_id.bits, limbs * sizeof(mp_limb_t))
     if tmp is not NULL:
         SC.gen_is_id.bits = tmp
     else:
@@ -959,9 +959,9 @@ cdef int SC_realloc_bitsets(StabilizerChain *SC, long size):
     SC.gen_is_id.limbs = limbs
     SC.gen_used.size = new_size
     SC.gen_is_id.size = new_size
-    SC.gen_used.bits[size_old >> index_shift] &= ((<unsigned long>1 << (size_old & offset_mask)) - 1)
+    SC.gen_used.bits[size_old >> index_shift] &= limb_lower_bits_down(size_old)
     memset(SC.gen_used.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long))
-    SC.gen_is_id.bits[size_old >> index_shift] &= ((<unsigned long>1 << (size_old & offset_mask)) - 1)
+    SC.gen_is_id.bits[size_old >> index_shift] &= limb_lower_bits_down(size_old)
     memset(SC.gen_is_id.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long))
     return 0
 

src/sage/libs/gmp/mpn.pxd

@@ -40,3 +40,15 @@ cdef extern from "gmp.h":
     unsigned long int mpn_popcount (mp_limb_t *s1p, mp_size_t n)
     unsigned long int mpn_hamdist (mp_limb_t *s1p, mp_limb_t *s2p, mp_size_t n)
     int mpn_perfect_square_p (mp_limb_t *s1p, mp_size_t n)
+    void mpn_and_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_andn_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_nand_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_ior_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_iorn_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_nior_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_xor_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_xnor_n(mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t n)
+    void mpn_com(mp_ptr rp, mp_srcptr sp, mp_size_t n)
+    void mpn_copyi(mp_ptr rp, mp_srcptr s1p, mp_size_t n)
+    void mpn_copyd(mp_ptr rp, mp_srcptr s1p, mp_size_t n)
+    void mpn_zero(mp_ptr rp, mp_size_t n)

src/sage/libs/gmp/types.pxd

@@ -2,7 +2,7 @@ from libc.stdio cimport FILE
 
 cdef extern from "gmp.h":
     # GMP's configuration of how many bits are stuffed into a limb
-    cdef int GMP_LIMB_BITS
+    cdef unsigned int GMP_LIMB_BITS
 
     ### Type Declarations ###
 
@@ -14,6 +14,7 @@ cdef extern from "gmp.h":
 
     ctypedef long mp_limb_t
     ctypedef mp_limb_t* mp_ptr
+    ctypedef mp_limb_t* mp_srcptr
 
     # This internal structure is not guaranteed to stay the same with
     # future releases of GMP or MPIR.

src/sage/misc/bitset.pxd

@@ -12,27 +12,25 @@ Cython bitset types
 
 # This file declares the bitset types. The implementation of the basic
 # Cython type bitset_t (inline functions) is in bitset.pxi, the
-# implementation of the Python later is in bitset.pyx. The latter file
+# implementation of the Python class is in bitset.pyx. The latter file
 # also contains all doctests.
 
 cdef extern from *:
-    # constant literals
-    int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)"
-    unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)"
-
     # Given an element index n in a set, (n >> index_shift) gives the
-    # corresponding limb number, while (n & offset_mask) gives the bit
-    # number inside of the limb.
+    # corresponding limb number.
+    int index_shift "(sizeof(mp_limb_t) == 8 ? 6 : 5)"
+
+from sage.libs.gmp.types cimport *
 
 cdef struct bitset_s:
     # The size of a bitset B counts the maximum number of bits that B can
     # hold. This size is independent of how many elements of B are toggled to
     # 1. For example, say B is the bitset 1001. Then B has size 4, with the
     # first and fourth elements toggled to 1, reading from left to right.
     # We can also think of the size of a bitset as its capacity.
-    long size
+    mp_bitcnt_t size
 
-    # A limb is that part of a bitset that can fit into an unsigned long
+    # A limb is that part of a bitset that can fit into an mp_limb_t
     # (typically, 32 bits on a 32-bit machine and 64 bits on a 64-bit
     # machine). This counts the number of limbs to represent a bitset.
     # If a bitset has size <= n, then the whole bitset fits into a limb
@@ -42,12 +40,12 @@ cdef struct bitset_s:
     # and the bitset has size 96 bits, then we require at most two limbs
     # to represent the bitset.
     #
-    # NOTE: we assume that a limb corresponds to an MPIR limb (this
-    # assumption is always true in practice).
-    long limbs
+    # NOTE: some code assumes that mp_limb_t is an unsigned long
+    # (this assumption is always true in practice).
+    mp_size_t limbs
 
     # The individual bits of a bitset.
-    unsigned long *bits
+    mp_limb_t* bits
 
 ctypedef bitset_s bitset_t[1]