Дальше будет рассказано, как вычисляется сомножитель.
Идея и случай для вещественных чисел
Например, если в коде 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 — формальное доказательство и ещё некоторые алгоритмы по теме.