bmfunc.h

ncbi源码

}

/*!
   \brief Checks if GAP block is all-one.
   \param buf - GAP buffer pointer.
   \param set_max - max possible bitset length
   \returns true if all-one.

   @ingroup gapfunc
*/
template<typename T> 
         bool gap_is_all_one(const T* buf, unsigned set_max)
{
    return ((*buf & 1)  && (*(++buf) == set_max - 1));
}

/*!
   \brief Returs GAP block length.
   \param buf - GAP buffer pointer.
   \returns GAP block length.

   @ingroup gapfunc
*/
template<typename T> unsigned gap_length(const T* buf)
{
    return (*buf >> 3) + 1;
}


/*!
   \brief Returs GAP block capacity.
   \param buf - GAP buffer pointer.
   \returns GAP block capacity.

   @ingroup gapfunc
*/
template<typename T> 
unsigned gap_capacity(const T* buf, const T* glevel_len)
{
    return glevel_len[(*buf >> 1) & 3];
}


/*!
   \brief Returs GAP block capacity limit.
   \param buf - GAP buffer pointer.
   \param glevel_len - GAP lengths table (gap_len_table)
   \returns GAP block limit.

   @ingroup gapfunc
*/
template<typename T> 
unsigned gap_limit(const T* buf, const T* glevel_len)
{
    return glevel_len[(*buf >> 1) & 3]-4;
}


/*!
   \brief Returs GAP blocks capacity level.
   \param buf - GAP buffer pointer.
   \returns GAP block capacity level.

   @ingroup gapfunc
*/
template<typename T> unsigned gap_level(const T* buf)
{
    return (*buf >> 1) & 3;
}


/*!
   \brief Sets GAP block capacity level.
   \param buf - GAP buffer pointer.
   \param level new GAP block capacity level.

   @ingroup gapfunc
*/
template<typename T> void set_gap_level(T* buf, 
                                        unsigned level)
{
    BM_ASSERT(level < bm::gap_levels);
    *buf = ((level & 3) << 1) | (*buf & 1) | (*buf & ~7); 
}


/*!
   \brief Calculates GAP block capacity level.
   \param len - GAP buffer length.
   \param glevel_len - GAP lengths table
   \return GAP block capacity level. 
            -1 if block does not fit any level.
   @ingroup gapfunc
*/
template<typename T>
inline int gap_calc_level(int len, const T* glevel_len)
{
    if (len <= (glevel_len[0]-4)) return 0;
    if (len <= (glevel_len[1]-4)) return 1;
    if (len <= (glevel_len[2]-4)) return 2;
    if (len <= (glevel_len[3]-4)) return 3;

    BM_ASSERT(bm::gap_levels == 4);
    return -1;
}

/*! @brief Returns number of free elements in GAP block array. 
    Difference between GAP block capacity on this level and actual GAP length.
    
    @param buf - GAP buffer pointer
    @parma glevel_len - GAP lengths table
    
    @return Number of free GAP elements
    @ingroup gapfunc
*/
template<typename T>
inline unsigned gap_free_elements(const T* buf, const T* glevel_len)
{
    unsigned len = gap_length(buf);
    unsigned capacity = gap_capacity(buf, glevel_len);
    return capacity - len;
}


/*! 
   \brief Lexicographical comparison of BIT buffers.
   \param buf1 - First buffer pointer.
   \param buf2 - Second buffer pointer.
   \param len - Buffer length in elements (T).
   \return  <0 - less, =0 - equal,  >0 - greater.

   @ingroup bitfunc 
*/
template<typename T> 
int bitcmp(const T* buf1, const T* buf2, unsigned len)
{
    BM_ASSERT(len);

    const T* pend1 = buf1 + len; 
    do
    {
        T w1 = *buf1++;
        T w2 = *buf2++;
        T diff = w1 ^ w2;
    
        if (diff)
        { 
            return (w1 & diff & -diff) ? 1 : -1;
        }

    } while (buf1 < pend1);

    return 0;
}


