Докладно про підрахунок одиничних бітів

Я хотів би подарувати спільноти Хабра статтю, в якій намагаюся дати досить повний опис підходів до алгоритмів підрахунку одиничних бітів у змінних розміром від 8 до 64 біт. Ці алгоритми належать до розділу так званої «бітової магії» або «бітової алхімії», яка заворожує своєю красою і неочевидностью багатьох програмістів. Я хочу показати, що в основах цієї алхімії немає нічого складного, і ви навіть зможете розробити власні методи підрахунку одиничних бітів, познайомившись з фундаментальними прийомами, складовими подібні алгоритми.



Перш ніж ми почнемо, я відразу хочу попередити, що це стаття не для новачків у програмуванні. Мені необхідно, щоб читач в загальних рисах уявляв собі найпростіші бітові операції (побітове «і», «або», зрушення), добре володів шістнадцятковій системою числення і досить впевнено користувався уявою, представляючи в ньому не завжди короткі бітові послідовності. По можливості, все буде супроводжуватися картинками, але самі розумієте, вони лише спрощують, але не замінюють повне уявлення.

Всі описані прийоми були реалізовані на мові Сі і протестовані в двох режимах: 32 і 64 біта. Таким чином, для більш повного розуміння статті буде краще, щоб ви хоча б приблизно розуміли мову Сі. Тестування проходило на процесорі Core 2 Duo E8400 @3GHz на 64-х бітної Windows 7. Вимір чистого часу роботи програм проводилося з допомогою утиліти runexe. Всі вихідні коди описуваних алгоритмів доступні архів на Яндекс диску, їх компіляція перевірена для компіляторів Visual C++, Intel C++, GCC і CLang, так що в принципі, проблем у вас бути не повинно, якщо хтось захоче перевірити результати у себе. Користувачі Linux, думаю, краще за мене знають, як їм тестувати час роботи програми у себе в системі, тому їм порад не даю.

Серед читачів, можливо, будуть такі, кому простіше подивитися все те ж саме на відео. Я записав таке відео (58 хвилин), в якому у форматі презентації викладено в точності все те ж саме, що буде нижче по тексту, але може трохи в іншому стилі, більш сухо і строго, тоді як текст я спробував трохи оживити. Тому вивчайте матеріал так, як кому зручніше.

відео

Зараз будуть послідовно описані алгоритми, породжувані тим чи іншим набором алхімічних прийомів, у кожному розділі буде таблиця порівняння часу роботи для змінних різного розміру, а в кінці буде зведена таблиця по всім алгоритмам. У всіх алгоритмах використовуються псевдоніми для чисел без знака від 8 до 64 біт.

typedef unsigned char u8;
typedef unsigned short int u16;
typedef unsigned int u32;
typedef unsigned long long u64;


Наївний підхід

Очевидно, що бітова алхімія застосовується зовсім не для того, щоб блищати на співбесіді, а з метою суттєвого прискорення програм. Прискорення по відношенню до чого? По відношенню до тривіальним прийомів, які можуть прийти в голову, коли немає часу більш детально вникнути в завдання. Таким прийомом і є наївний підхід до підрахунку бітів: ми просто «відкушуємо» від числа один біт за іншим і підсумовуємо їх, повторюючи процедуру до тих пір, поки не стане рівним нулю.

u8 CountOnes0 (u8 n) {
u8 res = 0;
while (n) {
res += n&1;
n >>= 1;
}
return res;
}

Я не бачу сенсу що-небудь коментувати в цьому тривіальне циклі. Неозброєним поглядом ясно, що якщо старший біт числа n дорівнює 1, то цикл змушений буде пройтися по всіх бітів числа, перш ніж дістанеться до старшого.

Змінюючи тип вхідного параметра u8 на u16, u32 і u64 ми отримаємо 4 різні функції. Давайте протестуємо кожну з них на потоці з 232 чисел, поданих в хаотичному порядку. Зрозуміло, що для u8 у нас 256 різних вхідних даних, але для одноманітності ми все одно проганяємо 232 випадкових чисел для всіх цих і всіх наступних функцій, причому завжди в одному і тому ж порядку (за подробицями можна звернутися до коду навчальної програми з архіву).

