Форум программистов, компьютерный форум, киберфорум
Наши страницы
Jin X
Войти
Регистрация
Восстановить пароль
Оценить эту запись

Реально быстрый алгоритм вычисления НОД

Запись от Jin X размещена 03.02.2018 в 02:24
Обновил(-а) Jin X 10.02.2018 в 20:59

Реально быстрый алгоритм вычисления НОД

Основными популярными алгоритмами вычисления наибольшего общего делителя (НОД) являются алгоритм Евклида и бинарный алгоритм. Первый очень простой и компактный, второй – якобы быстрый, поскольку в нём отсутствуют "долгие" операции деления, а присутствуют лишь операции сдвига. Однако на практике бинарный алгоритм зачастую работает раза в 2(!) медленнее алгоритма Евклида, поскольку в нём присутствует довольно много операций ветвления. Описанные в различных источниках бинарные алгоритмы вычисления НОД (в т.ч. в Википедии) подразумевают единичные операции на каждом цикле итерации и кучу проверок. По факту же за одну итерацию можно совершить сразу группу операций (сдвигов). Для тех, кто знает, как работает этот алгоритм, поясню, что сдвиг (деление пополам) можно производить безо всяких лишних проверок до тех пор, пока каждое из чисел не станет нечётным. После этого уже можно сделать сравнение значений с единицей, а далее (при отсутствии единицы) – нахождение разницы полученных нечётных значений с делением пополам и повтор цикла. Проверку же на 0 достаточно сделать 1 раз в самом начале.

А вот, собственно, и сам код (на fasm):
Assembler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
; (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
Приведу заодно и реализацию алгоритма Евклида:
Assembler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
; (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
(замечу, что при внесении описанных в комментариях изменений в алгоритм Евклида для работы с отрицательными числами результат может быть отрицательным при отрицательном значении одного или обеих параметров... в прикреплённом архиве есть signed-функция, которая выполняет приведение результата к абсолютному значению).

Скорость работы бинарного алгоритма выше алгоритма Евклида примерно на 35-40% (на моём компьютере с процессором Intel Core i5-2500K @ 3.7 ГГц). Для кого-то это может показаться несущественным, поэтому алгоритм Евклида вполне имеет право на жизнь . Кроме того, справедливости ради хочу сказать, что в алгоритме Евклида отсутствуют инструкции cmovcc, требующие процессора Pentium Pro+ (правда, сейчас трудно найти более слабый процессор), без которых скорость алгоритма значительно снизится. К тому же, бинарный алгоритм не может работать с отрицательными числами (даже если заменить shr на sar), хотя этот вопрос легко решается заменой первых инструкций (до метки .next) на следующие, приводящие оба числа к абсолютному (положительному) значению:
Assembler
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
        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|
При этом скорость работы алгоритма падает лишь на 1-2%

Прикрепляю к статье архив с исходниками и программами для теста корректности и скорости работы алгоритмов.
Наслаждайтесь

Эта и другие мои статьи на форуме.



Updated 03.02.2018:
  • Внёс небольшие изменения в код signed-функции, а также добавил архивчик GCDX.zip с немного изменёнными заголовками функций для обоих алгоритмов (стандартные прологи/эпилоги вырезаны и оптимизированы).
  • Через пару часов: добавил ещё одну версию (GCD386.zip), без инструкций cmovcc, т.е. работающую на любом 386+ процессоре. Этот код медленнее на 10-13%, но всё равно быстрее алгоритма Евклида.
Updated 09.02.2018:
  • Исправлен небольшой косяк в GCD.zip: вместо cinvoke использовался invoke.

p.s. Тестовые программы есть только в архиве GCD.zip !!!
Вложения
Тип файла: zip GCDX.zip (2.7 Кб, 81 просмотров)
Тип файла: zip GCD386.zip (1.6 Кб, 101 просмотров)
Тип файла: zip GCD.zip (14.3 Кб, 68 просмотров)
Просмотров 725 Комментарии 0
Всего комментариев 0
Комментарии
 
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru