1 Date:    Tue, 04 Jul 1995 13:10:58 -0400
2 From:    Chris_Okasaki@LOCH.MESS.CS.CMU.EDU
3 To:      simonpj@dcs.gla.ac.uk
4 Subject: Fibonacci Heaps
6 As I promised at the Haskell Workshop, here is a sample program
7 using encapsulated state.  I've translated this from SML, but
8 in doing so, I noticed that in fact accumArray is all the
9 encapsulated state you really need for this application.  In SML,
10 we are forced to use mutable arrays because we don't have such
11 fancy monolithic array "primitives" as accumArray.
13 I've written and tested this as a literate Gofer script because I've
14 never been able to get GHC to run under Mach. :-(
16 Let me know if you have any problems...
18 Chris
20 - ----------------------------------------------------------------------
22 FIBONACCI HEAPS
24 Fibonacci heaps are a priority queue data structure supporting the
25 following operations:
26         O(1) Insert
27         O(1) FindMin
28         O(1) Meld
29     O(log n) DeleteMin
31 (In an imperative settting, Fibonacci heaps also support
32         O(1) DecreaseKey (of an indicated element)
33     O(log n) Delete (an indicated element)
34 but these operations are problematic in a functional setting.)
36 There is one catch: for the DeleteMin operation, the bounds are
37 amortized instead of worst-case.  This means that the bounds are
38 only guaranteed if you use the data structure in a single-threaded manner.
39 Otherwise, you can take longer than expected by repeatedly going back
40 and operating on an "expensive" version of the data structure.
42 (Note: I am currently working on a paper with another student describing
43 a functional priority queue achieving the above bounds in the worst-case
44 instead of amortized.  This data structure may be freely used in a
45 non-single-threaded manner with no ill effects.)
47 To understand the implementation of Fibonacci heaps, it is helpful to
48 first understand binomial queues.  See, for example, David King's
49 "Functional Binomial Queues" from the last Glasgow workshop.
51 > -- partain
52 >module Main (main) where
53 >import ST
54 >import Array
55 >import System
57                          --------------------
59 Like binomial queues, Fibonacci heaps are based on heap-ordered
60 binomial trees.
62 >data Tree a = Node a [Tree a]
64 The degree of a binomial tree is equal to its number of children.
65 Every binomial tree of degree k has binomial trees of degrees
66 k-1...0 as children, in that order.  It is easy to show that
67 a binomial tree of degree k has size 2^k.
69 The fundamental operation on binomial trees is linking, which compares
70 the roots of two binomial trees and makes the larger a child of the
71 smaller (thus bumping its degree by one).  It is essential that this
72 only be called on binomial trees of equal degree.
74 >link (a @ (Node x as)) (b @ (Node y bs)) =
75 >  if x <= y then Node x (b:as) else Node y (a:bs)
77 It will also be useful to extract the minimum element from a tree.
79 >root (Node x _) = x
81 We will frequently need to tag trees with their degrees.
83 >type TaggedTree a = (Int,Tree a)
84 >
85 >degree (k, t) = k
86 >tree (k, t) = t
88 Given a tagged tree, extract and tag its children.
90 >getChildren (n, Node x ts) = zipWith (,) [n-1,n-2 .. ] ts
92 Extract the minimum element from a tagged tree.
94 >root' = root . tree
96                          --------------------
98 We also need a type for bags supporting constant time union.  The simple
99 representation given here is sufficient since we will always process bags
100 as a whole.  Note that for this application it is not necessary to
101 filter out occurences of EmptyBag.  Also, for this application order
102 is irrelevant.
104 >data Bag a = EmptyBag | ConsBag a (Bag a) | UnionBags (Bag a) (Bag a)
105 >
106 >bagToList b = flatten b []
107 >  where flatten EmptyBag xs = xs
108 >        flatten (ConsBag x b) xs = flatten b (x:xs)
109 >        flatten (UnionBags b1 b2) xs = flatten b1 (flatten b2 xs)
110 >
111 >applyToAll :: (a -> ST s ()) -> Bag a -> ST s ()
112 >applyToAll f EmptyBag = returnST ()
113 >applyToAll f (ConsBag x b) = f x `seqST` applyToAll f b
114 >applyToAll f (UnionBags b1 b2) = applyToAll f b1 `seqST` applyToAll f b2
116                          --------------------
118 Miscellaneous stuff.
120 >log2 1 = 0
121 >log2 n = 1 + log2 (n `div` 2)
123 >data MyMaybe a = Zero | One a
125                          --------------------
127 Since binomial trees only come in certain, fixed sizes, we need some
128 way to represent priority queues of other sizes.  We will do this
129 with a forest of trees summing to the correct size.
131 >type Forest a = Bag (TaggedTree a)
133 In binomial queues, this forest must be maintained in strictly increasing
134 order of degree.  For Fibonacci heaps, we adopt a more relaxed attitude:
135 degrees may be repeated and order does not matter.
137 To be able to find the minimum element quickly, we keep the tree with the
138 minimum root outside of the bag.  In addition, at the top level of each heap,
139 we store the total size of the heap.
141 >data FibHeap a = EmptyFH | FH Int (TaggedTree a) (Forest a)
143                          --------------------
145 Now, the following operations are trivial.
147 >emptyFH = EmptyFH
148 >
149 >isEmptyFH EmptyFH = True
150 >isEmptyFH (FH _ _ _) = False
151 >
152 >singleFH x = FH 1 (0, Node x []) EmptyBag
153 >
154 >insertFH x xs = meldFH (singleFH x) xs
155 >
156 >minFH EmptyFH = error "minFH EmptyFH"
157 >minFH (FH n tt f) = root' tt
159                          --------------------
161 Meld achieves its efficiency by simply unioning the two forests.
163 >meldFH EmptyFH xs = xs
164 >meldFH xs EmptyFH = xs
165 >meldFH (FH n1 tt1 f1) (FH n2 tt2 f2) =
166 >  if root' tt1 <= root' tt2 then
167 >      FH (n1+n2) tt1 (ConsBag tt2 (UnionBags f1 f2))
168 >  else
169 >      FH (n1+n2) tt2 (ConsBag tt1 (UnionBags f1 f2))
171                          --------------------
173 Finally, the only hard operation is deleteMin.  After throwing away the
174 minimum element, it repeatedly links trees of equal degree until
175 no such pairs are left.  The most efficient way to do this is with
176 an array.  I give two implementations, one using monadic arrays,
177 the other using accumArray.
179 In the first implementation, there are three steps.
180   1. Allocate an array indexed by degrees.
181   2. Insert every tree into this array.  If, when inserting a tree of
182      degree k, there already exists a tree of degree k, link the
183      two trees and reinsert the new larger tree.
184   3. Transfer the trees into a bag, keeping track of the minimum tree.
186 >deleteMinFH EmptyFH = error "deleteMinFH EmptyFH"
187 >deleteMinFH (FH 1 tt f) = EmptyFH
188 >deleteMinFH (FH n tt f) =
189 >  let
190 >    d = log2 (n-1) -- maximum possible degree
191 >
192 >    ins a (i, t) =
193 >        readArray a i `thenST` \e ->
194 >        case e of
195 >          Zero -> writeArray a i (One t)
196 >          One t2 -> writeArray a i Zero `seqST`
197 >                    ins a (i+1, link t t2)
199 Note that after inserting all the trees, the array contains trees
200 in the same pattern as the bits of n-1.  Since we know that the
201 highest order bit of n-1 is one, we know that there is a tree in
202 the highest slot of the array.
204 >    getMin a =
205 >        readArray a d `thenST` \e ->
206 >        case e of
207 >          Zero -> error "must be One" -- since array is filled as bits of n-1
208 >          One t -> getMin' a d t EmptyBag 0
209 >    getMin' a mini mint b i =
210 >        if i >= d then
211 >          returnST ((mini, mint),b)
212 >        else
213 >          readArray a i `thenST` \e ->
214 >          case e of
215 >            Zero -> getMin' a mini mint b (i+1)
216 >            One t -> if root mint <= root t then
217 >                       getMin' a mini mint (ConsBag (i, t) b) (i+1)
218 >                     else
219 >                       getMin' a i t (ConsBag (mini, mint) b) (i+1)
220 >
221 >  in
222 >    runST (newArray (0,d) Zero `thenST` \a ->
223 >           applyToAll (ins a) f `seqST`
224 >           listST (map (ins a) (getChildren tt)) `seqST`
225 >           getMin a `thenST` \ (tt,f) ->
226 >           returnST (FH (n-1) tt f))
228                          --------------------
230 The second version of deleteMin uses accumArray to group trees of like
231 size.  It then performs the linking and all remaining steps purely
232 functionally.
234 >deleteMinFH' EmptyFH = error "deleteMinFH EmptyFH"
235 >deleteMinFH' (FH 1 tt f) = EmptyFH
236 >deleteMinFH' (FH n tt f) =
237 >  let
238 >    d = log2 (n-1) -- maximum possible degree
239 >
240 >    a = accumArray (flip (:)) [] (0,d) (getChildren tt ++ bagToList f)
241 >
242 >    doLinks (ts:rest) = startup 0 ts rest
243 >      where startup i [] [] = []
244 >            startup i [] (ts:rest) = startup (i+1) ts rest
245 >            startup i ts [] = combine i ts [] []
246 >            startup i ts (next:rest) = combine i ts next rest
247 >
248 >            combine i [] next rest = startup (i+1) next rest
249 >            combine i [t] next rest = (i, t) : startup (i+1) next rest
250 >            combine i (t1:t2:ts) next rest =
251 >                combine i ts (link t1 t2 : next) rest
252 >
253 >    getMin (tt:rest) = foldl chooseMin (tt,EmptyBag) rest
254 >      where chooseMin (tt1,b) tt2 =
255 >                if root' tt1 <= root' tt2 then
256 >                    (tt1,ConsBag tt2 b)
257 >                else
258 >                    (tt2,ConsBag tt1 b)
259 >
260 >    (new_tt,new_f) = getMin (doLinks (elems a))
261 >  in
262 >    FH (n-1) new_tt new_f
264                          --------------------
266 Testing...
268 >fibToList :: Ord a => FibHeap a -> [a]
269 >fibToList xs = if isEmptyFH xs then []
270 >               else minFH xs : fibToList (deleteMinFH xs)
271 >
272 >fibToList' :: Ord a => FibHeap a -> [a]
273 >fibToList' xs = if isEmptyFH xs then []
274 >                else minFH xs : fibToList' (deleteMinFH' xs)
275 >
276 >makeFH :: Ord a => [a] -> FibHeap a
277 >makeFH xs = foldr insertFH emptyFH xs
278 >
279 >fibSort :: Ord a => [a] -> [a]
280 >fibSort = fibToList . makeFH
281 >
282 >fibSort' :: Ord a => [a] -> [a]
283 >fibSort' = fibToList' . makeFH
284 >
285 >randoms :: Int -> [Int]
286 >randoms n = take n (iterate (\seed-> (77*seed+1) `rem` 1024) 1967)
287 >
288 >test n = fibSort (randoms n) == fibSort' (randoms n)
290 >--partain
291 >main = getArgs >>= \ [n] ->
292 >       putStr (show (test (read n)))