Ticket #13352: trac_13352_bitset_len.patch

File trac_13352_bitset_len.patch, 10.8 KB (added by dcoudert, 9 years ago)
  • new file sage/misc/bitcount.pxi

    # HG changeset patch
    # User dcoudert <david.coudert@inria.fr>
    # Date 1344456406 -7200
    # Node ID c284c063b3da19e4b3cb535995d4412dc62dad3e
    # Parent  a3a5bdf1e8c8c9f17f141150bf07043c8115ec4b
    trac #13352 -- Running time improvement of the bitset_len method
    
    diff --git a/sage/misc/bitcount.pxi b/sage/misc/bitcount.pxi
    new file mode 100644
    - +  
     1"""
     2Fast methods for counting bits.
     3
     4Implements some well known methods for counting the number of 1-bits in a 32 and
     564-bits integer. Do your own test to know which is the most appropriate method
     6(i.e., the fastest one) for you.
     7
     8:meth:`parallel_bitcount_32 <parallel_bitcount_32>`
     9:meth:`parallel_bitcount_64 <parallel_bitcount_64>`
     10
     11:meth:`popcount_32 <popcount_32>`
     12:meth:`popcount_64 <popcount_64>`
     13"""
     14
     15#*****************************************************************************
     16#     Copyright (C) 2012 David Coudert <david.coudert@inria.fr>
     17#
     18#  Distributed under the terms of the GNU General Public License (GPL)
     19#
     20#    This code is distributed in the hope that it will be useful,
     21#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     22#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     23#    General Public License for more details.
     24#
     25#  The full text of the GPL is available at:
     26#
     27#                  http://www.gnu.org/licenses/
     28#*****************************************************************************
     29
     30# This is a .pxi file so that one can inline functions. Doctests in misc_c.pyx
     31
     32
     33from libc.stdint cimport uint32_t, uint64_t, UINT32_MAX, UINT64_MAX
     34
     35
     36##################
     37# Parallel Count #
     38##################
     39
     40cdef inline uint32_t __TWO32__(uint32_t c):
     41    return 0x1u << c
     42cdef inline uint32_t __MASK32__(uint32_t c):
     43    return UINT32_MAX / (__TWO32__(__TWO32__(c)) + 0x1u)
     44cdef inline uint32_t __COUNT32__(uint32_t x, uint32_t c):
     45    return ((x) & __MASK32__(c)) + (((x) >> (__TWO32__(c))) & __MASK32__(c))
     46
     47cdef inline uint64_t __TWO64__(uint64_t c):
     48    return 0x1u << c
     49cdef inline uint64_t __MASK64__(uint64_t c):
     50    return UINT64_MAX / (__TWO64__(__TWO64__(c)) + 0x1u)
     51cdef inline uint64_t __COUNT64__(uint64_t x, uint64_t c):
     52    return ((x) & __MASK64__(c)) + (((x) >> (__TWO64__(c))) & __MASK64__(c))
     53
     54cdef inline int parallel_bitcount_32(uint32_t n):
     55    """
     56    Returns the number of '1' bits in a 32-bits integer.
     57
     58    If sizeof(int) > 4, this function only returns the number of '1' bits
     59    in (i & ((1<<32) - 1)).
     60    """
     61    n = __COUNT32__(n, 0)
     62    n = __COUNT32__(n, 1)
     63    n = __COUNT32__(n, 2)
     64    n = __COUNT32__(n, 3)
     65    n = __COUNT32__(n, 4)
     66    return n
     67
     68cdef inline int parallel_bitcount_64(uint64_t n):
     69    """
     70    Returns the number of '1' bits in a 64-bits integer.
     71    """
     72    n = __COUNT64__(n, 0)
     73    n = __COUNT64__(n, 1)
     74    n = __COUNT64__(n, 2)
     75    n = __COUNT64__(n, 3)
     76    n = __COUNT64__(n, 4)
     77    n = __COUNT64__(n, 5)
     78    return n
     79
     80
     81##################
     82# MIT HAKMEM 169 #
     83##################
     84
     85cdef inline int popcount_32(uint32_t i):
     86    """
     87    Returns the number of '1' bits in a 32-bits integer.
     88
     89    If sizeof(int) > 4, this function only returns the number of '1' bits
     90    in (i & ((1<<32) - 1)).
     91    """
     92    i = i - ((i >> 1) & 0x55555555)
     93    i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
     94    return ((i + (i >> 4) & 0x0F0F0F0F) * 0x01010101) >> 24
     95
     96
     97cdef inline int popcount_64(uint64_t i):
     98    """
     99    Returns the number of '1' bits in a 64-bits integer.
     100    """
     101    i = i - ((i >> 1) & 0x5555555555555555ULL)
     102    i = (i & 0x3333333333333333ULL) + ((i >> 2) & 0x3333333333333333ULL)
     103    return ( ((i + (i >> 4)) & 0x0f0f0f0f0f0f0f0fULL) * 0x0101010101010101ULL ) >> 56
     104
  • sage/misc/bitset.pxi

    diff --git a/sage/misc/bitset.pxi b/sage/misc/bitset.pxi
    a b  
    402402    return -1
    403403
    404404
     405include 'bitcount.pxi'
    405406
    406407cdef inline long bitset_len(bitset_t bits):
    407408    """
    408409    Calculate the number of items in the set (i.e., the number of nonzero bits).
    409410    """
    410411    cdef long len = 0
    411     cdef long i = bitset_first(bits)
    412     while i>=0:
    413         len += 1
    414         i=bitset_next(bits, i+1)
     412    cdef long i
     413    for i from 0 <= i < bits.limbs:
     414        len += popcount_64(bits.bits[i])
    415415    return len
    416416
    417417cdef inline long bitset_hash(bitset_t bits):
  • sage/misc/misc_c.pyx

    diff --git a/sage/misc/misc_c.pyx b/sage/misc/misc_c.pyx
    a b  
    399399        k = (1+count) >> 1
    400400        return balanced_list_sum(L, offset, k, cutoff) + balanced_list_sum(L, offset+k, count-k, cutoff)
    401401
    402        
     402
    403403#############################################################################
    404404# Bitset Testing
    405405#############################################################################
     
    812812    bitset_free(a)
    813813
    814814
     815#############################################################################
     816# Testing bitcount.pxi
     817#############################################################################
     818
     819# We don't include it here since it is also included in bitset.pxi
     820# include "bitcount.pxi"
     821
     822def test_bitcount():
     823    """
     824    Test the Cython bit counting functions so we can have some relevant
     825    doctests.
     826
     827    TESTS:
     828
     829        sage: from sage.misc.misc_c import test_bitcount
     830        sage: test_bitcount()
     831        Compare all methods with random numbers in range 1..2^32 -1...
     832       
     833        Compare methods for 64-bits numbers with random numbers in range 1..2^64 -1...
     834       
     835        Methods for numbers on 32-bits are not able to count all bits of a 64-bits number
     836        Number 2^64-2^32= 18446744069414584320  has 32 bits set to 1
     837        parallel_bitcount_32: 0
     838        parallel_bitcount_64: 32
     839        popcount_32: 0
     840        popcount_64: 32
     841       
     842        Counting the number of 1 bits of the Integer 1234567890
     843        parallel_bitcount_32: 12
     844        parallel_bitcount_64: 12
     845        popcount_32: 12
     846        popcount_64: 12
     847       
     848        Counting the number of bits of the unsigned long long 12345678901234567890
     849        parallel_bitcount_32: 17
     850        parallel_bitcount_64: 32
     851        popcount_32: 17
     852        popcount_64: 32
     853       
     854        All methods fail with 2^80 = 1208925819614629174706176
     855        parallel_bitcount_32: 0
     856        parallel_bitcount_64: 0
     857        popcount_32: 0
     858        popcount_64: 0
     859    """
     860    from sage.rings.integer import Integer
     861    from sage.misc.prandom import randint
     862
     863    cdef uint64_t j, a,b,c,d
     864    cdef uint32_t i
     865    cdef unsigned long long xull
     866
     867    print 'Compare all methods with random numbers in range 1..2^32 -1...'
     868    for 0 <= i < 1000:
     869        j = randint(0,1<<32 -1)
     870        a = parallel_bitcount_32(j)
     871        b = parallel_bitcount_64(j)
     872        c = popcount_32(j)
     873        d = popcount_64(j)
     874        if a != b or a != c or a != d:
     875            print 'Error with int ',j,' = ',Integer(j).binary()
     876            print 'parallel_bitcount_32:',a
     877            print 'parallel_bitcount_64:',b
     878            print 'popcount_32:',c
     879            print 'popcount_64:',d
     880            break
     881
     882    print '\nCompare methods for 64-bits numbers with random numbers in range 1..2^64 -1...'
     883    for 0 <= i < 1000:
     884        j = randint(0,1<<64 -1)
     885        b = parallel_bitcount_64(j)
     886        d = popcount_64(j)
     887        if b != d:
     888            print 'Error with int ',j,' = ',Integer(j).binary()
     889            print 'parallel_bitcount_64:',b
     890            print 'popcount_64:',d
     891            break
     892
     893    print '\nMethods for numbers on 32-bits are not able to count all bits of a 64-bits number'
     894    j = <uint64_t>-1 - <uint32_t>-1
     895    print 'Number 2^64-2^32=',j,' has 32 bits set to 1'
     896    print 'parallel_bitcount_32:',parallel_bitcount_32(j)
     897    print 'parallel_bitcount_64:',parallel_bitcount_64(j)
     898    print 'popcount_32:',popcount_32(j)
     899    print 'popcount_64:',popcount_64(j)
     900
     901    print '\nCounting the number of 1 bits of the Integer 1234567890'
     902    print 'parallel_bitcount_32:',parallel_bitcount_32(Integer(1234567890))
     903    print 'parallel_bitcount_64:',parallel_bitcount_64(Integer(1234567890))
     904    print 'popcount_32:',popcount_32(Integer(1234567890))
     905    print 'popcount_64:',popcount_64(Integer(1234567890))
     906
     907    print '\nCounting the number of bits of the unsigned long long 12345678901234567890'
     908    xull =  12345678901234567890
     909    print 'parallel_bitcount_32:',parallel_bitcount_32(xull)
     910    print 'parallel_bitcount_64:',parallel_bitcount_64(xull)
     911    print 'popcount_32:',popcount_32(xull)
     912    print 'popcount_64:',popcount_64(xull)
     913
     914    print '\nAll methods fail with 2^80 =',2**80
     915    print 'parallel_bitcount_32:',parallel_bitcount_32(2**80)
     916    print 'parallel_bitcount_64:',parallel_bitcount_64(2**80)
     917    print 'popcount_32:',popcount_32(2**80)
     918    print 'popcount_64:',popcount_64(2**80)
     919
     920
     921def test_parallel_bitcount_32(sample=1000):
     922    """
     923    Returns the sum of the number of bits of the integers in 0..sample.
     924
     925    This method allows to evaluate the running time of the
     926    ``parallel_bitcount_32`` method.
     927
     928    EXAMPLE:
     929
     930        sage: from sage.misc.misc_c import test_parallel_bitcount_32
     931        sage: test_parallel_bitcount_32(1000)
     932        4932
     933    """
     934    cdef uint32_t i
     935    sum=0
     936    for 0 <= i < sample:
     937        sum += parallel_bitcount_32(i)
     938    return sum
     939
     940
     941def test_parallel_bitcount_64(sample=1000):
     942    """
     943    Returns the sum of the number of bits of the integers in 0..sample.
     944
     945    This method allows to evaluate the running time of the
     946    ``parallel_bitcount_64`` method.
     947
     948    EXAMPLE:
     949
     950        sage: from sage.misc.misc_c import test_parallel_bitcount_64
     951        sage: test_parallel_bitcount_64(1000)
     952        4932
     953    """
     954    cdef uint64_t i
     955    sum = 0
     956    for 0 <= i < sample:
     957        sum += parallel_bitcount_64(i)
     958    return sum
     959
     960
     961def test_popcount_32(sample=1000):
     962    """
     963    Returns the sum of the number of bits of the integers in 0..sample.
     964
     965    This method allows to evaluate the running time of the ``popcount_32``
     966    method.
     967
     968    EXAMPLE:
     969
     970        sage: from sage.misc.misc_c import test_popcount_32
     971        sage: test_popcount_32(1000)
     972        4932
     973    """
     974    cdef uint32_t i
     975    sum = 0
     976    for 0 <= i < sample:
     977        sum += popcount_32(i)
     978    return sum
     979
     980
     981def test_popcount_64(sample=1000):
     982    """
     983    Returns the sum of the number of bits of the integers in 0..sample.
     984
     985    This method allows to evaluate the running time of the ``popcount_64``
     986    method.
     987
     988    EXAMPLE:
     989
     990        sage: from sage.misc.misc_c import test_popcount_64
     991        sage: test_popcount_64(1000)
     992        4932
     993    """
     994    cdef uint64_t i
     995    sum = 0
     996    for 0 <= i < sample:
     997        sum += popcount_64(i)
     998    return sum
     999
     1000
    8151001#################################################################
    8161002# 32/64-bit computer?
    8171003#################################################################