1
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 m a-
ł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 skoordyn o-
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
i
ad
bc
bd
ac
i
d
c
i
b
a
)
(
)
)(
(
wymaga czterych mnożeń liczb rzeczywistaych, w istocie można go otrzymać, wykonyjąc zale-
dwie trzy mnożenia: ac ,
bd
oraz
)
)(
(
d
c
b
a
, gdyż
bd
ac
d
c
b
a
ad
bc
)
)(
(
.
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
2
/
n
bitów:
R
L
n
x
x
x
2
/
2
,
R
L
n
y
y
y
2
/
2
.
Na przykład, jeśli
2
10110110
x
(indeks dolny 2 oznacza zapis dwójkowy), to
2
1011
L
x
,
2
0110
R
x
oraz
2
4
2
0110
2
1011
x
. Iloczyn x i y może być zatem zapisany jako
R
R
L
R
R
L
n
R
L
n
R
L
n
R
L
n
y
x
y
x
y
x
x
x
y
y
x
x
xy
)
(
2
2
)
2
)(
2
(
2
/
2
/
2
/
.
2
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
2
/
n
-bitowe mnożenia:
L
L
y
x
,
R
L
y
x
,
L
R
y
x
oraz
R
R
y
x
, 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
2
/
n
-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
)
(n
O
. Oznaczając przez
)
(n
T
całkowity czas działania algorytmu na n-bitowych danych wej-
ściowych, dostajemy zależność rekurencyjną
)
(
)
2
/
(
4
)
(
n
O
n
T
n
T
.
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
)
(
2
n
O
, 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
2
/
n
-bitowych mnożeń, tak jak już było pokazane wcześniej,
wystarczą zaledwie trzy:
L
L
y
x
,
R
R
y
x
oraz
)
)(
(
R
L
R
L
y
y
x
x
, gdyż
R
R
L
L
R
L
R
L
L
R
R
L
y
x
y
x
y
y
x
x
y
x
y
x
)
)(
(
.
Czas działania wynikającego stąd algorytmu, przedstawionego na rysunku poniżej, jest krótszy:
)
(
)
2
/
(
3
)
(
n
O
n
T
n
T
.
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
)
(
)
(
59
.
1
3
log
n
O
n
O
.
funkcja multiply(x, y)
// wejście: dwie n-bitowe liczby naturalne x i y
// wyjście: ich iloczyn
if
1
n
then return
xy
L
x ,
R
x
lewe
2
/
n
, prawe
2
/
n
bitów x
L
y ,
R
y
lewe
2
/
n
, prawe
2
/
n
bitów y
P
1
= multiply(
L
x ,
L
y )
P
2
= multiply(
R
x ,
R
y )
P
3
= multiply(
R
L
x
x
,
R
L
y
y
)
return
2
2
/
2
1
3
1
2
)
(
2
P
P
P
P
P
n
n
Czas działania algorytmu można uzyskać po analize równania rekurencyjnego
cn
n
aT
n
T
)
2
/
(
)
(
.
Przy podstawianiu
k
n
2
otrzymamy
k
k
k
k
k
k
c
c
aT
a
c
aT
T
2
)
2
)
2
(
(
2
)
2
(
)
2
(
1
2
1
k
k
k
k
k
k
k
c
ac
c
aT
a
c
ac
T
a
2
2
)
2
)
2
(
(
2
2
)
2
(
1
2
3
2
1
2
2
k
k
k
k
c
ac
c
a
T
a
2
2
2
)
2
(
1
2
2
3
3
3
k
k
k
k
k
c
ac
c
a
c
a
T
a
2
2
2
2
)
2
(
1
2
2
1
1
0
a
a
c
a
a
a
a
a
c
a
a
k
k
k
k
k
k
/
2
1
)
/
2
(
1
2
)
)
/
2
(
...
)
/
2
(
/
2
1
(
2
1
1
2
1
2
2
2
2
1
2
2
2
2
2
1
2
2
2
log
a
cn
n
a
c
a
c
a
c
a
a
a
c
a
a
k
k
k
k
k
.
Dla szkolnej techniki mnożenia mamy
4
a
i wtedy
)
(
)
1
(
)
2
/
(
4
)
(
2
2
n
O
cn
n
c
cn
n
T
n
T
.
Zatem dla
3
a
)
(
2
)
2
1
(
)
2
/
(
3
)
(
3
log
3
log
n
O
cn
n
c
cn
n
T
n
T
,
w przybliżeniu jest równy
)
(
59
.
1
n
O
.
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 znaleźć 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
b
n /
każdy, a następnie łączenie tych rozwiązań w czasie
)
(
d
n
O
, dla pewnych
0
,
,
d
b
a
(w algo-
rytmie mnożenia
3
a
,
2
b
,
1
d
). Ich czas działania może być zatem wyrażony równaniem
)
(
)
/
(
)
(
d
n
O
b
n
aT
n
T
. Wyprowadzimy teraz zwarte rozwiązanie tej rekurencji, aby w przy-
szłości nie musieć rozwiązywać jej dla nowych danych.
Twerdzenie o rekurencji universalnej. Jeśli
)
(
)
/
(
)
(
d
n
O
b
n
aT
n
T
dla pewnych stałych
0
a
,
1
b
oraz
0
d
, to
a
d
gdy
n
O
a
d
gdy
n
n
O
a
d
gdy
n
O
n
T
b
a
b
d
b
d
b
log
)
(
log
)
log
(
log
),
(
)
(
log
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:
k
b
n
.
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 /
.
kd
d
k
k
kd
k
k
b
c
b
c
b
aT
a
b
c
b
aT
b
T
)
)
(
(
)
(
)
(
)
1
(
2
1
kd
d
k
d
k
k
kd
d
k
k
b
c
b
ac
b
c
b
aT
a
b
c
b
ac
b
T
a
)
1
(
)
2
(
3
2
)
1
(
2
2
)
)
(
(
)
(
kd
d
k
d
k
k
b
c
b
ac
b
c
a
b
T
a
)
1
(
)
2
(
2
3
3
)
(
4
kd
d
k
d
k
d
k
k
b
c
b
ac
b
c
a
b
c
a
b
T
a
)
1
(
2
2
1
0
)
(
1
)
1
(
2
)
2
(
2
2
1
1
)
1
(
k
d
k
k
d
k
d
d
d
k
k
a
b
a
b
a
b
a
b
b
c
a
T
a
.
Gdy
a
d
b
log
, otrzymujemy
a
b
b
a
d
b
log
, zatem
)
)
1
(
(
)
1
(
)
(
log
k
c
T
b
k
c
a
T
a
b
T
a
k
k
k
k
b
)
log
(
)
log
)
1
(
(
)
(
log
n
n
O
n
c
T
n
n
T
d
b
a
b
.
Gdy
a
d
b
log
i k zmierza się od 0 do
n
b
log
, wartoście w nawiasie tworzą ciąg geometryczny
o ilorazie
a
b
d
/ . Znaczenie sumy k wyrazów tego ciągu może być zapisane jako
a
b
a
b
d
k
kd
/
1
/
1
.
Zatem otrzymujemy
kd
d
d
d
d
k
d
k
kd
d
k
k
k
b
b
a
b
c
b
a
b
c
T
a
a
b
a
b
b
c
a
T
a
b
T
)
1
(
/
1
/
1
)
1
(
)
(
1
.
a
d
gdy
n
O
a
d
gdy
n
O
n
b
a
b
c
b
a
b
c
T
n
n
T
b
a
b
d
d
d
d
d
d
a
b
b
log
),
(
log
),
(
)
1
(
)
(
log
log
3. Mnożenie macierzy
Iloczynem dwóch macierzy
n
n
, X oraz Y, jest trzecia macierz
n
n
,
XY
Z
, której element o
współrzędnych
)
,
( j
i
jest dany wzorem
n
k
kj
ik
ij
Y
X
Z
1
.
Powyższy wzór implikuje algorytm mnożenia macierzy działający w czasie
)
(
3
n
O
; należy obli-
czyć wrtości
2
n elementów, a każde obliczenie zajmuje czas
)
(n
O
.
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
2
/
2
/
n
n
, to samo zróbmy dla Y:
D
C
B
A
X
,
H
G
F
E
Y
.
Wtedy ich iloczyn może być wyrażony za pomocą tych bloków, dokładnie tak, jakby te bloki
bały pojedynczymi elementami:
DH
CF
DG
CE
BH
AF
BG
AE
H
G
F
E
D
C
B
A
XY
.
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
2
/
n
każda, a następnie wykonujemy kilka dodawań w czasie
)
(
2
n
O
. Całkowity
czas działania jest opisany zależnością rekurencyjną:
5
)
(
)
2
/
(
8
)
(
2
n
O
n
T
n
T
.
To prowadzi do mało imponującego rozwiązania
)
(
3
n
O
, 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
2
/
2
/
n
n
:
7
3
5
1
4
3
2
1
6
2
4
5
P
P
P
P
P
P
P
P
P
P
P
P
XY
,
gdzie
)
(
1
H
F
A
P
,
)
)(
(
5
H
E
D
A
P
,
H
B
A
P
)
(
2
,
)
)(
(
6
H
G
D
B
P
,
E
D
C
P
)
(
3
,
)
)(
(
7
F
E
C
A
P
,
)
(
4
E
G
D
P
.
Nowy czas działania wynosi
)
(
)
2
/
(
7
)
(
2
n
O
n
T
n
T
,
który na podstawie twierdzenia o rekurencji uniwersalnej okazuje się być równy
)
(
)
(
81
.
2
7
log
n
O
n
O
.
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óźniejszego wyszukiwania elementów
sortowanego zbioru.
Niech U będzie zbiorem obiektów
n
a
a
a
,...,
,
2
1
Sortowanie polega na permutowaniu tych obiektów aż do chwili osiągnięcia uporządkowania
n
k
k
k
a
a
a
,...,
,
2
1
takiego, że dla zadanej funkcji porządkującej f zachodzi
)
(
...
)
(
)
(
2
1
n
k
k
k
a
f
a
f
a
f
.
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
n
a
a
a
,...,
,
2
1
znajdują się w tablicy
]
..
1
[ n
a
. 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
)
,
( j
i
Swap
pole-
gająca na zamianie elementów w tablicy a z pozycji i oraz j,
n
j
i
,
1
.
6
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
7
4.1. Poprawność algorytmu
Zauważmy, że elementy dodatkowe, dopisane do ciągów a i b, nie mogą się znaleźć w ciągu
wynikowym, gdyż pętla jest wykonywana tylko (n+m) razy. Gdyby, w którymś momencie
1
n
i
, tzn. a[i] =
, to przy porównywaniu tego elementu z b[j] dla
1
m
j
, wygra b[j] i
to b[j] zostanie zapisane w ciągu wynikowym. Analogicznie w przypadku, gdy
1
m
j
. Jeśli
zaś oba wskaź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
1
m
n
.
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]
8
// wyjście: tablica wejściwa po posortowaniu
if
1
n
then
return merge(mergesort(
2
/
..
1
[
n
a
),mergesort(
n
n
a
..
1
2
/
[
))
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
)
(n
O
. Zatem procedura scalania jest liniowa, a całkowity czas
potrzebny na wykonanie mergesort to
)
log
(
)
(
)
2
/
(
2
)
(
n
n
O
n
O
n
T
n
T
.
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ą
]
25
,
30
,
10
,
1
,
45
[
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 –
)
log
(
n
n
O
, 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 znaleźć ś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
1
k
, szukane jest minimum S, natomiast jeśli
2
/
|
| S
k
, 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 wyobraź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
L
S ,
v
S ,
R
S . Na przykład, jeśli
tablica
S :
2
36
5
21 8 13
11
20
5
4
1
jest podzielona dla
5
v
, trzy wygenerowane podtablice wyglądają następująco:
9
S
L
:
2
4
1
S
v
:
5
5 S
R
:
36
21
8
13
11
20
Przeszukiwanie może być natychmiast zawężone do jednej z tych podlist. Jeśli chcemy
znaleźć, powiedzmy, ósmy najmniejszy element S, wiemy, że będzie to trzeci najmniejszy ele-
ment
R
S , ponieważ
5
|
|
|
|
v
L
S
S
. Czyli selection(S,8) = selection(
R
S ,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:
|
|
|
|
|)
|
|
|
,
(
|
|
|
|
|
|
|
|
)
,
(
)
,
(
v
L
v
L
R
v
L
L
L
L
S
S
k
gdy
S
S
k
S
selection
S
S
k
S
gdy
v
S
k
gdy
k
S
selection
k
S
selection
Trzy podlisty
L
S ,
v
S ,
R
S 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{|
R
L
S
S
.
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ć
taka, by tablica została istotnie zmniejszona, sytuacja jest idealna dla
|
|
2
1
|
|
|,
|
S
S
S
R
L
. Gdyby-
śmy potrafili zawsze zagwarantować taką sytuację, dostalibyśmy czas działania
)
(
)
2
/
(
)
(
n
O
n
T
n
T
,
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ć
36
v
, następnie
21
v
itd. Ten czarny sce-
nariusz zmusiłby nasz algorytm selekcji do wykonania
)
(
2
...
)
2
(
)
1
(
2
n
n
n
n
n
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
)
(n
O
. Gdzie, w przedziale od
)
(n
O
do
)
(
2
n
, 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
L
S oraz
R
S 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
)
(n
T
oczekiwany czas działania dla
tablicy o rozmiarze n, dostajemy
)
(
)
4
/
3
(
)
(
n
O
n
T
n
T
.
10
Na podstawie tej rekurencji wnioskujemy, że
)
(
)
(
n
O
n
T
: dla dowolnych danych wejściowych
nasz algorytm zwraca poprawną odpowiedź 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
4
3
2
2
2
12
11
12
5
2
)
4
2
(
)
3
2
1
(
x
x
x
x
x
x
x
x
.
Ogólniej, jeśli
n
n
x
a
x
a
a
x
A
...
)
(
1
0
oraz
n
n
x
b
x
b
b
x
B
...
)
(
1
0
, to ich iloczyn
n
n
x
c
x
c
c
x
B
x
A
x
C
2
2
1
0
...
)
(
)
(
)
(
ma współczynniki
k
i
i
k
i
k
k
k
k
b
a
b
a
b
a
b
a
c
0
0
1
1
0
...
(dla
n
i
bierzemy
i
a oraz
i
b równe zero). Obliczenie
k
c na podstawie tego wzoru wymaga
)
(k
O
kroków, a znalezenie wszystkich
1
2
n
współczynników wydaje się wymagać czasu
)
(
2
n
. 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
Reprezentacja przez współczynniki wielomianu
n
j
j
j
x
a
x
A
0
)
(
stopnia n to wektor współczyn-
ników
)
,...,
,
(
1
0
n
a
a
a
a
. Reprezentacja za pomocą współczynników jest dogodna przy niektó-
rych operacjach na wielomianach. Na przykład operacja ewaluacji wielomianu
)
(x
A
w danym
punkcie
0
x polega na obliczeniu wartości
)
(
0
x
A
. Ewaluacje można wykonać w czasie
)
(n
,
korzystając t tzw. schematu Hornera:
))...))
(
(
...
(
(
)
(
0
1
0
2
0
1
0
0
0
n
n
a
x
a
x
a
x
a
x
a
x
A
.
Podobnie, dodawanie dwóch wielomianów reprezentowanych przez wektory współczynników
)
,...,
,
(
1
0
n
a
a
a
a
i
)
,...,
,
(
1
0
n
b
b
b
b
zajmuje czas
)
(n
: wynik stanowi wektor
)
,...,
,
(
1
0
n
c
c
c
c
,
gdzie
j
j
j
b
a
c
dla
n
j
,...,
1
,
0
.
Reprezentacja przez wartości w punktach wielomianu
)
(x
A
stopnia n to zbiór par
punkt-wartość
)}
,
(
...,
),
,
(
),
,
{(
1
1
0
0
n
n
y
x
y
x
y
x
taki, że wszystkie
k
x są parami różne oraz
)
(
k
k
x
A
y
dla
n
k
,...,
1
,
0
. 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
1
n
różnych punktów
0
x ,
1
x , ...,
n
x .
Obliczanie omawianej reprezentacji dla wielomianu danego w reprezentacji przez współ-
czynniki jest proste, poniaważ wystarczy w tym celu wybrać
1
n
różnych punktów
0
x ,
1
x , ...,
n
x
, a następnie obliczyć
)
(
k
x
A
dla
n
k
,...,
1
,
0
. Korzystając z metody Hornera, można zrobić
to w czasie
)
(
2
n
. Jak się później przekonamy, odpowiednio dobierając
k
x
, możemy zreduko-
wać czas obliczeń do
)
log
(
n
n
.
11
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
)}
,
(
...,
),
,
(
),
,
{(
1
1
0
0
n
n
y
x
y
x
y
x
złożonego z
1
n
par punkt-
wartość istnieje dokładnie jeden wielomian
)
(x
A
stopnia n taki, że
)
(
k
k
x
A
y
dla
n
k
,...,
1
,
0
.
Dowód. Dowód opiera się na istnieniu odwrotności pewnej macierzy. Równanie
)
(
k
k
x
A
y
jest
równoważnie z równaniem macierzowym
n
n
n
n
n
n
n
n
y
y
y
a
a
a
x
x
x
x
x
x
x
x
x
1
0
1
0
2
1
2
1
1
0
2
0
0
1
1
1
.
(1)
Macierz po lewej stronie oznaczamy jako
)
,...,
,
(
1
0
n
x
x
x
V
i nazywamy macierzą Vandermonde’a.
Wyznacznikiem tej macierzy jest
n
k
j
j
k
x
x
0
)
(
Zatem jeśli
k
x są różne, to jest ona odwracalna (to znaczy nieosobliwa). Współczynniki
j
a
można więc wyznaczyć jednoznacznie na podstawie reprezentacji przez wartości w punktach:
y
x
x
x
V
a
n
1
1
0
)
,...,
,
(
.
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
)
(
3
n
O
. Szybszy algorytm interpolacji w punktach opiera się na wzorze Lagrange’a:
n
k
k
j
j
k
k
j
j
k
x
x
x
x
y
x
A
0
)
(
)
(
)
(
.
(2)
Oczywiście prawa strona równania (2) jest wielomianem stopnia n, spełniającym równanie
k
k
y
x
A
)
(
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
1
n
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
)
(
2
n
.
Reprezentacja przez wartości w punktach jest dość wygodna do użycia przy wielu opera-
cjach na wielomianach. W dodawaniu, jeśli
)
(
)
(
)
(
x
B
x
A
x
C
, to
)
(
)
(
)
(
k
k
k
x
B
x
A
x
C
w
każdym punkcie
k
x
. Czas dodawania dwóch wielomianów stopnia n, zadanych przez wartości w
punktach, wynośi zatem
)
(n
.
Reprezentacja przez wartości w punktach jest równie dogodna do mnożenia wielomia-
nów. Jeśli
)
(
)
(
)
(
x
B
x
A
x
C
, to
)
(
)
(
)
(
k
k
k
x
B
x
A
x
C
w każdym punkcie
k
x
, 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ść.
12
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
)
(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? Odpowiedź 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
)
log
(
n
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
1
n
jest wielomianem stopnia
)
1
(
2
n
. Przed zmianą re-
prezentacji wejściowych wielomianów A i B podwajamy zatem najpierw ich stopni do wartości
)
1
(
2
n
, 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
n
2
.
Poniższa, korzystająca z FFT procedura mnoży dwa wielomiany
)
(x
A
i
)
(x
B
stopnia n
w czasie
)
log
(
n
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
)
(x
A
i
)
(x
B
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
)
(x
A
i
)
(x
B
, 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.
13
3. Mnożenie po współrzędnych: Oblicz reprezentację przez wartości w punktach wielo-
mianu
)
(
)
(
)
(
x
B
x
A
x
C
, 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
)
(n
, a 2 i 4 w czasie
)
log
(
n
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
)
log
(
n
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
taka, że
1
n
. Istnieje
dokladnie n zespolonych pierwiastków n-tego stopnia z jedności; są to liczby
n
ki
e
/
2
dla
1
,...,
1
,
0
n
k
. 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ść
n
i
n
e
/
2
.
Nazywamy glównym pierwiastkiem n-tego stopnia z jedności; wszystkie pozostałe zespolone
pierwiastki n-tego stopnia z jedności są potęgami
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
0
n
,
0
k
i
0
d
zachodzi
k
n
dk
dn
.
Lemat 2. (Lemat o redukcji) Jeśli
0
n
jest parzyste, to zbiór kwadratów n zespolonych pier-
wiastków n-tego stopnia z jedności to zarazem
2
/
n
zespolonych pierwiastków stopnia
2
/
n
z
jedności.
Lemat 3. (Lemat o sumowaniu) Dla dowolnej liczby całkowitej
1
n
i dowolnej nieujemnej
liczby całkowitej k niepodzielnej przez n zachodzi
14
0
)
(
1
0
n
j
j
k
n
.
6.4. Dyskretna transformata Fouriera (DFT)
Przypominamy, że chcemy dokonać ewaluacji wielomianu
1
0
)
(
n
j
j
j
x
a
x
A
stopnia n w punktach
0
n
,
1
n
,
2
n
, ...,
1
n
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:
)
,...,
,
(
1
1
0
n
a
a
a
a
. Zdefiniujmy wartości
k
y dla
1
,...,
1
,
0
n
k
wzorem
1
0
)
(
n
j
kj
n
j
k
n
k
a
A
y
.
Wektor
)
,...,
,
(
1
1
0
n
y
y
y
y
jest dyskretną transformatą Fouriera (DFT) wektora współczyn-
ników
)
,...,
,
(
1
1
0
n
a
a
a
a
. Piszemy także
)
(a
DFT
y
n
.
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ć
)
(a
DFT
n
w czasie
)
log
(
n
n
,
chociaż zwykła metoda wymaga czasu
)
(
2
n
.
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
2
/
n
, oznaczane
jako
)
(x
A
e
i
)
(x
A
o
:
1
2
/
2
2
4
2
0
...
)
(
n
n
e
x
a
x
a
x
a
a
x
A
,
1
2
/
1
2
5
3
1
...
)
(
n
n
o
x
a
x
a
x
a
a
x
A
.
Wielomian
)
(x
A
e
zawiera wszystkie współczynniki o parzystych indeksach w A (binarna repre-
zentacja indeksu kończy się zerem), a wielomian
)
(x
A
o
zawiera wszystkie współczynniki o nie-
parzystych indeksach (binarna reprezentacja indeksu kończy się jedynką). Wynika stąd, że
)
(
)
(
)
(
2
2
x
xA
x
A
x
A
o
e
,
zatem problem ewaluacji A(x) w punktach
0
n
,
1
n
,
2
n
, ...,
1
n
n
sprowadza się do:
i) ewaluacji wielomianów
)
(x
A
e
i
)
(x
A
o
stopnia
2
/
n
w punktach
2
0
n
,
2
1
n
,
2
2
n
, ...,
2
1
n
n
,
a następnie
ii) połączenia wyników zgodnie ze wzorem
)
(
)
(
)
(
2
2
x
xA
x
A
x
A
o
e
.
Na mocy lematu o redukcji, lista
2
0
n
,
2
1
n
,
2
2
n
, ...,
2
1
n
n
sklada się nie z n róż-
nych wartości, ale tylko z
2
/
n
zespolonych pierwiastków stopnia
2
/
n
z jedności, z których
każdy występuje dokładnie dwa razy. Dokonujemy zatem rekurencyjnie ewaluacji wielomianów
)
(x
A
e
i
)
(x
A
o
stopnia
2
/
n
we wszystkich
2
/
n
zespolonych pierwiastkach stopnia
2
/
n
z jed-
ności. Obudwa podproblemy są dokładnie tej samej postaci co problem pierwotny, ale dwukrot-
15
nie mniejszego rozmiaru. Udało nam się podzielić obliczenie n-elementowego
n
DFT na dwa ob-
liczenia
2
/
n
-elementowego
2
/
n
DFT
. Podział ten jest podstawą poniższego rekurencyjnego al-
gorytmu FFT, obliczającego DFT dla n-elementowego wektora
)
,...,
,
(
1
1
0
n
a
a
a
a
, 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
)
/
2
exp(
:
n
i
n
5
1
:
6
)
,...,
,
(
:
0
2
2
0
n
a
a
a
a
7
)
,...,
,
(
:
1
1
3
1
n
a
a
a
a
8
y0 := RECURSIVE-FFT(a0)
9
y1 := RECURSIVE-FFT(a1)
10
for k := 0 to
1
2
/
n
do
11
k
k
k
y
y
y
1
0
:
12
k
k
n
k
y
y
y
1
0
:
)
2
/
(
13
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
0
0
1
0
1
a
a
a
y
.
W wierszach 6-7 są definiowane wektory współczynników wielomianów
)
(x
A
e
i
)
(x
A
o
. Wier-
sze 4, 5 i 13 zapewniają poprawną aktualizację wartości
, dzięki czemu przy każdym wyko-
nywaniu instrukcji w wierszach 11-12 mamy
k
n
. (Aktualizacja wartości
zamiast obli-
czania
k
n
od początku w każdym przebiegu pętli for pozwala zaoszczędzić na czasie). W wier-
szach 8-9 obliczamy rekurencyjnie
2
/
n
DFT
, kładąc dla
1
2
/
,...,
1
,
0
n
k
)
(
0
2
/
k
n
e
k
A
y
,
)
(
1
2
/
k
n
o
k
A
y
lub
)
(
0
2k
n
e
k
A
y
,
)
(
1
2k
n
o
k
A
y
(ponieważ
k
n
k
n
2
2
/
na mocy lematu o skracaniu).
W wierszach 11-12 są lączone wyniki rekurencyjnych obliczeń
2
/
n
DFT
. Dla
0
y ,
1
y , ...,
1
2
/
n
y
w wierszu 11 obliczamy
)
(
)
(
)
(
1
0
2
2
k
n
k
n
o
k
n
k
n
e
k
k
n
k
k
A
A
A
y
y
y
,
gdzie ostatnia równość wynika ze wzoru
)
(
)
(
)
(
2
2
x
xA
x
A
x
A
o
e
. Dla
2
/
n
y
,
1
2
/
n
y
, ...,
1
n
y
,
przyjmując
1
2
/
,...,
1
,
0
n
k
, w wierszu 12 obliczmy
)
(
)
(
1
0
1
0
2
)
2
/
(
2
)
2
/
(
)
2
/
(
k
n
o
n
k
n
k
n
e
k
n
k
n
k
k
k
n
k
n
k
A
A
y
y
y
y
y
16
)
(
)
(
)
(
)
2
/
(
2
)
2
/
(
2
n
k
n
n
k
n
o
n
k
n
n
k
n
e
A
A
A
.
Druga równość wynika z pierwszej, ponieważ
k
n
n
k
n
)
2
/
(
. Czwarta równość wynika z trze-
ciej, bo z tego, że
1
n
n
, wynika, że
n
k
n
k
n
2
2
. Ostatnia równość wynika ze wzoru
)
(
)
(
)
(
2
2
x
xA
x
A
x
A
o
e
. 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
)
(n
, gdzie n jest długościa
wejściowego wektora. Równanie rekurencyjne na złożoność czasową wygląda następująco:
)
log
(
)
(
)
2
/
(
2
)
(
n
n
n
n
T
n
T
.
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
)
log
(
n
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
a
V
y
n
, gdzie
n
V jest
macierzą Vandermonde’a zawierającą odpowiednie potęgi
n
;
1
3
2
1
0
)
1
)(
1
(
)
1
(
3
)
1
(
2
1
)
1
(
3
9
6
3
)
1
(
2
6
4
2
1
3
2
1
3
2
1
0
1
1
1
1
1
1
1
1
1
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
a
a
a
a
a
y
y
y
y
y
.
Element na pozycji
)
,
(
j
k
w macierzy
n
V to
kj
n
dla
1
,...,
1
,
0
,
n
k
j
, ich wykładniki tworzą
więc tabliczkę mnożenia.
Operację odwrotną, którą zapisujemy jako
)
(
1
y
DFT
a
n
, wykonujemy, mnożąć y przez
1
n
V
(macierz odwrotną do
n
V ).
Twierdzenie. Dla
1
,...,
1
,
0
,
n
k
j
elementem na pozycji
)
,
( k
j
w macierzy
1
n
V
jest
n
jk
n
.
Znając postać macierzy odwrotnej
1
n
V
, wiemy, że
)
(
1
y
DFT
n
zadane jest wzorem
1
0
1
0
2
exp
1
1
n
j
k
n
k
jk
n
k
j
jk
n
i
y
n
y
n
a
dla
1
,...,
1
,
0
n
j
. 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
n
przez
1
n
i podzieleniu wartości każdej współrzędnej wyniku przez n. A więc
)
(
1
y
DFT
n
można
obliczyć również w czasie
)
log
(
n
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-
17
ści w punktach i z powrotem w czasie
)
log
(
n
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
)
(
)
(
2
2
1
2
b
DFT
a
DFT
DFT
b
a
n
n
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
s
m
^
2
:
5
)
/
2
exp(
:
m
i
m
6
1
:
7
for j := 0 to m/2-1 do
8
for k := j to n-1 step m do
9
]
2
/
[
:
m
k
A
t
10
]
[
:
k
A
u
11
t
u
k
A
:
]
[
12
t
u
m
k
A
:
]
2
/
[
od
13
m
:
od od
14
return A
BIT-REVERSE-COPY(a, A)
1
n := length(a)
2
for k := 0 to n - 1 do
3
k
a
k
rev
A
:
)]
(
[
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
)
(k
rev
: dla binarnej reprezentacji liczby k oblicza
się odwrotna kolejność bitów i znajduje odpowiedną liczbę w systemu dziesiątnym.