świstak.codes

O programowaniu, informatyce i matematyce przystępnym językiem

Podstawy algorytmiki: największy wspólny dzielnik

Obliczanie największego wspólnego dzielnika dwóch liczb to prawdopodobnie jeden z pierwszych algorytmów, które poznajemy podczas swojej edukacji. Zarazem jest to też bardzo proste do zapamiętania. Omówmy sobie to klasyczne podejście, a także różne inne podejścia.

Największy wspólny dzielnik

Zanim zaczniemy poznawać algorytmy znajdujące największy wspólny dzielnik, najpierw poznajmy definicję, co to jest. Jak sama nazwa wskazuje, jest to największa liczba naturalna dzieląca wskazane liczby całkowite. Funkcję, która znajdzie taką wartość dla liczb aa i bb, zapiszemy jako NWD(a,b)\text{NWD}(a,b). Dla ułatwienia sobie życia będę używać w tekście skrótu NWD zamiast pełnej nazwy. W literaturze angielskiej, która w świecie informatycznym dominuje, problem ten nazywa się greatest common divisor, a funkcję zapisujemy wówczas jako GCD(a,b)\text{GCD}(a,b) (czasami greatest common factor GCF\text{GCF}).

Z ciekawych własności funkcji NWD można wyróżnić to, że:

  • jest przemienna, czyli NWD(a,b)=NWD(b,a)\text{NWD}(a,b) = \text{NWD}(b,a)
  • jest łączna, czyli NWD(a,NWD(b,c))=NWD(NWD(a,b),c)\text{NWD}(a,\text{NWD}(b,c)) = \text{NWD}(\text{NWD}(a,b),c)
  • NWD(a,0)=a\text{NWD}(a,0) = \lvert a \rvert
  • dla nieujemnej liczby całkowitej mm mamy własność: NWD(ma,mb)=mNWD(a,b)\text{NWD}(m \cdot a, m \cdot b) = m \cdot \text{NWD}(a,b)

Największy wspólny dzielnik jest też powiązany z innym pojęciem matematycznym — najmniejszą wspólną wielokrotnością (NWW\text{NWW}). Znając sposób obliczania NWD, możemy ją obliczyć następująco:

NWW(a1,a2)=a1a2NWD(a1,a2)NWW(a1,a2,...,an)=NWW(a1,NWW(a2,a3,...,an))\text{NWW}(a_1,a_2) = \frac{a_1 \cdot a_2}{\text{NWD}(a_1,a_2)} \\ \text{NWW}(a_1,a_2, ..., a_n) = \text{NWW}(a_1,\text{NWW}(a_2, a_3, ..., a_n))

Znalezienie NWD przez rozkład na czynniki pierwsze

Zanim jednak omówimy to klasyczne podejście, które każdy zna, zacznijmy od podejścia pokazującego wprost, czym jest największy wspólny dzielnik. Mianowicie jest to podejście opierające się na rozkładzie liczb na czynniki pierwsze.

Rozkład liczby na czynniki pierwsze to znalezienie dzielników liczby, które są liczbami pierwszymi. Cała idea znalezienia największego wspólnego dzielnika z jego użyciem wygląda tak:

  1. Wykonujemy rozkład na czynniki pierwsze liczb, których NWD chcemy znaleźć.
  2. Znajdujemy część wspólną pośród czynników pierwszych liczb.
  3. Mnożymy wspólne czynniki ze sobą, w ten sposób obliczając NWD.

Dla przykładu znajdźmy NWD(80,180)\text{NWD}(80, 180). Pisemnie najłatwiej jest to zrobić, wykonując po kolei dzielenia przez najmniejszą liczbę pierwszą. Dla liczby 80 wygląda to następująco:

80/2=4040/2=2020/2=1010/2=55/5=180=2222580 / 2 = 40 \\ 40 / 2 = 20 \\ 20 / 2 = 10 \\ 10 / 2 = 5 \\ 5 / 5 = 1 \\ 80 = 2 \cdot 2 \cdot 2 \cdot 2 \cdot 5

Analogicznie postąpimy dla liczby 180:

180/2=9090/2=4545/3=1515/3=55/5=1180=22335180 / 2 = 90 \\ 90 / 2 = 45 \\ 45 / 3 = 15 \\ 15 / 3 = 5 \\ 5 / 5 = 1 \\ 180 = 2 \cdot 2 \cdot 3 \cdot 3 \cdot 5

Możemy zauważyć, że wspólne czynniki pierwsze to 2, 2 i 5. W takim razie:

NWD(80,180)=225=20\text{NWD}(80,180) = 2 \cdot 2 \cdot 5 = 20

