Ускорение деления

[B@R5uk]

Часто при написании программ приходится последовательно перебирать числа вида $m = a\cdot n$, где $n$ — это целочисленный индекс, изменяемый с единичным шагом, а $a$ — фиксированный множитель (как правило, тоже целочисленный). Разумеется, чтобы это считалось быстрее (при условии, что операция умножения значительно более затратная, чем сложение) такого рода вычисления лучше делать не в лоб:

Используется синтаксис Java

++n;
m = a * n;
 

а по-умному:

Используется синтаксис Java

++n;
m += a;
 

В связи с этим у меня возникла мысля: а нельзя ли точно так же ускорить вычисление результатов деления? Когда есть последовательные частные вида $a = q\cdot n + r$, где целочисленный делитель $n$ с каждой итерацией изменяется на единицу (или на другое какое-нибудь другое фиксированное число), $a$ — это фиксированное делимое, ну и частное $q$ с остатком $r$ подлежат вычислению на каждой итерации. То есть, можно ли вместо кода "в лоб":

Используется синтаксис Java

++n;
q = a / n;
r = a % n;
 

придумать какой-нибудь код, который используя знание о получившихся на прошлых итерациях частных и остатках, вычисляет величины на текущей итерации без необходимости применения собственно операции деления? Желательно, конечно, ещё бы обойтись и без умножения, но это не на столько критично (предполагается, что операция деления ещё более трудоёмка, чем умножение).

Проблема, естественно, в том, что в случае умножения мы двигаемся с каждой итерацией вдоль прямой, а при делении — вдоль гиперболы, у которой наклон меняется резко в начале координат, а в дальнейшем — всё медленнее и медленнее. Это надо будет как-то учитывать и корректировать.

[Dmitriy40]

Не знаю насчёт теории, на практике же часто применим известный способ замены вычислений выборкой из таблиц. Аппроксимируем нужную функцию (гиперболу) кусочно-линейно, храним лишь коэффициенты наклона (или приращения на каждом шаге). Если количество шагов мало и тем более заранее известно, то всё просто и даже без умножений. Если нужна большая точность (много шагов) или заранее неизвестно, то понадобится одно-два умножения (если сделать постоянный шаг невыгодно по соображениям размера или точности, то храним и наклон и величину шага в единицах минимального шага для каждого кусочка).

PS. Ну и ещё, если это не учебная задача, то стоит сначала выяснить действительно ли основные тормоза в операции деления. Просто потому что для современных процессоров деление конечно сильно медленнее умножения и тем более сложения/вычитания, но всего лишь раз в 20-30 (процессоры легко делают более сотни миллионов делений в секунду на ядро), и тормозить может например чтение памяти или тем более запись в память или слишком разветвлённый алгоритм (с плохим предсказанием переходов).

[Sender]

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

[photon]

Dmitriy40 в сообщении #1429957 писал(а):

Аппроксимируем нужную функцию (гиперболу) кусочно-линейно

Для столь простой функции не уверен, что дёрнуть нужные коэффициенты из памяти (а еще надо вычислить, откуда именно дёрнуть) и вычислить результат аппроксимации будет быстрее, чем посчитать в лоб. Конечно, для более сложных функций это работает и порой ускоряет в разы.

[mustitz]

Sender в сообщении #1429966 писал(а):

если речь об архитектурах x86, целочисленное деление с остатком реализуется одной машинной командой.

Вот только число тактов может быть различно. Какой-нить SYSENTER будет в сотни раз медленнее простой арифметики. У меня получилась разница приверно в три раза между сложением и умножением:

$ cat test1.c
#include <stdint.h>
#include <stdio.h>

#define N 10000000000ll

int main()
{
    const int64_t a = 7;
    int64_t dummy = 0;
    for (int64_t i=0; i<N; ++i) {
        const int64_t m = a * i;
        dummy ^= m;
    }
    printf("%ld\n", dummy);
    return 0;
}
$ gcc -Ofast test1.c -o test1
$ time ./test1
66780965888

real   0m4.639s
user   0m4.625s
sys   0m0.000s
$ time ./test1
66780965888

real   0m4.652s
user   0m4.635s
sys   0m0.000s
$ time ./test1
66780965888