Нижче зазначено в секундах. Для тестування програма запускалася тричі і вибиралося середнє час. Похибка не перевищує 0,1 секунди. Перший стовпець відображає режим компілятора (32-х бітовий вихідний код або 64-х бітовий), далі 4 стовпця відповідають за 4 варіанти вхідних даних.
Режим u8 u16 u32 u64
x86 38,18 72,00 130,49 384,76
x64 37,72 71,51 131,47 227,46
Як ми бачимо, швидкість роботи цілком закономірно зростає із зростанням розміру вхідного параметра. Трохи вибивається із загальної закономірності варіант, коли числа мають розмір 64 біта, а підрахунок йде режимі x86. Ясна річ, що процесор змушений робити чотирикратну роботу при подвоєнні вхідного параметра і навіть добре, що він справляється лише втричі повільніше.

Перша користь цього підходу в тому, що при його реалізації важко помилитися, тому написана таким чином програма може стати еталонною для перевірки більш складних алгоритмів (саме так і було зроблено в моєму випадку). Друга користь універсальності і відносно простий переносимості на числа будь-якого розміру.

Трюк з «відкушуванням» молодших одиничних бітів

Цей алхімічний прийом заснований на ідеї обнулення молодшого одиничного біта. Маючи число n, ми можемо вимовити заклинання n=n&(n-1), забираючи у числа n його молодшу одиничку. Картинка нижче для n=232 прояснить ситуацію для людей, які вперше дізналися про це трюк.


Код програми не сильно змінився.
u8 CountOnes1 (u8 n) {
u8 res = 0;
while (n) {
res ++;
n &= n-1; // Забираємо молодшу одиничку.
}
return res;
}

Тепер цикл виконається рівно стільки разів, скільки одиниць у числі n. Це не позбавляє від найгіршого випадку, коли всі біти в числі поодинокі, але значно скорочує середнє число ітерацій. Сильно даний підхід полегшить страждання процесора? Насправді не дуже, а для 8 біт буде навіть гірше. Нагадаю, що зведена таблиця результатів буде в кінці, а тут у кожному розділі буде своя таблиця.
Режим u8 u16 u32 u64
x86 44,73 55,88 72,02 300,78
x64 40,96 69,16 79,13 126,72


Предподсчет

Не будемо поспішати переходити до «жорстким» заклинанням, розглянемо останній простий прийом, який може врятувати навіть самого досвідченого мага. Даний варіант розв'язання задачі не відноситься безпосередньо до бітової алхімії, однак для повноти картини має бути розглянутий в обов'язковому порядку. Заведемо дві таблиці на 256 і 65536 значень, в яких заздалегідь пораховані відповіді для всіх можливих 1-байтових і 2-байтові величин відповідно.
u8 BitsSetTableFF[256]; // Тут всі відповіді для одного байта
u8 BitsSetTableFFFF[65536]; // Тут всі відповіді для двох байт

Тепер програма для 1 байта буде виглядати так:
u8 CountOnes2_FF (u8 n) {
return BitsSetTableFF[n];
}

Щоб розрахувати число біт у більш великих за розміром числах, їх потрібно розбити на байти. Наприклад, для u32 може бути ось такий код:
u8 CountOnes2_FF (u32 n) {
u8 *p = (u8*)&n;
n = BitsSetTableFF[p[0]]
+ BitsSetTableFF[p[1]]
+ BitsSetTableFF[p[2]]
+ BitsSetTableFF[p[3]];
return n;
}

Або такий, якщо ми застосовуємо таблицю предподсчета для 2-х байт:
u8 CountOnes2_FFFF (u32 n) {
u16 *p = (u16*)&n;
n = BitsSetTableFFFF[p[0]]
+ BitsSetTableFFFF[p[1]];
return n;
}

Ну а далі ви здогадалися, для кожного варіанту розміру вхідного параметра n (крім 8 біт) може існувати два варіанти предподсчета, в залежності від того, яку з двох таблиць ми застосовуємо. Думаю, читачеві зрозуміло, чому ми не можемо просто так взяти і завести таблицю BitsSetTableFFFFFFFF, однак цілком можуть існувати завдання, де і це буде виправданим.

