Wprowadzenie do Teorii Algorytmów
(Introduction to Algorithms Theory)
Prof. Dr hab. Alexander Prokopenya
Szkoła Główna Gospodarstwa Wiejskiego w Warszawie
Katedra Zastosowań Informatyki
Wykład 4-5. Algorytmy  dziel i zwyciężaj
Pry użyciu strategii  dziel i zwyciężaj (angl. Divide and conquer) problem roz-
wiązuje się, postępując zgodnie z następującym schematem:
1. Dzielimy problem na podproblemy, z których każdy jest mniejszym egzempla-
rzem problemu tego samego typu.
2. Kolejno rozwiązujemy te podproblemy.
3. Lączymy odpowiednio uzyskane rozwiązania.
Zadanie jest realizowane stopniowo, na trzech różnych etapach: pierwszy polega na
dzieleniu problemów na podproblemy; drugi ma miejsce, gdy podproblemy są już tak ma-
łe, że ich rozwiązania są natychmiastowe i nie wymagają stosowania rekursji; natomiast
na samym końcu rozwiązania częściowe skleja się w calość. Wszystko to jest skoordyno-
wane dzięki rekurencyjnej strukturze algorytmu.
Jako przykład wprowadzający zobaczymy, w jaki sposób, stosując tę technikę, uz y-
skamy nowy algorytm mnożenia liczb, który jest znacznie bardziej efektywny niż metoda
mnożenia, z jaką spotkaliśmy w skole.
1. Mnożenie
Pewnego razu matematyk Carl Friedrich Gauss (1777  1855) zauważył, że choć wydaje się, że
iloczyn dwóch liczb zespolonych
(a +� b��i)(c +� d ��i) =� ac -�bd +� (bc +� ad)��i
wymaga czterych mnożeń liczb rzeczywistaych, w istocie można go otrzymać, wykonyjąc zale-
dwie trzy mnożenia: ac , bd oraz (a +� b)(c +� d) , gdyż
bc +� ad =� (a +� b)(c +� d) -� ac -� bd .
Rozpatrując to punktu widzenia notacji O, zdawać by się mogło, że zredukowanie liczby
mnożeń z czterech do trzech jest grą niewartą świeczki. Jednak to niewielkie ulepszenie staje się
znaczące, gdy jest stosowane rekurencyjnie.
Porzućmy liczby zespolone i zobaczmy, jak to może pomoc w zwykłym mnożeniu. Niech
x i y będą n-bitowami liczbami natutalnymi. Założmy też dla wygody, że n jest potęga dwójki
(ogólny przypadek jest bardzo podobny). Pierwszym krokiem w stronę wymnożenia x oraz y bę-
dzie rozbicie obu liczb na ich lewe i prawe połowy, z których każda ma n / 2 bitów:
x =� 2n / 2 xL +� xR , y =� 2n / 2 yL +� yR .
Na przykład, jeśli x =� 101101102 (indeks dolny 2 oznacza zapis dwójkowy), to xL =� 10112 ,
xR =� 01102 oraz x =� 10112 ��24 +� 01102 . Iloczyn x i y może być zatem zapisany jako
xy =� (2n / 2 xL +� xR)(2n / 2 yL +� yR) =� 2n xLxR +� 2n / 2(xL yR +� xR yL) +� xR yR .
1
Obliczymy xy , korzystając z wyrażenia po prawej stronie. Czas wykonania dodawania
jest liniowy, podobnie jak mnożenie przez potęgi 2 (które jest po prostu przesunięciem w lewo).
Znaczącymi operacjami są cztery n / 2 -bitowe mnożenia: xL yL , xL yR , xR yL oraz xR yR , z któ-
rymi można sobie poradzić za pomocą czterych wywołań rekurencyjnych. Zatem nasza metoda
mnożenia n-bitowych liczb zaczyna się od wykonania wywołań rekurencyjnych w celu wymno-
żenia tych czterych par n / 2 -bitowych liczb (cztery podproblemy dla danych o połowę mniej-
szych), po czym następuje obliczenie wartości wyrażenia po prawej stronie równania w czasie
O(n) . Oznaczając przez T(n) całkowity czas działania algorytmu na n-bitowych danych wej-
ściowych, dostajemy zależność rekurencyjną
T(n) =� 4T(n/ 2) +� O(n) .
Wkrótce zobaczymy, jak wyglądają ogólne metody rozwiązywania takich równań. Tym-
czasem rozwiązanie tego konkretnego równania okazuje się być rzędu O(n2) , czyli czas działa-
nia jest taki sam jak dla tradycyjnej , szkolnej techniki mnożenia. Mamy zatem zupełnie nowy
algorytm, ale nie osiągneliśmy żadnego postępu w efektywności. W jaki sposób możemy jeszcze
usprawnić tę metodę?
W tym momencie przychodzi nam na myśl sztuczka Gaussa. Chociaż obliczenie iloczynu
xy zdaje się wymagać czterech n / 2 -bitowych mnożeń, tak jak już było pokazane wcześniej,
wystarczą zaledwie trzy: xL yL , xR yR oraz (xL +� xR)(yL +� yR) , gdyż
xL yR +� xR yL =� (xL +� xR)(yL +� yR) -� xL yL -� xR yR .
Czas działania wynikającego stąd algorytmu, przedstawionego na rysunku poniżej, jest krótszy:
T(n) =� 3T(n/ 2) +� O(n) .
Rzecz w tym, że poprawę stałego czynnika z 4 na 3 wykożystuje się na każdym poziomie rekur-
sji, a zwielokrotnienie tego wynika prowadzi do znacznego zmniejszenia ograniczenia czasu
działania  aż do O(nlog 3) =� O(n1.59) .
funkcja multiply(x, y)
// wejście: dwie n-bitowe liczby naturalne x i y
// wyjście: ich iloczyn
if n =�1 then return xy
xL , xR =� lewe / 2 prawe / 2 bitów x
��n ł�, ��n ��
yL , yR =� lewe / 2 prawe / 2 bitów y
��n ł�, ��n ��
P1 = multiply( xL , yL )
P2 = multiply( xR , yR )
P3 = multiply( xL +� xR , yL +� yR )
return P1 ��2n +� (P3 -� P1 -� P2)��2n / 2 +� P2
Czas działania algorytmu można uzyskać po analize równania rekurencyjnego
T(n) =� aT(n/ 2) +� cn .
Przy podstawianiu n =� 2k otrzymamy
T(2k ) =� aT(2k -�1) +� c �� 2k =� a(aT(2k -�2) +� c �� 2k -�1) +� c �� 2k =�
=� a2T(2k -�2) +� ac ��2k -�1 +� c��2k =� a2(aT(2k -�3) +� c��2k -�2) +� ac��2k -�1 +� c��2k =�
=� a3T(2k -�3) +� a2c �� 2k -�2 +� ac �� 2k -�1 +� c �� 2k =� L�=�
2
=� akT(20) +� ak -�1c��21 +� ak -�2c��22 +�L�+� ac ��2k -�1 +� c��2k =�
1-� (2/ a)k
=� ak +� ak -�1c��2(1+� 2/ a +� (2/ a)2 +�...+� (2/ a)k -�1) =� ak +� ak -�1c��2�� =�
1-� 2/ a
ak -� 2k 2c 2c 2c 2cn
ć�1+� ��
.
=� ak +� c��2�� =� ak ć�1+� �� -� �� 2k =� ��nlog a -�
�� �� �� ��
a -� 2 a -� 2 a -� 2 a -� 2 a -� 2
Ł� ł� Ł� ł�
Dla szkolnej techniki mnożenia mamy a =� 4 i wtedy
T(n) =� 4T(n/ 2) +� cn =� (1+� c)n2 -� cn =� O(n2) .
Zatem dla a =� 3
T(n) =� 3T(n / 2) +� cn =� (1+� 2c)nlog 3 -� 2cn =� O(nlog 3) ,
w przybliżeniu jest równy O(n1.59) .
W algorytmach  dziel i zwyciężaj liczba podproblemów odpowiada współczynnikowi
rozgałęzenia dzewa rekursji, niewielkie zmiany tego wspólczynnika mogą mieć wielki wpływ na
czas działania.
Uwaga praktyczna: w ogólnym przypadku nie ma potrzeby stosować rekursji aż do mo-
mentu osiągnięcia 1 bitu. Dla większości procesorów mnożenie 16- lub 32-bitowe jest pojedyn-
czą operacją, więc dopóki liczby nie wychodzą poza ten zakres, należy korzystać z wbudowa-
nych procedur.
Na koniec odwieczne pytanie: Czy możemy znalez
ć lepsze zorwiązanie? Okazuje się, że
istnieje szybszy algorytm mnożenia liczb oparty na innym ważnym algorytmie typu  dziel i
zwyciężaj : szybkiej transformacie Fouriera, o której pomówimy pózniej.
2. Zależności rekurencyjne
Algorytmy  dziel i zwyciężaj często działają według ogólnego schematu: radzą sobie z proble-
mem o rozmiarze n poprzez rekurencyjne rozwiązywanie a podproblemów o rozmiarze n / b
każdy, a następnie łączenie tych rozwiązań w czasie O(nd ) , dla pewnych a,b,d >� 0 (w algo-
rytmie mnożenia a =� 3, b =� 2, d =� 1). Ich czas działania może być zatem wyrażony równaniem
T(n) =� aT( +� O(nd ) . Wyprowadzimy teraz zwarte rozwiązanie tej rekurencji, aby w przy-
��n/bł�)
szłości nie musieć rozwiązywać jej dla nowych danych.
Twerdzenie o rekurencji universalnej. Jeśli T(n) =� aT( +� O(nd ) dla pewnych stałych
��n/bł�)
a >� 0 , b >� 1 oraz d ł� 0, to
��
O(nd ), gdy d >� logb a
��
d
T (n) =� log n) gdy d =� logb a
��O(n
��
b
O(nlog a ) gdy d <� logb a
��
To jedno twierdzenie mówi nam o czasach działania większości programów typu  dziel i
zwyciężaj , z którymi się spotkamy.
Dowód. Aby udowodnić tezę, dla wygody zacznijmy od założenia, że n jest potęga b: n =� bk .
Nie wpłynie to w istotny sposób na ostateczny wynik  w końcu n różni się od pewnej potęgi b
co najwyżej o stały czynnik  pomoże nam to natomiast zignorować zaokrąglanie / b
��n ł�.
T(bk ) =� aT(bk -�1) +� c��bkd =� a(aT(bk -�2) +� c��b(k -�1)d ) +� c��bkd =�
=� a2T(bk -�2) +� ac ��b(k -�1)d +� c ��bkd =� a2(aT(bk -�3) +� c ��b(k -�2)d ) +� ac ��b(k -�1)d +� c ��bkd =�
=� a3T(bk -�3) +� a2c ��b(k -�2)d +� ac ��b(k -�1)d +� c ��bkd =� L�=�
3
=� akT(b0) +� ak -�1c��bd +� ak -�2c��b2d +�L�+� ac��b(k -�1)d +� c��bkd =�
ć�
bd b2d b(k -�2)d b(k -�1)d ��
�� ��
=� akT(1) +� ak -�1c ��bd ��1+� +� +�L�+� +� .
a
a2 ak -�2 ak -�1 ��
Ł� ł�
Gdy d =� logb a , otrzymujemy bd =� blogb a =� a , zatem
T(bk ) =� akT(1) +� akc �� k =� bk logb a(T(1) +� c �� k)
T(n) =� nlogb a(T(1) +� c��logb n) =� O(nd log n) .
Gdy d ą� logb a i k zmierza się od 0 do logb n , wartoście w nawiasie tworzą ciąg geometryczny
o ilorazie bd / a . Znaczenie sumy k wyrazów tego ciągu może być zapisane jako
1-� bkd / ak
.
1-� bd / a
Zatem otrzymujemy
1-� bkd / ak ć� c ��bd �� c ��bd
�� ��
T(bk ) =� akT(1) +� ak -�1c ��bd �� =� ak ��T(1) +� -� ��bkd .
1-� bd / a a -� bd �� a -� bd
Ł� ł�
ć� ��
c ��bd �� c ��bd O(nd ), gdy d >� logb a
b
�� ��
T (n) =� nlog a��T(1) +� -� �� nd =�
��
a -� bd �� a -� bd ��O(nlogb a ), gdy d <� logb a
Ł� ł�
3. Mnożenie macierzy
Iloczynem dwóch macierzy n��n , X oraz Y, jest trzecia macierz n��n , Z =� XY , której element o
współrzędnych (i, j) jest dany wzorem
n
Zij =� XikYkj .
��
k =�1
Powyższy wzór implikuje algorytm mnożenia macierzy działający w czasie O(n3) ; należy obli-
czyć wrtości n2 elementów, a każde obliczenie zajmuje czas O(n) .
W 1969 roku niemieecki matematyk Volker Strassen ogłosił znacznie bardziej efektywny
algorytm, oparty na strategii  dziel i zwyciężaj . Mnożenie macierzy bardzo łatwo można po-
dzielić na podproblemy, ponieważ może być wykonywane blokowo. Aby zobaczyć, co to ozna-
cza, podzielmy X na cztery bloki o rozmiarach n / 2��n/ 2, to samo zróbmy dla Y:
A B E F
�� ł� �� ł�
X =� , .
ę�C Dś� Y =� ę�G H ś�
�� �� �� ��
Wtedy ich iloczyn może być wyrażony za pomocą tych bloków, dokładnie tak, jakby te bloki
bały pojedynczymi elementami:
A B E F AE +� BG AF +� BH
�� ł��� ł� �� ł�
XY =� =� .
ę�C Dś�ę�G H ś� ę�CE +� DG CF +� DHś�
�� ���� �� �� ��
Stosujemy teraz strategię  dziel i zwyciężaj : aby obliczyć iloczyn macierzy XY o roz-
miarze n, rekurencyjnie obliczamy osiem iloczynów AE, BG, CE, DG, AF, H, CF, DH macierzy
o rozmiarze n / 2 każda, a następnie wykonujemy kilka dodawań w czasie O(n2) . Całkowity
czas działania jest opisany zależnością rekurencyjną:
4
T(n) =� 8T(n / 2) +� O(n2) .
To prowadzi do mało imponującego rozwiązania O(n3) , takiego samego wynik jak przy
zwykłym algorytmie. Jednak efektywność może zostać ulepszona i podobnie jak przy mnożeniu
liczb całkowitych, kluczowe będzie sprytne wykorzystanie algebry. Okazuje się, że iloczyn XY
można obliczyć przy użyciu zaledwie siedmiu podproblemów o rozmiarach n / 2��n/ 2:
P5 +� P4 -� P2 +� P6 P1 +� P2
�� ł�
XY =� ,
ę�
P3 +� P4 P1 +� P5 -� P3 -� P7 ś�
�� ��
gdzie
P1 =� A(F -� H) , P5 =� (A+� D)(E +� H) ,
P2 =� (A+� B)H , P6 =� (B -� D)(G +� H) ,
P3 =� (C +� D)E , P7 =� (A-�C)(E +� F) ,
P4 =� D(G -� E) .
Nowy czas działania wynosi
T(n) =� 7T(n / 2) +� O(n2) ,
który na podstawie twierdzenia o rekurencji uniwersalnej okazuje się być równy
O(nlog 7) � O(n2.81) .
4. Sortowanie przez scalanie (merge sort)
Sortowanie to problem bardzo często rozwiązywany na komputerach. Jego popularność wiąże
się z faktem, że łatwiej jest korzystać ze zbiorów uporządkowanych niż nieuporządkowanych.
Sortowanie jest dobrym przykładem tego, że określone zadanie może być wykonane według
wielu różnych algorytmów. Każdy z algorytmów ma pewne zalety i wady, które trzeba przeana-
lizować dla konkretnego zastosowania.
Sortowaniem (ang. sorting) nazywamy proces ustawiania zbioru obiektów w określo-
nym porządku. Sortowanie stosuje się w celu ułatwienia póz
niejszego wyszukiwania elementów
sortowanego zbioru.
Niech U będzie zbiorem obiektów
a1, a2,..., an
Sortowanie polega na permutowaniu tych obiektów aż do chwili osiągnięcia uporządkowania
ak1, ak2 ,..., akn
takiego, że dla zadanej funkcji porządkującej f zachodzi
f (ak1 ) Ł� f (ak2 ) Ł� ... Ł� f (akn ) .
Wartość funkcji f nazywa się kluczem obiektu i służy do identyfikacji obiektów.
Zauważmy, że w sformułowanym powyżej problemie sortowania nic nie wiemy o natu-
rze elementów z U. Na U mogą składać się zarówno liczby całkowite lub rzeczywiste, jak i U
może być zbiorem rekordów, które należy posortować według ich kluczy. Przyjmujemy, że ele-
menty ciągu a1, a2,..., an znajdują się w tablicy a[1..n] . Jedynym sposobem ustalenia porządku w
tablicy a jest porównywanie jej elementów parami. Operacja porównania będzie operacją do-
minującą. Ponieważ będziemy chcieli ustalić wynik także w tablicy a, potrzebna nam jest jesz-
cze operacja zamiany dwóch elementów w tablicy. Operacją tą będzie operacja Swap(i, j) pole-
gająca na zamianie elementów w tablicy a z pozycji i oraz j, 1Ł� i, j Ł� n .
5
 Problem posortowania listy liczb prowadzi bezpośrednio do zastosowania strategii  dziel
i zwyciężaj : dzielimy listę na połowy, rekurencyjnie sortujemy każdą połowę, a następnie sca-
lamy (ang. merge) dwie posortowane podlisty (Rys. 1).
Rys. 1. Schemat sortowania przez scalanie
Przykład. Rozważmy ciąg 5,7,4,9,3,6,2,1. Potraktujmy, każdy z elementów tego ciągu jako jed-
noelementowy uporządkowany ciąg. Zastosujemy procedurę scalania do sąsiadujących ciągów,
otrzymując 4 dwuelementowe posortowane ciągi: {5,7}, {4,9}, {3,6}, {1,2}. Teraz ponownie
zastosujemy scalanie sąsiednich ciągów tworząc dwa czteroelementowe segmenty uporządko-
wane {4,5,7,9}, {1,2,3,6}. Wykonanie jeszcze jednego scalania pozwoli nam utworzyć ciąg upo-
rządkowany {1, 2, 3, 4, 5, 6, 7, 9}.
Niech a będzie ciągiem o długości n, a b ciągiem o długości m. W opisanym poniżej al-
gorytmie Merge, oba ciągi wejściowe a i b zostały rozszerzone o jeden element, oznaczony tu-
taj +� Ą� . Zakładamy, że jest on większy od wszystkich innych elementów obu ciągów. Taki za-
bieg pozwoli nam uprościć nieco strukturę algorytmu. Zadbamy jednak o to, by ciąg e otrzyma-
ny w wyniku miał dokładnie n+m elementów.
Algorytm Merge(a[1..n], b[1..m])
if n=0 return b[1..m]
if m=0 return a[1..n]
a[n+1] := +� �Ą�; b[m+1]:= +� �Ą�;
i := 1; j := 1; k :=1;
while (k �Ł� (n+m)) do
if (a[i] �Ł� b[j]) then
e[k] := a[i]; i := i+1
else
e[k] := b[j]; j := j+1
fi;
k := k+1
od;
return e
6
4.1. Poprawność algorytmu
Zauważmy, że elementy dodatkowe, dopisane do ciągów a i b, nie mogą się znalez
ć w ciągu
wynikowym, gdyż pętla jest wykonywana tylko (n+m) razy. Gdyby, w którymś momencie
i =� n +�1, tzn. a[i] = +� Ą� , to przy porównywaniu tego elementu z b[j] dla j <� m +�1, wygra b[j] i
to b[j] zostanie zapisane w ciągu wynikowym. Analogicznie w przypadku, gdy j =� m +�1. Jeśli
zaś oba wskaz
niki i, j wskazują element specjalny, to znaczy, że wszystkie poprzednie elementy
już zostały wpisane do tablicy wynikowej, a więc k musi być w tym momencie równe n +� m +�1.
Warunek pętli "while" nie będzie spełniony i zakończymy wykonywanie algorytmu.
Jest oczywiste, że algorytm zatrzymuje się dla dowolnych danych, gdyż zmienna kontro-
lująca pętlę zmienia się od 1 do n+m. Trzeba jeszcze pokazać, że otrzymany wynik jest ciągiem
uporządkowanym złożonym z elementów danych ciągów.
Przed pierwszym wejściem do pętli "while" jest trywialnie spełniony warunek e[1] Ł�
e[2] Ł�...Ł� e[k-1]. Załóżmy, że rozpoczynamy pewną iterację pętli "while" i spełnione są własno-
ści:
e[1] Ł� e[2] Ł�...Ł� e[k-1], i Ł� n+1, j Ł� m+1, k = (i + j -1), k Ł� (n+m+1), (*)
(ciąg e[1],...,e[k-1] jest permutacją ciągu a[1],...,a[i-1],b[1],...,b[j-1] ).
Formuły te mówią, że k-1 elementów w ciągu e tworzy ciąg uporządkowany, a ponadto i-1 ele-
mentów ciągu a oraz j-1 elementów ciągu b zostało już zapisanych na pozycjach od 1 do k-1
ciągu e.
Po wykonaniu instrukcji warunkowej "if (a[i] Ł� b[j])..." , na miejscu k-tym w ciągu e po-
jawi się mniejszy z elementów a[i], b[j] oraz, albo i, albo j wzrośnie o jeden. Ponieważ ciągi a i
b są uporządkowane niemalejąco zatem wstawiony element e[k] jest niemniejszy od wszystkich
elementów, które znajdują się na wcześniejszych pozycjach w ciągu e. Mamy więc w tablicy e
aż k elementów uporządkowanych, oraz i+j-1 = k+1.
e[1] Ł� ... Ł� e[k-1] Ł� e[k], i Ł� n+1, j Ł� m+1, (k+1) = (i + j -1), k Ł� (n+m+1),
(ciąg e[1],...,e[k] jest permutacją elementów ciągu a[1],...,a[i-1],b[1],...,b[j-1] ).
Po wykonaniu instrukcji przypisania "k:=k+1;" ponownie spełnione są własności (*). Wykazali-
śmy tym samym, że jest to niezmiennik pętli. Z twierdzenia o niezmienniku, ta sama formuła jest
spełniona w chwili wyjścia z pętli. Wtedy jednak k = n+m+1 i wszystkie elementy, zarówno cią-
gu a, jak i b, znalazły się już w ciągu e oraz e[1] Ł� ... Ł� e[n+m]. Spełniony jest zatem warunek
końcowy specyfikacji.
Twierdzenie. Algorytm Merge jest całkowicie poprawnym rozwiązaniem problemu scalania za
względu na specyfikację,
wp = {a[1] Ł� ... Ł� a[n], b[1] Ł� ... Ł� b[m], n>0, m>0}
wk = {ciąg e[1],...,e[n+m] jest uporządkowaną niemalejąco permutacją elementów a[1],..., a[n],
b[1],...,b[m] }
w każdej strukturze danych z liniowym porządkiem Ł�.
4.2. Koszt algorytmu
Ponieważ w pętli "while" wykonujemy w każdej iteracji tylko jedno porównanie, a liczba iteracji
wynosi dokładnie n+m, zatem koszt algorytmu jest liniowy w stosunku do długości ciągów sca-
lanych i wynosi
T(n) = O(n+m).
funkcja mergesort(a[1..n])
// wejście: tablica liczb a[1..n]
7
 // wyjście: tablica wejściwa po posortowaniu
if n >� 1 then
return merge(mergesort( a[1.. / 2 ),mergesort( a[ / 2 ))
��n �� ��n ��+�1..n
else return a
Poprawność tego algorytmu jest oczywista, o ile tylko określona jest poprawna procedura
merge. Funkcja merge wykonuje taką samą pracę przy każdym wywolaniu rekurencyjnym, a
jej całkowity czas działania to O(n) . Zatem procedura scalania jest liniowa, a całkowity czas
potrzebny na wykonanie mergesort to
T(n) =� 2T(n / 2) +� O(n) =� O(nlog n) .
Spoglądając raz jeszcze na mergesort, możemy zauważyć, że cała właściwa praca po-
lega na scalaniu, które jednak nie może się zacząć, dopoki rekurcja nie zejdzie do tablic jednoe-
lementowych. Takie tablice jednoelementowe są scalane parami, co prowadzi do powstania ta-
blic dwuelementowych. Następnie pary takich tablic dwuelementowych są scalane, tworząc czte-
roelementowe itd.
Takie podejście wskazuje również, w jaki sposób mergesort może być wykonywany
iteracyjne. W każdym momencie dany jest zbiór  aktywnych tablic  początkowo są to tablice
jednoelementowe  które scalane parami tworzą nowy zestaw aktywnych tablic. Tablice te mogą
być zorganizowane w kolejkę i przetwarzane przez sukcesywne usuwanie dwóch pierwszych
tablic, scalanie ich i umieszczenie otrzymanej tablicy na końcy kolejki.
5. Mediany
Medianą listy liczb nazywamy pięćdziesiąty percentyl listy: połowa liczb jest od niej większa, a
połowa mniejsza. Na przykład medianą [45,1,10,30,25] jest 25, poniewaz jest to środkowy ele-
ment uporządkowanej listy zlożonej z tych samych liczb. Jeśli lista ma parzystą długość, istnieją
dwie możliwe wartości dla elementu środkowego, wtedy, powiedzmy, wybieramy mniejszą z
nich.
Obliczenie mediany z n liczb jest proste: wystarczy je posortować. Wadą jest czas trwa-
nia algorytmu  O(nlog n) , podczas gdy zdecydowanie wolelibyśmy czas liniowy. Nie powinni-
śmy jednak tracić nadziei, ponieważ podczas sortowania wykonujemy znacznie więcej pracy niż
w rzeczywistości potrzebujemy  chcemy tylko znalez
ć środkowy element i nie interesuje nas
właściwe uporządkowanie pozostalych.
Szukając rozwiązania rekurencyjnego, często paradoksalnie łatwiej jet pracować z ogól-
niejszą wersją problemu  z tego prostego powodu, że możemy się oprzeć na silniejszej rekursji.
W naszym przypadku rozważanym uogólnieniem będzie selekcja.
Selekcja.
Wejście: lista liczb S; liczba naturalna k,
Wyjście: k-ty najmniejszy element S.
Na przykład, jeśli k =� 1, szukane jest minimum S, natomiast jeśli k =� S | / 2
��| ��, to szu-
kamy mediany.
5.1. Randomizowany algorytm  dziel i zwyciężaj dla selekcji
Oto podejście  dziel i zwyciężaj do problemu selekcji. Dla dowolnej liczby v wyobraz
my sobie
podzielenie listy S na trzy kategorie elementów: elementy mniejsze od v, równe v (elementy mo-
gą się powtarać) oraz większe od v. Nazwijmy je odpowiednio SL , Sv , SR . Na przykład, jeśli
tablica
S : 2 36 5 21 8 13 11 20 5 4 1
jest podzielona dla v =� 5, trzy wygenerowane podtablice wyglądają następująco:
8
SL : 2 4 1 Sv : 5 5 SR : 36 21 8 13 11 20
Przeszukiwanie może być natychmiast zawężone do jednej z tych podlist. Jeśli chcemy
znalez
ć, powiedzmy, ósmy najmniejszy element S, wiemy, że będzie to trzeci najmniejszy ele-
ment SR , ponieważ | SL | +� | Sv |=� 5. Czyli selection(S,8) = selection( SR ,3). Ogólnie rzecz bio-
rąc, porównując k z rozmiarem podtablic, możemy szybko stwierdzić, która z nich zawiera po-
szukiwany element:
selection (SL,k) gdy k Ł�| SL |
��
��
selection (S,k) =� v gdy | SL |<� k Ł�| SL | +� | Sv |
��
��selection (SR,k-� | SL | -� | Sv |) gdy k >�| SL | +� | Sv |
��
Trzy podlisty SL , Sv , SR można uzyskać z S w czasie liniowym; w rzeczywistości takie
obliczenie może być nawet wykonane bez przydzielania dodatkowej pamięci. Następnie wyko-
nuje się obliczania rekurencyjnie na właściwej liście. Wynikiem podziału tablicy jest zatem
zmniejszenie liczby elementów z | S | do co najwyżej max{| SL |,| SR |}.
Nasz algorytm  dziel i zwyciężaj dla selekcji jest teraz w pełni opisany, z wyjątkiem
kluczowego szczegółu, jakim jest wybór v. Wartość v musi być wybrana szybko i powinna być
1
taka, by tablica została istotnie zmniejszona, sytuacja jest idealna dla | SL |,| SR |� | S |. Gdyby-
2
śmy potrafili zawsze zagwarantować taką sytuację, dostalibyśmy czas działania
T(n) =� T(n / 2) +� O(n) ,
który zgodnie z naszym życzeniem jest liniowy. To jednak wymaga wybrania na v mediany, co
przecież jest naszym ostatecznym celem. Zamiast tego posłużymy się znacznie prostszym roz-
wiązaniem: wybieramy v z S losowo.
5.2. Analiza efektywności
Oczywiście czas działania naszego algorytmu zależy od losowych wyborów v. Może się zdarzyć,
że każdy wybór będzie pechowy i wartością v będzie największy (lub najmniejszy) element ta-
blicy, a zatem tableca będzie się zmniejszać w każdym kroku tylko o jeden element. We wcze-
śniejszym przykładzie mogliśmy najpierw wybrać v =� 36, następnie v =� 21 itd. Ten czarny sce-
nariusz zmusiłby nasz algorytm selekcji do wykonania
n
n +� (n -�1) +� (n -� 2) +�...+� =� Q�(n2)
2
operacji (podczas obliczania mediany), niemniej takie zdarzenie jest wyjątkowo mało prawdo-
podobne. Równie mało prawdopodobny jest najlepszy przypadek przedyskutowany wcześniej, w
którym każde losowo wybrane v okazuje się rozdzielać tablicę idealnie na pół, skutkując czasem
działania O(n) . Gdzie, w przedziale od O(n) do Q�(n2) , znajduje się średni czas działania? Na
szczęście znajduje się on bardzo blisko czasu działania w najlepszym przypadku.
Aby rozróżnic szczęsliwe wybory v od nieszczęśliwych, powiemy, że v jest dobre, jeśli
znajduje się pomiędzy 25 a 75 percentylem tablicy, z której jest wybierany. Takie wybory v się
nam podobają, ponieważ gwarantują, że podlisty SL oraz SR mają rozmiar równy co najwyżej
trzy czwarte S, a zatem tablica znacznie się zmniejsza. Na szczęście dobre wartości v są liczne:
połowa elementów dowolnej listy musi leżeć między 25 a 75 percentylem.
Zatem średnio po dwóch operacjach podziału tablica zmniejszy się do co najwyżej trzech
czwartych swego wyjściowego rozmiaru. Oznaczając przez T(n) oczekiwany czas działania dla
tablicy o rozmiarze n, dostajemy
T(n) Ł� T(3n/ 4) +� O(n) .
9
Na podstawie tej rekurencji wnioskujemy, że T(n) =� O(n) : dla dowolnych danych wejściowych
nasz algorytm zwraca poprawną odpowiedz
po  średnio  liniowej liczbie kroków.
6. Szybka transformata Fouriera
Widzieliśmy do tej pory, w jaki sposób stosując strategię  dziel i zwyciężaj , można uzyskać
szybkie algorytmy mnożenia liczb całkowitych i macierzy. Naszym kolejnym celem są wielo-
miany. Iloczynem dwóch wielomianów stopnia n jest wielomian stopnia 2n, na przykład
(1+� 2x +� 3x2)��(2 +� x +� 4x2) =� 2 +� 5x +�12x2 +�11x3 +�12x4 .
Ogólniej, jeśli A(x) =� a0 +� a1x +�...+� anxn oraz B(x) =� b0 +� b1x +�...+� bnxn , to ich iloczyn
C(x) =� A(x)�� B(x) =� c0 +� c1x +�...+� c2nx2n ma współczynniki
k
ck =� a0bk +� a1bk -�1 +�...+� akb0 =�
��aibk -�i
i=�0
(dla i >� n bierzemy ai oraz bi równe zero). Obliczenie ck na podstawie tego wzoru wymaga
O(k) kroków, a znalezenie wszystkich 2n +�1 współczynników wydaje się wymagać czasu
Q�(n2) . Czy możemy szybciej mnożyć wielomiany?
Rozwiązanie, którym będziemy się teraz zajmować, szybka transformata Fouriera, zre-
wolucjonizowalo  a właściwe stworzyło  dziedzinę przetwarzania sygnalów. Z uwagi na
ogromne znaczenie i bogactwo zastosowań w różnych dziedzinach badawczych, podejdziemy do
tego zagadnienia nieco dokładniejniż zazwyczaj.
6.1. Alternatywne reprezentacje wielomianów
n
Reprezentacja przez współczynniki wielomianu A(x) =� x stopnia n to wektor współczyn-
��a j j
j=�0
ników a =� (a0,a1,...,an) . Reprezentacja za pomocą współczynników jest dogodna przy niektó-
rych operacjach na wielomianach. Na przykład operacja ewaluacji wielomianu A(x) w danym
punkcie x0 polega na obliczeniu wartości A(x0) . Ewaluacje można wykonać w czasie Q�(n) ,
korzystając t tzw. schematu Hornera:
A(x0) =� a0 +� x0(a1 +� x0(a2 +�...+� x0(an-�1 +� x0(an))...)).
Podobnie, dodawanie dwóch wielomianów reprezentowanych przez wektory współczynników
a =� (a0,a1,...,an) i b =� (b0,b1,...,bn) zajmuje czas Q�(n) : wynik stanowi wektor c =� (c0,c1,...,cn) ,
gdzie cj =� aj +� bj dla j =� 0, 1,..., n.
Reprezentacja przez wartości w punktach wielomianu A(x) stopnia n to zbiór par
punkt-wartość
{(x0, y0), (x1, y1), ..., (xn, yn)}
taki, że wszystkie xk są parami różne oraz yk =� A(xk ) dla k =� 0, 1,..., n . Wielomian może mieć
wiele różnych reprezentacji przez wartości w punktach, ponieważ jako podstawy tej reprezenta-
cji można użyć dowolnego zbioru n +�1 różnych punktów x0 , x1, ..., xn .
Obliczanie omawianej reprezentacji dla wielomianu danego w reprezentacji przez współ-
czynniki jest proste, poniaważ wystarczy w tym celu wybrać n +�1 różnych punktów x0 , x1, ...,
xn , a następnie obliczyć A(xk ) dla k =� 0, 1,..., n . Korzystając z metody Hornera, można zrobić
to w czasie Q�(n2) . Jak się póz
niej przekonamy, odpowiednio dobierając xk , możemy zreduko-
wać czas obliczeń do Q�(nlog n).
10
 Zadanie odwrotne do ewaluacji  wyznaczanie współczynników wielomianu na podsta-
wie reprezentacji przez wartości w punktach  nosi nazwę interpolacji.
Twierdzenie. Dla dowolnego zbioru {(x0, y0), (x1, y1), ..., (xn, yn)} złożonego z n +�1 par punkt-
wartość istnieje dokładnie jeden wielomian A(x) stopnia n taki, że yk =� A(xk ) dla
k =� 0, 1,..., n .
Dowód. Dowód opiera się na istnieniu odwrotności pewnej macierzy. Równanie yk =� A(xk ) jest
równoważnie z równaniem macierzowym
2 n
ć�1 x0 x0 L� x0 �� a0 y0
ć� �� ć� ��
�� ��
�� �� �� ��
2 n
a1 y1
��1 x1 x1 L� x1 ��
�� �� �� ��
. (1)
��M� M� M� O� M� ���� =�
�� �� �� ��
M� M�
�� ��
��
�� �� ��
2 n ��a �� �� yn ��
��1 xn xn L� xn ��
Ł� n ł� Ł� ł�
Ł� ł�
Macierz po lewej stronie oznaczamy jako V (x0, x1,..., xn) i nazywamy macierzą Vandermonde a.
Wyznacznikiem tej macierzy jest
��(x -� x )
k j
0Ł� j<�kŁ�n
Zatem jeśli xk są różne, to jest ona odwracalna (to znaczy nieosobliwa). Współczynniki a
j
można więc wyznaczyć jednoznacznie na podstawie reprezentacji przez wartości w punktach:
a =� V (x0, x1,..., xn)-�1 y .
Dowód twierdzenia opisuje algorytm interpolacji polegający na rozwiązaniu układu (1)
równań liniowych. Korzystając z metody eleminacji Gaussa, ten układ możemy rozwiązać w
czasie O(n3) . Szybszy algorytm interpolacji w punktach opiera się na wzorze Lagrange a:
��(x -� xj )
n
A(x) =� yk ją�k . (2)
��
k =�0 ��(x -� xj )
k
ją�k
Oczywiście prawa strona równania (2) jest wielomianem stopnia n, spełniającym równanie
A(xk ) =� yk dla każdego k.
Reprezentacje wielomianów za pomocą współczynników i wartości w punktach są w
pewnym sensie równoważne; to znaczy, wielomian representowany przez wartości w punktach
ma swój wyznaczony jednoznacznie odpowiednik w reprezentacji przez współczynniki. Tak
więc ewaluacja i interpolacja w n +�1 punktach są dobrze zdefiniowanymi, wzajemnie odwrot-
nymi operacjami, realizującymi przejście między reprezentacją wielomianu przez współczynniki
a reprezentacją przez wartości w punktach. Opisane powyżej algorytmy dla tych problemów
działają w czasie Q�(n2) .
Reprezentacja przez wartości w punktach jest dość wygodna do użycia przy wielu opera-
cjach na wielomianach. W dodawaniu, jeśli C(x) =� A(x) +� B(x) , to C(xk ) =� A(xk ) +� B(xk ) w
każdym punkcie xk . Czas dodawania dwóch wielomianów stopnia n, zadanych przez wartości w
punktach, wynośi zatem Q�(n) .
Reprezentacja przez wartości w punktach jest równie dogodna do mnożenia wielomia-
nów. Jeśli C(x) =� A(x)�� B(x) , to C(xk ) =� A(xk )B(xk ) w każdym punkcie xk , więc w celu
otrzymania reprezentacji przez wartości w punktach wielomianu C możemy w każdym punkcie z
osobna przemnożyć wartość wielomianu A przez wartość wielomianu B. Ponieważ jednak sto-
pień wielomianu C jest równy 2n, do reprezentowania C potrzebujemy 2n par punkt-wartość.
11
Musimy zatem wyjść od  rozszerzonych reprezentacji przez wartości w punktach dla A i B, z
których każda będzie się składać z 2n par punkt-wartość. Widać stąd, że dla dwóch wejściowych
wielomianów w rozszerzonej reprezentacji przez wartości w punktach czas potrzebny na obli-
czanie reprezentacji przez wartości w punktach ich iloczynu wynosi Q�(n) , a więc znacznie
mniej niż w przypadku reprezentacji przez wspołczynniki.
Rozważmy na koniec problem ewaluacji w nowym punkcie wielomianu zadanego przez
wartości w punktach. Nie widać tu żadnego prostszego sposobu niż przekształcenie wielomianu
do reprezentacji przez współczynniki, a następnie obliczenie jego wartości w nowym punkcie.
6.2. Szybkie mnożenie wielomianów reprezentowanych przez współczynniki
Czy możemy skorzystać a działającej w czasie liniowym metody mnożenia wielomianów repre-
zentowanych przez wartości w punktach, żeby przyspieszyć mnożenie wielomianów w reprezen-
tacji przez współczynniki? Odpowiedz
na to pytanie zależy od tego, czy umiemy szybko wyko-
nywać przekształcenie wielomianu z reprezentacji przez współczynniki do reprezentacji przez
wartości w punktach (ewaluacja) i na odwrót (interpolacja).
Do obliczania wartości możemy użyć zupełnie dowolnych punktów, ale jeśli wybierzemy
je odpowiednio, będziemy mogli dokonywać konwersji między obiema reprezentacjami w czasie
Q�(nlog n). Strategię tę ilustruje rysunek poniżej.
Należy jeszcze poruszyć drobną kwestię związaną z ograniczeniami stopni wielomianów.
Iloczyn dwóch wielomianów stopnia n -�1 jest wielomianem stopnia 2(n -�1) . Przed zmianą re-
prezentacji wejściowych wielomianów A i B podwajamy zatem najpierw ich stopni do wartości
2(n -�1) , dodając n zerowych współczynników przy najwyższych potęgach. Ponieważ wektory
mają po 2n elementów, korzystamy z  zespolonych pierwiastków stopnia 2n z jedności , ozna-
czonych na rysunku symbolami w�2n .
Poniższa, korzystająca z FFT procedura mnoży dwa wielomiany A(x) i B(x) stopnia n
w czasie Q�(nlog n), przy czym wielomiany wejściowe i wyjściowy są reprezentowane przez
współczynniki. Zakładamy, że n jest potęga dwójki; warunku tego można zawsze dotrzymać,
dostawiając zerowe współczynniki przy najwyższych potęgach x.
1. Podwojenie stopnia wielomianów: Rozszerz reprezentacje przez współczynniki wielo-
mianów A(x) i B(x) do wartości stopnia 2n, dodając do każdej po n zerowych współczynników
przy najwyższych potęgach.
2. Ewaluacja: Oblicz reprezentacje przez wartości w punktach dla wielomianów A(x) i
B(x) , stosując dwukrotnie FFT rzędu 2n. Reprezentacje te składają się z wartości wielomianów
dla pierwiastków stopnia 2n z jedności.
12
3. Mnożenie po współrzędnych: Oblicz reprezentację przez wartości w punktach wielo-
mianu C(x) =� A(x)B(x) , wymnażając odpowiadające sobie wartości. Reprezentacja ta składa się
z wartości C(x) we wszystkich pierwiastkach stopnia 2n z jedności.
4. Interpolacja: Utwórz reprezentację przez współczynniki wielomianu C(x), stosując
jednokrotnie FFT do wektora 2n wartości w celu oblicznia odwrotnej DFT.
Kroki 1 i 3 realizuje się w czasie Q�(n) , a 2 i 4 w czasie Q�(nlog n). Jeśli więc pokażemy,
jak wykonywać FFT, udowodnimy następujące twierdzenie.
Twierdzenie. Iloczyn dwóch wielomianów stopnia n można obliczyć w czasie Q�(nlog n) , przy
czym wielomiany wejściowe i wyjściowy są reprezentowane przez współczynniki.
6.3. Zespolone pierwiastki z jedności
Zespolony pierwiastek n-tego stopnia z jedności to liczba zespolona w� taka, że w�n =� 1. Istnieje
dokladnie n zespolonych pierwiastków n-tego stopnia z jedności; są to liczby e2p�ki/ n dla
k =� 0, 1,..., n -�1. Na rysunku poniżej widać, że n zespolonych pierwiastków z jedności jest roz-
mieszczonych w równych odstępach na okręgu o promieniu jednostkowym i środku w początku
układu współrzędnych na plaszczyznie zespolonej. Wartość
w�n =� e2p� ��i / n .
Nazywamy glównym pierwiastkiem n-tego stopnia z jedności; wszystkie pozostałe zespolone
pierwiastki n-tego stopnia z jedności są potęgami w�n .
Podstawowe własności zespolonych pierwiastków n-tego stopnia z jedności opisują po-
niższe lematy.
Lemat 1. (Lemat o skracaniu) Dla dowolnych liczb całkowitych n ł� 0 , k ł� 0 i d >� 0 zachodzi
dk k
w�dn =� w�n .
Lemat 2. (Lemat o redukcji) Jeśli n >� 0 jest parzyste, to zbiór kwadratów n zespolonych pier-
wiastków n-tego stopnia z jedności to zarazem n / 2 zespolonych pierwiastków stopnia n / 2 z
jedności.
Lemat 3. (Lemat o sumowaniu) Dla dowolnej liczby całkowitej n ł� 1 i dowolnej nieujemnej
liczby całkowitej k niepodzielnej przez n zachodzi
13
n-�1
k
��(w�n ) j =� 0 .
j=�0
6.4. Dyskretna transformata Fouriera (DFT)
Przypominamy, że chcemy dokonać ewaluacji wielomianu
n-�1
j
A(x) =� ajx
��
j=�0
0 1 2 n-�1
stopnia n w punktach w�n , w�n , w�n , ..., w�n (to znaczy w n zespolonych pierwiastkach n-tego
stopnia z jedności). Bez straty ogólności możemy zalożyć, że n jest potęga 2, ponieważ dane
ograniczenie stopnia można zawsze powiększyć  zawsze możemy w miarę potrzeb dodawać
zerowe współczynniki przy najwyższych potęgach. Zakładamy, że wielomian A jest zadany
przez współczynniki: a =� (a0, a1,..., an-�1) . Zdefiniujmy wartości yk dla k =� 0, 1,..., n -�1 wzorem
n-�1
k kj
yk =� A(w�n ) =� w�n .
��aj
j =�0
Wektor y =� (y0, y1,..., yn-�1) jest dyskretną transformatą Fouriera (DFT) wektora współczyn-
ników a =� (a0, a1,..., an-�1) . Piszemy także y =� DFTn(a) .
6.5. Szybkie przekształcenie Fouriera (FFT)
Stosując metodę znaną jako szybkie przekształcenie Fouriera, korzystającą ze szczególnych wła-
sności zespolonych pierwiastków z jedności, możemy obliczyć DFTn(a) w czasie Q�(nlog n),
chociaż zwykła metoda wymaga czasu Q�(n2) .
Metoda FFT opiera się na strategii  dziel i zwyciężaj , wyodrębniając współczynniki o
parzystych i nieparzystych indeksach i definiując dwa nowe wielomiany stopnia n / 2 , oznaczane
jako Ae(x) i Ao(x) :
Ae(x) =� a0 +� a2x +� a4x2 +�...+� an-�2xn / 2-�1,
Ao(x) =� a1 +� a3x +� a5x2 +�...+� an-�1xn / 2-�1.
Wielomian Ae(x) zawiera wszystkie współczynniki o parzystych indeksach w A (binarna repre-
zentacja indeksu kończy się zerem), a wielomian Ao(x) zawiera wszystkie współczynniki o nie-
parzystych indeksach (binarna reprezentacja indeksu kończy się jedynką). Wynika stąd, że
A(x) =� Ae(x2) +� xAo(x2) ,
0 1 2 n-�1
zatem problem ewaluacji A(x) w punktach w�n , w�n , w�n , ..., w�n sprowadza się do:
2 2 2
0 1 2
i) ewaluacji wielomianów Ae(x) i Ao(x) stopnia n / 2 w punktach (�w�n)� , (�w�n)� , (�w�n)� , ...,
2
n-�1
(�w�n )� ,
a następnie
ii) połączenia wyników zgodnie ze wzorem A(x) =� Ae(x2) +� xAo(x2) .
2 2 2 2
0 1 2 n-�1
Na mocy lematu o redukcji, lista (�w�n)� , (�w�n)� , (�w�n)� , ..., (�w�n )� sklada się nie z n róż-
nych wartości, ale tylko z n / 2 zespolonych pierwiastków stopnia n / 2 z jedności, z których
każdy występuje dokładnie dwa razy. Dokonujemy zatem rekurencyjnie ewaluacji wielomianów
Ae(x) i Ao(x) stopnia n / 2 we wszystkich n / 2 zespolonych pierwiastkach stopnia n / 2 z jed-
ności. Obudwa podproblemy są dokładnie tej samej postaci co problem pierwotny, ale dwukrot-
14
nie mniejszego rozmiaru. Udało nam się podzielić obliczenie n-elementowego DFTn na dwa ob-
liczenia n / 2 -elementowego DFTn / 2 . Podział ten jest podstawą poniższego rekurencyjnego al-
gorytmu FFT, obliczającego DFT dla n-elementowego wektora a =� (a0,a1,...,an-�1) , gdzie n jest
potęga dwójki.
RECURSIVE-FFT(a)
1 n := length(a) // n jest potęga 2
2 if n = 1
3 then return a
4
w�n :=� exp(2p�i / n)
5 w� :=� 1
6
a0 :=� (a0, a2,..., an-�2)
7
a1:=� (a1, a3,..., an-�1)
8 y0 := RECURSIVE-FFT(a0)
9 y1 := RECURSIVE-FFT(a1)
10 for k := 0 to n / 2 -�1 do
11
yk :=� y0k +�w� �� y1k
12
yk +�(n / 2) :=� y0k -�w� �� y1k
13
w� :=� w� ��w�n
14 return y od // y jest wektorem kolumnowym
Procedura RECURSIVE-FFT działa następująco. Wiersze 2-3 odpowiadają największemu za-
głębieniu rekursji; wartość DFT pojedynczego elementu jest równa jemu samemu, bo wówczas
0
y =� a0w�1 =� a0 ��1 =� a0 .
W wierszach 6-7 są definiowane wektory współczynników wielomianów Ae(x) i Ao(x) . Wier-
sze 4, 5 i 13 zapewniają poprawną aktualizację wartości w� , dzięki czemu przy każdym wyko-
k
nywaniu instrukcji w wierszach 11-12 mamy w� =� w�n . (Aktualizacja wartości w� zamiast obli-
k
czania w�n od początku w każdym przebiegu pętli for pozwala zaoszczędzić na czasie). W wier-
szach 8-9 obliczamy rekurencyjnie DFTn / 2 , kładąc dla k =� 0, 1,..., n / 2 -�1
k k
y0k =� Ae(w�n / 2) , y1k =� Ao(w�n / 2)
lub
2k 2k
y0k =� Ae(w�n ) , y1k =� Ao(w�n )
k 2k
(ponieważ w�n / 2 =� w�n na mocy lematu o skracaniu).
W wierszach 11-12 są lączone wyniki rekurencyjnych obliczeń DFTn / 2 . Dla y0 , y1 , ...,
yn / 2-�1 w wierszu 11 obliczamy
k 2k k 2k k
yk =� y0k +�w�n �� y1k =� Ae(w�n ) +�w�n Ao(w�n ) =� A(w�n ) ,
gdzie ostatnia równość wynika ze wzoru A(x) =� Ae(x2) +� xAo(x2) . Dla yn / 2 , yn / 2+�1 , ..., yn-�1 ,
przyjmując k =� 0, 1,..., n / 2 -�1, w wierszu 12 obliczmy
k k 2k k 2k
yk +�(n / 2) =� y0k -�w�n �� y1k =� y0k +�w�n +�(n / 2) �� y1k =� Ae(w�n ) +�w�n +�(n / 2)Ao(w�n ) =�
15
2k k 2k k
=� Ae(w�n +�n) +�w�n +�(n / 2)Ao(w�n +�n) =� A(w�n +�(n / 2)) .
k k
Druga równość wynika z pierwszej, ponieważ w�n +�(n / 2) =� -�w�n . Czwarta równość wynika z trze-
n 2k 2k
ciej, bo z tego, że w�n =� 1, wynika, że w�n =� w�n +�n . Ostatnia równość wynika ze wzoru
A(x) =� Ae(x2) +� xAo(x2) . Wektor y obliczany przez procedurę RECURSIVE-FFT jest zatem
rzeczywiście dyskretną transformatą Fouriera (DFR) wejściowego wektora a.
W celu oszacowania czasu działania procedury RECURSIVE-FFT zauważmy, że poza
wywołaniami rekurencyjnymi wykonanie procedury zajmuje czas Q�(n) , gdzie n jest długościa
wejściowego wektora. Równanie rekurencyjne na złożoność czasową wygląda następująco:
T(n) =� 2T(n/ 2) +� Q�(n) =� Q�(nlog n) .
Korzystając z szybkiego przekształcenia Fouriera, możemy zatem dokonywać ewaluacji wielo-
mianu stopnia n w zespolonych pierwiastkach n-tego stopnia z jedności w czasie Q�(nlog n).
6.6. Interpolacja w zespolonych pierwiastkach z jedności
Dla przejścia od reprezentacji wielomianu przez wartości w punktach z powrotem do reprezenta-
cji przez współczynniki musimy obliczyć interpolację w zespolonych pierwiastkach z jedności.
Wzór na interpolację wyprowadzamy zapisując DFT jako równanie macierzowe, a następnie
analizyjąc postać macierzy odwrotnej.
Obliczenie DFT możemy przedstawić jako mnożenie macierzy y =� Vn ��a , gdzie Vn jest
macierzą Vandermonde a zawierającą odpowiednie potęgi w�n ;
y0 1 1 1 1 L� 1 a0
ć� �� ć� �� ć� ��
�� �� �� �� �� ��
2 3 n-�1
y1 a1
�� �� ��1 w�n w�n w�n L� w�n �� �� ��
2 4 6 2(n-�1)
�� �� ��1 w�n w�n w�n L� w�n �� �� ��
y2 a2
�� �� �� ������ ��
=� .
3 6 9 3
y3 1 w�n w�n w�n L� w�n(n-�1) a3
�� �� �� �� �� ��
�� �� ��M� M� M� �� �� ��
M� M� O� M� M�
�� �� �� �� �� ��
n-�1 2(n-�1) 3 (
�� �� ��1 w�n w�n w�n(n-�1) L� w�nn-�1)(n-�1) �� ��a ��
yn-�1
Ł� ł� Ł� ł� Ł� n-�1 ł�
kj
Element na pozycji (k, j) w macierzy Vn to w�n dla j, k =� 0, 1,..., n -�1, ich wykładniki tworzą
więc tabliczkę mnożenia.
Operację odwrotną, którą zapisujemy jako a =� DFTn-�1(y) , wykonujemy, mnożąć y przez
Vn-�1 (macierz odwrotną do Vn ).
-�
w�n jk
Twierdzenie. Dla j, k =� 0, 1,..., n -�1 elementem na pozycji ( j,k) w macierzy Vn-�1 jest .
n
Znając postać macierzy odwrotnej Vn-�1, wiemy, że DFTn-�1(y) zadane jest wzorem
n-�1 n-�1
1 1 2p�i
��
-�
aj =� ykw�n jk =� yk expć�-� jk
�� �� �� ��
n n n
k =�0 j=�0 Ł� ł�
dla j =� 0, 1,..., n -�1. Widzimy, że transformatę odwrotną do DFT można obliczyć za pomocą
modyfikacji algorytmu FFT, polegającej na zamianie rolami wektorów a i y, zastępieniu w�n
-�1
przez w�n i podzieleniu wartości każdej współrzędnej wyniku przez n. A więc DFTn-�1(y) można
obliczyć również w czasie Q�(nlog n) .
Korzystając z FFT i odwrotnego przekształcenia FFT, możemy zatem dla danego wielo-
mianu stopnia n przechodzić od reprezentacji przez współczynniki do reprezentacji przez warto-
16
ści w punktach i z powrotem w czasie Q�(nlog n) . W kontekście mnożenia wielomianów udo-
wodniliśmy następujące twierdzenie:
Twierdzenie o splocie. Dla dowolnej pary wektorów a i b długości n, gdzie n jest potęgą
dwójki, mamy
a ��b =� DFT2-�1(�DFT2n(a)�� DFT2n(b))�,
n
gdzie wektory a i b są uzupełnone zerami do długości 2n, a  �� oznacza mnożenie po współrzęd-
nych dwóch 2n-elementowych wektorów.
6.7. Iteracyjna wersja algorytmu FFT
Przedstawimy teraz wersję iteracyjną algorytmu FFT, korzystającą z pomocniczej procedury
BIT-REVERSE-COPY(a, A) w celu skopiowania elementów wektora a do tablicy A w odpo-
wiednim porządku.
ITERATIVE-FFT(a)
1 BIT-REVERSE-COPY(a, A)
2 n := length(a) // n jest potęga 2
3 for s := 1 to log n do
4 m :=� 2^ s
5
w�m :=� exp(2p�i / m)
6 w� :=� 1
7 for j := 0 to m/2-1 do
8 for k := j to n-1 step m do
9 t :=� w� �� A[k +� m/ 2]
10 u :=� A[k]
11 A[k]:=� u +� t
12 A[k +� m/ 2]:=� u -�t od
13
w� :=� w� ��w�m od od
14 return A
BIT-REVERSE-COPY(a, A)
1 n := length(a)
2 for k := 0 to n - 1 do
3
A[rev(k)]:=� ak
W jaki sposób procedura BIT-REVERSE-COPY(a, A) wstawia elementy wejściowego wektora
a na właściwe miejsca w tablicy A? Funkcja rev(k) : dla binarnej reprezentacji liczby k oblicza
się odwrotna kolejność bitów i znajduje odpowiedną liczbę w systemu dziesiątnym.
17