Создание списка целых чисел с заданным количеством битовых наборов и суммой битовых индексов

Question

Я хотел бы эффективно сгенерировать список целых чисел (желательно упорядоченных) со следующими определяющими свойствами:

1. Все целые числа имеют одинаковое количество битов N .

2. Все целые числа имеют одинаковую сумму битовых индексов K.

Чтобы быть определенным, для целого числа I его двоичный представление:

$I=\sum_{j=0}^M c_j 2^j$ where $c_j=0$ or $1$

Количество битовых наборов:

$N(I)=\sum_{j=0}^M c_j$

Сумма битовых индексов:

$K(I)=\sum_{j=0}^M j c_j$

У меня есть неэффективный способ сгенерируйте список следующим образом: сделайте do / for l oop над целыми числами с приращением с помощью функции "snoob" - наименьшее следующее целое число с таким же набором битов и при каждом приращении проверяя, имеет ли оно правильное значение K

это крайне неэффективно, потому что в целом, начиная с целого числа с правильным значением N и K, целое число snoob от I не имеет правильного K, и приходится выполнять много вычислений snoob чтобы получить следующий intege r с обоими N и K, равными выбранным значениям. Использование snoob дает упорядоченный список, который удобен для дихотоми c поиска, но не является абсолютно обязательным.

Подсчет количества элементов в этом списке легко выполняется с помощью рекурсии, если рассматривать его как подсчет номеров разделов. вот рекурсивная функция в fortran 90, выполняющая эту работу:

=======================================================================
recursive function BoundedPartitionNumberQ(N, M, D)  result (res)
implicit none

  ! number of partitions of N into M distinct integers, bounded by D
  ! appropriate for Fermi counting rules

   integer(8) :: N, M, D, Nmin
   integer(8) :: res
    
    Nmin = M*(M+1)/2       ! the Fermi sea
    
    if(N < Nmin) then
        res = 0

    else if((N == Nmin) .and. (D >= M)) then
        res = 1

    else if(D < M) then
       res = 0

    else if(D == M)  then
       if(N == Nmin) then
              res = 1
       else 
              res = 0  
       endif

    else if(M == 0) then
       res = 0

     else

     res = BoundedPartitionNumberQ(N-M,M-1,D-1)+BoundedPartitionNumberQ(N-M,M,D-1)

     endif

    end function BoundedPartitionNumberQ
========================================================================================

Мое нынешнее решение неэффективно, когда я хочу сгенерировать списки с несколькими $10^7$ элементами. В конечном итоге я хочу оставаться в сфере C / C ++ / Fortran и достигать списков длиной до нескольких $10^9$

мой текущий код f90 следующий:

program test
implicit none

integer(8) :: Nparticles
integer(8) :: Nmax, TmpL, CheckL, Nphi
integer(8) :: i, k, counter
integer(8) :: NextOne

Nphi = 31        ! word size is Nphi+1
Nparticles = 16  ! number of bit set

print*,Nparticles,Nphi

Nmax = ishft(1_8, Nphi + 1) - ishft(1_8, Nphi + 1 - Nparticles)

i = ishft(1, Nparticles) - 1

counter = 0

! integer CheckL is the sum of bit indices

CheckL = Nparticles*Nphi/2  ! the value of the sum giving the largest list

do while(i .le. Nmax)   ! we increment the integer

    TmpL = 0

    do k=0,Nphi
        if (btest(i,k)) TmpL = TmpL + k
    end do

    if (TmpL == CheckL) then    ! we check whether the sum of bit indices is OK

        counter = counter + 1

    end if

    i = NextOne(i)   ! a version of "snoob" described below

end do

print*,counter

end program

!==========================================================================
function NextOne (state)
implicit none

integer(8) :: bit    
integer(8) :: counter 
integer(8) :: NextOne,state,pstate

bit     =  1
counter = -1
  
!  find first one bit 

do  while (iand(bit,state) == 0)

    bit = ishft(bit,1)

end do

!  find next zero bit 

do  while (iand(bit,state) /= 0)
    
    counter = counter + 1
    bit = ishft(bit,1)

end do

if (bit == 0) then 

    print*,'overflow in NextOne'
    NextOne = not(0)
  
else 

    state = iand(state,not(bit-1))  ! clear lower bits i &= (~(bit-1));

    pstate = ishft(1_8,counter)-1 ! needed by IBM/Zahir compiler

 !  state = ior(state,ior(bit,ishft(1,counter)-1)) ! short version OK with gcc

    state = ior(state,ior(bit,pstate))

    NextOne = state

end if

end function NextOne

Answer 1