Швидко працює предподсчет? Все сильно залежить від розміру, дивіться таблиці нижче. Перша для однобайтового предподсчета, а друга для двобайтового.
Режим u8 u16 u32 u64
x86 0,01 1,83 21,07 36,25
x64 0,01 1,44 24,79 26,84
Цікавий момент: для режиму x64 предподсчет для u64 працює помітно швидше, можливо, це особливості оптимізації, хоча подібне не проявляється у другому випадку.
Режим u8 u16 u32 u64
x86 --- 0,05 7,95 13,01
x64 --- 0,07 8,49 13,01
Важливе зауваження: цей алгоритм з використанням предподсчета виявляється вигідним лише при дотриманні наступних двох умов: (1) у вас є зайва пам'ять, (2) вам необхідно виконувати розрахунок числа одиничних бітів набагато більше, ніж розмір самої таблиці, тобто є можливість «відіграти» час, витрачений на попереднє заповнення таблиці якимось з простих алгоритмів. Мабуть, можна також мати на увазі екзотичне умова, яка на практиці завжди виконано. Ви повинні гарантувати, що звернення до пам'яті саме по собі швидке і не уповільнює роботу інших функцій системи. Справа в тому, що звернення до таблиці може викинути з пам'яті те, що там було спочатку і уповільнити таким чином якийсь інший ділянку коду. Косяк це ви навряд чи знайдете швидко, однак подібні жахливі оптимізації навряд чи комусь знадобляться на практиці при реалізації звичайних програм.

Множення і залишок від ділення

Візьмемо ж нарешті більш сильні зілля з нашої алхімічної полиці. За допомогою множення і залишку від ділення на степінь двійки без одиниці можна робити досить цікаві речі. Почнемо творити заклинання з одного байта. Для зручності позначимо всі біти одного байта латинськими літерами від «a» до «h». Наше число n прийме вигляд:


Множення n' = n⋅0x010101 (так, через префікс «0x», я позначаю шістнадцяткові числа) робить ні що інше як «тиражування» одного байта в трьох примірниках:

Тепер розіб'ємо подумки наші 24 біта на 8 блоків по 3 біти в кожному (див. нижченаведену картинку, перший рядок таблички). Потім за допомогою двійкового «і» з маскою 0x249249 (другий рядок таблички) обнулим у кожному блоці два старших біти.


Третій рядок таблиці пояснює шістнадцяткову запис маски. В останньому рядку показаний результат, якого ми домагалися: всі біти початкового байта міститися кожен у своєму трехбитовом блоці, але в іншому порядку (порядок нам і не важливий).

Тепер увага: ми повинні скласти ці 8 блоків – і отримаємо суму наших біт!

Виявляється, що залишок від ділення деякого числа A на 2k-1 якраз дає суму k-бітових блоків числа A, теж взяту по модулю 2k-1.
ПруфРозіб'ємо число A (в двійковій запису) на блоки по k біт у кожному (при необхідності можна доповнити самий останній, старший, блок нулями). Позначимо через Aі i-й блок. Тепер запишемо значення числа A через суму цих блоків, помножених на відповідний ступінь двійки:
A= A0⋅20⋅k+ A1⋅21⋅k+...+ AN-1⋅2(N-1)⋅k,
де N – число блоків.
Тепер розрахуємо A mod (2k-1).
A mod (2k-1)= (A0⋅20⋅k+ A1⋅21⋅k+...+ AN-1⋅2(N-1)⋅k) mod (2k-1) = (A0+ A1+...+ AN-1) mod (2k-1).
Все завдяки тому, що 2i⋅k = 1 (mod (2k-1)) для будь-якого цілого ненегативного i. (Тут, правда, важливо, що трюк має сенс, коли k>1, інакше не зовсім зрозуміло як нам інтерпретувати модуль 1). Ось ми і отримали суму блоків по модулю 2k-1.