/*! 
   \brief Converts bit block to GAP. 
   \param dest - Destinatio GAP buffer.
   \param src - Source bitblock buffer.
   \param bits - Number of bits to convert.
   \param dest_len - length of the dest. buffer.
   \return  New ength of GAP block or 0 if conversion failed 
   (insufficicent space).

   @ingroup gapfunc
*/
template<typename T> 
    unsigned bit_convert_to_gap(T* BMRESTRICT dest, 
                                const unsigned* BMRESTRICT src, 
                                bm::id_t bits, 
                                unsigned dest_len)
{
    register T* BMRESTRICT pcurr = dest;
    T* BMRESTRICT end = dest + dest_len; 
    register int bitval = (*src) & 1;
    *pcurr |= bitval;

    ++pcurr;
    *pcurr = 0;
    register unsigned bit_idx = 0;
    register int bitval_next;

    unsigned val = *src;

    do
    {
        // We can fast pace if *src == 0 or *src = 0xffffffff

        while (val == 0 || val == 0xffffffff)
        {
           bitval_next = val ? 1 : 0;
           if (bitval != bitval_next)
           {
               *pcurr++ = bit_idx-1; 
               BM_ASSERT((pcurr-1) == (dest+1) || *(pcurr-1) > *(pcurr-2));
               if (pcurr >= end)
               {
                   return 0; // OUT of memory
               }
               bitval = bitval_next;
           }
           bit_idx += sizeof(*src) * 8;
           if (bit_idx >= bits)
           {
               goto complete;
           }
           ++src;
           val = *src;
        }


        register unsigned mask = 1;
        while (mask)
        {
            // Now plain bitshifting. Optimization wanted.

            bitval_next = val & mask ? 1 : 0;
            if (bitval != bitval_next)
            {
                *pcurr++ = bit_idx-1;
                BM_ASSERT((pcurr-1) == (dest+1) || *(pcurr-1) > *(pcurr-2));
                bitval = bitval_next;
                if (pcurr >= end)
                {
                    return 0; // OUT of memory
                }
            }

            mask <<= 1;
            ++bit_idx;

        } // while mask

        if (bit_idx >= bits)
        {
            goto complete;
        }

        ++src;
        val = *src;

    } while(1);

complete:
    *pcurr = bit_idx-1;
    unsigned len = (unsigned)(pcurr - dest);
    *dest = (*dest & 7) + (len << 3);
    return len;
}

/*!
    @ingroup bitfunc 
*/
template<typename T> T bit_convert_to_arr(T* BMRESTRICT dest, 
                                          const unsigned* BMRESTRICT src, 
                                          bm::id_t bits, 
                                          unsigned dest_len)
{
    register T* BMRESTRICT pcurr = dest;
    T* BMRESTRICT end = dest + dest_len; 
    register unsigned bit_idx = 0;

    do
    {
        register unsigned val = *src;
        // We can skip if *src == 0 

        while (val == 0)
        {
            bit_idx += sizeof(*src) * 8;
            if (bit_idx >= bits)
            {
               return (T)(pcurr - dest);
            }
            val = *(++src);
        }

        if (pcurr + sizeof(val)*8 > end) // insufficient space
        {
            return 0;
        }

        for (int i = 0; i < 32; i+=4)
        {
            if (val & 1)
                *pcurr++ = bit_idx;
            val >>= 1; ++bit_idx;
            if (val & 1)
                *pcurr++ = bit_idx;
            val >>= 1; ++bit_idx;
            if (val & 1)
                *pcurr++ = bit_idx;
            val >>= 1; ++bit_idx;
            if (val & 1)
                *pcurr++ = bit_idx;
            val >>= 1; ++bit_idx;
        }
        if (bits <= bit_idx)
            break;

        val = *(++src);

    } while (1);

    return (T)(pcurr - dest);
}



/*! 
    @brief Bitcount for bit string
    
	Function calculates number of 1 bits in the given array of words.
    Make sure the addresses are aligned.

    @ingroup bitfunc 
*/
inline 
bm::id_t bit_block_calc_count(const bm::word_t* block, 
							  const bm::word_t* block_end)
{
    BM_ASSERT(block < block_end);
	bm::id_t count = 0;

#ifdef BM64OPT

    // 64-bit optimized algorithm.

    const bm::id64_t* b1 = (bm::id64_t*) block;
    const bm::id64_t* b2 = (bm::id64_t*) block_end;

    bm::id64_t  acc = *b1++;  // accumulator (sparse vectors optimization)

    do
    {
        bm::id64_t in = *b1++;
        bm::id64_t acc_prev = acc;
        acc |= in;

        if (acc_prev &= in)  // counting bits in accumulator
        {
            acc = (acc & 0x5555555555555555LU) + (acc >> 1 & 0x5555555555555555LU);
            acc = (acc & 0x3333333333333333LU) + (acc >> 2 & 0x3333333333333333LU);
            acc = acc + (acc >> 4) & 0x0F0F0F0F0F0F0F0FLU;
            acc = acc + (acc >> 8);
            acc = acc + (acc >> 16);
            acc = acc + (acc >> 32) & 0x0000007F;
            count += (unsigned)acc;

            acc = acc_prev;
        }
    } while (b1 < b2);
    count += word_bitcount64(acc);  // count-in remaining accumulator 

#else
    // For 32 bit code the fastest method is
    // to use bitcount table for each byte in the block.
    // As optimization for sparse bitsets used bits accumulator
    // to collect ON bits using bitwise OR. 
    bm::word_t  acc = *block++;
    do
    {
        bm::word_t in = *block++;
        bm::word_t acc_prev = acc;
        acc |= in;
        if (acc_prev &= in)  // accumulator miss: counting bits
        {
            BM_INCWORD_BITCOUNT(count, acc);
            acc = acc_prev;
        }
    } while (block < block_end);

    BM_INCWORD_BITCOUNT(count, acc); // count-in remaining accumulator 

#endif
	
    return count;
}

