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

[Jin X]

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

Основными популярными алгоритмами вычисления наибольшего общего делителя (НОД) являются алгоритм Евклида и бинарный алгоритм. Первый очень простой и компактный, второй – якобы быстрый, поскольку в нём отсутствуют "долгие" операции деления, а присутствуют лишь операции сдвига. Однако на практике бинарный алгоритм зачастую работает раза в 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

(замечу, что при внесении описанных в комментариях изменений в алгоритм Евклида для работы с отрицательными числами результат может быть отрицательным при отрицательном значении одного или обеих параметров... в прикреплённом архиве есть signed-функция, которая выполняет приведение результата к абсолютному значению).

Скорость работы бинарного алгоритма выше алгоритма Евклида примерно на 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|

При этом скорость работы алгоритма падает лишь на 1-2%

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

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)

[Pacific]

Цитата Jin X @ 03.02.18, 09:26

поскольку лично я с трудом представляю, как на каком-либо языке высокого уровня можно организовать работу инструкций 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. Хотя нет, код не совсем как у тебя. В моем варианте внутренний цикл вообще полностью линейный, ни одного ветвления.

Сообщение отредактировано: Pacific - 03.02.18, 10:32

[Славян]

Jin X, ваша ссылка на википедию у меня несколько странно выглядит. Точно всё верно?

[Jin X]

Цитата Pacific @ 03.02.18, 10:26

Никаких проблем, современные компиляторы имеют интринсики для bsf и умеют заменять некоторые if на cmov

Ок, убрал это из статьи

Цитата Pacific @ 03.02.18, 10:26

В моем варианте внутренний цикл вообще полностью линейный, ни одного ветвления.

do {} while ?
А вообще, чуть позже проанализирую твой код (вижу пару прикольных моментов)

Цитата Славян @ 03.02.18, 10:49

ваша ссылка на википедию у меня несколько странно выглядит

Fixed, thanks!

Добавлено 03.02.18, 11:09
Добавил GCD386.zip, кстати

[Pacific]

Цитата Jin X @ 03.02.18, 11:03

do {} while ?

Ну, одно ветвление всегда есть, т.к. это цикл. А вот внутри самого цикла ветвлений нет.

[Федосеев Павел]

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.

[Qraizer]

Кто-нибудь может объяснить сей феномен?

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
}

Для Оптимизирующий компилятор Microsoft ® C/C++ версии 19.00.24215.1 для x86 результаты в целом ожидаемы.

[Pacific]

Цитата Федосеев Павел @ 03.02.18, 12:17

Остаётся ощущение, что нет обработки шага 2.

В явном виде писать в коде ее не обязательно, все равно на выходе будет 1. Хотя если во входных данных одно из чисел часто равно 1, то имеет смысл добавить явную проверку.

Добавлено 03.02.18, 15:40

Цитата Qraizer @ 03.02.18, 13:09

Кто-нибудь может объяснить сей феномен?

Надо смотреть сгенерированный asm-код.

[Qraizer]

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, понятное дело, без изменений.

[Pacific]

Qraizer
Сдается мне, что Интеловский компилятор слишком умный и перехитрил код, замеряющий время. Но не смог сделать этого с asm-кодом.

[Qraizer]

В натуре слишком умный. Только дело в другом оказалось. Я предусмотрительно вынес определения тестируемых функций в другую единицу трансляции, но по привычке включил многофайловую мультипроцедурную оптимизацию, так что как будто и не выносил. В результате этот крендель навстраивал вызовы testSpeed(), отпромоушил указатели на fun(), opt() и cls(), попытался их тоже встроить и обнаружил отсутствие сайдэффектов от их вызовов. В результате в циклах остались только беготня по векторам и замеры времени.
Вот такие результаты без этой хрени:

no intrinsic    13717
with intrinsic  7964
classic:        12770
assembly:       8167

[Славян]

А если в последнем варианте (cls) заменить строки:
1. "test ecx, ecx ;80.15" на "test edx, edx";
2. пару "mov ecx, edx", "test edx, edx" поменять местами - "test edx, edx", "mov ecx, edx",
то это незаметно ускорит, или совершенное крохоборство?