real   0m4.393s
user   0m4.382s
sys   0m0.004s
$ cat test2.c
#include <stdint.h>
#include <stdio.h>

#define N 10000000000ll

int main()
{
    const int64_t a = 7;
    int64_t dummy = 0;
    int64_t m = 0;
    for (int64_t i=0; i<N; ++i) {
        dummy ^= m;
        m += a;
    }
    printf("%ld\n", dummy);
    return 0;
}
$ gcc -Ofast test2.c -o test2
$ time ./test2
66780965888

real   0m1.297s
user   0m1.249s
sys   0m0.036s
$ time ./test2
66780965888

real   0m1.300s
user   0m1.296s
sys   0m0.000s
$ time ./test2
66780965888

real   0m1.306s
user   0m1.302s
sys   0m0.000s

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

[wrest]

B@R5uk в сообщении #1429956 писал(а):

Разумеется, чтобы это считалось быстрее (при условии, что операция умножения значительно более затратная, чем сложение) такого рода вычисления лучше делать не в лоб:

А вы экспрементировали? Точно будет быстрее?
Ничего не делаем:
? a=3;s=0;for(i=1,10^7,s=a)
time = 15,781 ms.
Добавляем константу в аккумулятор:
?a=3;s=0;for(i=1,10^7,s=s+a)
time = 18,220 ms.
Умножаем константу на итератор:
? a=3;s=0;for(i=1,10^7,s=i*a)
time = 19,180 ms.
Делим (целочисленно) константу на итератор:
? a=3;s=0;for(i=1,10^7,s=a\i)
time = 18,380 ms.

Может, с большими числами будет по-другому?

? a=2^1333;s=0;for(i=1,10^7,s=a)
time = 17,000 ms.
? a=2^1333;s=0;for(i=1,10^7,s=s+a)
time = 21,390 ms.
? a=2^1333;s=0;for(i=1,10^7,s=i*a)
time = 22,190 ms.
? a=2^1333;s=0;for(i=1,10^7,s=a\i)
time = 28,461 ms.

Помойму затраты на организацию цикла тут заметно больше, чем на вычисления -- хоть с умножением хоть со сложением хоть с делением...

[mihaild]

mustitz, тут еще вопрос оптимизаций от компилятора. У меня ваш вариант со сложением на gcc отрабатывает за секунду, с умножением за 3.3. На clang - за 0.65 и 0.85 соответственно.
wrest, это же интерпретируемый язык? На нем производительность таких штук смотреть бессмысленно.

[mustitz]

mihaild в сообщении #1429973 писал(а):

mustitz, тут еще вопрос оптимизаций от компилятора. У меня ваш вариант со сложением на gcc отрабатывает за секунду, с умножением за 3.3. На clang - за 0.65 и 0.85 соответственно.

Да, тут надо смотреть ассемблер. Код мой нейсколько неудачен, потому как компилятор заменил 7*a = 8*a - a = (a<< 3) - a. Плюс в таком простом коде компилятор начинает интенсивно использовать SIMD-ы:

Используется синтаксис Text

$ gdb -q ./test1
Reading symbols from ./test1...done.
(gdb) b main
Breakpoint 1 at 0x570
(gdb) r
Starting program: /home/avs/sandbox/opt/test1 

Breakpoint 1, 0x0000555555554570 in main ()
(gdb) x/19 $pc
=> 0x555555554570 <main>:       sub    $0x8,%rsp
   0x555555554574 <main+4>:     pxor   %xmm2,%xmm2
   0x555555554578 <main+8>:     movdqa 0x230(%rip),%xmm0        # 0x5555555547b0
   0x555555554580 <main+16>:    xor    %eax,%eax
   0x555555554582 <main+18>:    movabs $0x12a05f200,%rdx
   0x55555555458c <main+28>:    movdqa 0x22c(%rip),%xmm3        # 0x5555555547c0
   0x555555554594 <main+36>:    nopl   0x0(%rax)
   0x555555554598 <main+40>:    movdqa %xmm0,%xmm1
   0x55555555459c <main+44>:    add    $0x1,%rax
   0x5555555545a0 <main+48>:    cmp    %rdx,%rax
   0x5555555545a3 <main+51>:    psllq  $0x3,%xmm1       # !!! Умножаем на 8 !!!
   0x5555555545a8 <main+56>:    psubq  %xmm0,%xmm1  # !!! Вычитаем !!!
   0x5555555545ac <main+60>:    paddq  %xmm3,%xmm0
   0x5555555545b0 <main+64>:    pxor   %xmm1,%xmm2
   0x5555555545b4 <main+68>:    jne    0x555555554598 <main+40>
   0x5555555545b6 <main+70>:    movdqa %xmm2,%xmm0
   0x5555555545ba <main+74>:    lea    0x1e3(%rip),%rsi        # 0x5555555547a4
   0x5555555545c1 <main+81>:    mov    $0x1,%edi
   0x5555555545c6 <main+86>:    xor    %eax,%eax
 

