Быстрое целочисленное деление на константу

На всех CPU операция деления выполняется сравнительно медленно, с этим ничего поделать нельзя. Но если делитель константа, то деление можно заменить на умножение на какую-то другую константу (обратное число, которое вычисляется во время компиляции). Тогда код будет чуть быстрее работать (и потреблять меньше энергии). Такую оптимизацию делают многие компиляторы (gcc, MSVC), но оказывается, многие разработчики не знают, как вычисляется сомножитель, а это не тривиально.

Дальше будет рассказано, как вычисляется сомножитель.

Идея и случай для вещественных чисел

Например, если в коде double x = a / 2.5, то это можно заменить на x = a * (1.0 / 2.5) или x = a * 0.4. Код будет быстрее, но иногда, менее точным. Для целых чисел такой трюк в лоб не получится, тк. сомножитель будет меньше единицы (те ноль). Можно использовать числа с фиксированной точкой (мантисса N): int x = a / b = a * B >> N, где B = (1 << N) / b. Проблема в том, что иногда так рассчитанный х будет на 1 отличатся от правильного ответа, что для целых чисел неприемлемо.

Алгоритм

Для упрощения будем работать только с неотрицательными числами. Математическая формулировка:

Для известного целого b найти такую пару целых M, B чтобы для всех целых a: 0 ≤ a < (1 << N ) выполнялось x = [a / b] = [a * B >> M], где [ ] – целая часть, а >> M это разделить a на 2 в степени M (и << умножить).

Запишем 1.0 / b в двоичном виде (иногда бесконечно длинное число), первые M бит которого обозначим как целое C, остальные биты (хвост-остаток) как D, те 1.0 / b = C >> M + D, D < (1 >> M).

Посмотрим, что будет, если B = C и a * D ≤ 1 (те N ≤ M):

x = [a / b] = [a * (C >> M + D)]= [a * C >> M + a*D] ≥ [a * C >> M] + [a*D] = [a * C >> M]

Получается, что если B это первые M бит числа 1/b, а остальные биты отбросить (заменить нулями), тогда будем получать иногда правильные значения, а иногда на 1 меньше чем x. Это происходит, потому что отбрасывается a*D, которое хоть и меньше единицы, но всё равно может изменять целую часть.

Попробуем увеличить B на единичку (те B = C + 1), тогда:

[a*(C + 1) >> M] = [a*C >> M + a * (1 >> M)] ≥ [a * C >> M + a*D] = [a * (C >> M + D)] = [a / b] = x

Теперь получается, что будем получать иногда правильный ответ x, а иногда на 1 больше.

Может показаться, что точно вычислять x через умножение невозможно тк при одном приближении будут иногда получатся числа на 1 меньше чем надо, при другом на 1 больше. Но это не так. Ведь a – это целые числа, а максимальная дробная часть a/b равна в точности (b-1) / b, и при вычислениях x оно увеличивается на a * ((1 >> M) – D). Ага, значит если a * ((1>>M) – D) < 1 / b тогда неравенство превращается в равенство!!! Пускай число L такое, что (1 << L) ≥ b, тогда для равенства достаточно (1 >> M) ≤ (1 >> (L + N)) или M ≥ L + N.

Пример на Java

В примере делитель равен 127 (L = 7). Результат умножения на JVM должен помещаться в 32-а бита, выбираем N = (32 – 7) / 2 = 12. Пример:

    final static int N = 12;

    private static int divideBy127(int a) {
        final int b = 127;
        final int L = 7;
        final int B = ((1 << (N + L)) / b) + 1;
        return (a * B) >>> (N + L);
    }

    public static void main(String[] args) {
        final int size = 20000;
        final int times = 10001;
        final int[] data = new int[size];
        for (int i = 0; i < size; ++i)  
            data[i] = (i * i) & ((1 << N) - 1);     // fill by any random data

        for (int iteration = 0; iteration < 10; ++iteration) {
            long time0 = System.currentTimeMillis();
            
            // Test 1: Using div
            int v0 = 0;
            for (int j = 0; j < times; ++j)
                for (int i = 0; i < size; ++i)
                    v0 ^= data[i] / 127;     // keep result
            long time1 = System.currentTimeMillis();
            
            // Test 2: Using mul
            int v1 = 0;
            for (int j = 0; j < times; ++j)
                for (int i = 0; i < size; ++i)
                    v1 ^= divideBy127(data[i]);     // keep result
            long time2 = System.currentTimeMillis();

            System.out.println((time1 - time0) + " vs " + (time2 - time1) + " " + (v0 - v1));
        }
    }

На моем Intel Core-i5 2500, JVM 1.7.0u5 с JIT и опцией -XX:+AggressiveOpts тест с делением работает в 2 раза медленнее, чем с умножением.

Дополнение

Остальное по теме не рассмотренное в этой статье:

• делимое может быть отрицательным;
• как избежать переполнения при умножении;
• результат умножения целых на x86 занимает 64-бита и при определенных условиях можно убрать сдвиг вправо;
• если известно, что числа делятся нацело, то можно использовать другой алгоритм;
• в некоторых случаях умножение можно заменить на сдвиги влево и сложения;
• частные случаи для некоторых делителей с ещё большей скоростью выполнения.

Литература

«Division by Invariant Integers using Multiplication», Torbjorn Granlund, Peter L. Montgomery — формальное доказательство и ещё некоторые алгоритмы по теме.