Pokazane powyżej działania możemy zapisać następująco językiem matematyki, tym samym tworząc formalną definicję funkcji NWD:

a=2a23a35a57a711a11...=p pierwszepapNWD(a,b)=p pierwszepmin(ap,bp)a = 2^{a_2}3^{a_3}5^{a_5}7^{a_7}11^{a_{11}}... = \prod_{p \text{ pierwsze}} p^{a_p} \\ \text{NWD}(a,b) = \prod_{p \text{ pierwsze}} p^{\min(a_p,b_p)}

Pierwsza linia to definicja rozkładu liczby na czynniki pierwsze, natomiast druga to bazująca na niej definicja NWD.

Warto jednak zauważyć, że obliczanie NWD w ten sposób jest niepraktyczne, dlatego też nie będziemy go implementować.

Algorytm Euklidesa

Przejdźmy więc do tej właściwej, klasycznej metody znanej jako algorytm Euklidesa. Jego autorstwo przypisuje się Euklidesowi, ponieważ najwcześniejszy znany zapis pochodzi z jego dzieła Elementy z IV wieku p.n.e. Prawdopodobnie jednak Euklides jedynie skompilował wówczas znaną wiedzę w jedno dzieło, więc algorytm może być jeszcze starszy.

Wersja z odejmowaniem

Oryginalny algorytm Euklidesa działał tylko dla dodatnich liczb całkowitych i można go rozpisać jak poniżej.

Na wejściu algorytm otrzymuje liczby aa i bb, natomiast wyjściem jest ich największy wspólny dzielnik.

  1. Tak długo, jak aba \neq b:
    1. Jeśli a>ba > b, przypisz aa wartość aba - b.
    2. W przeciwnym wypadku przypisz bb wartość bab - a.
  2. Zwróć aa.

Oznacza to, że w kodzie zapiszemy go w taki oto krótki sposób:

function gcd(a, b) {
  while (a !== b) {
    if (a > b) {
      a = a - b;
    } else {
      b = b - a;
    }
  }
  return a;
}

Działanie możesz przetestować w serwisie repl.it.

Algorytm opiera się na tym, że największy wspólny dzielnik dwóch liczb nie zmienia się, gdy od większej odejmujemy mniejszą. Na wcześniej pokazanym przykładzie możemy łatwo obliczyć, że NWD(80,180)=NWD(80,100)=NWD(80,20)=NWD(40,20)=NWD(20,20)=20\text{NWD}(80,180) = \text{NWD}(80,100) = \text{NWD}(80,20) = \text{NWD}(40,20) = \text{NWD}(20,20) = 20. Wynikiem jest 20, bo oczywiście największym dzielnikiem liczby jest ona sama.

Wersja z modulo

Poprzednia wersja ma dwie wady. Po pierwsze wymaga dużo obliczeń, a po drugie liczby muszą być dodatnie całkowite. Chcielibyśmy, aby nasz algorytm nie wymagał tylu operacji, a po drugie, żeby operował na całym przedziale liczb całkowitych.

W celu zmniejszenia liczby operacji spójrzmy jeszcze raz na to, jak działa poprzednio opisany algorytm. Warto zauważyć w nim, że gdy a>ba > b i odejmujemy bb odpowiednio długo, to tak naprawdę otrzymujemy resztę z dzielenia aa przez bb, czyli a mod ba \text{ mod } b. Zmieniając działanie na modulo, zyskujemy także obsługę liczb ujemnych. Dzięki temu algorytm można zapisać następująco (zakładając to samo wejście i wyjście co poprzednio):

  1. Tak długo, jak b0b \neq 0:
    1. Zapamiętaj dotychczasową wartość bb.
    2. Przypisz bb wynik działania a mod ba \text{ mod } b.
    3. Przypisz aa poprzednią wartość bb.
  2. Zwróć aa.

Zanim przejdziemy do implementacji, wyjaśnijmy sobie różnice. Zacznijmy najpierw od tego, dlaczego nie sprawdzamy, która liczba jest większa. Otóż w przypadku jeśli a<ba < b, to a mod b=aa \text{ mod } b = a, a w następnej iteracji wartości zostaną zamienione miejscami, więc wyjdzie na to samo. Musimy wykonać jedną iterację więcej, ale z drugiej strony nie musimy robić warunków komplikujących algorytm. Następną rzeczą jest zmiana warunku pętli będąca naturalnym następstwem tego, że gdy aa jest podzielne przez bb, to reszta z dzielenia wynosi 0. Gdybyśmy w poprzedniej wersji z odejmowaniem w momencie zrównania liczb wykonali odejmowanie, to też otrzymalibyśmy 0, więc algorytmy są tożsame.

