Сложение очень-очень длинных чисел

Складываем два длинных-предлинных числа X и Y (Y<X).

Вычисляем X=X+Y.

Основной цикл сложения выглядит так:



cikl: mov eax,[esi]

lea esi,[esi+4]

adc [edi],eax

lea edi,[edi+4]

dec ecx

jnz cikl



Как бы сделать это еще быстрее?

 

А в какой системе счисления записаны числа? В 2³²-ричной? Изврат...

afaik add работает быстрее, чем lea - команда попроще, микроинструкций меньше



Исправлено

Начет системы - извините, лох. Двоичная. Только работать с такими "очень длинными" числами все равно изврат imho...

 

Arvensis, числа хранятся двойными словами, но можно формат хранения и поменять, если это позволит складывать быстрее. Как предлагаешь?

Может add и быстрее, но он влияет на флаг переполнения,

поэтому lea здесь придется заменять на pushf / add / popf,

что уже будет медленнее...

 

add повлияет на adc , и помойму 1 такт будет что lea , что add

 

А числа ооочень большие? Может юзать SSE/SSE2?

 

Gray >

На сколько эти числа длинные и что с ними потом делать нужно?



По поводу цикла: IMHO тут особо не разгонишся, т.к. зависимости по результату из-за переноса..



Интереснго, если это складывать как-нибудь двумя потоками через DWORD, то как потом переносы учитывать?

 

Хотел предложить использовать lodsd, stosd, loop для уменьшения размера... хорошо, что не поленился проверить, прежде чем предлагать



Оказалось, что в случае с использованием lodsd и stosd скорость падает на 30-40%, а использование loop на скорость не влияет.



Правда, есть один нюанс. При размере чисел в 8 мегабайт все вычисление (Seleron 2,6) занимает порядка 100-150 миллисекунд. А при размере в 16 мегабайт винда начинает свопить массивы, и в результате ни о каких измерениях скорости речи быть не может - результат колеблется от 2 до 15 секунд.



Для желающих принять участие в тестировании привожу код:

Код (Text):

1. #include <windows.h>
2.  
3. #define Size 0x800000
4.  
5. DWORD x1[Size];
6. DWORD x2[Size];
7. DWORD y[Size];
8.  
9. void add1()
10. {
11.     __asm
12.     {
13.         mov edi,offset x1
14.         mov esi,offset y
15.         mov ecx,Size
16.         cld
17. cikl1:  lodsd
18.         adc eax,[edi]
19.         stosd
20.         loop cikl1
21.     }
22. }
23.  
24. void add2()
25. {
26.     __asm
27.     {
28.         mov edi,offset x2
29.         mov esi,offset y
30.         mov ecx,Size
31. cikl2:  mov eax,[esi]
32.         lea esi,[esi+4]
33.         adc [edi],eax
34.         lea edi,[edi+4]
35.         dec ecx
36.         jnz cikl2
37.     }
38. }
39.  
40. DWORD my_rand()
41. {
42.     return rand()%0x100 + ((rand()%0x100)<<8) + ((rand()%0x100)<<16) + ((rand()%0x100)<<24);
43. }
44.  
45. #pragma comment(linker, "/entry:main")
46. void main()
47. {
48. //  srand(0);
49.     srand(GetTickCount());
50.     for (DWORD i=0;i<Size;i++)
51.     {
52.         DWORD tmp=my_rand();
53.         x1[i]=tmp;
54.         x2[i]=tmp;
55.         y[i]=my_rand();
56.     }
57.     DWORD t0=GetTickCount();
58.     add1();
59.     DWORD t1=GetTickCount();
60.     add2();
61.     DWORD t2=GetTickCount();
62.    
63.     char buf[128];
64.     wsprintf(buf,"add1:%i, add2:%i",t1-t0,t2-t1);
65.     MessageBox(0,buf,"Finished",0);
66.     ExitProcess(0);
67. }

 

1. Использовано ADC и в сложении первых (младших) 32 бит числа. Что если на входе будет CF? Будет неверный результат. Если хочется экономить размер кода - сбрось CF