Тобто від отриманого нами числа потрібно взяти залишок від ділення на 23-1 (сім) — і ми отримуємо суму наших 8-і блоків за модулем 7. Біда в тому, що сума біт може бути дорівнює 7 або 8, в такому випадку алгоритм видасть 0 і 1 відповідно. Але давайте подивимося: в якому випадку ми можемо отримати відповідь 8? Тільки коли n=255. А в якому випадку можемо отримати 0? Тільки коли n=0. Тому якщо алгоритм після взяття залишку на 7 дасть 0, то або ми на вході отримали n=0, або в числі рівно 7 одиничних біт. Підсумовуючи цю міркування, отримуємо наступний код:
u8 CountOnes3 (u8 n) { 
if (n == 0) return 0; // Єдиний випадок, коли відповідь 0.
if (n == 0xFF) return 8; // Єдиний випадок, коли відповідь 8.
n = (0x010101*n & 0x249249) % 7; // Вважаємо число біт за модулем 7.
if (n == 0) return 7; // Гарантовано маємо 7 одиничних бітів.
return n; // Випадок, коли в числі від 1 до 6 одиничних бітів.
} 

У випадку коли n має розмір 16 біт можна розбити його на дві частини по 8 біт. Наприклад, так:
u8 CountOnes3 (u16 n) { 
return CountOnes3 (u8(n&0xFF)) + CountOnes3 (u8(n>>8));

Для випадку 32-х і 64-х розрядних подібне розділення не має сенсу вже навіть у теорії, множення і залишок від ділення з трьома ветвлениями будуть дуже дорого коштувати, якщо виконувати їх 4 або 8 разів підряд. Але я залишив для вас порожні місця в нижченаведеній таблиці, тому якщо ви мені не вірите – заповніть їх, будь ласка, самі. Там швидше всього будуть результати, порівнянні з процедурою CountBits1, якщо у вас схожий процесор (я не кажу про те, що тут можливі оптимізації з допомогою SSE, це вже буде інша розмова).
Режим u8 u16 u32 u64
x86 12,42 30,57 --- ---
x64 13,88 33,88 --- ---
Даний трюк, звичайно, можна зробити і без розгалужень, але тоді нам потрібно, щоб при розбитті числа на блоки в блок вмістилися всі числа від 0 до 8, а цього можна домогтися лише у випадку 4-бітових блоків (і більше). Щоб виконати підсумовування 4-бітових блоків, потрібно підібрати множник, який дозволить правильно «розтиражувати» число і взяти залишок від ділення на 24-1=15, щоб скласти утворені блоки. Досвідчений алхімік (який знає математику) легко підбере такий множник: 0x08040201. Чому він обраний таким?


Справа в тому, що нам необхідно, щоб всі біти початкового числа зайняли правильні позиції в своїх 4-бітових блоках (зображення вище), а коль скоро 8 і 4 не є взаємно простими числами, звичайне копіювання 8 бітів 4 рази не дасть правильного розташування потрібних бітів. Нам доведеться додати на нашому байту один нулик, тобто тиражувати 9 бітів, так як 9 взаємно просто з 4. Так ми отримаємо число, що має розмір 36 біт, але в якому всі біти початкового байта стоять на молодших позиціях 4-бітових блоків. Залишилося тільки взяти побітове «і» з числом 0x111111111 (другий рядок на картинці вище), щоб обнулити по три старших біта в кожному блоці. Потім блоки потрібно скласти.

При такому підході програма підрахунку одиничних бітів у байті буде гранично простий:
u8 CountOnes3_x64 (u8 n) {
return ((u64)0x08040201*n & 0x111111111) % 15;
}

Недолік програми очевидний: потрібно вихід 64-бітову арифметику зі всіма витікаючими звідси наслідками. Можна помітити, що в дійсності ця програма використовує лише 33 біта з 64-х (старші 3 біта обнуляються), і в принципі можна зрозуміти, як перенести дані обчислення в 32-х бітову арифметику, але розповіді про подібні оптимізацію не входять в тему цього керівництва. Давайте поки просто вивчати прийоми, а оптимізувати їх вам доведеться самим вже під конкретну задачу.

Відповімо на питання про те, якого розміру може бути змінна n, щоб даний трюк правильно працював для неї. Коль скоро ми беремо залишок від ділення на 15, така змінна не може мати розмір більше 14 біт, в іншому випадку доведеться застосувати розгалуження, як ми робили це раніше. Але для 14 біт прийом працює, якщо додати до 14-ти біт один нулик, щоб всі біти встали на свої позиції. Тепер я буду вважати, що ви в цілому засвоїли суть прийому і зможете самі без праці підібрати множник для тиражування і маску для обнулення непотрібних бітів. Покажу відразу готовий результат.
u8 CountOnes3_x64 (u14 n) { // Це не помилка (і не повинно працювати)!
return (n*0x200040008001llu & 0x111111111111111llu) % 15;
}

Ця програма вище показує, як міг би виглядати код, якби у вас мінлива розміром 14 бітів без знака. Цей код буде працювати зі змінною в 15 біт, але за умови, що максимум лише 14 з них дорівнює одиниці, або якщо випадок, коли n=0x7FFF ми розберемо окремо. Це все потрібно розуміти для того, щоб написати правильний код для змінної типу u16. Ідея в тому, щоб спочатку «відкусити» молодший біт, порахувати біти в останньому 15-ти бітовому числі, а потім назад додати «відкушений» біт.
static u8 CountOnes3_x64 (u16 n) { 
u8 leastBit = n&1; // Беремо молодший біт.
n >>= 1; // Залишаємо тільки 15 біт вихідного числа.
if (n == 0) return leastBit; // Якщо вийшов 0, значить відповіддю буде молодший біт.
if (n == 0x7FFF) return leastBit + 15; // Єдиний випадок, коли відповідь 15+молодший біт.
return leastBit + (n*0x200040008001llu & 0x111111111111111llu) % 15; // Відповідь (максимум 14+молодший біт).
}

Для n розміром 32 біта доводиться чаклувати вже з більш серйозним обличчям… По-перше, відповідь влазить тільки в 6 біт, але можна розглянути окремий випадок, коли n=232-1 і спокійно робити розрахунки в полях розміром 5 біт. Це економить місце і дозволяє нам розбити 32-бітове поле числа n на 3 частини по 12 біт у кожній (так, останнє поле буде не повним). Оскільки 12⋅5=60, ми можемо тиражувати одне 12-бітове поле 5 разів, підібравши множник, а потім скласти 5-бітові блоки, взявши залишок від ділення на 31. Зробити це треба 3 рази для кожного поля. Підбиваючи підсумок, отримуємо такий код:
static u8 CountOnes3_x64 ( u32 n ) { 
if (n == 0) return 0;
if (n + 1 == 0) return 32;
u64 res = (n&0xFFF)*0x1001001001001llu & 0x84210842108421llu;
res += ((n&0xFFF000)>>12)*0x1001001001001llu & 0x84210842108421llu;
res += (n>>24)*0x1001001001001llu & 0x84210842108421llu;
res %= 0x1F;
if (res == 0) return 31;
return res;
}

Тут точно також можна було замість трьох розгалужень взяти 3 залишку від ділення, але я вибрав гіллястий варіант, на моєму процесорі він буде працювати краще.

Для n розміром 64 біта мені не вдалося придумати відповідного заклинання, в якому було б не так багато множень і додавань. Виходило або 6, або 7, а це надто багато для такого завдання. Інший варіант — вихід до 128-бітову арифметику, а це вже не зрозумій яким «відкатом» для нас обернеться, непідготовленого мага може і до стінки відкинути :)