Implementacja wygląda następująco:

function gcd(a, b) {
  while (b !== 0) {
    const temp = b;
    b = a % b;
    a = temp;
  }
  return Math.max(a, -a);
}

Działanie możesz przetestować w serwisie repl.it.

Obserwując implementację, możesz się zastanawiać, skąd nagle wzięło się branie wartości maksymalnej z aa i a-a. Otóż wynika ono z definicji reszty z dzielenia. W matematycznej (euklidesowej) definicji operacja modulo zawsze zwraca liczbę dodatnią, niezależnie od znaku liczb biorących udział w operacji. Niestety, większość języków programowania tego nie przestrzega, co opisałem szczegółowo w artykule Dziwny przypadek reszty z dzielenia.

Wersja rekurencyjna

Połączmy teraz w jedno to, co pokazałem przy obu wersjach, aby uzyskać jeszcze inną wersję — rekurencyjną. Jak pokazałem przy wersji z odejmowaniem, wersja ta działa dlatego, że odejmując liczbę mniejszą od większej, wartość NWD pozostaje taka sama. Skoro metoda rekurencyjna jest uproszczeniem metody z odejmowaniem, to zachodzi tutaj dosłownie ta sama własność. Powtarzając ten sam przykład, mamy teraz sytuację następującą:

NWD(80,180)=NWD(180,80)=NWD(80,20)=NWD(20,0)=20\text{NWD}(80,180) = \text{NWD}(180,80) = \text{NWD}(80,20) = \text{NWD}(20,0) = 20

Ostatnia część wynika oczywiście z własności NWD, o której wspomniałem na początku artykułu. Natomiast oprócz tego zauważmy, że operacje wyżej pokazane możemy zapisać jako:

NWD(a,b)=NWD(b,a mod b)\text{NWD}(a, b) = \text{NWD}(b, a \text{ mod } b)

Dopisując do tego własność związaną z liczeniem NWD, gdy jedna z liczb jest zerem, możemy ułożyć bardzo prosty wzór rekurencyjny:

