Обстоятельно о подсчете единичных битов

Тема в разделе "WASM.ARTICLES", создана пользователем Mikl___, 11 мар 2023.

  1. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.707
    Статья @Zealint унесена мною с хабра Было бы жалко, если такая прелесть потеряется :)

    Статья, в которой дано достаточно полное описание подходов к алгоритмам подсчета единичных битов в переменных размером от 8 до 64 битов. Алгоритмы относятся к так называемой «битовой магии» или «битовой алхимии», которая завораживает своей красотой и неочевидностью. В этой алхимии нет ничего сложного, и вы сможете разработать собственные методы подсчета единичных битов, познакомившись с фундаментальными приемами, составляющими подобные алгоритмы.
    Это статья не для новичков в программировании. Необходимо, чтобы читатель в общих чертах представлял себе простейшие битовые операции («и», «или», сдвиг), владел шестнадцатеричной системой счисления и пользовался воображением, представляя не всегда короткие битовые последовательности. По возможности, всё будет сопровождаться таблицами, но они лишь упрощают, но не заменяют полное представление.
    Описанные приемы были реализованы на языке Си и протестированы в режимах 32 и 64 бита. Для более полного понимания статьи необходимо, чтобы вы хотя бы приблизительно понимали язык Си. Тестирование проходило на процессоре Core 2 Duo E8400 @3GHz на 64-х битовой Windows 7. Измерение чистого времени работы программ проводилось с помощью утилиты runexe. Все исходные коды описываемых алгоритмов были доступны в архиве (https://yadi.sk/d/XiHe6MNtomAeh) на Яндекс диске, их компиляция проверена для компиляторов Visual C++, Intel C++, GCC и CLang (на данный момент архив не доступен).
    Возможно, будут такие, кому проще посмотреть всё то же самое на видео. В 58 минутном видео изложено то же самое, но более сухо и строго.

    Будут последовательно описаны алгоритмы, порождаемые тем или иным набором приемов, в каждом разделе будет таблица сравнения времени работы для переменных разного размера, в конце будет сводная таблица по всем алгоритмам. Во всех алгоритмах используются псевдонимы для чисел без знака от 8 до 64 бит.
    unsigned charu8
    unsigned short intu16
    unsigned intu32
    unsigned long longu64

    Наивный подход

    Битовая алхимия применяется с целью ускорения программ. Ускорения по отношению к чему? По отношению к тривиальным приемам, которые могут прийти в голову, когда нет времени более детально вникнуть в задачу. Таким приемом и является наивный подход к подсчету битов — «откусываем» от числа один бит за другим и суммируем их, повторяя процедуру до тех пор, пока оставшееся количество бит не станет равным нулю.
    Код (C):
    1. u8 CountOnes0 (u8 n) {
    2.   u8 res = 0;
    3.   while (n) {
    4.     res += n&1;
    5.     n >>= 1;
    6.   }
    7.   return res;
    8. }
    Нет смысла что-либо комментировать в этом цикле. Если старший бит числа n равен 1, то цикл вынужден будет пройтись по всем битам числа, прежде чем доберется до старшего.
    Меняя тип входного параметра u8 на u16, u32 и u64 получим 4 различные функции. Протестируем каждую из них на потоке из 232 чисел, подаваемых в хаотичном порядке. Для u8 имеем 256 различных входных данных, но для единообразия прогоняем 232 случайных чисел для этих и всех последующих функций, всегда в одном и том же порядке (за подробностями можно обратиться к коду учебной программы из архива).
    Время в таблице указано в секундах. Для тестирования программа запускалась трижды и выбиралось среднее время. Погрешность не превышает 0,1 секунды. Первый столбец отражает режим компилятора (32-х битовый исходный код или 64-х битовый), 4 столбца отвечают за 4 варианта входных данных.
    Режимu8u16u32u64
    x8638,1872,00130,49384,76
    x6437,7271,51131,47227,46
    Скорость работы возрастает с ростом размера входного параметра. Из общей закономерности выбивается вариант, когда числа имеют размер 64 бита, а подсчет идет режиме x86. Процессор делает четырехкратную работу при удвоении входного параметра, но справляется втрое медленнее.
    Польза этого подхода в том, что при его реализации трудно ошибиться, поэтому написанная таким образом программа может стать эталонной для проверки более сложных алгоритмов. Вторая польза в универсальности и относительно простой переносимости на числа любого размера.

    Трюк с «откусыванием» младших единичных битов

    Этот прием основан на идее обнуления младшего единичного бита. Имея число n, можем произвести [math]n=n\&(n-1)[/math], забирая у числа n его младшую единичку. Таблица ниже для [math]n[/math] = 0E8h = 11101000b прояснит ситуацию для людей, впервые узнавших об этом трюке.
    binhex
    n11101000
    E8​
    n-111100111
    E7​
    n&(n-1)11100000
    E0​
    Код программы не сильно изменился.
    Код (C):
    1. u8 CountOnes1 (u8 n) {
    2.  u8 res = 0;
    3.  while (n) {
    4.  res ++;
    5.  n &= n-1; // Забираем младшую единичку.
    6.  }
    7.  return res;
    8. }
    Цикл выполнится ровно столько раз, сколько единиц в числе n. Это не избавляет от худшего случая, когда все биты в числе единичные, но значительно сокращает среднее число итераций. Сильно ли данный подход ускоряет работу процессора? Не очень, а для 8 бит будет даже хуже.
    Режимu8u16u32u64
    x8644,7355,8872,02300,78
    x6440,9669,1679,13126,72
    Код программы не сильно изменился.
    Код (C):
    1. u8 CountOnes1 (u8 n) {
    2.  u8 res = 0;
    3.  while (n) {
    4.  res ++;
    5.  n &= n-1; // Забираем младшую единичку.
    6.  }
    7.  return res;
    8. }
    Режимu8u16u32u64
    x8644,7355,8872,02300,78
    x6440,9669,1679,13126,72

    Предподсчет

    Рассмотрим простой прием, который может спасти даже самого неопытного мага. Заведем две таблицы на 256 (=28) и 65536 (=216) значений, в которых заранее посчитаны ответы для всех возможных 1-байтовых и 2-байтовых величин соответственно.
    Код (C):
    1. u8 BitsSetTableFF[256];     // Здесь все ответы для одного байта
    2.   u8 BitsSetTableFFFF[65536]; // Здесь все ответы для двух байт
    Программа для 1 байта выглядит так
    Код (C):
    1. u8 CountOnes2_FF (u8 n) {
    2.     return BitsSetTableFF[n];
    3.   }
    Чтобы рассчитать число бит в более крупном по размеру числе, нужно разбить его на байты. Для u32 может быть вот такой код:
    Код (C):
    1. u8 CountOnes2_FF (u32 n) {
    2.     u8 *p = (u8*)&n;
    3.     n = BitsSetTableFF[p[0]]
    4.       + BitsSetTableFF[p[1]]
    5.       + BitsSetTableFF[p[2]]
    6.       + BitsSetTableFF[p[3]];
    7.     return n;
    8.   }
    Или такой, если применить таблицу предподсчета для 2-х байт:
    Код (C):
    1.   u8 CountOnes2_FFFF (u32 n) {
    2.     u16 *p = (u16*)&n;
    3.     n = BitsSetTableFFFF[p[0]]
    4.       + BitsSetTableFFFF[p[1]];
    5.     return n;
    6.   }
    Для каждого варианта размера входного параметра n (кроме 8 бит) может существовать два варианта предподсчета, в зависимости от того, которую из двух таблиц применяют. Понятно, почему не можем завести таблицу BitsSetTableFFFFFFFF, вполне могут существовать задачи, где это будет оправданным.
    Быстро ли работает предподсчет? Это сильно зависит от размера числа. Первая таблица для 1 байтового предподсчета, вторая для 2-х байтового.
    Режимu8u16u32u64
    x860,011,8321,0736,25
    x640,011,4424,7926,84
    Для режима x64 предподсчет для u64 работает быстрее, возможно, это особенности оптимизации.
    Режимu8u16u32u64
    x860,057,9513,01
    x640,078,4913,01
    Данный алгоритм с использованием предподсчета выгоден при соблюдении трех условий:
    1. у вас есть лишняя память,
    2. вам требуется выполнять расчет числа единичных битов намного больше раз, чем размер самой таблицы, то есть имеется возможность «отыграть» время, потраченное на предварительное заполнение таблицы каким-то из нехитрых алгоритмов.
    3. экзотическое условие, которое на практике всегда выполнено. Вы должны гарантировать, что обращение к памяти само по себе быстрое и не замедляет работу других функций системы. Дело в том, что обращение к таблице может выбросить из кэша то, что там было изначально и замедлить таким образом какой-то другой участок кода. Косяк это вы вряд ли найдете быстро, однако подобные чудовищные оптимизации едва ли кому-то понадобятся на практике при реализации обычных программ.

    Умножение и остаток от деления

    С помощью умножения и остатка от деления на степень двойки без единицы можно делать довольно интересные вещи. Обозначим все биты одного байта латинскими буквами от «a» до «h». Число n примет вид: [math]n=a\;b\;c\;d\;e\;f\;g\;h[/math] Умножение [math]n' = n\cdot 10101h[/math] делает «тиражирование» одного байта в трех экземплярах:
    +abcdefgh
    abcdefgh
    abcdefgh
    =abcdefghabcdefghabcdefgh
    Разобьем 24 бита на 8 блоков по 3 бита в каждом (нижеследующая таблица, первая строка). С помощью побитового «и» с маской 249249h (вторая строка таблицы) обнулим в каждом блоке два старших бита.
    n⋅10101habcdefghabcdefghabcdefgh
    249249hbin001001001001001001001001
    hex249249
    n' & 249249h00c00f00a00d00g00b00e00h
    Третья строка таблицы поясняет шестнадцатеричную запись маски. В последней строке — результат: все биты исходного байта содержаться каждый в своем 3-х битовом блоке, но в ином порядке.
    Мы должны сложить эти 8 блоков — и получить сумму бит в числе.
    Остаток от деления некоторого числа A на 2k-1 дает сумму k-битовых блоков числа A, тоже взятую по модулю 2k-1.
    Разобьем число A (в двоичной записи) на блоки по k бит в каждом (при необходимости можно дополнить самый последний, старший, блок нулями). Обозначим через Ai 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).
    2ik = 1 (mod (2k-1)) для любого целого неотрицательного i. (Трюк имеет смысл, когда [math]k>1[/math], иначе не совсем понятно как интерпретировать модуль 1). Мы получили сумму блоков по модулю 2k-1.
    От полученного числа нужно взять остаток от деления на 23-1 (семь) — и получаем сумму наших 8-и блоков по модулю 7. Сумма бит может быть равна 7 или 8, в таком случае алгоритм выдаст 0 и 1 соответственно. В каком случае можем получить 8? Только когда [math]n=255[/math]. В каком случае можем получить 0? Только когда [math]n=0[/math]. Если алгоритм после взятия остатка на 7 даст 0, тогда на входе получили [math]n=0[/math], либо в числе ровно 7 единичных бит. Получаем следующий код:
    Код (C):
    1. u8 CountOnes3 (u8 n) {
    2.  if (n == 0) return 0; // Единственный случай, когда ответ 0.
    3.  if (n == 0xFF) return 8; // Единственный случай, когда ответ 8.
    4.  n = (0x010101*n & 0x249249) % 7; // Считаем число бит по модулю 7.
    5.  if (n == 0) return 7; // Гарантированно имеем 7 единичных битов.
    6.  return n; // Случай, когда в числе от 1 до 6 единичных битов.
    7. }
    В случае когда у n размер 16 бит можно разбить его на две части по 8 бит.
    Код (C):
    1. u8 CountOnes3 (u16 n) {
    2.  return CountOnes3 (u8(n&0xFF)) + CountOnes3 (u8(n>>8));
    Для случая 32-х и 64-х бит подобное разбиение не имеет смысла, умножение и остаток от деления с тремя ветвлениями будут слишком дорого стоить, если выполнять их 4 или 8 раз подряд. В следующей таблице, оставлены пустые места, если вы не верите — заполните их сами. Скорее всего, будут результаты, сравнимые с процедурой CountBits1, если у вас похожий процессор (возможны оптимизации с помощью SSE).
    Режимu8u16u32u64
    x8612,4230,57
    x6413,8833,88
    Данный трюк можно сделать без ветвлений, но тогда нужно, чтобы при разбиении числа на блоки в блок вместились все числа от 0 до 8, этого можно добиться в случае 4-битовых блоков (и больше). Чтобы выполнить суммирование 4-битовых блоков, нужно подобрать множитель, который позволит правильно «растиражировать» число и взять остаток от деления на 24-1=15, чтобы сложить получившиеся блоки. Легко подобрать такой множитель: 8040201h. Почему?
    Необходимо, чтобы все биты исходного числа заняли правильные позиции в 4-битовых блоках (таблица выше), 8 и 4 не являются взаимно простыми числами, копирование 8 битов 4 раза не даст правильного расположения нужных битов. Придется добавить к байту ноль, то есть тиражировать 9 битов, так как 9 взаимно просто с 4. Получим число, имеющее размер 36 бит, в котором все биты исходного байта стоят на младших позициях 4-битовых блоков. Возьмем побитовое «и» с числом 111111111h, чтобы обнулить по три старших бита в каждом блоке. Далее блоки нужно сложить.
    0abcdefgh0abcdefgh0abcdefgh0abcdefgh
    000100010001000100010001000100010001
    000c000g000b000f000a000e0000000d000h
    Программа подсчета единичных битов в байте предельно простая:
    Код (C):
    1. u8 CountOnes3_x64 (u8 n) {
    2.  return ((u64)0x08040201*n & 0x111111111) % 15;
    3.  }
    Недостаток программы — требуется 64-битовая арифметика. Программа задействует 33 бита из 64-х (старшие 3 бита обнуляются).
    Какого размера может быть переменная n, чтобы данный трюк правильно работал? Берем остаток от деления на 15, такая переменная не может иметь размер больше 14 бит, иначе придется применить ветвление. Для 14 бит прием работает, если добавить к 14-ти битам один ноль, чтобы все биты встали на свои позиции. Подберем множитель для тиражирования и маску для обнуления ненужных битов.
    Код (C):
    1. u8 CountOnes3_x64 (u14 n) { // Это не опечатка (и не должно работать)!
    2.  return (n*0x200040008001llu & 0x111111111111111llu) % 15;
    3.  }
    Программа показывает, как мог бы выглядеть код, если бы переменная была размером 14 бит без знака. Этот код будет работать с переменной в 15 бит, при условии, что максимум 14 бит равны единице, либо если [math]n=7FFFh[/math] (этот случай разберем отдельно). Сначала «откусим» младший бит, посчитаем биты в оставшемся 15-ти битовом числе, далее обратно прибавим «откушенный» бит.
    Код (C):
    1. u8 CountOnes3_x64 (u16 n) {
    2.  u8 leastBit = n&1; // Берём младший бит.
    3.  n >>= 1; // Оставляем только 15 бит исходного числа.
    4.  if (n == 0) return leastBit;
    5. // Если получился 0, значит ответом будет младший бит.
    6.  if (n == 0x7FFF) return leastBit + 15;
    7. // Единственный случай, когда ответ 15+младший бит.
    8.  return leastBit + (n*0x200040008001llu & 0x111111111111111llu) % 15;
    9. // Ответ (максимум 14+младший бит).
    10.  }
    Для n размером 32 бита ответ влезает только в 6 бит, можно рассмотреть отдельный случай, когда n=232-1 и сделать расчеты в полях размером 5 бит. Это экономит место и позволяет разбить 32-битовое поле числа n на 3 части по 12 бит в каждой (последнее поле будет не полным). 12×5=60, можно тиражировать 12-битовое поле 5 раз, подобрав множитель, сложить 5-битовые блоки, взяв остаток от деления на 31. Это нужно делать 3 раза для каждого поля. Получаем:
    Код (C):
    1. u8 CountOnes3_x64 ( u32 n ) {
    2.  if (n == 0) return 0;
    3.  if (n + 1 == 0) return 32;
    4.  u64 res = (n&0xFFF)*0x1001001001001llu & 0x84210842108421llu;
    5.  res += ((n&0xFFF000)>>12)*0x1001001001001llu & 0x84210842108421llu;
    6.  res += (n>>24)*0x1001001001001llu & 0x84210842108421llu;
    7.  res %= 0x1F; if (res == 0) return 31;
    8.  return res;
    9.  }
    Можно вместо трех ветвлений взять 3 остатка от деления.
    Для n размером 64 бита не удалось придумать подходящее заклинание, в котором не так много умножений и сложений. Получалось либо 6, либо 7. Другой вариант — выход в 128-битовую арифметику.
    Время работы.
    Режимu8u16u32u64
    x8639,7860,48146,78
    x646,7812,2831,12
    Вывод из этой таблицы — 64-х битовая арифметика плохо воспринимается в 32-х битовом режиме исполнения. Из скорости алгоритма предподсчета в режиме x64 для однобайтовой таблицы для случая u32 (24,79 с), получается, что данный алгоритм отстает на 25%.

    Замена взятия остатка на умножение и сдвиг

    Недостаток операции взятия остатка — деление это долго. Современные компиляторы заменяют деление на умножение со сдвигом, для получения остатка, вычитают из делимого частное, умноженное на делитель. Можно суммировать k-битовые блоки не взятием остатка от деления, а еще одним умножением на маску, с помощью которой обнуляют лишние биты в блоках. Для n размером в 1 байт тиражируем байт трижды и удаляем по два старших бита у каждого 3-битового блока с помощью операций [math]10101h⋅n \& 249249h[/math].
    00c00f00a00d00g00b00e00h
    CFADGBEH
    Каждый 3-х битовый блок для обозначен заглавной латинской буквой. Умножаем полученный результат на маску 249249h. Маска содержит единичный бит в каждой 3-й позиции, поэтому умножение эквивалентно сложению числа самого с собой 8 раз со сдвигом на 3 бита:
    CFADGBEH
    +​
    CFADGBEH
    +​
    CFADGBEH
    +​
    CFADGBEH
    +​
    CFADGBEH
    +​
    CFADGBEH
    +​
    CFADGBEH
    +CFADGBEH
    Биты с 21 по 23 дают нужную сумму. Переполнение в каком-либо из блоков справа невозможно, так как ни в одном блоке нет числа больше 7. Если сумма равна 8 —получаем 0, этот единственный случай можно рассмотреть отдельно.
    Код (C):
    1. u8 CountOnes3_x64_m (u8 n) {
    2.   if (n == 0xFF)  return 8;
    3.   return ((u64(0x010101*n & 0x249249) * 0x249249) >> 21) & 0x7;
    4. }
    Взят код предыдущего раздела и в нем заменили взятие остатка от деления на 7 на умножение, сдвиг и побитовое «И». Вместо 3-х ветвлений осталось одно.
    Чтобы составить аналогичную программу для 16 бит, нужно взять код предыдущего раздела, в котором взятие остатка от деления на 15 заменяют на умножение.
    Код (C):
    1. u8 CountOnes3_x64_m (u16 n) {
    2.  u8 leastBit = n&1; // Берём младший бит
    3.  n >>= 1; // Оставляем только 15 бит.
    4.  return leastBit + (( (n*0x200040008001llu & 0x111111111111111llu)*0x111111111111111llu >> 56) & 0xF);
    5. // Ответ (максимум 15 + младший бит).
    6.  }
    Для 32-х бит делаем то же самое: берем код предыдущего раздела, соображаем, каким будет сдвиг, если заменить остаток на умножение.
    Код (C):
    1. u8 CountOnes3_x64_m ( u32 n ) {
    2.  if (n+1 == 0) return 32;
    3.  u64 res = (n&0xFFF)*0x1001001001001llu & 0x84210842108421llu;
    4.  res += ((n&0xFFF000)>>12)*0x1001001001001llu & 0x84210842108421llu;
    5.  res += (n>>24)*0x1001001001001llu & 0x84210842108421llu;
    6.  return (res*0x84210842108421llu >> 55) & 0x1F;
    7.  }
    Режимu8u16u32u64
    x8612,6642,3799,90
    x643,544,5118,35
    Удивили результаты для режима x64. Мы обогнали предподсчет с однобайтовой таблицей для случая u32. Можно ли обогнать предподсчет?

    Параллельное суммирование

    Это самый распространённый трюк, который очень повторяют, не понимая, как он точно работает.
    Начнем с 1 байта. Байт — 4-х поля по 2 бита, суммируем биты в этих полях:
    Код (C):
    1. n = (n>>1)&0x55 + (n&0x55);
    Обозначаем биты одного байта первыми латинскими буквами):
    na bc de fg h
    x=n&55h0 b0 d0 f0 h
    y=(n>>1)&55h0 a0 c0 e0 g
    x+ya+bc+de+fg+h
    Первое «И» оставляет младшие биты каждого 2-х битового блока, второе оставляет старшие биты и сдвигает их на позиции, соответствующие младшим битам. В результате суммирования получаем сумму смежных битов в каждом 2-х битовом блоке (последняя строка на таблице выше).
    Теперь сложим парами числа, находящиеся в двухбитовых полях, помещая результат в 2 четырёхбитовых поля:
    Код (C):
    1. n = (n>>2)&0x33 + (n&0x33);
     
    Marylin, Intro и mantissa нравится это.
  2. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.707
    Нижеследующая таблица поясняет результат:
    na+bc+de+fg+h
    x=n&0x330c+d0g+h
    y=(n>>2)&0x330a+b0e+f
    x+ya+b+c+de+f+g+h
    Сложим два числа в 4-х битовых полях:
    Код (C):
    1. n = (n>>4)&0x0F + (n&0x0F);
    na+b+c+de+f+g+h
    x=n&0x0F0e+f+g+h
    y=(n>>4)&0x0F0a+b+c+d
    x+ya+b+c+d+e+f+g+h
    Действуя по аналогии, можно распространить прием на любое число бит, равное степени двойки. Число строк равно двоичному логарифму от числа бит. Взгляните на 4 функции, записанных ниже, чтобы убедиться в правильности своего понимания.
    Код (C):
    1. u8 CountOnes4 (u8 n) {
    2.  n = ((n>>1) & 0x55) + (n & 0x55);
    3.  n = ((n>>2) & 0x33) + (n & 0x33);
    4.  n = ((n>>4) & 0x0F) + (n & 0x0F);
    5.  return n;
    6.  }
    7. u8 CountOnes4 (u16 n) {
    8.  n = ((n>>1) & 0x5555) + (n & 0x5555);
    9.  n = ((n>>2) & 0x3333) + (n & 0x3333);
    10.  n = ((n>>4) & 0x0F0F) + (n & 0x0F0F);
    11.  n = ((n>>8) & 0x00FF) + (n & 0x00FF);
    12.  return n;
    13.  }
    14. u8 CountOnes4 (u32 n) {
    15.  n = ((n>>1) & 0x55555555) + (n & 0x55555555);
    16.  n = ((n>>2) & 0x33333333) + (n & 0x33333333);
    17.  n = ((n>>4) & 0x0F0F0F0F) + (n & 0x0F0F0F0F);
    18.  n = ((n>>8) & 0x00FF00FF) + (n & 0x00FF00FF);
    19.  n = ((n>>16) & 0x0000FFFF) + (n & 0x0000FFFF);
    20.  return n;
    21.  }
    22. u8 CountOnes4 (u64 n) {
    23.  n = ((n>>1) & 0x5555555555555555llu) + (n & 0x5555555555555555llu);
    24.  n = ((n>>2) & 0x3333333333333333llu) + (n & 0x3333333333333333llu);
    25.  n = ((n>>4) & 0x0F0F0F0F0F0F0F0Fllu) + (n & 0x0F0F0F0F0F0F0F0Fllu);
    26.  n = ((n>>8) & 0x00FF00FF00FF00FFllu) + (n & 0x00FF00FF00FF00FFllu);
    27.  n = ((n>>16) & 0x0000FFFF0000FFFFllu) + (n & 0x0000FFFF0000FFFFllu);
    28.  n = ((n>>32) & 0x00000000FFFFFFFFllu) + (n & 0x00000000FFFFFFFFllu);
    29.  return n;
    30. }
    На этом параллельное суммирование не заканчивается. Развить идею позволяет то, что в каждой строчке дважды используется одна и та же битовая маска, что наводит на мысль «а нельзя ли как-нибудь только один раз выполнить побитовое «И»?». Можно, но не сразу. Вот что можно сделать, если взять в качестве примера код для u32 (смотрите комментарии).
    Код (C):
    1. u8 CountOnes4 (u32 n) {
    2.  n = ((n>>1) & 0x55555555) + (n & 0x55555555); // Можно заменить на разность
    3.  n = ((n>>2) & 0x33333333) + (n & 0x33333333); // Нельзя исправить
    4.  n = ((n>>4) & 0x0F0F0F0F) + (n & 0x0F0F0F0F); // Можно вынести & за скобку
    5.  n = ((n>>8) & 0x00FF00FF) + (n & 0x00FF00FF); // Можно вынести & за скобку
    6.  n = ((n>>16) & 0x0000FFFF) + (n & 0x0000FFFF); // Можно вообще убрать &
    7.  return n; // Неявное обрезание по 8-и младшим битам.
    8. }
    В качестве упражнения предлагается доказать самостоятельно, почему нижеследующий код будет точным отображением предыдущего. Для первой строки дается подсказка: 2-х битовый блок ab имеет точное значение 2a+b, значит вычитание сделает его равным…?
    Код (C):
    1. u8 CountOnes4_opt (u32 n) {
    2.  n -= (n>>1) & 0x55555555;
    3.  n = ((n>>2) & 0x33333333 ) + (n & 0x33333333);
    4.  n = ((n>>4) + n) & 0x0F0F0F0F;
    5.  n = ((n>>8) + n) & 0x00FF00FF;
    6.  n = ((n>>16) + n);
    7.  return n;
    8. }
    Аналогичные варианты оптимизации возможны и для остальных типов данных.
    Ниже приводятся две таблицы: одна для обычного параллельного суммирования, а вторая для оптимизированного.
    Режимu8u16u32u64
    x867,5214,1021,1262,70
    x648,0611,8921,3022,59
    Режимu8u16u32u64
    x867,1811,8918,8665,00
    x648,0910,2719,2019,20
    Оптимизированный алгоритм работает хорошо, но проигрывает обычному в режиме x86 для u64.

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

    Наилучшие варианты подсчета единичных битов – параллельный метод (с оптимизацией) и метод тиражирования с умножением для подсчета суммы блоков. Можно объединить оба метода, получая комбинированный алгоритм.
    Первое, что нужно сделать — выполнить первые три строки параллельного алгоритма. Это даст точную сумму битов в каждом байте числа. Для u32 выполним:
    Код (C):
    1. n -= (n>>1) & 0x55555555;
    2.  n = ((n>>2) & 0x33333333) + (n & 0x33333333);
    3.  n = (((n>>4) + n) & 0x0F0F0F0F );
    Число n состоит из 4 байт, которые рассматриваем как 4 числа, сумму которых мы ищем:
    n=ABCD
    Можно найти сумму 4-х байт, если умножить число n на 1010101h. Для удобства определения позиции, в которой будет находиться ответ, приводится таблица:
    ABCD
    +​
    ABCD
    +​
    ABCD
    +ABCD
    Ответ находится в 3-байте (если считать от 0). Комбинированный прием для u16 выглядит так:
    Код (C):
    1. u8 CountOnes5 (u16 n) {
    2.  n -= (n>>1) & 0x5555;
    3.  n = ((n>>2) & 0x3333) + (n & 0x3333);
    4.  n = ((((n>>4) + n) & 0x0F0F) * 0x0101) >> 8;
    5.  return n; // Здесь происходит неявное обрезание по 8 младшим битам.
    6. }
    Он же для u32:
    Код (C):
    1. u8 CountOnes5 ( u32 n ) {
    2.  n -= (n>>1) & 0x55555555;
    3.  n = ((n>>2) & 0x33333333 ) + (n & 0x33333333);
    4.  n = ((((n>>4) + n) & 0x0F0F0F0F) * 0x01010101) >> 24;
    5.  return n; // Здесь происходит неявное обрезание по 8 младшим битам.
    6. }
    Он же для u64:
    Код (C):
    1. u8 CountOnes5 ( u64 n ) {
    2.  n -= (n>>1) & 0x5555555555555555llu;
    3.  n = ((n>>2) & 0x3333333333333333llu ) + (n & 0x3333333333333333llu);
    4.  n = ((((n>>4) + n) & 0x0F0F0F0F0F0F0F0Fllu) * 0x0101010101010101) >> 56;
    5.  return n; // Здесь происходит неявное обрезание по 8 младшим битам.
    6. }
    Скорость работы этого метода можно посмотреть в итоговой таблице.

    Итоговое сравнение

    Читателю предлагается самостоятельно сделать выводы, из двух нижеследующих таблиц.
    Итоговое сравнения для режима компиляции x86.
    ФункцияТип n
    u8u16u32u64
    Наивный метод38,1872,00130,49384,76
    Удаление младшей единички44,7355,8872,02300,78
    Предподсчет 1 байт0,011,8321,0736,25
    Предподсчет 2 байта0,057,9513,01
    Умножение и остаток (32 бита)12,4230,57
    Умножение и остаток (64 бита)39,7860,48146,78
    Умножение и умножение (64 бита)12,6642,3799,90
    Параллельное суммирование7,5214,1021,1262,70
    Параллельное суммирование (оптимизация)7,1811,8918,8665,00
    Комбинированный метод17,6513,6052,65
    Итоговое сравнения для режима компиляции x64.
    ФункцияТип n
    u8u16u32u64
    Наивный метод37,7271,51131,47227,46
    Удаление младшей единички40,9669,1679,13126,72
    Предподсчет 1 байт0,011,4424,7926,84
    Предподсчет 2 байта0,078,4913,01
    Умножение и остаток (32 бита)13,8833,88
    Умножение и остаток (64 бита)6,7812,2813,12
    Умножение и умножение (64 бита)3,544,5118,35
    Параллельное суммирование8,0611,8921,3022,59
    Параллельное суммирование (оптимизация)8,0910,2719,2019,20
    Комбинированный метод10,5313,609,41

    Замечание

    Ни в коем случае не рассматривайте итоговую таблицу как доказательство в пользу того или иного подхода. Поверьте, что на вашем процессоре и с вашим компилятором некоторые числа в такой таблице будут совершенно иными. К сожалению, нельзя точно сказать, который из алгоритмов окажется лучше в том или ином случае. Под каждую задачу нужно затачивать конкретный метод, а универсального быстрого алгоритма не существует.
    В статье изложены те идеи, о которых автор знал, но это лишь идеи, конкретные реализации которых в разных комбинациях могут быть очень разными. Объединяя эти идеи разными способами, можно получать огромное количество разных алгоритмов подсчета единичных битов, каждый из которых вполне может оказаться хорошим в каком-то своем случае.
     
    aa_dav и Marylin нравится это.
  3. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    559
    Прикольно, но для современных камней, уже не очень актуально, SSE4.2 и SSE4A уже давно поддерживаются(от атлон/феном 2 и начиная i'шаки 1-го поколения Nehalem), а значит POPCNT присутствует, является простой командой, выполняется за 1 такт на всех портах(конвейерах), это 3 на Атлон/Феном II, и 4 на райзенах, на бульдозерах только 2, т.к. пессимизированая архитектура. Я её использовал в задаче 19 ферзей, очень хорошо ускорила задачу, ещё используется в XRay. Так же подобная команда присутствует в процессорах альфа, и в некоторых АРМ.
    --- Сообщение объединено, 14 мар 2023 ---
    Может ли эта инструкция выполнятся на всех портах intel'ах я точно не знаю, но как минимум на 4-х портах точно. Но это надо проверить. Команда реализуется каскадом простых сумматоров.
     
    Последнее редактирование: 14 мар 2023
  4. alex_dz

    alex_dz Active Member

    Публикаций:
    0
    Регистрация:
    26 июл 2006
    Сообщения:
    324
    Однако...
    POPCNT was introduced at the same time as SSE4.2, but is not a part of it. It has its own CPUID bit.
     
  5. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    559
    alex_dz, ну это правильно, т.к. SSE4.2 и SSE4A это не одно и тоже, в них только эта одна команда совпадает, так что проще проверить конкретно эту команду. У меня конечно проверяется, Кора 2 дуба всё ещё используется, так что даже в 23, нельзя 100% гарантировать что эта команда есть в целевой машине.