Давайте краще подивимося на час роботи.
Режим u8 u16 u32 u64
x86 39,78 60,48 146,78 ---
x64 6,78 12,28 31,12 ---
Очевидним висновком з цієї таблиці буде те, що 64-х бітова арифметика погано сприймається в 32-х бітовому режимі виконання, хоча загалом-то алгоритм непоганий. Якщо згадати швидкість алгоритму предподсчета в режимі x64 для однобайтовим таблиці для випадку u32 (24,79 с), то отримаємо, що даний алгоритм відстає лише на 25%, а це привід до змагання, втіленому в наступному розділі.

Заміна взяття залишку на множення і зсув

Недолік операції взяття залишку всім очевидний. Це ділення, а розподіл – це довго. Зрозуміло, сучасні компілятори знають алхімію і вміють замінювати поділ на множення зі зсувом, а щоб отримати залишок, потрібно відняти з діленого приватне, помножене на дільник. Тим не менш, це все одно довго! Виявляється, що у стародавніх сувоях заклинателів коду зберігся один цікавий спосіб оптимізації попереднього алгоритму. Ми можемо підсумувати k-бітові блоки не взяттям залишку від ділення, а ще одним множенням на маску, за допомогою якої обнуляли зайві біти в блоках. Ось як це виглядає для n розміром в 1 байт.

