Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[3.141.41.187] |
|
Сообщ.
#1
,
|
|
|
Реально быстрый алгоритм вычисления НОД Основными популярными алгоритмами вычисления наибольшего общего делителя (НОД) являются алгоритм Евклида и бинарный алгоритм. Первый очень простой и компактный, второй – якобы быстрый, поскольку в нём отсутствуют "долгие" операции деления, а присутствуют лишь операции сдвига. Однако на практике бинарный алгоритм зачастую работает раза в 2(!) медленнее алгоритма Евклида, поскольку в нём присутствует довольно много операций ветвления. Описанные в различных источниках бинарные алгоритмы вычисления НОД (в т.ч. в Википедии) подразумевают единичные операции на каждом цикле итерации и кучу проверок. По факту же за одну итерацию можно совершить сразу группу операций (сдвигов). Для тех, кто знает, как работает этот алгоритм, поясню, что сдвиг (деление пополам) можно производить безо всяких лишних проверок до тех пор, пока каждое из чисел не станет нечётным. После этого уже можно сделать сравнение значений с единицей, а далее (при отсутствии единицы) – нахождение разницы полученных нечётных значений с делением пополам и повтор цикла. Проверку же на 0 достаточно сделать 1 раз в самом начале. А вот, собственно, и сам код (на fasm): ; (c) 2018 by Jin X ; Return GCD of val1 and val2 in EAX (binary algorithm, no recursion calls) ; val1 and val2 are treated as unsigned proc FastGCD uses ebx esi edi, val1, val2 mov edx,[val2] mov eax,[val1] test edx,edx jz .val1noshift ; if val2 = 0 return val1 xor ebx,ebx ; bl = left shift count before return (result must be multiplied by 2^bl) test eax,eax jz .val2 ; if val1 = 0 return val2 .next: bsf esi,eax ; esi = count of lower zero bits in val1 bsf edi,edx ; edi = count of lower zero bits in val2 mov ecx,esi shr eax,cl ; shift val1 right until it will be odd mov ecx,edi shr edx,cl ; shift val2 right until it will be odd cmp esi,edi cmova esi,edi ; esi = min(esi,edi) add ebx,esi ; add this number to ebx (bl) cmp eax,1 je .val1 ; if val1 = 1 return 1 (val1) cmp edx,1 je .val2 ; if val2 = 1 return 1 (val2) cmp eax,edx je .val1 ; if val1 = val2 return val1 mov ecx,eax cmovb eax,edx cmovb edx,ecx ; swap val1 and val2 if val1 < val2 sub eax,edx shr eax,1 ; val1 = (val1-val2)/2 jmp .next ; go to next iteration .val2: mov eax,edx ; eax = val2 .val1: mov ecx,ebx shl eax,cl ; multiple result (eax) by 2^bl .val1noshift: ret endp Приведу заодно и реализацию алгоритма Евклида: ; (c) 2018 by Jin X ; Return GCD of val1 and val2 in EAX (Euclid algorithm, no recursion calls) ; val1 and val2 are treated as unsigned proc EuclidGCD val1, val2 mov eax,[val1] mov ecx,[val2] @@: jecxz @F ; return val1 if val2 = 0 xor edx,edx div ecx ; eax = val1 / val2; edx = val1 % val2 (replace this instruction by 'idiv' and previous one by 'cdq' to treat val1 and val2 as signed) mov eax,ecx ; eax = val2 mov ecx,edx ; ecx = edx = val1 % val2 jmp @B @@: ret endp Скорость работы бинарного алгоритма выше алгоритма Евклида примерно на 35-40% (на моём компьютере с процессором Intel Core i5-2500K @ 3.7 ГГц). Для кого-то это может показаться несущественным, поэтому алгоритм Евклида вполне имеет право на жизнь . Кроме того, справедливости ради хочу сказать, что в алгоритме Евклида отсутствуют инструкции cmovcc, требующие процессора Pentium Pro+ (правда, сейчас трудно найти более слабый процессор), без которых скорость алгоритма значительно снизится. К тому же, бинарный алгоритм не может работать с отрицательными числами (даже если заменить shr на sar), хотя этот вопрос легко решается заменой первых инструкций (до метки .next) на следующие, приводящие оба числа к абсолютному (положительному) значению: mov eax,[val2] cdq xor eax,edx sub eax,edx ; eax = abs(eax) mov ecx,eax xor ebx,ebx ; bl = left shift count before return (result must be multiplied by 2^bl) mov eax,[val1] cdq xor eax,edx sub eax,edx ; eax = abs(eax) = |val1| mov edx,ecx ; edx = |val2| jz .val2 ; if val1 = 0 return |val2| test edx,edx jz .val1noshift ; if val2 = 0 return |val1| Прикрепляю к статье архив с исходниками и программами для теста корректности и скорости работы алгоритмов. Наслаждайтесь p.s. Прикрепляю также архивчик GCDX.zip с немного изменёнными заголовками функций для обоих алгоритмов (стандартные прологи/эпилоги вырезаны и оптимизированы). Updated через пару часов: добавил ещё одну версию (GCD386.zip), без инструкций cmovcc, т.е. работающую на любом 386+ процессоре. Этот код медленнее на 10-13%, но всё равно быстрее алгоритма Евклида. Updated 09.02.2018: исправлен небольшой косячок в тестовых программах GCD.zip: вместо cinvoke использовался invoke. asm7x Прикреплённый файлGCD.zip (14,27 Кбайт, скачиваний: 177) – тестовые программы только тут! Прикреплённый файлGCDX.zip (2,73 Кбайт, скачиваний: 92) Прикреплённый файлGCD386.zip (1,63 Кбайт, скачиваний: 92) |
Сообщ.
#2
,
|
|
|
Цитата Jin X @ поскольку лично я с трудом представляю, как на каком-либо языке высокого уровня можно организовать работу инструкций bsf и cmovcc без ассемблерных вставок Никаких проблем, современные компиляторы имеют интринсики для bsf и умеют заменять некоторые if на cmov. Вот мой код на C++, в котором студия оба оператора ?: заменяет на cmov: FORCEINLINE num64 GCD(num64 a, num64 b) { if (a == 0) return b; if (b == 0) return a; unsigned long shift; _BitScanForward64(&shift, a | b); unsigned long index_a; _BitScanForward64(&index_a, a); a >>= index_a; do { unsigned long index_b; _BitScanForward64(&index_b, b); b >>= index_b; const num64 a1 = a; const num64 b1 = b; a = (a1 > b1) ? b1 : a1; b = (a1 > b1) ? a1 : b1; b -= a; } while (b); return (a << shift); } Генерируемый asm код почти 1 в 1, как у тебя, и так же на 40-50% быстрее классической реализации с делениями. P.S. Хотя нет, код не совсем как у тебя. В моем варианте внутренний цикл вообще полностью линейный, ни одного ветвления. |
Сообщ.
#3
,
|
|
|
Jin X, ваша ссылка на википедию у меня несколько странно выглядит. Точно всё верно?
|
Сообщ.
#4
,
|
|
|
Цитата Pacific @ Ок, убрал это из статьи Никаких проблем, современные компиляторы имеют интринсики для bsf и умеют заменять некоторые if на cmov Цитата Pacific @ do {} while ?В моем варианте внутренний цикл вообще полностью линейный, ни одного ветвления. А вообще, чуть позже проанализирую твой код (вижу пару прикольных моментов) Цитата Славян @ Fixed, thanks! ваша ссылка на википедию у меня несколько странно выглядит Добавлено Добавил GCD386.zip, кстати |
Сообщ.
#5
,
|
|
|
Цитата Jin X @ do {} while ? Ну, одно ветвление всегда есть, т.к. это цикл. А вот внутри самого цикла ветвлений нет. |
Сообщ.
#6
,
|
|
|
Pacific, хороший анализ алгоритма! Наверное меньше уже и не получить.
Сначала даже не понял, где шаги 4-7. Цитата 1. НОД(0, n) = n; НОД(m, 0) = m; НОД(m, m) = m; 2. НОД(1, n) = 1; НОД(m, 1) = 1; 3. Если m, n чётные, то НОД(m, n) = 2*НОД(m/2, n/2); 4. Если m чётное, n нечётное, то НОД(m, n) = НОД(m/2, n); 5. Если n чётное, m нечётное, то НОД(m, n) = НОД(m, n/2); 6. Если m, n нечётные и n > m, то НОД(m, n) = НОД((n-m)/2, m); 7. Если m, n нечётные и n < m, то НОД(m, n) = НОД((m-n)/2, n); Остаётся ощущение, что нет обработки шага 2. |
Сообщ.
#7
,
|
|
|
Кто-нибудь может объяснить сей феномен?
no intrinsic 4218 with intrinsic 4083 classic: 4094 assembly: 8705 Компилятор Intel® C++ Intel® 64 Compiler for applications running on IA-32, Version 16.0.3.207 Build 20160415 Код на ревью #include <iostream> #include <vector> #include <chrono> #if defined _WIN64 typedef unsigned long long dataType; #define BitScanForward _BitScanForward64 #else typedef unsigned long dataType; #define BitScanForward _BitScanForward #endif inline dataType BSF(dataType val) { int i; for (i = 0; (val & 1) == 0; ++i) val >>=1; return i; } dataType fun(dataType m, dataType n) { if (m == 0) return n; if (n == 0) return m; unsigned long shift = BSF(m | n); unsigned long index_a=BSF(m); m >>= index_a; do { unsigned long index_b = BSF(n); n >>= index_b; const dataType a = m; const dataType b = n; m = (a > b) ? b : a; n = (a > b) ? a : b; n -= m; }while (n); return m << shift; } dataType opt(dataType m, dataType n) { if (m == 0) return n; if (n == 0) return m; unsigned long shift; BitScanForward(&shift, m | n); unsigned long index_a; BitScanForward(&index_a, m); m >>= index_a; do { unsigned long index_b; BitScanForward(&index_b, n); n >>= index_b; const dataType a = m; const dataType b = n; m = (a > b) ? b : a; n = (a > b) ? a : b; n -= m; } while (n); return (m << shift); } dataType cls(dataType m, dataType n) { while (n != 0) { const dataType a = m % n; m = n; n = a; } return m; } #ifndef _WIN64 __declspec(naked) dataType asmed(dataType, dataType) { __asm { FastGCD: ; val1, val2 (parameters) mov edx,[esp+4] mov eax,[esp+8] test eax,eax jz val2noshift ; if val1 = 0 return val2 test edx,edx jz val1noshift ; if val2 = 0 return val1 push ebx push esi push edi xor ebx,ebx ; bl = left shift count before return (result must be multiplied by 2^bl) next: bsf esi,eax ; esi = count of lower zero bits in val1 bsf edi,edx ; edi = count of lower zero bits in val2 mov ecx,esi shr eax,cl ; shift val1 right until it will be odd mov ecx,edi shr edx,cl ; shift val2 right until it will be odd cmp esi,edi cmova esi,edi ; esi = min(esi,edi) add ebx,esi ; add this number to ebx (bl) cmp eax,1 je val1 ; if val1 = 1 return 1 (val1) cmp edx,1 je val2 ; if val2 = 1 return 1 (val2) cmp eax,edx je val1 ; if val1 = val2 return val1 mov ecx,eax cmovb eax,edx cmovb edx,ecx ; swap val1 and val2 if val1 < val2 sub eax,edx shr eax,1 ; val1 = (val1-val2)/2 jmp next ; go to next iteration val2: mov eax,edx ; eax = val2 val1: mov ecx,ebx shl eax,cl ; multiple result (eax) by 2^bl pop edi pop esi pop ebx val1noshift: ret val2noshift: mov eax,edx ; eax = val2 ret } } #endif template <typename Ch, typename Tr> std::basic_ostream<Ch, Tr>& operator<<(std::basic_ostream<Ch, Tr>& ostream, std::chrono::system_clock::duration time) { return ostream << std::chrono::duration_cast<std::chrono::milliseconds>(time).count(); } void testSpeed(const dataType *m, const dataType *n, int count, dataType (*func)(dataType, dataType), const char *msg) { std::chrono::duration<std::chrono::system_clock::rep, std::chrono::system_clock::period> res; res = res.zero(); for (int i = 0; i < count; ++i) { auto start= std::chrono::system_clock::now(); func(m[i], n[i]); auto diff = std::chrono::system_clock::now() - start; res += diff; } std::cout << msg << res << std::endl; } int main() { const int COUNT = 100000000; std::vector<dataType> m(COUNT), n(COUNT); for (auto &x : m) x = rand(); for (auto &x : n) x = rand(); testSpeed(&m[0], &n[0], COUNT, fun, "no intrinsic\t"); testSpeed(&m[0], &n[0], COUNT, opt, "with intrinsic\t"); testSpeed(&m[0], &n[0], COUNT, cls, "classic:\t"); #ifndef _WIN64 testSpeed(&m[0], &n[0], COUNT,asmed,"assembly:\t"); #endif } |
Сообщ.
#8
,
|
|
|
Цитата Федосеев Павел @ Остаётся ощущение, что нет обработки шага 2. В явном виде писать в коде ее не обязательно, все равно на выходе будет 1. Хотя если во входных данных одно из чисел часто равно 1, то имеет смысл добавить явную проверку. Добавлено Цитата Qraizer @ Кто-нибудь может объяснить сей феномен? Надо смотреть сгенерированный asm-код. |
Сообщ.
#9
,
|
|
|
fun
_TEXT SEGMENT PARA PUBLIC FLAT 'CODE' ; COMDAT ?fun@@YAKKK@Z ; mark_begin; ALIGN 16 PUBLIC ?fun@@YAKKK@Z ; --- fun(dataType, dataType) ?fun@@YAKKK@Z PROC NEAR ; parameter 1: eax ; parameter 2: edx .B9.1: ; Preds .B9.0 L9:: ;22.1 mov eax, DWORD PTR [4+esp] ;22.1 mov edx, DWORD PTR [8+esp] ;22.1 PUBLIC ?fun.@@YAKKK@Z ?fun.@@YAKKK@Z:: push ebx ;22.1 push ebp ;22.1 mov ebx, edx ;22.1 test eax, eax ;23.12 je .B9.19 ; Prob 28% ;23.12 ; LOE eax ebx esi edi .B9.2: ; Preds .B9.1 test ebx, ebx ;24.12 jne .B9.4 ; Prob 42% ;24.12 ; LOE eax ebx esi edi .B9.3: ; Preds .B9.2 pop ebp ;24.22 pop ebx ;24.22 ret ;24.22 ; LOE .B9.4: ; Preds .B9.2 mov ebp, eax ;26.33 xor edx, edx ;26.25 or ebp, ebx ;26.33 test ebp, 1 ;26.25 jne .B9.8 ; Prob 10% ;26.25 ; LOE eax edx ebx ebp esi edi .B9.6: ; Preds .B9.4 .B9.6 shr ebp, 1 ;26.25 inc edx ;26.25 test ebp, 1 ;26.25 je .B9.6 ; Prob 82% ;26.25 ; LOE eax edx ebx ebp esi edi .B9.8: ; Preds .B9.6 .B9.4 mov ebp, eax ;27.25 xor ecx, ecx ;27.25 test al, 1 ;27.25 jne .B9.12 ; Prob 10% ;27.25 ; LOE eax edx ecx ebx ebp esi edi .B9.10: ; Preds .B9.8 .B9.10 shr ebp, 1 ;27.25 inc ecx ;27.25 test ebp, 1 ;27.25 je .B9.10 ; Prob 82% ;27.25 ; LOE eax edx ecx ebx ebp esi edi .B9.12: ; Preds .B9.10 .B9.8 shr eax, cl ;29.3 ; LOE eax edx ebx esi edi .B9.13: ; Preds .B9.17 .B9.12 mov ebp, ebx ;33.29 xor ecx, ecx ;33.29 test bl, 1 ;33.29 jne .B9.17 ; Prob 10% ;33.29 ; LOE eax edx ecx ebx ebp esi edi .B9.15: ; Preds .B9.13 .B9.15 shr ebp, 1 ;33.29 inc ecx ;33.29 test ebp, 1 ;33.29 je .B9.15 ; Prob 82% ;33.29 ; LOE eax edx ecx ebx ebp esi edi .B9.17: ; Preds .B9.15 .B9.13 shr ebx, cl ;34.5 mov ebp, eax ;36.24 cmp ebx, eax ;39.5 cmovb eax, ebx ;39.5 cmp ebp, ebx ;42.5 cmova ebx, ebp ;42.5 sub ebx, eax ;40.5 jne .B9.13 ; Prob 82% ;43.11 ; LOE eax edx ebx esi edi .B9.18: ; Preds .B9.17 mov ecx, edx ;45.15 shl eax, cl ;45.15 pop ebp ;45.15 pop ebx ;45.15 ret ;45.15 ; LOE .B9.19: ; Preds .B9.1 mov eax, ebx ;23.22 pop ebp ;23.22 pop ebx ;23.22 ret ;23.22 ALIGN 16 ; LOE ; mark_end; ?fun@@YAKKK@Z ENDP ;?fun@@YAKKK@Z ENDS opt _TEXT SEGMENT PARA PUBLIC FLAT 'CODE' ; COMDAT ?opt@@YAKKK@Z ; mark_begin; ALIGN 16 PUBLIC ?opt@@YAKKK@Z ; --- opt(dataType, dataType) ?opt@@YAKKK@Z PROC NEAR ; parameter 1: eax ; parameter 2: edx .B7.1: ; Preds .B7.0 L7:: ;49.1 mov eax, DWORD PTR [4+esp] ;49.1 mov edx, DWORD PTR [8+esp] ;49.1 PUBLIC ?opt.@@YAKKK@Z ?opt.@@YAKKK@Z:: push ebp ;49.1 test eax, eax ;50.12 je .B7.6 ; Prob 28% ;50.12 ; LOE eax edx ebx esi edi .B7.2: ; Preds .B7.1 test edx, edx ;51.12 je .B7.7 ; Prob 4% ;51.12 ; LOE eax edx ebx esi edi .B7.3: ; Preds .B7.2 bsf ecx, eax ;57.3 mov ebp, eax ;54.3 or ebp, edx ;54.3 bsf ebp, ebp ;54.3 shr eax, cl ;59.3 ; LOE eax edx ebx ebp esi edi .B7.4: ; Preds .B7.4 .B7.3 bsf ecx, edx ;63.7 shr edx, cl ;64.7 mov ecx, eax ;66.26 cmp edx, eax ;69.7 cmovb eax, edx ;69.7 cmp ecx, edx ;72.7 cmova edx, ecx ;72.7 sub edx, eax ;70.7 jne .B7.4 ; Prob 82% ;73.12 ; LOE eax edx ebx ebp esi edi .B7.5: ; Preds .B7.4 mov ecx, ebp ;75.16 shl eax, cl ;75.16 pop ebp ;75.16 ret ;75.16 ; LOE .B7.6: ; Preds .B7.1 mov eax, edx ;50.22 pop ebp ;50.22 ret ;50.22 ; LOE .B7.7: ; Preds .B7.2 ; Infreq pop ebp ;51.22 ret ;51.22 ALIGN 16 ; LOE ; mark_end; ?opt@@YAKKK@Z ENDP ;?opt@@YAKKK@Z ENDS _TEXT ENDS cls _TEXT SEGMENT PARA PUBLIC FLAT 'CODE' ; COMDAT ?cls@@YAKKK@Z ; mark_begin; ALIGN 16 PUBLIC ?cls@@YAKKK@Z ; --- cls(dataType, dataType) ?cls@@YAKKK@Z PROC NEAR ; parameter 1: eax ; parameter 2: edx .B6.1: ; Preds .B6.0 L6:: ;79.1 mov eax, DWORD PTR [4+esp] ;79.1 mov edx, DWORD PTR [8+esp] ;79.1 PUBLIC ?cls.@@YAKKK@Z ?cls.@@YAKKK@Z:: mov ecx, edx ;79.1 test ecx, ecx ;80.15 je .B6.5 ; Prob 10% ;80.15 ; LOE eax ecx ebx ebp esi edi .B6.3: ; Preds .B6.1 .B6.3 xor edx, edx ;82.28 div ecx ;82.28 mov eax, ecx ;84.5 mov ecx, edx ;85.5 test edx, edx ;80.15 jne .B6.3 ; Prob 82% ;80.15 ; LOE eax ecx ebx ebp esi edi .B6.5: ; Preds .B6.3 .B6.1 ret ;87.10 ALIGN 16 ; LOE ; mark_end; ?cls@@YAKKK@Z ENDP ;?cls@@YAKKK@Z ENDS _TEXT ENDS asmed, понятное дело, без изменений. |
Сообщ.
#10
,
|
|
|
Qraizer
Сдается мне, что Интеловский компилятор слишком умный и перехитрил код, замеряющий время. Но не смог сделать этого с asm-кодом. |
Сообщ.
#11
,
|
|
|
В натуре слишком умный. Только дело в другом оказалось. Я предусмотрительно вынес определения тестируемых функций в другую единицу трансляции, но по привычке включил многофайловую мультипроцедурную оптимизацию, так что как будто и не выносил. В результате этот крендель навстраивал вызовы testSpeed(), отпромоушил указатели на fun(), opt() и cls(), попытался их тоже встроить и обнаружил отсутствие сайдэффектов от их вызовов. В результате в циклах остались только беготня по векторам и замеры времени.
Вот такие результаты без этой хрени: no intrinsic 13717 with intrinsic 7964 classic: 12770 assembly: 8167 |
Сообщ.
#12
,
|
|
|
А если в последнем варианте (cls) заменить строки:
1. "test ecx, ecx ;80.15" на "test edx, edx"; 2. пару "mov ecx, edx", "test edx, edx" поменять местами - "test edx, edx", "mov ecx, edx", то это незаметно ускорит, или совершенное крохоборство? |