Полировка глюкал

[spamsink]

Рассмотрим целочисленное деление, точнее, вычисление остатка. Эта операция долгая, многотактная, поэтому при каждой удобной возможности - например, при вычислении остатка от деления на константу - ее нужно оптимизировать. Тривиальный случай остатка от деления на степень двойки исключаем; а деление 32-битного целого на произвольное наперед заданное число M делается просто: умножаем на "обратную величину" - число 2^N/M операцией, дающей в результате 64-битное число, и берем старшие разряды (N и количество старших разрядов определяются общим количеством бит M и количеством младших нулевых бит в нем, вручную это делать не нужно, компилятор уже сам умеет). Если нужно найти остаток, то умножаем результат на M и вычитаем из исходного числа. К примеру,

unsigned mod(unsigned x) {
    return x % 2037795097;
}

с помощью весьма нового GCC 4.7.2 превращается в

	movl	4(%esp), %ecx
	movl	$-2031890883, %edx
	movl	%ecx, %eax
	mull	%edx                    ; умножаем на "обратную величину"
	shrl	$30, %edx               ; берем старшие разряды
	imull	$2037795097, %edx, %edx ; умножаем частное на аргумент
	subl	%edx, %ecx              ; вычитаем
	movl	%ecx, %eax
	ret

Нетрудно проверить, что -2031890883 - это то же самое, что и 2³²-2031890883 = 2263076413, а 2⁶²/2037795097 = 2263076412.94.
Казалось бы, круто. Миллиард таких вычислений остатка делается около 4 с четвертью секунд на моем ноутбуке трехлетней давности под виртуальной машиной, чего еще желать?

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

unsigned mod(unsigned x) {
	unsigned a;
	a = x - 2037795097;
	if (a < x) x = a;
	a = x - 2037795097;
	if (a < x) x = a;
	return x;
}

что выходит из-под пресса как

	movl	4(%esp), %eax
	leal	-2037795097(%eax), %edx
	cmpl	%eax, %edx
	cmovbe	%edx, %eax
	leal	-2037795097(%eax), %edx
	cmpl	%eax, %edx
	cmovbe	%edx, %eax
	ret

Ни единого разрыва потока команд, прошу заметить. Три секунды ровно.

До этой оптимизации, которая годится для нахождения остатка от деления на любое число, помещающееся в разрядной сетке не более 3 раз, гнутые товарищи не дотумкали.

Comments

[fatoff.livejournal.com]

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

[spamsink]

Точнее, на однопоточном коде без системных вызовов и paging.

[lider.livejournal.com]

Вообще ничего вьїчитать не нужно

unsigned mod(unsigned x, unsigned int d) {
if(d>x)
return x;
return x % d;
}

[spamsink]

Откуда уверенность, что это будет быстрее? Условный переход нынче дорог.

[lider.livejournal.com]

Две секундьї?

[spamsink]

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

[lider.livejournal.com]

На том же где и 4 сек, 3 сек. ;)


unsigned mod(unsigned x) {
if(2037795097>х)
return x;
return x % 2037795097;
}


а потом вот такое:
#define mod(х) ((2037795097>х)?х:х%2037795097)

[spamsink]

См. коммент и добавление к нему ниже.

[spamsink]

Факты налицо (на более быстрой машине):

% foreach i (VARMOD CONSTMOD SUBTRACT)
foreach? gcc -O5 -D$i -o mod mod.c
foreach? time ./mod > /dev/null
foreach? end
4.275u ...
1.589u ...
0.789u ...

Печатают все три варианта, естественно, одно и то же: 326015744

Текст программы:

#include <stdio.h>

#if defined(CONSTMOD)
unsigned mod(unsigned x) {
    return x % 2037795097;
}
#elif defined(SUBTRACT)
unsigned mod(unsigned x) {
        unsigned a;
        a = x - 2037795097;
        if (a < x) x = a;
        a = x - 2037795097;
        if (a < x) x = a;
        return x;
}
#elif defined(VARMOD)
unsigned y = 2037795097;
unsigned mod(unsigned x) {
    return x % y;
}
#endif

main() {
        unsigned i, x = 0;
        static unsigned base = 100000000;
        for (i = base; i < base+1000000000; ++i) {
                x += mod(i);
        }
        printf("%d\n", x);
}

[spamsink]

Слегка меняем тест, для большей рандомизации и упрощения:

        unsigned i, x = 1;
        for (i = 0; i < 1000000000; ++i) {
                x += mod(x);
        }
        printf("%d\n", x);

Добавляем

#elif defined(CHECK)
unsigned mod(unsigned x) {
    return x < 2037795097 ? x : x % 2037795097;
}

Получаем, на foreach i ( VARMOD CONSTMOD SUBTRACT CHECK )

7.906u
3.908u
2.043u
3.612u

Архитектуру не проведёшь.

[izard.livejournal.com]

Если бы % встречалась в SpecINT и занимала сколько-нибудь времени, то обязательно бы додумали.

[spamsink]

Как в школе: чего не будет на экзамене, то и не изучают.

[archaicos.livejournal.com]

Но таки оцени прогресс в компиляторах и компьтерах за последние 25 лет. А пропуски всегда какие-то будут.

[spamsink]

Идея использовать умножение на обратное для деления на константу мне известна примерно те же 25 лет.

[dvv.livejournal.com]

И такая подмена была ещё в ричевском компиляторе для PDP-11. Если склероз мне не изменяет. А если изменяет — то в джонсоновском точно была.

[spamsink]

Нет, в джонсоновском в BSD 2.9 нет. Проверяется за считаные минуты: образ рута BSD 2.9, плюс apout (собирать с -m32); компилируем
unsigned mod(a) unsigned a; { return a / 3; }
с помощью
apout bin/cc -S -O2 mod.c
получаем

_mod:
~~mod:
jsr     r5,csv
~a=4
mov     4(r5),r1
clr     r0
div     $3,r0
jmp     cret

[dvv.livejournal.com]

Не, столько минут у меня нет. Особенно рич(и)евский компилер искать.

(Альцгеймер, правда, подсказывает, что это могло быть в run time библиотеке и/или эмуляторе отсутствующих команд…)

[spamsink]

Ричевский компилятор искать совсем просто :) http://pdp11.aiju.de/

[dvv.livejournal.com]

Дык. Всегда можешь им патч засабмитить.