Największy Wspólny Dzielnik - algorytm Euklidesa

Dryobates

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ę...

12 komentarzy

rafcio88_put: Przecież NWD(12, 24, 6) = NWD(12, NWD(24, 6)), a to już prosto zaimplementować.

  1. Pobierz tablicę A
  2. Ustaw wartość zmiennej I na długość tablicy A
  3. Póki I > 0 (1 w Delphi)
    a. ustaw wartość element tablicy A o numerze I-1 na NWD(wartość elementu tablicy A o numerze I, wartość elementu tablicy A o numerze I-1)
    b. usuń element tablicy a o numerze I
    c. zmniejsz I o 1
  4. Zwróć wartość pierwszego elementu tablicy A

Może zaciekawi was taka oto ultaszybka wersja gcd:

template <long x, long y>
struct gcd
{
static const long value = gcd<y,x%y>::value;
};

template<long x="x">
struct gcd<x,0>
{
static const long value = x;
};

int main(int argc, char* argv[])
{
long y = gcd<1000, 2132424>::value;
return 0;
}

bati87: no nie rozśmieszaj, przecież to prawie jak pseudokod, tylko pomyśleć trzeba trochę...

1) function GcdR(u, v : Cardinal):Cardinal zamień na np.: int GcdR(int u, int v)
2) begin i end zamień na { i }
3) if v = 0 then zamień na if(v == 0)
4) Result := GcdR(v, u mod v); zamień na return GcdR(v, u % v);

Reszta jest analogicznie.

A może ktoś w javie napisać?

A czy ktoś jest tak dobry, że potrafi napisac program (w Delphi) do wyznaczania NWD n liczb naturalnych? Użytkownik podaje liczbę naturalną n < 20, oraz zbiór n liczb naturalnych, dla których znajdowany będzie NWD.

[quote]function nwd (a, b : integer) : integer;
begin
while not (a = b) do
if a > b then
a := a - b
else
b := b - a;
nwd := a;
end;
// tez dziala ale nie wiem jak szybko sie wykonuje[/quote]

to jest to samo co na poczatu, tyle, ze bez rekurencji

function nwd (a, b : integer) : integer;
begin
while not (a = b) do
if a > b then
a := a - b
else
b := b - a;
nwd := a;
end;
// tez dziala ale nie wiem jak szybko sie wykonuje

a wspólny przez u? :)

Ciekawe, przyda się.

A Ty Johny rzeczywiście jesteś Wielbłond, a raczej Wielkibłąd. Od kiedy słówko rzecz pisze się przez ż??

Największy wspulny dzielnik/wielokrotność to najgłupsza żecz w matematyce jaką kiedykolwiek widziałem.
Szkoda że nie interesowałem się programowaniem w czasie gdy to przerabiałem w szkole.
Ale itak mi sie to prztyda.
Thanks

Hmmm przyda sie w szkole;)