Największy Wspólny Dzielnik - algorytm Euklidesa

Dawno, dawno temu ...
Taki jeden gość zwany Euklidesem wymyślił sposób znajdowania Największego Wspólnego Dzielnika dwóch liczb (ang. the Greatest Common Divisor).

1. Mamy dwie liczby u i v
2. Jeżeli liczby są równe to NWD(u, v) = u
3. Jeżeli u > v to u := u - v i idziemy do punktu 2.
4. Jeżeli u < v to v := v - u i idziemy do punktu 2.

Jeżeli to takie proste to piszemy funkcję która nam to obliczy, wykorzystując rekurencję:

function GcdOd(u, v: Cardinal): Cardinal;
begin
 if v > u then
   Result := GcdOd(u, v-u);
 if v < u then
   Result := GcdOd(u-v, v);
 if u = v then
   Result := v;
end;

GcdOd(1000, 2132424);
Średni czas wykonywania: 40457 cykli zegara

Fajnie, ale trochę długo oblicza (wywołajcie to kilkadziesiąt tysięcy razy)

Jak się temu dobrze przyjrzeć to odejmowanie możemy zastąpić resztą z dzielenia jednej liczby przez drugą i wyowływać na przemian funkcję. Jeżeli u < v to u-v = u mod v, jeżeli u > v to u mod v = u i nie ma zmiany.

function GcdR(u, v: Cardinal): Cardinal;
begin
 if v = 0 then
   Result := u
 else
   Result := GcdR(v, u mod v);
end;

GcdR(1000, 2132424);
Średni czas wykonywania: 399 cykli zegara

O ta funkcja działa znacznie szybciej (101 razy szybciej)

Ale czy nie można tego jeszcze trochę przyspieszyć?
Jak wiadomo każdą funkcję rekurencyjną (wywołującą samą siebie) można przerobić na funkcję iteracyjną (wykorzystującą pętle). Funkcje iteracyjne są zazwyczaj szybsze, ale w porównaniu do rekurencyjnych zazwyczaj są mniej czytelne i zużywają więcej pamięci.

Spróbujmy:

function GcdI(u, v: Cardinal): Cardinal;
var
  t: Integer;
begin
 while v <> 0 do
 begin
   t := u mod v;
   u := v;
   v := t;
 end;
 Result := u;
end;

GcdI(1000, 2132424);
Średni czas wykonywania: 360 cykli zegara

Całkiem nieźle. Zyskaliśmy jakieś 10% szybkości.

Jak jeszcze można to przyspieszyć? Może napisać wszystko w assemblerze? Choć kompilatory zazwyczaj doskonale optymalizują kod, to nie ma to jak wkład człowieka:

function GcdASM(u, v: Cardinal): Cardinal;
asm
  {w edx jest v, a w eax jest u}
  push ebx //jakby coś tam było to odkładamy na stos
  mov ecx, edx  //przenosimy v do ecx, bo będzie nam potrzebne edx
  //jeżeli v=0 to koniec
  test ecx, ecx
  jz @koniec
  @poczatek:
  xor edx, edx //tu będzie reszta, więc zerujemy rejestr
  div ecx //dzielimy przez v. Dzielona liczba jest w eax, dzielnik w ecx
  //t := u mod v. t to jest nasze edx
  mov eax, ecx //u := v
  mov ecx, edx //v := t
  //jeżeli v<>0 to poczatek
  test ecx, ecx
  jnz @poczatek
  @koniec:
 //Result := u; czyli eax := u;
  pop ebx //ściągamy pierwotną wartość ze stosu
end;

GcdOd(1000, 2132424);
Średni czas wykonywania: 348 cykli zegara (i jeszcze 1%)

Jak widać algorytm jest bardzo pomysłowy, ale można go zepsuć źle pisząc funkcję...