При программировании на старых процессорах, на которых операции умножения и деления чисел выполнялись медленно, программисты прибегали к трюкам, позволявшим ускорить вычисления. Так, битовый трюк, позволяющий получить остаток от деления на число, равное точной степени двойки, остаётся актуальным и сейчас. Операция типа a&((1<<s)−1) всё ещё работает быстрее обычного деления (в том случае, когда компилятор не имеет возможности выполнить соответствующую оптимизацию). Но с тех времён забытым остался трюк, позволяющий похожим набором операций заменить вычисление остатка от деления на число, на единицу меньшее степени двойки. Рассмотрим, как он работает.

Давайте договоримся, что s>1 и целое. Требуется вычислить а mod (2s−1).

Теория

Любое неотрицательное целое число a можно разложить на сумму

Представление числа в виде суммы

Поэтому а mod (2s−1) можно представить в виде

Рекурсивное вычисление модуля

Последнее равенство возможно в силу того, что 2s ≡ 1 mod (2s−1). Полученная формула даёт рекурсивный алгоритм вычисление остатка: чтобы вычислить а mod (2s−1), нужно вычислить остаток для числа a, сдвинутого вправо на s бит, а затем прибавить младшие s бит к результату, и от суммы затем снова взять остаток от деления на 2s−1. Эту операцию необходимо делать до тех пор, пока получившаяся сумма не станет меньше либо равной 2s−1.

Можно представить алгоритм иначе: разбиваем битовое представление числа на части по s бит и суммируем их как числа, затем к полученной сумме применяем ту же операцию до тех пор, пока эта сумма превосходит делитель.

Например, пусть s=3. Посчитаем остаток от деления числа 100(10) на 7. Для этого изобразим число 100 в двоичном коде, разбив по битам на три части:

Число 100 в двоичной записи

Вычисляем: 100 mod 7 = ((12 mod 7) + 4) mod 7 = (1 + 4 + 4) mod 7 = 9 mod 7 =[т.к. 9=1|001(2)]= 1+1=2.

Здесь есть небольшая проблема: когда исходное число a равно 2s−1, то алгоритм не работает, и нужно этот случай разбирать отдельно. То есть если бы мы захотели посчитать 7 mod 7, то получили бы 7, поэтому такие результаты нужно обнулять вручную.

Практика

Приведенная теория позволяет состряпать несложную программку:
P = ( 1 << s ) - 1; for ( ; x > P; x = z )
  for ( z = 0; x; x >>= s )
    z += x & P;
if ( x == P )
  x = 0;

Число x — это исходное число, оно же будет результатом. Число z — вспомогательное, оно после подсчётов не потребуется.

Когда смотришь на этот код, то сразу понятно, что он бесполезен на практике. На современных процессорах операция деления уже не такая медленная, чтобы проигрывать такому вложенному циклу. По крайней мере, в тех задачах, что мне приходилось решать, данная процедура проигрывала обычной операции деления. А если ещё учесть, что современные компиляторы умеют заменять деление на 2s−1 на битовые операции, хотя и делают это не всегда хорошо, то лучше не сопротивляться им. Можете проверить сами, с таким кодом Вы проиграете.

Другой разговор начинается, когда число s известно заранее (а так обычно и бывает). Ещё лучше, когда оно является степенью двойки. В этом случае внутренний цикл нашей процедуры можно заменить на серию битовых сдвигов, операций И и сложений. Существует известный трюк, предназначенный для подсчёта количества единичных битов в двоичном числе. Я предполагаю, что этот способ подсчета Вам хорошо известен. Так вот, небольшая его модификация позволят считать не биты, а сумму s-битовых чисел, входящих в исходное число a. Так, если s=8, то внутренний цикл уходит, остаётся только внешний:
for ( z = x; x > P; x = z ) {
  z = ( z & 0x00FF00FFu ) + ( ( z >> 8 ) & 0x00FF00FFu );
  z = ( z & 0x0000FFFFu ) + ( z >> 16 );
}
if ( x == P )
  x = 0;

Можно пойти дальше и доказать, что количество итераций этого цикла будет заведомо меньше какой-нибудь границы (зависит решаемой задачи). Это позволит развернуть цикл нужное число раз.

В конкурсе по вычислению обратной по модулю матрицы число s=31, однако заранее известно, что исходные числа не превышают 262, что делает возможным также удалить внутренний цикл, заменив его на операции с младшим тридцатью одним битом и старшим тридцатью одним битом. При этом внешний цикл совершает не более двух итераций, вторая из которых оперирует уже 32-х битовыми числами, а не 64-х битовыми.

Может быть кто-то соберётся делать масштабные эксперименты, чтобы проверить эффективность метода на практике? Тогда после компилирования программы убедитесь, что компилятор оставил инструкцию div на том месте, где ей положено быть, а также не применил к Вашему коду инструкции SSE, иначе сравнение будет не совсем корректным. Как я уже сказал выше, современные компиляторы уже встраивают подобный битовый трюк для деления на степень двойки без единицы, если заранее знают, что делитель имеет именно такой вид и «решили», что в данном случае это будет быстрее других методов вычисления остатка (например, через деление чисел типа float). То есть при тестировании убедитесь, что Вы тестируете именно то, что хотите.

Я таких экспериментов не проводил, но могу сказать, что трюк полезен, так как реально ускоряет программы. Например, в вышеупомянутом конкурсе версия алгоритма Гаусса «в лоб» работала около 2400 с, а если заменить операцию взятия остатка на функцию
int Mod ( int64 x ) {
  int32u z;
  z = int32u ( x & P ) + int32u ( x >> 31 );
  z = ( z & P ) + ( z >> 31 );
  if ( z == P )
    return 0;
  return int ( z );
}
то время работы сразу уменьшается до 1200 с.

Конечно, максимальный эффект достигается переписыванием на ассемблере, но это уже другая история. В этой статье предполагалось, что нет доступа к низкоуровневому программированию.

Обсуждение статьи на форуме в этой теме.