Уже в таком варианте честно используется pmuludq, и это уже 4.7 секунд.

Используется синтаксис C

#include <stdint.h>
#include <stdio.h>

#define N 10000000000ll

int main()
{
    int64_t a;
    sscanf("7\n", "%ld", &a);
    
    int64_t dummy = 0;
    for (int64_t i=0; i<N; ++i) {
        const int64_t m = a * i;
        dummy ^= m;
    }
    printf("%ld\n", dummy);
    return 0;
}
 

Ну а clang, я думаю, провёл еще более агрессивную оптимизацию с арифметикой.

[B@R5uk]

Dmitriy40 в сообщении #1429957 писал(а):

Не знаю насчёт теории, на практике же часто применим известный способ замены вычислений выборкой из таблиц.

В моём случае это исключено. Делимое может имеет порядок, например, $10^{30}$, а делитель пробегает значения от единицы до корня из этого значения. Кроме того, может проверяться множество различных делимых, поэтому ни у какого суперкопьютера на это дело памяти не хватит. Плюс здесь целочисленная арифметика, а не приблизительные вычисления, поэтому важна каждая единица в частном и остатке.

-- 13.12.2019, 14:19 --

mihaild в сообщении #1429973 писал(а):

wrest, это же интерпретируемый язык? На нем производительность таких штук смотреть бессмысленно.

Почему это бессмысленно? Я вот, например, пользуюсь Java. Поэтому для меня вопрос актуален не только в применении к какой-то единственной архитектуре, а вообще в широком смысле. Хотя, вопрос в первую очередь практический.

-- 13.12.2019, 14:24 --

wrest в сообщении #1429972 писал(а):

А вы экспрементировали? Точно будет быстрее?

Да, пробовал. Разницу во времени выполнения сложения и умножения засечь не получилось, а вот при добавлении в цикл деления быстродействие падает камнем вниз. К сожалению, я не знаю точной методологии замера времени работы подобных вещей, так как Java, во-первых, интерпретируемый язык, а во-вторых, компилятор в псевдокод выполняет определённые оптимизации. Но, опять же, меня интересует алгоритм, а не вопрос его применимости.

[wrest]

mihaild в сообщении #1429973 писал(а):

это же интерпретируемый язык? На нем производительность таких штук смотреть бессмысленно.

Смотреть надо там где будет применяться. Компиляция (в случае pari/gp) даёт какой-то выигрыш, но всё равно, пропорциональрый, т.е. пропорции можно оценить и так.
Конечно, какие-то специальные случаи (числа целиком помещаются в регистры, нет контролей переполнения и т.п.) дадут другие результаты, но практически, как пишет ТС, выполняется джавовский байт-код, так что все равно на организацию цикла затрачивается львиная часть времени, а если в цикле помимо вычисления m=a*n делаются какие-то ещё, то на них.

-- 13.12.2019, 14:58 --

B@R5uk в сообщении #1429982 писал(а):

Делимое может имеет порядок, например, $10^{30}$

Ну то есть в 64-битный регистр делимое не влезает.

[mihaild]

mustitz в сообщении #1429978 писал(а):

Уже в таком варианте честно используется pmuludq

Зависит от контекста, что честно, а что нет. Вполне может быть что и в реальном коде после разворачивания цикла почти всё можно будет сделать SIMD.

B@R5uk в сообщении #1429982 писал(а):

Почему это бессмысленно?

Потому что они делают что-то своё и очень странное.

B@R5uk в сообщении #1429982 писал(а):