Для початку знову тиражуємо байт тричі і видаляємо по два старших біти у кожного 3-бітового блоку за допомогою вже розглянутою вище формули 0x010101⋅n & 0x249249.


Кожен трехбитовый блок я для зручності позначив великою латинською літерою. Тепер отриманий результат множимо на ту ж саму маску 0x249249. Маска містить одиничний біт у кожної 3-ї позиції, тому таке множення еквівалентно складання числа самого з собою 8 разів, кожен раз зі зсувом на 3 біти:


Що ми бачимо? Біти з 21 по 23 і дають нам потрібну суму! При цьому переповнення в якому-небудь з блоків праворуч неможливо, так як там в одному блоці не буде числа, більшого 7. Проблема лише в тому, що якщо наша сума дорівнює 8, ми отримаємо 0, але це не страшно, адже це єдиний випадок можна розглянути окремо.
u8 CountOnes3_x64_m (u8 n) { 
if (n == 0xFF) return 8;
return ((u64(0x010101*n & 0x249249) * 0x249249) >> 21) & 0x7;
}

По суті, ми взяли код з попереднього розділу і замінили в ньому взяття залишку від ділення на 7 множення, зсув і побітове «І» в кінці. При цьому замість 3-х розгалужень залишилося лише одне.

Щоб скласти аналогічну програму для 16 біт, нам потрібно взяти код з попереднього розділу, в якому показано, як це робиться за допомогою взяття залишку від ділення на 15 і замінити цю процедуру множенням. При цьому неважко помітити те, які умови можна прибрати з коду.
u8 CountOnes3_x64_m (u16 n) { 
u8 leastBit = n&1; // Беремо молодший біт
n >>= 1; // Залишаємо тільки 15 біт.
return leastBit
+ (( (n*0x200040008001llu & 0x111111111111111llu)*0x111111111111111llu 
>> 56) & 0xF); // Відповідь (максимум 15 + молодший біт).
}

Для 32-х біт ми робимо те ж саме: беремо код з попереднього розділу і, порисовав трохи на папері, розуміємо, яким буде зсув, якщо замінити залишок на множення.
u8 CountOnes3_x64_m ( u32 n ) { 
if (n+1 == 0) return 32;
u64 res = (n&0xFFF)*0x1001001001001llu & 0x84210842108421llu;
res += ((n&0xFFF000)>>12)*0x1001001001001llu & 0x84210842108421llu;
res += (n>>24)*0x1001001001001llu & 0x84210842108421llu;
return (res*0x84210842108421llu >> 55) & 0x1F;
}

Для 64-х розрядних я теж не зміг придумати чогось такого, щоб не змушувало б мій процесор виконувати роль грубки.
Режим u8 u16 u32 u64
x86 12,66 42,37 99,90 ---
x64 3,54 4,51 18,35 ---
Приємно здивували результати для режиму x64. Як і очікувалося, ми обігнали предподсчет з однобайтовим таблицею для випадку u32. Можна взагалі обігнати предподсчет? Хороше питання :)

Паралельне підсумовування

Мабуть, це самий поширений трюк, який дуже часто повторюють один за одним не цілком досвідчені заклинателі, не розуміючи, як він точно працює.

Почнемо з 1 байта. Байт складається з 4-х полів по 2 біта, спочатку просумуємо біти в цих полях, сказавши щось на кшталт:
n = (n>>1)&0x55 + (n&0x55);

Ось пояснювальна картинка до цієї операції (раніше, позначаємо біти одного байта першими латинськими літерами):

Одне з двійковими «І» залишає тільки молодші біти кожного двухбитового блоку, друге залишає старші біти, але зсуває їх на позиції, що відповідають молодшим бітам. У результаті підсумовування отримуємо суму суміжних бітів в кожному двухбитовом блоці (останній рядок на картинці вище).

Тепер складемо парами числа, що знаходяться в двухбитовых полях, поміщаючи результат в 2 четырехбитовых поля:
n = (n>>2)&0x33 + (n&0x33);

Нижченаведена картинка пояснює результат. Наводжу її тепер без зайвих слів:

