Оптимизация длинной арифметики на C++

  • Tutorial

С Новым годом! Опишу классический сюжет — оптимизацию длинной арифметики в C++ при помощи ассемблерных вставок. Однако, на Хабре его еще не было, поэтому после некоторых колебаний решил запостить сюда, вы уж простите, если сами когда-то писали то же самое и продвинулись дальше меня :-)


Пусть в классе BigInt, представляющем длинное целое, объявлены поля void *words — массов слов, int wc — количество слов в массиве (длина слова в байтах задается на этапе компиляции). Тогда реализация сложения на чистом C++ может быть такой:
...
typedef unsigned int uInt;
typedef unsigned long long int uLong;

#define MAX_uInt ((uInt)(-1))
#define MUL ((uLong)MAX_uInt) // MUL: Max_uInt in Unsigned Long

int BigInt::WS = sizeof(uInt);

BigInt & BigInt::operator+=(const BigInt &rhs)
{
	...	// увеличить приемник, если необходимо
	uLong carry = 0;
	for (   uInt *th = (uInt*)words, *oz = (uInt*)rhs.words, 
			*limit = oz + rhs.wc, *thl = th + wc;
			oz < limit || carry && th < thl;
			oz++, th++)
	{
			uLong t = carry + *th;
			t += oz < limit ? *oz : 0;
			if (t > MUL)
			{
				carry = 1;
				t &= MUL;
			}
			else carry = 0;
			*th = (uInt)t;
	}
	if (carry) {
		...	//увеличить приемник на одно слово
	}
	return *this;
}

Видно, что Си не очень подходит для решения нашей задачи: нельзя складывать сразу по 8 байт, потому что не остается места для единственного бита переполнения.

Будем считать, сколько тактов процессор тратит на сложение каждых 4 байтов из массива слов длинного целого. Примерно так:
struct timespec t1, t2;
BigInt *a = ..., *b = ...;	// длина в байтах - 400
clock_gettime(CLOCK_MONOTONIC_RAW, &t1);
for (int j = 0; j < 10000000; j++) {
	*a += *b;
	*a -= *b;
}	
clock_gettime(CLOCK_MONOTONIC_RAW, &t2);
printf(...);	// пошаманить с t1 и t2

Вот прогресс по ключам оптимизации GCC:
-O1 -O2 -O3
7 6 6


Перепишем тело метода сложения с использованием ассемблерной вставки:
// uInt это int ИЛИ long long int
BigInt & BigInt::operator+=(const BigInt &rhs)
{
	...	// увеличить приемник, если необходимо
	uInt *th = (uInt*)words, *oz = (uInt*)rhs.words;
	
	asm goto (
		"clc\n"
		"l1:\t"
		"mov\t(%[oz]), %[ax]\n\t"
		"adc\t%[ax], (%[th])\n\t"
		"lahf\n\t"
		"add\t%[ws], %[th]\n\t"
		"add\t%[ws], %[oz]\n\t"
		"sahf\n\t"
		"loop\tl1\n\t"
		"jnc\t%l[nocarry]\n"
		:
		: [th] "r" (th), [oz] "r" (oz),
		  [ws] "I" (sizeof(uInt)), [ax] "a" ((uInt)0), "c" (wc)
		:
		: nocarry
	);

	...	//увеличить приемник на одно слово

	nocarry:	return *this;
}


В дальнейшем сразу буду писать результат компиляции (для uInt ≡ long long int):
...
	clc
lo: mov		(%rbx), %rax
	adc		%rax, (%rdx)
	lahf
	add		$8, %rdx
	add		$8, %rbx
	sahf
	loop	lo
	jnc		nocarry
...

Вкратце о том, что тут происходит: в регистрах bx и dx лежат указатели на текущее слово, adc — очень удобная инструкция, которая складывает источник, приемник и флаг переполнения, то что нужно! lahf и sahf сохраняют и загружают базовые флаги — эти инструкции нужны, потому что add сбрасывает флаг переполнения. loop итерирует столько раз, сколько лежит в регистре cx.

Это именно то (честно), что я сначала написал. Данный вариант выполняется за 7 тактов в обоих вариантах, соответственно за 3,5 такта на каждые четыре байта в «длинном варианте». GCC пока на высоте, но получившийся ассемблер еще оптимизировать и оптимизировать.

.

Дурные слухи ходят вокруг инструкции loop. Надо попробовать без нее, тем более, что это почти ничего не стоит:
...
	clc
lo: mov		(%rbx), %rax
	adc		%rax, (%rdx)
	lahf
	add		$8, %rdx
	add		$8, %rbx
	sahf
	dec		%rcx
	jnz		lo
	jnc		nocarry
...


Результат — 5 (2,5) тактов! Никогда не используйте loop — он экономит строчку, зато тратит целых 2 жирных такта на каждой итерации.

.

Хотелось бы разобраться с увеличением указателей — нет ли способа не задевать CF (carry flag — флаг «переноса», переполнения)? К счастью, есть — инструкция lea вычисляет ссылку на память и кладет ее в регистр-приемник, не трогая флаги. Вот как можно использовать ее для инкремента указателей:
...
	clc