Поскольку вы упомянули C / C ++ / Fortran, я попытался сохранить этот язык относительно независимым / легко переносимым , но также добавил быстрее встроенные альтернативы где это применимо.

Все целые числа имеют одинаковое количество установленных битов N

Тогда мы также можем сказать, все допустимые целые числа будут перестановками N установить биты.

Сначала мы должны сгенерировать начальную / минимальную перестановку:

uint32_t firstPermutation(uint32_t n){
    // Fill the first n bits (on the right)
    return (1 << n) -1;
}

Затем мы должны установить конечную / максимальную перестановку - указав «точку остановки»:

uint32_t lastPermutation(uint32_t n){
    // Fill the last n bits (on the left)
    return (0xFFFFFFFF >> n) ^ 0xFFFFFFFF;
}

Наконец, нам нужен способ получить следующую перестановку.

uint32_t nextPermutation(uint32_t n){
    uint32_t t = (n | (n - 1)) + 1;
    return t | ((((t & -t) / (n & -n)) >> 1) - 1);
}

// or with builtins:
uint32_t nextPermutation(uint32_t &p){
    uint32_t t = (p | (p - 1));
    return (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(p) + 1));
}

Все целые числа имеют одинаковую сумму битовых индексов K

Предполагая, что это целые числа (32 бита), вы можете использовать эту последовательность ДеБрюйна для быстрого определения индекса первого установленного бита - fsb . Аналогичные последовательности существуют для других типов / битовых чисел, например этот один может быть адаптирован для использования.

Удалив текущий fsb , мы можем применить вышеупомянутую технику к определить индекс следующего fsb и т. д.

int sumIndices(uint32_t n){
    const int MultiplyDeBruijnBitPosition[32] = {
      0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
      31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
    };

    int sum = 0;
    // Get fsb idx
    do sum += MultiplyDeBruijnBitPosition[((uint32_t)((n & -n) * 0x077CB531U)) >> 27];        
    // strip fsb
    while (n &= n-1);   

    return sum;
}

// or with builtin
int sumIndices(uint32_t n){
    int sum = 0;
    do sum += __builtin_ctz(n);
    while (n &= n-1);
    return sum;
}

Наконец, мы можем перебирать каждую перестановку, проверяя, соответствует ли сумма всех индексов указанному значению K.

p = firstPermutation(n);
lp = lastPermutation(n);

do {
    p = nextPermutation(p);
    if (sumIndices(p) == k){
        std::cout << "p:" << p << std::endl;
    }
} while(p != lp);

Вы можете легко изменить код «обработчика», чтобы сделать что-то подобное, начиная с заданное целое число - используя его значения N и K.

тест встроенного c по сравнению с самореализацией

Answer 2

Базовая c рекурсивная реализация может быть:

void listIntegersWithWeight(int currentBitCount, int currentWeight, uint32_t pattern, int index, int n, int k, std::vector<uint32_t> &res)
{
    if (currentBitCount > n ||
        currentWeight > k)
        return;

    if (index < 0)
    {
        if (currentBitCount == n && currentWeight == k)
            res.push_back(pattern);
    }
    else
    {
        listIntegersWithWeight(currentBitCount, currentWeight, pattern, index - 1, n, k, res);
        listIntegersWithWeight(currentBitCount + 1, currentWeight + index, pattern | (1u << index), index - 1, n, k, res);
    }
}

Это не мое предложение, это просто отправная точка. На моем P C для n = 16, k = 248 и эта версия, и итеративная версия занимают почти (но не совсем) 9 секунд. Почти столько же времени, но это просто совпадение. Можно выполнить дополнительную обрезку:

• currentBitCount + index + 1 < n если количество установленных битов не может достигнуть n с количеством оставшихся незаполненных позиций, продолжать бессмысленно.
• currentWeight + (index * (index + 1) / 2) < k если сумма позиций не может достичь k, продолжать бессмысленно.

Вместе:

void listIntegersWithWeight(int currentBitCount, int currentWeight, uint32_t pattern, int index, int n, int k, std::vector<uint32_t> &res)
{
    if (currentBitCount > n || 
        currentWeight > k ||
        currentBitCount + index + 1 < n ||
        currentWeight + (index * (index + 1) / 2) < k)
        return;

    if (index < 0)
    {
        if (currentBitCount == n && currentWeight == k)
            res.push_back(pattern);
    }
    else
    {
        listIntegersWithWeight(currentBitCount, currentWeight, pattern, index - 1, n, k, res);
        listIntegersWithWeight(currentBitCount + 1, currentWeight + index, pattern | (1u << index), index - 1, n, k, res);
    }
}

На моем P C с теми же параметрами это только занимает полсекунды. Вероятно, его можно улучшить и дальше.