Нарешті, складемо два числа в четырехбитовых полях:
n = (n>>4)&0x0F + (n&0x0F);


Діючи за аналогією, можна поширити прийом на будь-яке число біт, що дорівнює ступеню двійки. Число рядків заклинання одно двійковому логарифму від кількості бітів. Вловивши ідею, погляньте побіжно на 4 функції, записаних нижче, щоб переконатися у правильності свого розуміння.
u8 CountOnes4 (u8 n) {
n = ((n>>1) & 0x55) + (n & 0x55);
n = ((n>>2) & 0x33) + (n & 0x33);
n = ((n>>4) & 0x0F) + (n & 0x0F);
return n;
}

u8 CountOnes4 (u16 n) {
n = ((n>>1) & 0x5555) + (n & 0x5555);
n = ((n>>2) & 0x3333) + (n & 0x3333);
n = ((n>>4) & 0x0F0F) + (n & 0x0F0F);
n = ((n>>8) & 0x00FF) + (n & 0x00FF);
return n;
}

u8 CountOnes4 (u32 n) {
n = ((n>>1) & 0x55555555) + (n & 0x55555555);
n = ((n>>2) & 0x33333333) + (n & 0x33333333);
n = ((n>>4) & 0x0F0F0F0F) + (n & 0x0F0F0F0F);
n = ((n>>8) & 0x00FF00FF) + (n & 0x00FF00FF);
n = ((n>>16) & 0x0000FFFF) + (n & 0x0000FFFF);
return n;
}

u8 CountOnes4 (u64 n) {
n = ((n>>1) & 0x5555555555555555llu) + (n & 0x5555555555555555llu);
n = ((n>>2) & 0x3333333333333333llu) + (n & 0x3333333333333333llu);
n = ((n>>4) & 0x0F0F0F0F0F0F0F0Fllu) + (n & 0x0F0F0F0F0F0F0F0Fllu);
n = ((n>>8) & 0x00FF00FF00FF00FFllu) + (n & 0x00FF00FF00FF00FFllu);
n = ((n>>16) & 0x0000FFFF0000FFFFllu) + (n & 0x0000FFFF0000FFFFllu);
n = ((n>>32) & 0x00000000FFFFFFFFllu) + (n & 0x00000000FFFFFFFFllu);
return n;
}

На цьому паралельне підсумовування не закінчується. Розвинути ідею дозволяє те спостереження, що в кожній сходинці двічі використовується одна і та ж бітова маска, що ніби наводить на думка «а чи не можна як-небудь тільки один раз виконати побітове «І»?». Можна, але не одразу. Ось що можна зробити, якщо взяти в якості прикладу код для u32 (дивіться коментарі).
u8 CountOnes4 (u32 n) {
n = ((n>>1) & 0x55555555) + (n & 0x55555555); // Можна замінити на різницю
n = ((n>>2) & 0x33333333) + (n & 0x33333333); // не Можна виправити
n = ((n>>4) & 0x0F0F0F0F) + (n & 0x0F0F0F0F); // Можна винести & за дужку
n = ((n>>8) & 0x00FF00FF) + (n & 0x00FF00FF); // Можна винести & за дужку
n = ((n>>16) & 0x0000FFFF) + (n & 0x0000FFFF); // Можна взагалі прибрати &
return n; // Неявне обрізання на 8-й молодшим бітам.
}

В якості вправи я б хотів запропонувати довести самостійно те, чому нижченаведений код буде точним відображенням попереднього. Для першого рядка я даю підказку, але не дивіться в неї відразу:
ПідказкаДвухбитовый блок ab має точне значення 2a+b, значить віднімання зробить його рівним...?

u8 CountOnes4_opt (u32 n) {
n - = n>>1) & 0x55555555;
n = ((n>>2) & 0x33333333 ) + (n & 0x33333333);
n = ((n>>4) + n) & 0x0F0F0F0F;
n = ((n>>8) + n) & 0x00FF00FF;
n = ((n>>16) + n);
return n;
}

Аналогічні варіанти оптимізації можливі і для інших типів даних.