Java, во-первых, интерпретируемый язык

Всё-таки нет. Там довольно честная JIT компиляция (проблемы из-за автоматического управления памятью конечно остаются, но вряд ли они тут влияют).

wrest в сообщении #1429985 писал(а):

т.е. пропорции можно оценить и так

Нельзя, потому что истинность утверждения

wrest в сообщении #1429985 писал(а):

на организацию цикла затрачивается львиная часть времени

зависит от языка (и компилятора). И элементарно разворачивание цикла (в тех языках, где оно возможно и компилятор его поддерживает) автоматически снижает накладные расходы на него.

[Dmitriy40]

wrest в сообщении #1429985 писал(а):

B@R5uk в сообщении #1429982 писал(а):

Делимое может имеет порядок, например, $10^{30}$

Ну то есть в 64-битный регистр делимое не влезает.

Главное что при x64 делитель влезает в регистр и деление выполняется одной-двумя командами, суммарно 100-200 тактов (одной командой можно обойтись при условии что делитель больше сдвинутого вправо на 64 бита делимого, что для делимого $10^{30}$ выполняется для делителей больше $5{,}4\cdot10^{10}$, т.е. для 99.995% всего диапазона). На x86 делитель в регистр не влезает и возможно быстрее будет классический цикл на 50 сдвигов и вычитаний, частично развёрнутый в условные пересылки командами cmov.

А ещё можно сделать финт ушами: использовать FDIV, которая выдаст примерный результат за 24 такта, который парой-тройкой умножений или полдесятком сложений/вычитаний (и то и то за пару десятков тактов) довести до точного. А в AVX есть команда VDIVPD, которая за 35 тактов выдаст сразу 4 итерации, которые тоже тактов за 30-40 наверняка можно досчитать до точных значений, сразу все 4, итого в среднем 15-20 тактов на итерацию.

Как всё это объяснить компилятору ЯВУ - не спрашивайте!

-- 13.12.2019, 18:28 --

Ну и чисто алгоритмически.

B@R5uk в сообщении #1429982 писал(а):

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

Можно не вычислять на каждой итерации деление, а вычислять скажем с шагом 16-1024, промежуточные же получать линейной интерполяцией, что потребует лишь трёх умножений (третье - для проверки и коррекции результата если ошиблись на единицу) и одного сдвига вправо. Ускорение в десятки раз, до скорости трёх умножений (два из которых кстати выполнятся параллельно) вместо одного деления. И работает для достаточно больших делителей, например при делимом $10^{30}$ и интерполяции 1:16 делитель должен быть больше $10^{11}$, что покрывает 99.99% диапазона. Для интерполяции 1:1024 делитель должен быть больше $10^{12}$, но и это 99.9% всего диапазона.
Примерно так (алгоритм, не исходный код):

Используется синтаксис C

w=10**30;//w=qn+r, w=10^30, n=10^11..10^15
for (n=10**11; n<10**15; n+=16) {
        a=w\n; b=w\(n+16);
        for (i=0; i<16; i++) {
                ni=n+i;
                q=(a*(16-i)+b*i)\16;//Линейная интерполяция
                r=w-q*ni;
                while (r<0) {q--; r+=ni}//Коррекция если результат завышен, выполняется 0-1 раз
                while (r>=ni) {q++; r-=ni}//Коррекция если результат занижен, выполняется 0-1 раз
                ...
        }
}

Hint: оба умножения в линейной интерполяции тоже можно заменить на сложение/вычитание, что ускорит вообще до одного умножения и сложений/вычитаний.

И да, пороги я бы не считал заранее, а сделал бы автовыбор алгоритма, или прямое деление если делитель маленький, или линейную интерполяцию если он уже достаточно большой. И переходил бы от первого ко второму проверяя точность интерполяции одного (или пары) из чисел внутри диапазона интерполяции. При 1:16 это прилично времени займёт, но можно же и сразу 1:256 или более сделать, где это время растворится. Да и в любом случае это будет не более чем 0.01% от всего времени.

[B@R5uk]

