Алгоритмы манипуляций с битами
- воскресенье, 2 марта 2025 г. в 00:00:10
TL; DR в статье приведены алгоритмы обработки коротких битовых строк, обычно вмещающихся в машинное слово, в большей степени эти алгоритмы предназначены для обработки строк длины 32 или 64, но многие из них можно применять для SIMD инструкций или даже GPU.
Суть в двух словах
В общем то большинство рассматриваемых алгоритмов основаны на двух идеях:
Для строк длины 8, а иногда даже 16 можно просто подсчитать результат операций.
Операция над машинным словом не менее эффективна операции над отдельными битами, в связи с чем например можно применить параллельные версии алгоритмов не потеряв при этом эффективности.
Где применяется?
Вот пару примеров что нашел в open source:
Обработка позиций в шахматных движках Stockfish
Для реализации succint структур данных, в частности rank-select словарей SDSL-lite
Реализация нестандартной арифметики, например полей Галуа: вот например довольно популярный C репозиторий использует альтернативную реализацию деления на 2^k-1.
Поиск по репозиторию Linux выдаёт 11 мест применения __builtin_popcount
Преамбула
Давайте для начала зафиксируем несколько небольших вещей для лучшей структурированности и наглядности:
Все алгоритмы в статьи работают с 32-разрядными числами, я буду использовать тип uint32_t
, в большинстве случаев они очевидным образом обобщаются для длины любой степени двойки, в случае наличия нюансов я об этом явно напишу.
Если некоторая величина представляет собой позицию или что-то количественное, то она для неё будет использоваться тип size_t
, единственное исключение - байтовые таблицы, где размер будет влиять на производительность.
Я буду намеренно избегать и использовать альтернативные методы выделению меньшего значимого бита с помощью x & -x из-за того, что он предполагает конкретное представление знакового целого числа, на данный момент не являющегося единым стандартом.
Включаю в сравнение методы __builtin_popcount
, __builtin_ctz
, но понятия не имею что внутри происходит. Точно могу сказать, что в случае их доступности надо использовать их. Пока готовил статью внезапно для себя обнаружил, что они могут и не быть самыми быстрыми. UPD спасибо @slonopotamusза наводку, первичные замеры были сделаны без использования процессорных инструкций, c -march=native
__builtin функции заметно быстрее почти всех используемых в статье методов
Изначально я подготовил материал в markdown и надеялся его просто скопировать на хабр, но не тут то было. Формулы в хабравских статьях -- это боль, поэтому я оставил упрощенную запись, если встретите 2^k в тексте, то это "два в степени k", а вот в коде ^
будет обозначать битовый ксор.
Дано число, хотим развернуть порядок его бит, т.е. чтобы 31-ый бит поменялся местами с 0-ым, 30-ый с 1-ым и т.д.
Очевидный алгоритм
Для начала очевидный алгоритм выглядит вот так
uint32_t _reverse_naive(uint32_t x) {
uint32_t result = 0;
for (size_t i = 0; i < 32; ++i) {
result = (result << 1) | (x & 1);
x >>= 1;
}
return result;
}
Псевдопараллельный алгоритм
Если обозначить как R(S)
развернутую строку S
, то для конкатенации выполняется простое свойство R(AB)=R(B)R(A).
Из этого сразу следует простой рекурсивный алгоритм вида "разделяй и властвуй": разделить число на две части, развернуть половинки рекурсивно и записать их в обратном порядке. Оказывается удобно объединить все операции на определенном уровне рекурсии в одну и превратить в цикл длины 5.
Для лучшего понимания происходящего посмотрим двоичную запись используемых констант
0x55555555 -> 01010101 01010101 01010101 01010101
0xAAAAAAAA -> 10101010 10101010 10101010 10101010
0x33333333 -> 00110011 00110011 00110011 00110011
0xCCCCCCCC -> 11001100 11001100 11001100 11001100
0x0F0F0F0F -> 00001111 00001111 00001111 00001111
0xF0F0F0F0 -> 11110000 11110000 11110000 11110000
0x00FF00FF -> 00000000 11111111 00000000 11111111
0xFF00FF00 -> 11111111 00000000 11111111 00000000
0x0000FFFF -> 00000000 00000000 11111111 11111111
0xFFFF0000 -> 11111111 11111111 00000000 00000000
По сути вышеописанный код выполняет рекурсию в обратном порядке: на первой итерации меняются все четные и нечетные биты, на втором -- блоки по два бита, и т.д. Код для генерации таких констант для чисел длины 2^n выглядит так
std::vector<uint32_t> left_mask;
std::vector<uint32_t> right_mask;
for (size_t j = 0; j < n; ++j) {
left_mask.emplace_back(0);
for (size_t i = 0; i < (1 << n); ++i) {
left_mask[j] |= (1ULL << i) * (((i / (1 << j)) & 1) == 0);
}
right_mask.emplace_back(~left_mask[j]);
}
В общем виде такой алгоритм имел бы сложность O(n log n
).
Табличный метод
Результат разворота для всех возможных значений байтов можно предподсчитать
static const unsigned char BitReverseTable256[256] = {
# define R2(n) n, n + 2*64, n + 1*64, n + 3*64
# define R4(n) R2(n), R2(n + 2*16), R2(n + 1*16), R2(n + 3*16)
# define R6(n) R4(n), R4(n + 2*4 ), R4(n + 1*4 ), R4(n + 3*4 )
R6(0), R6(2), R6(1), R6(3)
};
Остаётся только разделить 32-разрядное число на 4 байта, применить по отдельности и собрать в обратном порядке
uint32_t _reverse_lookup(uint32_t x) {
return (BitReverseTable256[x & 0xff] << 24) |
(BitReverseTable256[(x >> 8) & 0xff] << 16) |
(BitReverseTable256[(x >> 16) & 0xff] << 8) |
(BitReverseTable256[(x >> 24) & 0xff]);
}
Есть число, мы хотим посчитать количество единичных бит в нём.
Очевидные методы
Самый простой способ - перебрать все биты:
size_t _count_bits_set_naive(uint32_t x) {
size_t result = 0;
while (x > 0) {
result += x & 1;
x >>= 1;
}
return result;
}
Можно сделать аналогичный перебор только единичных битов, здесь мы впервые встречаемся с трюком, связанным с выделением наименьшего бита числа. Идея в том, что у чисел x
и x-1
все биты после младшего единичного бита x
совпадают, младший единичный бит x
превращается в ноль у x-1
, а биты до него являются нулями у x
и единичками у x-1
x 10010101101000
x-1 10010101100111
Таким образом Выражение x ^ (x & (x-1)) оставляет только младший единичный бит x
. То же самое можно получить более простым выражением x & -x в случае если одно и то же битовое представление знакового и беззнаковых типов совпадает по модулю 2^{32}.
Так или иначе получаем алгоритм для подсчета единичным бит последовательным их перебором.
size_t _count_bits_set_iterate_set_bits(uint32_t x) {
size_t result = 0;
while (x > 0) {
++result;
x = x & (x - 1);
}
return result;
}
Псевдопараллельный метод
Общая идея схожа с разворотом: разделяем по уровням, на i-ом уровне что-то делаем с блоками длины 2^i. В случае с разворотом мы брали два блока и меняли их местами, в случае с подсчетом битов будем записывать в блок 2^{i+1} сумму битов в этом блоке предполагая, что на предыдущем шаге мы записали в каждый блок 2^i его сумму. Ну и опять же все это объединяется в одну операцию с числом
size_t _count_bits_set_parallel(uint32_t x) {
uint32_t result = x - ((x >> 1) & 0x55555555);
result = ((result >> 2) & 0x33333333) + (result & 0x33333333);
result = ((result >> 4) + result) & 0x0F0F0F0F;
result = ((result >> 8) + result) & 0x00FF00FF;
result = ((result >> 16) + result) & 0x0000FFFF;
return result;
}
Операции 3, 4, 5 идентичны. Операции 1, 2 немного изменены для избежания проблем с переполнением блока.
Табличный метод
Как и с разворотом числа можно предподсчитать значения для байтов
static const unsigned char BitsSetTable256[256] = {
# define B2(n) n, n+1, n+1, n+2
# define B4(n) B2(n), B2(n+1), B2(n+1), B2(n+2)
# define B6(n) B4(n), B4(n+1), B4(n+1), B4(n+2)
B6(0), B6(1), B6(1), B6(2)
};
и использовать соответственно так
size_t _count_bits_set_lookup(uint32_t x) {
return BitsSetTable256[x & 0xff] +
BitsSetTable256[(x >> 8) & 0xff] +
BitsSetTable256[(x >> 16) & 0xff] + BitsSetTable256[x >> 24];
}
Подсчет до заданной позиции
Для rank-select словарей возникает задача подсчёта количества единичных битов до некоторой позиции pos
, довольно легко это свести обнулив старшие биты маской (1 << pos) - 1 с небольшой модификацией для учета максимального значения для pos
size_t count_prefix_bits_set(uint64_t x, size_t pos) {
uint32_t mask = (pos > 0) * (1 + (((1 << (pos - 1)) - 1) << 1));
return count_bits_set(x & mask);
}
Задача нахождения позиции самого младшего единичного бита (least significant bit, bit forward scan, count trailing zeros).
LSB тривиальный подход
Самое простое, что можно сделать для нахождения LSB - просто пройтись циклом
size_t _least_significant_bit_naive(uint32_t x) {
for (size_t i = 0; i < 32; ++i) {
if (x & (1 << i)) {
return i;
}
}
return 0;
}
LSB бинарный поиск
В общем то ничего интересного, развернутый бинарный поиск длины 5
size_t _least_significant_bit_binary_search(uint32_t x) {
size_t result = 1;
if ((x & 0xFFFF) == 0) {
x >>= 16;
result += 16;
}
if ((x & 0xFF) == 0) {
x >>= 8;
result += 8;
}
if ((x & 0xF) == 0) {
x >>= 4;
result += 4;
}
if ((x & 0x3) == 0) {
x >>= 2;
result += 2;
}
result -= x & 0x1;
return result;
}
LSB параллельный вариант
Очень схожий вариант с бинарным поиском, но без необходимости постоянной модификации исходного числа за счёт псевдопараллельности c использованием констант, которые мы уже видели ранее.
size_t _least_significant_bit_parallel(uint32_t x) {
size_t result = 31;
x ^= x & (x - 1);
if (x & 0x0000FFFF)
result -= 16;
if (x & 0x00FF00FF)
result -= 8;
if (x & 0x0F0F0F0F)
result -= 4;
if (x & 0x33333333)
result -= 2;
if (x & 0x55555555)
result -= 1;
return result;
}
Одну модификацию исходного числа все-таки нужно произвести - изоляция LSB второй операцией функции.
LSB через конвертацию в float
Довольно интересный вариант, полезен в случае наличия эффективной машинной ковертации целочисленного типа в float. Идея довольно простая: стандартное представление float IEEE 754 как известно выглядит как (-1)^s 2^{e-127} 1.m , где s - знаковый бит, e - экспонента длины 8, m - мантиса длины 23. Если перевести число 2^i к такому виду, то мы получим s=m=0, e-127=i, остаётся только выделить биты экспоненты (c 24-ого по 31-ый).
size_t _least_significant_bit_via_float(uint32_t x) {
float y = x & ~(x - 1);
return (*(uint32_t *)&y >> 23) - 0x7F;
}
LSB с помощью байтовой таблицы
Я думаю, суть понятна из подзаголовка, так что просто приведу код генерации таблицы
static const uint8_t LeastSignificantBitTable256[256] = {
#define LSBT0 0, 1, 0
#define LSBT1 LSBT0, 2, LSBT0
#define LSBT2 LSBT1, 3, LSBT1
#define LSBT3 LSBT2, 4, LSBT2
#define LSBT4 LSBT3, 5, LSBT3
#define LSBT5 LSBT4, 6, LSBT4
0, LSBT5, 7, LSBT5};
Использование прямолинейно: разбиваем число на байты и проверяем по очереди есть ли хотя бы одна единичка, если есть - берем результат из таблицы.
size_t _least_significant_bit_lookup(uint32_t x) {
size_t result = 0;
if (x & 0xFF)
return result + LeastSignificantBitTable256[x & 0xFF];
x >>= 8;
result += 8;
if (x & 0xFF)
return result + LeastSignificantBitTable256[x & 0xFF];
x >>= 8;
result += 8;
if (x & 0xFF)
return result + LeastSignificantBitTable256[x & 0xFF];
x >>= 8;
result += 8;
if (x)
return result + LeastSignificantBitTable256[x];
return 32;
}
LSB с помощью компактной таблицы и модульной арифметики
В нашем диапазоне есть всего 32 степени двойки, по идее все эти 32 значения можно предподсчитать в компактную таблицу, но не совсем понятно как это эффективно сделать. Один из способов - найти такой модуль, остатки по которому от этих 32 степеней двойки будут различными. Минимальный такой модуль - 37.
static const uint8_t LeastSignificantBitMod37[] = {
32, 0, 1, 26, 2, 23, 27, 0, 3, 16, 24, 30, 28, 11, 0, 13, 4, 7, 17,
0, 25, 22, 31, 15, 29, 10, 12, 6, 0, 21, 14, 9, 5, 20, 8, 19, 18};
size_t _least_significant_bit_mod_lookup(uint32_t x) {
return LeastSignificantBitMod37[(x & ~(x - 1)) % 37];
}
В этой таблице есть несколько неиспользованных значений, используется операция деления с остатком, которая обычно более вычислительно затратная, чем другие арифметические действия, впрочем очень похоже, что в современных компиляторах это не так.
LSB через последовательность Де Брюина
Последний и самый эффективный табличный способ получения LSB. Идея следующая: умножение числа на степень двойки соответствует битовому сдвигу, существует всего 32 битовых последовательности длины 5 и если бы мы нашли такое 32-разрядное число, у которого любые 5 подряд идущих битов были бы разные, то домножая на такое числа степень двойки мы бы получали различные битовые строки на конкретных позициях. Таким числом например является
0x077CB531U = 00000111 01111100 10110101 00110001
и оно может быть использовано для получения LSB следующим образом
static const uint8_t LeastSignificantBitDeBruijn[] = {
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};
size_t _least_significant_bit_de_bruijn(uint32_t x) {
return LeastSignificantBitDeBruijn[((x & ~(x - 1)) * 0x077CB531U) >> 27];
}
Вообще говоря, найти такую константу для 32 можно просто перебором, но для более длинных чисел это будет очень долго. Быстрый способ получить константу для чисел длины 2^k - построить эйлеров цикл на графе с вершинами 0, 1, ..., 2^{k-1} и переходами x -> 2x mod 2^{k-1}, x -> 2x+1 mod 2^{k-1}. Первый переход соответствует дописывание нуля в начало и отбрасывание k-ого бита, второй - дописывание единички.
В прошлом разделе мы обсуждали алгоритмы нахождения наименьшего единичного бита, а что если нам нужно находить позицию i-ого единичного бита? Эта операция является базовой для rank-select словарей и succint структур данных в целом, обычно обозначается select.
Наивная реализация алгоритма
size_t _select_naive(uint32_t x, size_t rank) {
if (rank == 0) {
return 32;
}
size_t result = 0;
while (x > 0 && rank > 0) {
if (x & 1)
--rank;
x >>= 1;
++result;
}
if (x == 0 && rank > 0) {
return 32;
}
return result - 1;
}
Улучшенная реализация -- комбинация параллельного подсчета и параллельного бинарного поиска:
Можно использовать бинарный поиск для поиска целевой позиции: на каждой итерации подсчитываем количество единичных бит до позиции, запрашиваемой на этой итерации.
Вместо того, чтобы делать вычисления количества единичных бит на каждой итерации, можно сделать это один раз в параллельном формате с использованием уже знакомых ранее
0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF, 0x0000FFFF
size_t _select_parallel(uint32_t x, size_t rank) {
// parallel rank
uint32_t s1 = ((x >> 1) & BitSetParallelMasks[0]) + (x & BitSetParallelMasks[0]);
uint32_t s2 = ((s1 >> 2) & BitSetParallelMasks[1]) + (s1 & BitSetParallelMasks[1]);
uint32_t s3 = ((s2 >> 4) & BitSetParallelMasks[2]) + (s2 & BitSetParallelMasks[2]);
uint32_t s4 = ((s3 >> 8) & BitSetParallelMasks[3]) + (s3 & BitSetParallelMasks[3]);
uint32_t s5 = ((s4 >> 16) & BitSetParallelMasks[4]) + (s4 & BitSetParallelMasks[4]);
// the next code selects rank most significant bit with 0-based indexing
// redefine rank to find least significant bit.
if (rank > s5) {
return 32;
}
rank = s5 - rank + 1;
// parallel binary search
size_t t = (s4 >> 16);
size_t s = 32;
if (rank > t) { s -= 16; rank -= t; }
t = (s3 >> (s - 8)) & 0xf;
if (rank > t) { s -= 8; rank -= t; }
t = (s2 >> (s - 4)) & 0x7;
if (rank > t) { s -= 4; rank -= t; }
t = (s1 >> (s - 2)) & 0x3;
if (rank > t) { s -= 2; rank -= t; }
t = (x >> (s - 1)) & 0x1;
if (rank > t) s--;
return s - 1;
}
В этом разделе пытаемся быстро делить на 2^k-1. Не секрет, что деление с остатком на 2^k выражается через бинарные операции
x % (1 << k) == x & ((1 << k) - 1)
x / (1 << k) == x >> k
деление на 2^k-1 так красиво не получается, но всё-таки можно что-то сделать. Такое деление актуально например для арифметики в полях Галуа.
Наивная реализация
Для полноты картины (и бенчмарков) привожу очевидное
uint32_t _mod_by_power_of_two_minus_one_naive(uint32_t x, size_t k) {
return x % ((1 << k) - 1);
}
Две реализации с помощью последовательного деления на 2^k
Битовым делением на 2^k можно воспользоваться с помощью следующего равенства
Соответственно вычислить остаток при делении на 2^k-1 можно вычислить последовательным делением на 2^k
uint32_t _mod_by_power_of_two_minus_one_no_division_v1(uint32_t x, size_t k) {
uint32_t m = (1 << k) - 1;
while (x >= m) {
x -= m;
x = (x >> k) + (x & m);
}
return x;
}
Можно заметить, что на каждой итерации x отбрасывается хотя бы k бит и соответственно количество итераций этого цикла не превосходит (log x)/k.
Эта идея в целом схожа со свойством, что остаток при делении числа на 9 такой же, как и у суммы его десятичных цифр. Перейдя от числу к сумме его цифр несколько раз пока цифра не станет одна мы получаем остаток при делении на 9. Примерно тоже самое происходит в следующем коде, только по основанию 2^k
uint32_t _mod_by_power_of_two_minus_one_no_division_v2(uint32_t x, size_t k) {
uint32_t result = x, d = (1 << k) - 1;
for (; x > d; x = result) {
for (result = 0; x; x >>= k) {
result += x & d;
}
}
return result == d ? 0 : result;
}
Табличный параллельный метод
Наконец какой-то странный параллельный метод взятый отсюда, который я плохо понимаю, судя по всему работает не очень быстро. Привожу для полноты картины
Проверял в трех сценариях
* Один аргумент, все $32$-битные инты
* Один аргумент, все $32$-битные инты, у которых только один единичный бит
* Два аргумента, первый -- все 24-битные инты, второй -- число в диапазоне 0, 1, ..., 23.
Весь код и замеры можно найти в репозитории
1: -------------------------------------------------------------------------------------------
1: Benchmark Time CPU Iterations
1: -------------------------------------------------------------------------------------------
1: All 32-bit uint/Count bits set naive 45524 ms 45522 ms 1
1: All 32-bit uint/Count bits set iterate set bits 27579 ms 27578 ms 1
1: All 32-bit uint/Count bits set in parallel 6770 ms 6770 ms 1
1: All 32-bit uint/Count bits set with lookup table 4029 ms 4029 ms 1
1: All 32-bit uint/Count bits set with builtin 13350 ms 13349 ms 1
1: All 32-bit uint/Reverse naive 69577 ms 69573 ms 1
1: All 32-bit uint/Reverse parallel 6012 ms 6011 ms 1
1: All 32-bit uint/Reverse lookup 4011 ms 4011 ms 1
1: All 24-bit uint/Select naive 8310 ms 8310 ms 1
1: All 24-bit uint/Select parallel 973 ms 973 ms 1
1: All 32-bit uint/LSB naive 2763 ms 2763 ms 1
1: All 32-bit uint/LSB parallel 7654 ms 7654 ms 1
1: All 32-bit uint/LSB float 2767 ms 2767 ms 1
1: All 32-bit uint/LSB lookup 2767 ms 2767 ms 1
1: All 32-bit uint/LSB mod lookup 5015 ms 5014 ms 1
1: All 32-bit uint/LSB De Bruijn 1338 ms 1338 ms 1
1: All 32-bit uint/LSB builtin ctz 1337 ms 1337 ms 1
1: All 32-bit single bit uint/LSB naive 179 ns 179 ns 3893997
1: All 32-bit single bituint/LSB parallel 45.0 ns 45.0 ns 15452914
1: All 32-bit single bituint/LSB float 30.8 ns 30.8 ns 22331270
1: All 32-bit single bit uint/LSB lookup 29.6 ns 29.6 ns 23498091
1: All 32-bit single bituint/LSB mod lookup 41.0 ns 41.0 ns 17077605
1: All 32-bit single bituint/LSB De Bruijn 20.3 ns 20.3 ns 34526862
1: All 32-bit single bituint/LSB builtin ctz 20.3 ns 20.3 ns 34507293
1: All 24-bit uint/Division by 2^k-1 naive 969 ms 969 ms 1
1: All 24-bit uint/Division by 2^k-1 cycle v1 1161 ms 1161 ms 1
1: All 24-bit uint/Division by 2^k-1 cycle v2 1876 ms 1876 ms 1
1: All 24-bit uint/Division by 2^k-1 lookup 3278 ms 3278 ms 1
UPD замеры с "-march=native", __builtin функции становятся в разы быстрее. Забавно, что весь код _count_bits_set_iterate_set_bits
в этом случае оптимизируется в popcnt, а вот De Brujin оказывается таким же эффективным как ctz
1: -------------------------------------------------------------------------------------------
1: Benchmark Time CPU Iterations
1: -------------------------------------------------------------------------------------------
1: All 32-bit uint/Count bits set naive 46854 ms 46853 ms 1
1: All 32-bit uint/Count bits set iterate set bits 1335 ms 1335 ms 1
1: All 32-bit uint/Count bits set in parallel 6514 ms 6513 ms 1
1: All 32-bit uint/Count bits set with lookup table 4030 ms 4030 ms 1
1: All 32-bit uint/Count bits set with builtin 1336 ms 1336 ms 1
1: All 32-bit uint/Reverse naive 69520 ms 69518 ms 1
1: All 32-bit uint/Reverse parallel 6229 ms 6229 ms 1
1: All 32-bit uint/Reverse lookup 4018 ms 4018 ms 1
1: All 24-bit uint/Select naive 8137 ms 8136 ms 1
1: All 24-bit uint/Select parallel 2000 ms 2000 ms 1
1: All 32-bit uint/LSB naive 2568 ms 2568 ms 1
1: All 32-bit uint/LSB parallel 7366 ms 7366 ms 1
1: All 32-bit uint/LSB float 2954 ms 2954 ms 1
1: All 32-bit uint/LSB lookup 2761 ms 2761 ms 1
1: All 32-bit uint/LSB mod lookup 8002 ms 8002 ms 1
1: All 32-bit uint/LSB De Bruijn 1346 ms 1346 ms 1
1: All 32-bit uint/LSB builtin ctz 1338 ms 1338 ms 1
1: All 32-bit single bit uint/LSB naive 176 ns 176 ns 3998695
1: All 32-bit single bituint/LSB parallel 49.3 ns 49.3 ns 14127259
1: All 32-bit single bituint/LSB float 27.3 ns 27.3 ns 25716520
1: All 32-bit single bit uint/LSB lookup 31.5 ns 31.5 ns 22298392
1: All 32-bit single bituint/LSB mod lookup 57.8 ns 57.8 ns 12121144
1: All 32-bit single bituint/LSB De Bruijn 17.4 ns 17.4 ns 40198935
1: All 32-bit single bituint/LSB builtin ctz 17.5 ns 17.5 ns 40201840
1: All 24-bit uint/Division by 2^k-1 naive 969 ms 969 ms 1
1: All 24-bit uint/Division by 2^k-1 cycle v1 1145 ms 1145 ms 1
1: All 24-bit uint/Division by 2^k-1 cycle v2 1454 ms 1454 ms 1
1: All 24-bit uint/Division by 2^k-1 lookup 3261 ms 3261 ms 1
1/1 Test #1: Benchmarks ....................... Passed 189.91 sec
http://www-graphics.stanford.edu/~seander/bithacks.html
Donald E. Knuth's "The Art of Computer Programming", Volume 4A, Chapter 7.1.3