Нижче наводяться дві таблиці: одна для звичайного паралельного підсумовування, а друга для оптимізованого.
Режим u8 u16 u32 u64
x86 7,52 14,10 21,12 62,70
x64 8,06 11,89 21,30 22,59
Режим u8 u16 u32 u64
x86 7,18 11,89 18,86 65,00
x64 8,09 10,27 19,20 19,20


У цілому ми бачимо, що оптимізований алгоритм працює добре, але програє в звичайному режимі x86 для u64.

Комбінований метод

Ми бачимо, що найкращі варіанти підрахунку одиничних бітів – це паралельний метод (з оптимізацією) і метод тиражування з множенням для підрахунку суми блоків. Ми можемо об'єднати обидва методу, отримуючи комбінований алгоритм.

Перше, що потрібно зробити — виконати перші три рядки паралельного алгоритму. Це дасть нам точну суму бітів в кожному байті числа. Наприклад, для u32 виконаємо наступне:
n= (n>>1) & 0x55555555;
n = ((n>>2) & 0x33333333) + (n & 0x33333333);
n = (((n>>4) + n) & 0x0F0F0F0F );

Тепер наше число n складається з 4 байт, які слід розглядати як 4 числа, суму яких ми шукаємо:


Ми можемо знайти суму цих 4-х байт, якщо помножимо число n на 0x01010101. Ви добре розумієте, що означає таке множення, для зручності визначення позиції, в якій буде знаходитися відповідь, наводжу картинку:


Відповідь знаходиться в 3-байті (якщо вважати їх від 0). Таким чином, комбінований прийом для u32 буде виглядати так:
u8 CountOnes5 ( u32 n ) {
n - = n>>1) & 0x55555555;
n = ((n>>2) & 0x33333333 ) + (n & 0x33333333);
n = (((( n>>4) + n) & 0x0F0F0F0F) * 0x01010101) >> 24;
return n; // Тут відбувається неявне обрізання на 8 молодшим бітам.
}

Він же для u16:
static u8 CountOnes5 (u16 n) {
n - = n>>1) & 0x5555;
n = ((n>>2) & 0x3333) + (n & 0x3333);
n = (((( n>>4) + n) & 0x0F0F) * 0x0101) >> 8;
return n; // Тут відбувається неявне обрізання на 8 молодшим бітам.
}

Він же для u64:
u8 CountOnes5 ( u64 n ) {
n - = n>>1) & 0x5555555555555555llu;
n = ((n>>2) & 0x3333333333333333llu ) + (n & 0x3333333333333333llu);
n = (((( n>>4) + n) & 0x0F0F0F0F0F0F0F0Fllu)
* 0x0101010101010101) >> 56;
return n; // Тут відбувається неявне обрізання на 8 молодшим бітам.
}

Швидкість роботи цього методу ви можете подивитися відразу в підсумковій таблиці.

Підсумкове порівняння

Я пропоную читачеві самостійно зробити цікавлять його висновки, вивчивши дві нижченаведені таблиці. В них я позначив назва методів, програми до яких ми реалізували, а також позначив прямокутною рамкою ті підходи, які я вважаю найкращими в кожному конкретному випадку. Тих, хто думав, що предподсчет завжди виграє, чекає невеликий сюрприз для режиму x64.

Підсумкове порівняння для режиму компіляції x86.


Підсумкове порівняння для режиму компіляції x64.


Примітка

Ні в якому разі не розглядайте підсумкову таблицю як доказ на користь того або іншого підходу. Повірте, що на вашому процесорі і з вашим компілятором деякі числа в такій таблиці будуть зовсім іншими. На жаль, ми ніколи не можемо точно сказати, який з алгоритмів виявиться краще в тому чи іншому випадку. Під кожну задачу потрібно заточувати конкретний метод, а універсального швидкого алгоритму, на жаль, не існує.

Я виклав ті ідеї, про які знаю сам, але це лише ідеї, конкретні реалізації яких у різних комбінаціях можуть бути дуже різними. Об'єднуючи ці ідеї різними способами, ви можете отримувати величезну кількість різних алгоритмів підрахунку одиничних бітів, кожен з яких цілком може виявитися гарним в якомусь своєму випадку.

Дякую за увагу. До нових зустрічей!

UPD: Інструкція POPCNT з SSE4.2 не включена в список тестування, тому що у мене немає процесора, який підтримує SSE4.2.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.