lo: mov		(%rbx), %rax
	adc		%rax, (%rdx)
	lea		8(%rdx), %rdx
	lea		8(%rbx), %rbx
	dec		%rcx
	jnz		lo
	jnc		nocarry
...


По сути
lea 8(%rdx), %rdx
делает то же самое, что и
add $8, %rdx
— прибавляет к регистру размер слова.

На этот вариант процессор тратит 2 такта (то есть 1 такт на 4 байта). Признаться, сам не ожидал, что lahf/sahf занимают так много времени.

.

Что же делать дальше? SSE, к сожалению, тут не помощники: пока используется 64-битная архитектура, несуществующая инструкция padc xmm, xmm генерировала бы по сути тот же поток микроинструкций, что и последовательность из двух обычных сложений, распараллелить сложение с переполнением нельзя. Точно, значит надо и написать два сложения подряд, развернуть ассемблерный цикл. Золотой прием оптимизации.
// WS = 16
...
	clc
lo: mov		(%rbx), %rax
	adc		%rax, (%rdx)
	mov		8(%rbx), %rax
	adc		%rax, 8(%rdx)
	lea		16(%rbx), %rbx
	lea		16(%rdx), %rdx
	dec		%rcx
	jnz		lo
	jnc		nocarry
...


Теперь одна итерация выполняется за 3 такта, что соответствует 0.75 такта на 4 байта.

Ура! GCC -O2 сделан в 8 раз! Дальше думать некогда! Еще раз всех с наступающим!
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама
Комментарии 22
  • +8
    Самый новогодний программерский пост! :)
    • +1
      А вот эти «такты» — это окргугление попугаев по среднему? Просто слабо представляется работа clock_gettime на таких промежутках времени в многозадачной ОС (ведь многозадачной?)

      А так спасибо за пост.
      • +1
        Такты = наносекунды *3.2
        У меня много ядер. Вообще, есть предпосылки считать, что это именно такты. Большинство значений получились круглыми + 0,04 (остаток от временной доли обвязки после деления). Я эти 0,04 отовсюду вычел.
        • 0
          Во-первых спасибо за хороший пост, во-вторых — оптимизация за счет инлайнов — это уже какой-то прошлый век, ИМХО.

          Вот про оптимизацию циклов чуть по методу Loop-invariant Code Motion + про инлайны в комментариях:

          blogerator.ru/page/zametki-ob-optimizacii-programm-1-loop-invariant-code-motion-licm-singleton-i-optimizacija-funkcij

          • 0
            Не очень понял, к чему тут про LICM, в сишном цикле из поста как раз все что можно оптимизировано.

            По секрету, в посте действительно неправильно почти все что можно, зря спешил. Напишу «опровержение».
      • 0
        Ухх, было дело, помню-помню, такой же велосипед писал.
        • +2
          Интересно, что написано внутри GMP. К сожалению, не хватило времени туда залезть.
          • +2
            Да, мне тоже всегда было интересно. Версии GMP всегда меня по производительности уделывали, но оно и понятно — в документации сказано, что библиотека при компиляции умеет определять множество разных процессоров и соответствующим образом выбирать исходник.
            Кстати, система Mathematica использует какую-то свою версию GMP, потому, что она, в свою очередь, уделывает собранный мной GMP. Что, кстати, удивительно, ведь я-то собираю на целевой машине, а WR вынуждены собирать универсально.
            • 0
              Я его ее не померил. Будет чем первого заняться :)
              • 0
                * ее даже
                • +1
                  Буду ждать более подробную статью со сравнением с GMP.
                  Кстати, неплохо было бы еще разобраться с умножением и делением. Вот тут настоящий челендж вам, а сложение так, фигня :)
                  • +1
                    По поводу умножения — у меня помнится на x86 катастрофически нехватало регистра. Там же, если делать столбиком, необходимо хранить 2 счетчика цикла. Ух, как я мучался, помню… Хотя это наверное потому, что ассемблер я довольно ознакомительно знаю.
                    • +1
                      Умножение делается очень быстро свертками Фурье. Перемножаешь свертки множителей, получаешь свертку произведения; обратная БПФ — и вот тебе результат.

                      Я этот велосипед на ассемблере писал году эдак в 2002, вот мне делать было нехер тогда. :-)
            • +1
              Кстати, по поводу gcc версии. Вы, наверное, зря циклы переделывали из индексной формы в хитрые операции с указателями — компилятору сложнее оптимизировать код в этом случае. Возможно, если бы вы переписали код более высокоуровнево, получилось бы даже быстрее.
              • 0
                А там проще и не напишешь, оба лимита используются и в теле, и в условии цикла. Так писалось само, индексно — может быть уже целью.
              • 0
                а можно на код полностью взглянуть, пожалуйста?
              • 0
                Кстати, это походу последний пост в 2011 году по местному времени.
                С наступившим всех.
                • +2
                  Респект за умелое использование constraints в inline assembly.
                  • 0
                    меня почему-то тянет использовать esi/edi для хранения указателей на исходное число и память, куда записывается сумма.
                    • 0
                      bx и dx для логичности написал, в реальности от компиляции к компиляции GCC выбирает из dx, di и r8.

                      Меня тоже тянуло, попробовал один раз навести порядок:
                      : [th] "D" (th), [oz] "S" (oz),
                      Итерация замедлилась на 0,15 такта!

                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.