если важнее скорость - перед входом в цикл сделай сложение первых двойных слов (счётчик понятно в цикле должен быть уменьшен на 1 (т.е. (колл_байт_в_числах/4)-1)

2. Далее - можно использовать ecx одновременно как счётчик и индекс. Для этого нужно указать в счётчике вместо x -x

а в указатели загрузить предельный (первый невходящий) адрес.

т.е. если адресс a, а колличество двойных слов x

то адресс а можно представить как a+4x-4x.

т.е. если ecx это x, а скажем edi это a

то [edi]= [edi+ecx*4][4*-ecx]

далее вместо уменьшения наоброт умеличиваем счётчик.

т.е. пусть ecx - счётчик





lea edi,[edi][ecx*4]

lea esi,[esi][ecx*4]

neg ecx

clc



@loop:

mov eax,[esi][ecx*4]

adc [edi][ecx*4],eax

inc ecx

jne @loop

 

Их вообще нет смысла юзать - они работают с floating point, а не с целыми



А вот насчет MMX я пытался думать, и ничего хорошего из этого не вышло: если я все правильно понял, то PADDB/PADDW/PADDD/PADDQ суммируют массивы чисел, а перенос между элементами массивов не делают.

 

RobinFood

Их вообще нет смысла юзать - они работают с floating point, а не с целыми



PADDx могут работать и с XMM регистрами для SSE2



А вот насчет MMX я пытался думать, и ничего хорошего из этого не вышло

я тоже думал

как вариант, можно распаковывать длинные числа не по 32 бита, а по 31 бит.

тогда 31-й бит (считая от 0) можно рассматривать как заем/перенос при сложении/вычитании.

 

RobinFood >

Что-то не то с виндой



AXP @1666:



---------------------------

Size 0x800000

---------------------------

add1:109, add2:110

---------------------------





---------------------------

Size (2*0x800000)

---------------------------

add1:219, add2:203

---------------------------





Способ The Svin'а даёт прирост до 16%



---------------------------

Size (2*0x800000)

---------------------------

add1:188, add2:218

---------------------------





ЗЫ

Как сюда MMX/SSE прикрутить? Что с переносом-то делать

 

RobinFood

Че? SSE2 работает с целыми данными, а вот что делать с переносами действительно не понятно... Может хранить не 64 бита, а 63 и последий юзать как детект переноса? Хотя действительно вряд-ли выигрыш даст - скорее наоборот.

 

Когда мы задаем мегабайтные числа, существенную роль играет скорость обмена с памятью (в этом случае не мешало бы prefetch использовать). Если же все данные лежат в кэше, то результаты получаются, мягко говоря, не столь однозначные. Например, такой эксперимент - берем size = 4096, а для увеличения статистики используем 10^4 - 10^5 повторений запуска процедуры. У меня на P4 1.8GHz получилось, что из приведенных методов исходный add2 самый быстрый, метод TheSwin на 27% хуже, а lodsd\stosd - ни в какие ворота не лезет (почти в 2 раза хуже).

Возможно для S_T_A_S_ это очередное подтвержение отстойности и тупиковости P4 NetBurst ?



А известным методом повышения быстродействия циклов является увеличение числа операций на цикл. Если в исходном методе в цикле складывать по 4 dword-а, то получим реальный выигрыш (по крайней мере хуже никогда не будет). Например, в приведенном выше тесте на P4 получился выигрыш на 20%.

Код (Text):

1.     shr  ecx,2
2.     clc
3.   @@loop:
4.     mov eax,[esi]
5.     adc [edi],eax
6.     mov eax,[esi+4]
7.     adc [edi+4],eax
8.     mov eax,[esi+8]
9.     adc [edi+8],eax
10.     mov eax,[esi+12]
11.     adc [edi+12],eax
12.     lea esi,[esi+16]
13.     lea edi,[edi+16]
14.     dec ecx
15.     jnz @@loop

(Здесь, кстати, перестановки операций и задействование доп.регистров для уменьшения write-after-read delay ничего не дают - на P4)

 

TheSvin, очень изящное решение, сэр. Быстрее на 10-15%.



RobinFood, leo, S_T_A_S_ то что lods/stos вариант работает медленнее - проявление простой закономерности - чаще всего новые инструкции работают быстрее, чем старые,

и то, что было оптимальным для 8086 для Пентиумов уже таковым не является. Причина, IMHO, в том, что разработчики ленятся (и опасаются) переделывать то, что уже было сделано в кристалле предыдущих процессоров.



В реальной задаче цифры 1000-2000 бит длинной. Так, что кэш работает

При сложениии есть еще один трюк. При вычисленни b=b+a, если длинна а (обозначим ее A) меньше чем длинна b, то достаточно сделать А сложений, а затем если CF=1 cкладывать с 0. Код привожу ниже. (случай А>B пока не поддерживается)



void vadd(char* a,char* b) // b=b+a; a<b

{

_asm{

mov edi,b

mov esi,a

mov edx,edi // save pointer to b

movzx ecx, word ptr [esi]

movzx ebx, word ptr [edi]

sub ebx,ecx

lea edi,[edi+4*ecx+4]

lea esi,[esi+4*ecx+4]

neg ecx

clc

cikl3: mov eax,[esi+4*ecx] // основной цикл

adc [edi+4*ecx],eax

inc ecx

jne cikl3

jnc end_vadd

mov eax,ecx

cikl0: adc [edi+4*ecx],eax // eax=0

inc ecx

jc cikl0

end_vadd3: lea eax,[edi+4*ecx-4] // корректируем размер b

sub eax,edx

shr eax,2

mov [edx],ax

end_vadd:

}

}

 

leo >

IMHO это как раз не тот случай. Скорость чтения данных порядка 80Mb/sec. О кеше нужно думать, когда скорость напорядок больше.

И пока есть зависимость - учёт переноса, никуда от такой низкой скорости не деться.



>

А это вроде бы (полу)документировано intel'ом - из-за умножения ecx на 4. Если вынести уножение за цикл, то должно быть лучше.. только вот inc на add так просто не поменяешь.





Gray



Поясню, зачем я спрашивал о том, как эти числа используются.

Если 80% времени происходит сложение, а потом числа преобразуются, например, в ASCII, то можно попробовать распараллелить сложение (см. пост Max), а переносы учитывать именно при преобразовании - оно и так медленное.

 

leo

> „на P4 1.8GHz получилось, что из приведенных методов исходный add2 самый быстрый, метод TheSwin на 27% хуже“



Странно, а у меня метод TheSwin показал на 10-15% лучший результат. Проверял на трех машинах с P4.



S_T_A_S_

Поясню как эти числа используются:

n=<очень-очень большое число>;

y=z=0; d=1;

for(i=0;i<10000000000000000000000000;i++)

{

y+=z; z+=d; // вычисляем y=i^2 без умножений )

if(простой_и_быстрый_тест(n,y))){

{

сложный_и_долгий_тест(n,y) // чаще всего до этого теста дело не оходит

}

}

 

Если брать во внимание , что значительная часть тактов тратиться на jxx , то код leo должен давать ускорение т.к. обрабатывает несколько двордов за один проход цикла . У меня это подтвердилось , но может я неправильно тестирую , вот результаты (кол. тактов основного цикла на celeron 667mhz , где-то после 8-го прохода) :

Код (Text):

1. ============================================
2.           || Gray   ||  The Svin  ||  Leo
3. ============================================
4. 4096 byte || 9222   ||  8884      ||  5891
5. ============================================
6. 256  byte || 582    ||  564       ||  371
7. ============================================

1767312314__test.asm

 

Распараллеливание не противоречит тому что я написал.

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

Пока же это не сделаны - сравнения абсолютно некорректны и нечего тут в перегонки играть пока маршрут неясен.

Всё что я сделал

1. Поскольку автор не сказал чему кратно колличество байт в "ооочень больших числах" мне пришлось индуктивно предположить из его собственного кода что кратно 4м.

Разворачивание же цикла не может быть приведено в сравнение пока либо эта кратность неизвестна либо не внесены в код корректные пролог\эпилог которые вычисляют остаток от деления размера части числа ушедшего на интерацию в цикле и обрабатывают этот остаток.

2. Приведённый мною кусок кода всего лишь показывает как использовать счётчик в качестве индекса чтобы полности избежать команд на изменение источника и приемника.

Всё. Это идея а не процедура, один лишь винтик, не более.

3. Я отметил что CF должен быть сброшен в оригинальном тексте иначе процедура будет работать некорректно, либо первое слоежение должно быть с помощью add и лишь остальные - adc.



Теперь если рассматривать в тестировании мой кусочек как процедуру которая должна выполняться на P4.

1. Начало цикла итераций ОБЯЗАТЕЛЬНО должно быть выравнено по границе 16 байт. Т.е. нужно вставить ALIGN 16 перед меткой лупа.

2. P4 не любят inc \ dec (долго разъяснять - связано с архитектурой внутренних регистор) нужно здесь искать эквиваленты, пока мысли которые мне пришли в голову не кажутся заслуживающими внимания.

3. Команды на начале цикла зависимы - вот тут (если либо выравненность на 8 либо есть пролог\эпилог для остатка) очень хорошо можно распараллелить команды при разворачивании цикла.



Грубая прикидка (в предположении что длина чисел кратна 8и байтам)

Пусть в edi,esi адреса первого непренадлежащего числу байта лежащего сверху числа.

ecx - минус (число учетверённых слов в числе)

тогда цикл может выглядеть так:

clc

align 16

@@loop:

mov eax,[edi][ecx*8]

mov edx,[esi][ecx*8]

mov ebx,[esi][ecx*8][4]

adc eax,edx

mov edx,[edi][ecx*8][4]

mov [edi][ecx],eax

adc edx,ebx

inc ecx

mov [edi][ecx*8][4],edx

jne @@loop

 

The Swin >"Давайте поточнее определимся с задачей, тогда можно будет делать корректные сравнения..."



Прелесть задачи в том, что формат чисел можно выбрать тот, который обеспечит максимальную скорость вычислений. Хоть байты двух чисел поочереди ставить (а1,b1,a2,b2...). Ведь преобразовать исходный формат в тот который нужен достаточно один раз, а потом все делаем в "быстром" формате крутясь в циклах.



P.S. В реальной задаче цифры 1000-2000 бит длинной.

 

Gray

Кратность длины суперчисел (о которой спрашивает The Swin) тоже можно выбирать любую ?

Т.е. выбирать кратность дворду (1024,1056,1088,...,1984 бит) и выбирать кратность определённому количеству двордов . Или там могут встречаться такие длины -1008,1040,...,1992,2000 бит ?

В RSA например , ведь не даром используют конкретные 512,1024,2048 бит .



P.S. Я пока не оптимизатор , но хочу поучиться на реальном примере и опыте других