{NWD(a,0)=aNWD(a,b)=NWD(b,a mod b)\begin{cases} \text{NWD}(a, 0) = a \\ \text{NWD}(a, b) = \text{NWD}(b, a \text{ mod } b) \end{cases}

W kodzie będzie to wyglądać następująco:

function gcd(a, b) {
  if (b === 0) {
    return Math.max(a, -a);
  } else {
    return gcd(b, a % b);
  }
}

Działanie możesz przetestować w serwisie repl.it.

Binarne NWD

Innym podejściem do obliczania NWD jest algorytm binarne NWD (ang. binary GCD), znany też jako algorytm Steina. W swojej aktualnej formie został opublikowany przez J. Steina w 1967 r., aczkolwiek prawdopodobnie podobne podejście było znane już w II w. p.n.e. w Chinach (jest spór, jak należy rozumieć niedokładny opis zachowany do naszych czasów).

Algorytm bazuje na wykorzystaniu poniższych obserwacji (zakładamy, że wejściem są liczby aa i bb):

  1. Jeśli zarówno aa, jak i bb są parzyste, wówczas: NWD(a,b)=2NWD(a/2,b/2)\text{NWD}(a,b) = 2\cdot\text{NWD}(a/2,b/2)
  2. Jeśli aa jest parzyste, a bb nieparzyste, to: NWD(a,b)=NWD(a/2,b)\text{NWD}(a,b) = \text{NWD}(a/2,b)
  3. Jeśli obie liczby są nieparzyste, to: NWD(a,b)=NWD(ab/2,b)\text{NWD}(a,b) = \text{NWD}(|a-b|/2,b)

Jak widać powyżej, większość z nich opiera się na mnożeniu lub dzieleniu przez 2, co w systemie binarnym możemy bardzo łatwo zaimplementować za pomocą przesunięć bitowych, co opisałem w artykule 1 0 0 0? 0 1 0 1! 1 0 0 1 – czyli matematyka zero-jedynkowa. Zresztą od tego przystosowania pod system binarny wywodzi się nazwa algorytmu.

Zanim przejdziemy do implementacji, warto zaznaczyć, że algorytm działa tylko dla dodatnich liczb całkowitych.

Podstawowa wersja

Kroki podstawowej wersji algorytmu są następujące (znowu wejście i wyjście są takie same):

  1. Ustaw wartość gg na 1. Wartość ta oznacza, przez ile należy pomnożyć wynik, co wynika z obserwacji nr 1.
  2. Tak długo, jak aa i bb są parzyste:
    1. Przypisz a=a/2a = a/2.
    2. Przypisz b=b/2b = b/2.
    3. Przypisz g=g2g = g \cdot 2.
  3. Tak długo, jak a>0a > 0:
    1. Jeśli aa jest parzyste, przypisz a=a/2a = a/2.
    2. W przeciwnym wypadku, jeśli bb jest parzyste, przypisz b=b/2b = b/2.
    3. W przeciwnym wypadku (aa i bb nieparzyste) wykonaj:
      1. Jeśli a<ba < b, przypisz b=ab/2b = |a-b|/2.
      2. W przeciwnym wypadku przypisz a=ab/2a = |a-b|/2.
  4. Zwróć wynik mnożenia bgb \cdot g.

Implementacja w JavaScript wygląda następująco:

function gcd(a, b) {
  let g = 1;
  while (a % 2 === 0 && b % 2 === 0) {
    a = a / 2;
    b = b / 2;
    g = g * 2;
  }
  while (a > 0) {
    if (a % 2 === 0) {
      a = a / 2;
    } else if (b % 2 === 0) {
      b = b / 2;
    } else {
      const temp = Math.abs(a - b) / 2;
      if (a < b) {
        b = temp;
      } else {
        a = temp;
      }
    }
  }
  return b * g;
}

Działanie możesz przetestować w serwisie repl.it.

Wersja z przesunięciami bitowymi

Jak zaznaczyłem wcześniej, to, czym wyróżnia się ten algorytm, to użycie dzieleń i mnożeń przez 2, co można w systemie binarnym bardzo łatwo zaimplementować przesunięciami bitowymi. Jak wiemy, liczby w komputerze są właśnie w ten sposób reprezentowane, więc właśnie tak możemy niskopoziomowo zoptymalizować algorytm.

Nie będę podawać od nowa listy kroków, a jedynie zmiany, które należy zrobić:

  • W punkcie 2.3 zamiast mnożenia będziemy zwiększać gg o 1 co iterację. W ten sposób zliczymy, przez którą potęgę liczby 2 musimy pomnożyć wynik. Z tego też powodu w punkcie 1 inicjujemy zmienną zerem.
  • Wszystkie dzielenia przez 2 zastępujemy przesunięciem w prawo o 1 (>>1>> 1).
  • W punkcie 4 zamiast bgb \cdot g wykonujemy b<<gb << g (jest to binarny odpowiednik działania b2gb \cdot 2^{g}).

Implementacja wygląda jak poniżej. Co prawda została zrobiona w JavaScript, gdzie takie niskopoziomowe optymalizacje nie mają sensu, ale analogicznie można zaprogramować to w dowolnym innym języku.

function gcd(a, b) {
  let g = 0;
  while ((a & 1) === 0 && (b & 1) === 0) {
    a = a >> 1;
    b = b >> 1;
    g++;
  }
  while (a > 0) {
    if ((a & 1) === 0) {
      a = a >> 1;
    } else if ((b & 1) === 0) {
      b = b >> 1;
    } else {
      const temp = Math.abs(a - b) >> 1;
      if (a < b) {
        b = temp;
      } else {
        a = temp;
      }
    }
  }
  return b << g;
}

Działanie możesz przetestować w serwisie repl.it. W implementacji oprócz wyżej opisanej optymalizacji związanej z wykorzystaniem przesunięć bitowych zmieniłem też sposób sprawdzania parzystości na sprawdzanie wartości ostatniego bitu (liczby parzyste zawsze mają 0).

Podsumowanie

W artykule pokazałem trzy różne sposoby, jak możemy obliczać największy wspólny dzielnik. Są to bardzo podstawowe i proste algorytmy, ale nawet takie warto pisać na własną rękę dla poszerzenia swojej wiedzy i rozwinięcia umiejętności programistycznych. Warto też dodać, że na bazie algorytmu Euklidesa bazuje rozszerzony algorytm Euklidesa, który swoje zastosowanie znalazł w kryptografii.

Literatura

  • Cormen, T. H.; Leiserson, C. E.; Rivest R. L.; Stein C. “Greatest common divisor” w Introduction to algorithms, 3rd ed.. MIT Press, Cambridge, MA, U.S.A., s. 933-939.
  • Knuth, D. E. “The Greatest Common Divisor” w The art of computer programming: Volume 2.. Addison-Wesley, 2011, s. 333-356
  • Paul E. Black, "binary GCD", in Dictionary of Algorithms and Data Structures [online], Paul E. Black, ed. 2 November 2020. (ostatni dostęp: 22.11.2022) Dostępne na: https://www.nist.gov/dads/HTML/binaryGCD.html
Zdjęcie na okładce wygenerowane przez DALL-E. Oryginał znajduje się tutaj.