Dmitriy40, спасибо. Вы предложили действительно интересный подход. И при этом работающий. Но в его основе лежат аналитические соображения, а мне хотелось бы из алгебраических соображений построить (и ещё желательно доказать верность) такую рекуррентную последовательность, чтобы она давала правильные частные и остатки без деления вообще (и желательно без умножения). Я понимаю, что это уже не программирование, и мне важна сама задача, а не её применимость (хотя корни её лежат в программировании), так что эта тема, наверное, больше для раздела "Помогите решить/разобраться". Но тем не менее. Это и есть суть моего вопроса.

Для иллюстрации сказанного, привожу следующий код. Я его угадал методом научного тыка, и совершенно без понятия, почему он работает. Код рассчитывает последовательные частные от корня из делимого до половины делимого (верхняя граница найдена экспериментально). То есть, если делимое — $x$, то код верно работает на отрезке $\left[ \sqrt{x},{x}/{2}\right]$.

Используется синтаксис Java

import java.util.*;

public class StudyDivisionReplacement {
        
        private static final Random rng = new Random();
        
        public static void main(String[] args) {
                long dividend, divisor, startDivisor, quotient, remainder,
                                                        quotientRef, remainderRef, remainderShift;
                dividend = 100000000l + rng.nextInt(1000000);
                startDivisor = (long) Math.sqrt(dividend);
                divisor = startDivisor;
                quotient = dividend / divisor;
                remainder = dividend % divisor;
                if (divisor == quotient) {
                        remainderShift = 1;
                } else {
                        remainderShift = 0;
                }
                while (divisor > remainderShift) {
                        quotientRef = dividend / divisor;
                        remainderRef = dividend % divisor;
                        if (quotientRef != quotient || remainderRef != remainder) {
                                System.out.println(String.format("Error: %d / %d = %d (%d) not %d (%d)",
                                                dividend, divisor, quotientRef, remainderRef, quotient, remainder));
                                ++divisor;
                                break;
                        }
                        ++divisor;
                        remainder += remainderShift;
                        if (remainder >= divisor) {
                                remainder -= divisor;
                                ++remainderShift;
                        } else {
                                --quotient;
                                remainderShift += 2;
                        }
                }
                --divisor;
                System.out.println(String.format("Dividend: %d\nStart: %d\nFinish: %d",
                                                                                dividend, startDivisor, divisor));
        }
}
 

Как видно из примера, весь расчёт одной итерации производится четырьмя сложениями/вычитаниями и одной проверкой. К сожалению меня больше интересует возможность рекуррентного расчёта на отрезке $\left[2, \sqrt{x}\right]$, но приведённый пример даёт надежду на существование такого алгоритма.

[Dmitriy40]

B@R5uk в сообщении #1430011 писал(а):

и совершенно без понятия, почему он работает.

Особенно не вникал, но очень похоже он работает потому что для достаточно больших делителей частные меняются относительно медленно, гипербола аппроксимируется длинными горизонтальными участками, вот их Вы и проходите. Да, если бы сразу сказали что нужен такой диапазон (делитель сильно больше корня из частного), нечто похожее и предложил бы.
Но для малых делителей это работать не будет т.к. аппроксимирующие прямые будут сильно вертикальны, т.е. при изменении делителя на единицу частное будет меняться на много единиц, что и ломает этот алгоритм (собственно он ломается при производной меньше -1, т.е. точно до точки корня из делимого). А при x/2 он ломается (если ломается, не разбирался) лишь из-за огрехов реализации. Завернув внутренность в цикл while можно расширить область работы вниз от корня из делимого, раз в 2-3-5, но ценой соответствующего уменьшения скорости.
Линейная же интерполяция будет работать и в этом случае тоже (правда со всё возрастающей погрешностью с уменьшением делителя).
Насчёт рекуррентной формулы я сомневаюсь, но подождём вдруг кто из математиков заинтересуется.

[B@R5uk]

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

Так же, ключевой момент в моём алгоритме заключается в квадратичном росте остатка, а не в линейном убывании делимого. Обратите внимание: к остатку добавляется линейно растущая величина (в начале интервала, она равна 2, в конце — 1, а в середине — эти величины чередуются).

А ещё, я сейчас заметил, что частные лежат как раз на интересующем меня отрезке $\left[2, \sqrt{x}\right]$, только проходят его от конца к началу с повторениями. Как бы так исхитриться и вывернуть алгоритм наизнанку, чтобы делители лежали на этом отрезке?