/*!
    Function calculates number of times when bit value changed 
    (1-0 or 0-1).
    
    For 001 result is 2
        010 - 3
        011 - 2
        111 - 1
    
    @ingroup bitfunc 
*/

inline 
bm::id_t bit_count_change(bm::word_t w)
{
    unsigned count = 1;
    w ^= (w >> 1);

    BM_INCWORD_BITCOUNT(count, w);
    count -= (w >> ((sizeof(w) * 8) - 1));
    return count;
}


/*!
    Function calculates number of times when bit value changed 
    (1-0 or 0-1) in the bit block.
        
    @ingroup bitfunc 
*/
inline 
bm::id_t bit_block_calc_count_change(const bm::word_t* block, 
							         const bm::word_t* block_end)
{
    BM_ASSERT(block < block_end);
    bm::id_t count = 1;
    
#ifdef BM64OPT

    // 64-bit optimized algorithm.

    const bm::id64_t* b1 = (bm::id64_t*) block;
    const bm::id64_t* b2 = (bm::id64_t*) block_end;

    bm::id64_t w, w0, w_prev, w_l;
    w = w0 = *b1;
    const int w_shift = sizeof(w) * 8 - 1;
    w ^= (w >> 1);
    count += word_bitcount64(w);
    count -= (w_prev = (w0 >> w_shift)); // negative value correction
    
    for (++b1 ;b1 < b2; ++b1)
    {
        w = w0 = *b1;
        ++count;
        
        if (!w)
        {
            count -= !w_prev;
            w_prev = 0;
        }
        else
        {
            w ^= (w >> 1);
            count += word_bitcount64(w);
            
            w_l = w0 & 1;
            count -= (w0 >> w_shift);  // negative value correction
            count -= !(w_prev ^ w_l);  // word border correction
            
            w_prev = (w0 >> w_shift);
        }
    } // for

#else
    
    bm::word_t  w, w0, w_prev, w_l; 
    
    w = w0 = *block;
    const int w_shift = sizeof(w) * 8 - 1;    
    w ^= (w >> 1);
    BM_INCWORD_BITCOUNT(count, w);
    count -= (w_prev = (w0 >> w_shift)); // negative value correction

    for (++block ;block < block_end; ++block)
    {
        w = w0 = *block;
        ++count;

        if (!w)
        {       
            count -= !w_prev;
            w_prev = 0;
        }
        else
        {
            w ^= (w >> 1);
            BM_INCWORD_BITCOUNT(count, w);
            
            w_l = w0 & 1;
            count -= (w0 >> w_shift);  // negative value correction
            count -= !(w_prev ^ w_l);  // word border correction
            
            w_prev = (w0 >> w_shift);
        }
    } // for
#endif
    return count;
}


/*!
	Function calculates number of 1 bits in the given array of words in
    the range between left anf right bits (borders included)
    Make sure the addresses are aligned.

    @ingroup bitfunc
*/
inline 
bm::id_t bit_block_calc_count_range(const bm::word_t* block,
                                    bm::word_t left,
                                    bm::word_t right)
{
    BM_ASSERT(left <= right);
    
	bm::id_t count = 0;

    unsigned nbit  = left; // unsigned(left & bm::set_block_mask);
    unsigned nword = unsigned(nbit >> bm::set_word_shift);
    nbit &= bm::set_word_mask;

    const bm::word_t* word = block + nword;

    if (left == right)  // special case (only 1 bit to check)
    {
        return (*word >> nbit) & 1;
    }
    unsigned acc;
    unsigned bitcount = right - left + 1;

    if (nbit) // starting position is not aligned
    {
        unsigned right_margin = nbit + (right - left);

        if (right_margin < 32)
        {
            unsigned mask =
                block_set_table<true>::_right[nbit] &
                block_set_table<true>::_left[right_margin];
            acc = *word & mask;
            
            BM_INCWORD_BITCOUNT(count, acc);
            return count;
        }
        else
        {
            acc = *word & block_set_table<true>::_right[nbit];
            BM_INCWORD_BITCOUNT(count, acc);
            bitcount -= 32 - nbit;
        }
        ++word;
    }

    // now when we are word aligned, we can count bits the usual way
    for ( ;bitcount >= 32; bitcount -= 32)