Быстрые тесты на делимость (на 2,3,4,5,.., 16)?
Каковы самые быстрые тесты делимости? Скажем, учитывая архитектуру с прямым порядком байтов и 32-разрядное целое число со знаком: как очень быстро вычислить, что число делится на 2,3,4,5,... до 16?
ВНИМАНИЕ: данный код только ПРИМЕР. Каждая строка независима! Просто очевидное решение, использующее работу по модулю, является медленным на многих процессорах, которые не имеют аппаратного обеспечения DIV (как и многие ARM). Некоторые компиляторы также не могут выполнять такие оптимизации (скажем, если делитель является аргументом функции или зависит от чего-либо).
Divisible_by_1 = do();
Divisible_by_2 = if (!(number & 1)) do();
Divisible_by_3 = ?
Divisible_by_4 = ?
Divisible_by_5 = ?
Divisible_by_6 = ?
Divisible_by_7 = ?
Divisible_by_8 = ?
Divisible_by_9 = ?
Divisible_by_10 = ?
Divisible_by_11 = ?
Divisible_by_12 = ?
Divisible_by_13 = ?
Divisible_by_14 = ?
Divisible_by_15 = ?
Divisible_by_16 = if(!number & 0x0000000F) do();
и особые случаи:
Divisible_by_2k = if(number & (tk-1)) do(); //tk=2**k=(2*2*2*...) k times
16 ответов
В общем, неплохая идея - найти альтернативы инструкциям деления (включая модуль по x86/x64), потому что они очень медленные. Медленнее (или даже намного медленнее), чем думает большинство людей. Те, кто предлагает "% n", где n - это переменная, дают глупый совет, потому что это неизменно приведет к использованию инструкции деления. С другой стороны, "% c" (где c - константа) позволит компилятору определить лучший алгоритм, доступный в его репертуаре. Иногда это будет инструкция деления, но в большинстве случаев это не так.
В этом документе Torbjörn Granlund показывает, что соотношение тактов, необходимых для 32-разрядных битов без знака:divs, составляет 4:26 (6,5x) для Sandybridge и 3:45 (15x) для K10. для 64-битных соответствующие соотношения составляют 4:92 (23x) и 5:77 (14,4x).
Столбцы "L" обозначают задержку. "T" столбцы обозначают пропускную способность. Это связано со способностью процессора обрабатывать несколько команд параллельно. Sandybridge может выдавать одно 32-битное умножение каждый цикл или один 64-битный цикл. Для K10 соответствующая пропускная способность меняется на противоположную. Для делений K10 должен завершить всю последовательность, прежде чем он может начать другую. Я подозреваю, что это то же самое для Sandybridge.
Использование K10 в качестве примера означает, что во время циклов, требуемых для 32-разрядного деления (45), может быть выдано одинаковое число (45) умножений, и последним и последним из последних будет выполнено одно и два тактовые циклы после завершения деления. МНОГО работы может быть выполнено в 45 умножениях.
Также интересно отметить, что div стали менее эффективными с переходом от K8-K9 к K10: от 39 до 45 и от 71 до 77 тактов для 32- и 64-разрядных.
Страница Гранлунда на gmplib.org и в Королевском технологическом институте в Стокгольме содержит больше полезностей, некоторые из которых были включены в компилятор gcc.
В каждом случае (в том числе делится на 2):
if (number % n == 0) do();
Обработка маски младших разрядов - это просто запутывание, а современный компилятор не будет быстрее, чем написание кода в удобочитаемой форме.
Если вам нужно протестировать все случаи, вы можете улучшить производительность, поместив некоторые случаи в if
для другого: нет смысла проверять делимость на 4, если, например, делимость на 2 уже не удалась.
Как упомянул @James, позвольте компилятору упростить его для вас. Если n
является константой, любой спускаемый компилятор может распознать шаблон и изменить его на более эффективный эквивалент.
Например, код
#include <stdio.h>
int main() {
size_t x;
scanf("%u\n", &x);
__asm__ volatile ("nop;nop;nop;nop;nop;");
const char* volatile foo = (x%3 == 0) ? "yes" : "no";
__asm__ volatile ("nop;nop;nop;nop;nop;");
printf("%s\n", foo);
return 0;
}
составлен с g++-4.5 -O3, соответствующая часть x%3 == 0
станет
mov rcx,QWORD PTR [rbp-0x8] # rbp-0x8 = &x
mov rdx,0xaaaaaaaaaaaaaaab
mov rax,rcx
mul rdx
lea rax,"yes"
shr rdx,1
lea rdx,[rdx+rdx*2]
cmp rcx,rdx
lea rdx,"no"
cmovne rax,rdx
mov QWORD PTR [rbp-0x10],rax
что в переводе на C-код означает
(hi64bit(x * 0xaaaaaaaaaaaaaaab) / 2) * 3 == x ? "yes" : "no"
// equivalatent to: x % 3 == 0 ? "yes" : "no"
здесь нет разделения (Обратите внимание, что 0xaaaaaaaaaaaaaaab == 0x20000000000000001L/3
)
Редактировать:
- Магическая постоянная 0xaaaaaaaaaaaaaaab может быть вычислена в http://www.hackersdelight.org/magic.htm
- Для делителей вида 2n - 1 проверьте http://graphics.stanford.edu/~seander/bithacks.html
Предполагать number
является unsigned
(32 бит). Затем следующие очень быстрые способы вычисления делимости до 16. (Я не измерял, но код сборки указывает на это.)
bool divisible_by_2 = number % 2 == 0;
bool divisible_by_3 = number * 2863311531u <= 1431655765u;
bool divisible_by_4 = number % 4 == 0;
bool divisible_by_5 = number * 3435973837u <= 858993459u;
bool divisible_by_6 = divisible_by_2 && divisible_by_3;
bool divisible_by_7 = number * 3067833783u <= 613566756u;
bool divisible_by_8 = number % 8 == 0;
bool divisible_by_9 = number * 954437177u <= 477218588u;
bool divisible_by_10 = divisible_by_2 && divisible_by_5;
bool divisible_by_11 = number * 3123612579u <= 390451572u;
bool divisible_by_12 = divisible_by_3 && divisible_by_4;
bool divisible_by_13 = number * 3303820997u <= 330382099u;
bool divisible_by_14 = divisible_by_2 && divisible_by_7;
bool divisible_by_15 = number * 4008636143u <= 286331153u;
bool divisible_by_16 = number % 16 == 0;
Относительно делимости на d
действуют следующие правила:
когда
d
это степень 2:Как отметил James Kanze, вы можете использовать
is_divisible_by_d = (number % d == 0)
, Компиляторы достаточно умны, чтобы заменить это(number & (d - 1)) == 0
который более эффективен, но также и запутан.Однако когда
d
это не степень 2, похоже, что запутывание, показанное выше, более эффективно, чем то, что делают современные компиляторы. (Подробнее об этом позже).когда
d
странныйТехника принимает форму
is_divisible_by_d = number * a <= b
гдеa
а такжеb
умно полученные константы. Обратите внимание, что все, что нам нужно, это 1 умножение и 1 сравнение:когда
d
четное, но не степень 2:Затем написать
d = p * q
гдеp
это степень 2 иq
странно и использовать "язык в щеке", предложенный unpythonic, то естьis_divisible_by_d = is_divisible_by_p && is_divisible_by_q
, Опять же, только 1 умножение (в расчетеis_divisible_by_q
) выполняется.
Многие компиляторы (я тестировал clang 5.0.0, gcc 7.3, icc 18 и msvc 19, используя godbolt) заменяют number % d == 0
от (number / d) * d == number
, Они используют умную технику (см. Ссылки в Olof Forshell), чтобы заменить деление умножением и сдвигом битов. Они заканчивают тем, что делали 2 умножения. В отличие от вышеописанных методик выполняют только 1 умножение.
Обновление 01 октября 2018
Похоже, что приведенный выше алгоритм скоро поступит в GCC (уже в транке):
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82853
Внедрение GCC кажется еще более эффективным. Действительно, реализация выше имеет три части: 1) делимость на четную часть делителя; 2) делимость на нечетную часть делителя; 3) &&
связать результаты двух предыдущих шагов. Используя инструкцию на ассемблере, которая не всегда доступна в стандарте C++ (ror
), GCC объединяет три части в одну, что очень похоже на разделение на нечетную часть. Качественный товар! Имея эту реализацию доступной, лучше (для ясности и производительности) вернуться к %
все время.
Немного в щеке, но при условии, что вы получите остальные ответы:
Divisible_by_6 = Divisible_by_3 && Divisible_by_2;
Divisible_by_10 = Divisible_by_5 && Divisible_by_2;
Divisible_by_12 = Divisible_by_4 && Divisible_by_3;
Divisible_by_14 = Divisible_by_7 && Divisible_by_2;
Divisible_by_15 = Divisible_by_5 && Divisible_by_3;
Прежде всего, напомню, что число в виде bn...b2b1b0 в двоичном виде имеет значение:
number = bn*2^n+...+b2*4+b1*2+b0
Теперь, когда вы говорите номер%3, у вас есть:
number%3 =3= bn*(2^n % 3)+...+b2*1+b1*2+b0
(Я использовал =3= для обозначения конгруэнтности по модулю 3). Обратите внимание, что b1*2 =3= -b1*1
Теперь я напишу все 16 делений, используя + и - и, возможно, умножение (обратите внимание, что умножение может быть записано как сдвиг или сумма одного и того же значения сдвинута в разные места. Например, 5*x
средства x+(x<<2)
в котором вы вычисляете x
только однажды)
Давайте назовем номер n
и скажем Divisible_by_i
является логическим значением В качестве промежуточного значения представьте Congruence_by_i
это значение, соответствующее n
по модулю i
,
Кроме того, скажем n0
означает бит ноль n, n1
означает бит 1 и т. д., то есть
ni = (n >> i) & 1;
Congruence_by_1 = 0
Congruence_by_2 = n&0x1
Congruence_by_3 = n0-n1+n2-n3+n4-n5+n6-n7+n8-n9+n10-n11+n12-n13+n14-n15+n16-n17+n18-n19+n20-n21+n22-n23+n24-n25+n26-n27+n28-n29+n30-n31
Congruence_by_4 = n&0x3
Congruence_by_5 = n0+2*n1-n2-2*n3+n4+2*n5-n6-2*n7+n8+2*n9-n10-2*n11+n12+2*n13-n14-2*n15+n16+2*n17-n18-2*n19+n20+2*n21-n22-2*n23+n24+2*n25-n26-2*n27+n28+2*n29-n30-2*n31
Congruence_by_7 = n0+2*n1+4*n2+n3+2*n4+4*n5+n6+2*n7+4*n8+n9+2*n10+4*n11+n12+2*n13+4*n14+n15+2*n16+4*n17+n18+2*n19+4*n20+n21+2*n22+4*n23+n24+2*n25+4*n26+n27+2*n28+4*n29+n30+2*n31
Congruence_by_8 = n&0x7
Congruence_by_9 = n0+2*n1+4*n2-n3-2*n4-4*n5+n6+2*n7+4*n8-n9-2*n10-4*n11+n12+2*n13+4*n14-n15-2*n16-4*n17+n18+2*n19+4*n20-n21-2*n22-4*n23+n24+2*n25+4*n26-n27-2*n28-4*n29+n30+2*n31
Congruence_by_11 = n0+2*n1+4*n2+8*n3+5*n4-n5-2*n6-4*n7-8*n8-5*n9+n10+2*n11+4*n12+8*n13+5*n14-n15-2*n16-4*n17-8*n18-5*n19+n20+2*n21+4*n22+8*n23+5*n24-n25-2*n26-4*n27-8*n28-5*n29+n30+2*n31
Congruence_by_13 = n0+2*n1+4*n2+8*n3+3*n4+6*n5-n6-2*n7-4*n8-8*n9-3*n10-6*n11+n12+2*n13+4*n14+8*n15+3*n16+6*n17-n18-2*n19-4*n20-8*n21-3*n22-6*n3+n24+2*n25+4*n26+8*n27+3*n28+6*n29-n30-2*n31
Congruence_by_16 = n&0xF
Или когда факторизовано:
Congruence_by_1 = 0
Congruence_by_2 = n&0x1
Congruence_by_3 = (n0+n2+n4+n6+n8+n10+n12+n14+n16+n18+n20+n22+n24+n26+n28+n30)-(n1+n3+n5+n7+n9+n11+n13+n15+n17+n19+n21+n23+n25+n27+n29+n31)
Congruence_by_4 = n&0x3
Congruence_by_5 = n0+n4+n8+n12+n16+n20+n24+n28-(n2+n6+n10+n14+n18+n22+n26+n30)+2*(n1+n5+n9+n13+n17+n21+n25+n29-(n3+n7+n11+n15+n19+n23+n27+n31))
Congruence_by_7 = n0+n3+n6+n9+n12+n15+n18+n21+n24+n27+n30+2*(n1+n4+n7+n10+n13+n16+n19+n22+n25+n28+n31)+4*(n2+n5+n8+n11+n14+n17+n20+n23+n26+n29)
Congruence_by_8 = n&0x7
Congruence_by_9 = n0+n6+n12+n18+n24+n30-(n3+n9+n15+n21+n27)+2*(n1+n7+n13+n19+n25+n31-(n4+n10+n16+n22+n28))+4*(n2+n8+n14+n20+n26-(n5+n11+n17+n23+n29))
// and so on
Если эти значения оказываются отрицательными, добавьте i
пока они не станут положительными.
Теперь то, что вы должны сделать, это рекурсивно передать эти значения через тот же процесс, который мы только что сделали, пока Congruence_by_i
становится меньше чем i
(и, очевидно, >= 0
). Это похоже на то, что мы делаем, когда хотим найти остаток числа на 3 или 9, помните? Суммируйте цифры, если в них было больше одной цифры, несколько раз суммируйте цифры результата, пока не получите только одну цифру.
Теперь для i = 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 16
:
Divisible_by_i = (Congruence_by_i == 0);
А в остальном:
Divisible_by_6 = Divisible_by_3 && Divisible_by_2;
Divisible_by_10 = Divisible_by_5 && Divisible_by_2;
Divisible_by_12 = Divisible_by_4 && Divisible_by_3;
Divisible_by_14 = Divisible_by_7 && Divisible_by_2;
Divisible_by_15 = Divisible_by_5 && Divisible_by_3;
Изменить: Обратите внимание, что некоторые из дополнений можно было избежать с самого начала. Например n0+2*n1+4*n2
такой же как n&0x7
, так же n3+2*n4+4*n5
является (n>>3)&0x7
и поэтому с каждой формулой вам не нужно получать каждый бит по отдельности, я написал это так для ясности и схожести в работе. Чтобы оптимизировать каждую из формул, вы должны поработать над ней самостоятельно; групповые операнды и факторизация операций.
LCM этих чисел, по-видимому, равен 720720. Он довольно мал, так что вы можете выполнить одну операцию модуля и использовать остаток в качестве индекса в предварительно вычисленном LUT.
Вы должны просто использовать (i % N) == 0 в качестве теста.
Мой компилятор (довольно старая версия gcc) генерировал хороший код для всех случаев, которые я пробовал. Там, где битовые тесты были уместны, он делал это. Где N была константой, она не создавала очевидного "деления" для любого случая, она всегда использовала некоторый "трюк".
Просто позвольте компилятору сгенерировать код для вас, он почти наверняка будет знать больше об архитектуре машины, чем вы:) И это простая оптимизация, когда вы вряд ли придумаете что-то лучше, чем это делает компилятор.
Это интересный вопрос. Я не могу перечислить приемы, используемые компилятором для каждой константы, так как мне приходится компилировать на другом компьютере... Но я обновлю этот ответ позже, если меня никто не побьет:)
Это, вероятно, не поможет вам в коде, но в некоторых случаях есть полезная хитрость, которая может помочь вам сделать это:
Для деления на 3: Для числа, представленного в десятичном виде, вы можете сложить все цифры и проверить, делится ли сумма на 3.
Пример: 12345 => 1+2+3+4+5 = 15 => 1+5 = 6
, который делится на 3 (3 x 4115 = 12345)
,
Более интересно, что та же самая техника работает для всех факторов X-1, где X - основа, в которой представлено число. Таким образом, для десятичного числа вы можете проверить деление на 3 или 9. Для шестнадцатеричного вы можете проверить деление на 3,5 или 15. А для восьмеричных чисел вы можете проверить деление на 7.
В предыдущем вопросе я показал быстрый алгоритм проверки базы N для делителей, которые являются факторами N-1. Базовые преобразования между различными степенями 2 тривиальны; это просто группировка.
Таким образом, проверка на 3 легко в базе 4; проверка на 5 легка в базе 16, и проверка на 7 (и 9) легка в базе 64.
Непростые делители тривиальны, поэтому только 11 и 13 - сложные случаи. Для 11 вы можете использовать основание 1024, но в этот момент оно не очень эффективно для маленьких целых чисел.
Метод, который может помочь по модулю уменьшения всех целочисленных значений, использует нарезку битов и popcount.
mod3 = pop(x & 0x55555555) + pop(x & 0xaaaaaaaa) << 1; // <- one term is shared!
mod5 = pop(x & 0x99999999) + pop(x & 0xaaaaaaaa) << 1 + pop(x & 0x44444444) << 2;
mod7 = pop(x & 0x49249249) + pop(x & 0x92492492) << 1 + pop(x & 0x24924924) << 2;
modB = pop(x & 0x5d1745d1) + pop(x & 0xba2e8ba2) << 1 +
pop(x & 0x294a5294) << 2 + pop(x & 0x0681a068) << 3;
modD = pop(x & 0x91b91b91) + pop(x & 0xb2cb2cb2) << 1 +
pop(x & 0x64a64a64) << 2 + pop(x & 0xc85c85c8) << 3;
Максимальные значения для этих переменных - 48, 80, 73, 168 и 203, которые вписываются в 8-битные переменные. Второй раунд может проводиться параллельно (или может быть применен некоторый метод LUT)
mod3 mod3 mod5 mod5 mod5 mod7 mod7 mod7 modB modB modB modB modD modD modD modD
mask 0x55 0xaa 0x99 0xaa 0x44 0x49 0x92 0x24 0xd1 0xa2 0x94 0x68 0x91 0xb2 0x64 0xc8
shift *1 *2 *1 *2 *4 *1 *2 *4 *1 *2 *4 *8 *1 *2 *4 *8
sum <-------> <------------> <-----------> <-----------------> <----------------->
Стоит учесть одну вещь: поскольку вы заботитесь только о делимости до 16, вам действительно нужно проверять делимость только на простые числа до 16. Это 2, 3, 5, 7, 11 и 13.
Разделите ваш номер на каждое из простых чисел, следя за логическим значением (например, div2 = true). Числа два и три являются особыми случаями. Если div3 имеет значение true, попробуйте снова разделить на 3, установив div9. Два и его возможности очень просты (обратите внимание: "&" - это одна из самых быстрых возможностей процессора):
if n & 1 == 0:
div2 = true
if n & 3 == 0:
div4 = true
if n & 7 == 0:
div8 = true
if n & 15 == 0:
div16 = true
Теперь у вас есть логические значения div2, div3, div4, div5, div7, div8, div9, div11, div13 и div16. Все остальные числа являются комбинациями; например, div6 - это то же самое, что и (div2 && div3)
Таким образом, вам нужно только сделать 5 или 6 фактических делений (6, только если ваше число делится на 3).
Для себя я бы, вероятно, использовал биты в одном регистре для моих логических значений; например, bit_0 означает div2. Затем я могу использовать маски:
if (flags & (div2+div3)) == (div2 + div3): do_6()
обратите внимание, что div2+div3 может быть предварительно вычисленной константой. Если div2- это bit0, а div3 - это bit1, то div2+div3 == 3. Это позволяет оптимизировать приведенное выше выражение "if" так:
if (flags & 3) == 3: do_6()
Так что теперь... мод без дележа
def mod(n,m):
i = 0
while m < n:
m <<= 1
i += 1
while i > 0:
m >>= 1
if m <= n: n -= m
i -= 1
return n
div3 = mod(n,3) == 0
...
Кстати: наихудший случай для вышеупомянутого кода - 31 раз через любой цикл для 32-битного числа
К вашему сведению: только что посмотрел на пост Msalter выше. Его техника может использоваться вместо мода (...) для некоторых простых чисел.
Вы можете заменить деление на константу, не являющуюся степенью двойки, умножением, по сути, умножением на обратную величину вашего делителя. Детали, чтобы получить точный результат этим методом, сложны.
Восторг Хакера обсуждает это подробно в главе 10 (к сожалению, недоступно онлайн).
Из частного вы можете получить модуль другим умножением и вычитанием.
Быстрые тесты на делимость сильно зависят от базы, в которой представлено число. В случае, когда base равен 2, я думаю, что вы можете выполнять "быстрые тесты" на делимость только на степени 2. Бинарное число делится на 2n, если последние n двоичных цифр этого числа равны 0. Для других тестов я не Я думаю, что вы можете найти что-нибудь быстрее, чем %
,
Немного злого, запутанного, немного перепутанного, вы можете получить делимость на 15.
Для 32-разрядного числа без знака:
def mod_15ish(unsigned int x) {
// returns a number between 0 and 21 that is either x % 15
// or 15 + (x % 15), and returns 0 only for x == 0
x = (x & 0xF0F0F0F) + ((x >> 4) & 0xF0F0F0F);
x = (x & 0xFF00FF) + ((x >> 8) & 0xFF00FF);
x = (x & 0xFFFF) + ((x >> 16) & 0xFFFF);
// *1
x = (x & 0xF) + ((x >> 4) & 0xF);
return x;
}
def Divisible_by_15(unsigned int x) {
return ((x == 0) || (mod_15ish(x) == 15));
}
Вы можете построить аналогичные процедуры делимости для 3
а также 5
основанный на mod_15ish
,
Если у вас есть 64-битные целые числа без знака, расширьте каждую константу выше *1
линии очевидным образом, и добавьте строку над *1
линия, чтобы сделать правильный сдвиг на 32 бита с маской 0xFFFFFFFF
, (Последние две строки могут остаться прежними) mod_15ish
затем соблюдает тот же основной договор, но возвращаемое значение теперь между 0
а также 31
, (так что утверждается, что x % 15
== mod_15ish(x) % 15
)
Вот несколько советов, которые я еще никому не предлагал:
Одна идея состоит в том, чтобы использовать switch
оператор или предварительно вычислить некоторый массив. Тогда любой порядочный оптимизатор может просто индексировать каждый случай напрямую. Например:
// tests for (2,3,4,5,6,7)
switch (n % 8)
{
case 0: break;
case 1: break;
case 2: do(2); break;
case 3: do(3); break;
case 4: do(2); do(4) break;
case 5: do(5); break;
case 6: do(2); do(3); do(4); break;
case 7: do(7); break;
}
Ваша заявка немного двусмысленна, но вам может потребоваться только проверить простые числа меньше, чем n=16. Это потому, что все числа являются факторами текущих или предыдущих простых чисел. Таким образом, для n = 16 вы можете обойтись только проверкой 2, 3, 5, 7, 11, 13
как-то. Просто мысль.