Date:    Tue, 04 Jul 1995 13:10:58 -0400
From:    Chris_Okasaki@LOCH.MESS.CS.CMU.EDU
To:      simonpj@dcs.gla.ac.uk
Subject: Fibonacci Heaps

As I promised at the Haskell Workshop, here is a sample program
using encapsulated state.  I've translated this from SML, but 
in doing so, I noticed that in fact accumArray is all the
encapsulated state you really need for this application.  In SML,
we are forced to use mutable arrays because we don't have such
fancy monolithic array "primitives" as accumArray.

I've written and tested this as a literate Gofer script because I've
never been able to get GHC to run under Mach. :-(

Let me know if you have any problems...

Chris


- ----------------------------------------------------------------------

FIBONACCI HEAPS

Fibonacci heaps are a priority queue data structure supporting the
following operations:
        O(1) Insert
        O(1) FindMin
        O(1) Meld
    O(log n) DeleteMin

(In an imperative settting, Fibonacci heaps also support
        O(1) DecreaseKey (of an indicated element)
    O(log n) Delete (an indicated element)
but these operations are problematic in a functional setting.)

There is one catch: for the DeleteMin operation, the bounds are
amortized instead of worst-case.  This means that the bounds are
only guaranteed if you use the data structure in a single-threaded manner.
Otherwise, you can take longer than expected by repeatedly going back
and operating on an "expensive" version of the data structure.

(Note: I am currently working on a paper with another student describing
a functional priority queue achieving the above bounds in the worst-case
instead of amortized.  This data structure may be freely used in a
non-single-threaded manner with no ill effects.)

To understand the implementation of Fibonacci heaps, it is helpful to
first understand binomial queues.  See, for example, David King's
"Functional Binomial Queues" from the last Glasgow workshop.

> -- partain
>module Main (main) where
>import ST
>import Array
>import MutableArray
>import System

                         --------------------

Like binomial queues, Fibonacci heaps are based on heap-ordered
binomial trees.

>data Tree a = Node a [Tree a]

The degree of a binomial tree is equal to its number of children.  
Every binomial tree of degree k has binomial trees of degrees
k-1...0 as children, in that order.  It is easy to show that
a binomial tree of degree k has size 2^k.

The fundamental operation on binomial trees is linking, which compares
the roots of two binomial trees and makes the larger a child of the 
smaller (thus bumping its degree by one).  It is essential that this
only be called on binomial trees of equal degree.

>link (a @ (Node x as)) (b @ (Node y bs)) =
>  if x <= y then Node x (b:as) else Node y (a:bs)

It will also be useful to extract the minimum element from a tree.

>root (Node x _) = x

We will frequently need to tag trees with their degrees.

>type TaggedTree a = (Int,Tree a)
>
>degree (k, t) = k
>tree (k, t) = t

Given a tagged tree, extract and tag its children.

>getChildren (n, Node x ts) = zipWith (,) [n-1,n-2 .. ] ts

Extract the minimum element from a tagged tree.

>root' = root . tree

                         --------------------

We also need a type for bags supporting constant time union.  The simple
representation given here is sufficient since we will always process bags
as a whole.  Note that for this application it is not necessary to
filter out occurences of EmptyBag.  Also, for this application order
is irrelevant.

>data Bag a = EmptyBag | ConsBag a (Bag a) | UnionBags (Bag a) (Bag a)
>
>bagToList b = flatten b []
>  where flatten EmptyBag xs = xs
>        flatten (ConsBag x b) xs = flatten b (x:xs)
>        flatten (UnionBags b1 b2) xs = flatten b1 (flatten b2 xs)
>
>applyToAll :: (a -> ST s ()) -> Bag a -> ST s ()
>applyToAll f EmptyBag = return ()
>applyToAll f (ConsBag x b) = f x >> applyToAll f b
>applyToAll f (UnionBags b1 b2) = applyToAll f b1 >> applyToAll f b2

                         --------------------

Miscellaneous stuff.

>log2 1 = 0
>log2 n = 1 + log2 (n `div` 2)

>data MyMaybe a = Zero | One a

                         --------------------

Since binomial trees only come in certain, fixed sizes, we need some
way to represent priority queues of other sizes.  We will do this
with a forest of trees summing to the correct size.

>type Forest a = Bag (TaggedTree a)

In binomial queues, this forest must be maintained in strictly increasing
order of degree.  For Fibonacci heaps, we adopt a more relaxed attitude: 
degrees may be repeated and order does not matter.

To be able to find the minimum element quickly, we keep the tree with the 
minimum root outside of the bag.  In addition, at the top level of each heap,
we store the total size of the heap.

>data FibHeap a = EmptyFH | FH Int (TaggedTree a) (Forest a)

                         --------------------

Now, the following operations are trivial.

>emptyFH = EmptyFH
>
>isEmptyFH EmptyFH = True
>isEmptyFH (FH _ _ _) = False
>
>singleFH x = FH 1 (0, Node x []) EmptyBag
>
>insertFH x xs = meldFH (singleFH x) xs
>
>minFH EmptyFH = error "minFH EmptyFH"
>minFH (FH n tt f) = root' tt

                         --------------------

Meld achieves its efficiency by simply unioning the two forests.

>meldFH EmptyFH xs = xs
>meldFH xs EmptyFH = xs
>meldFH (FH n1 tt1 f1) (FH n2 tt2 f2) =
>  if root' tt1 <= root' tt2 then
>      FH (n1+n2) tt1 (ConsBag tt2 (UnionBags f1 f2))
>  else
>      FH (n1+n2) tt2 (ConsBag tt1 (UnionBags f1 f2))

                         --------------------

Finally, the only hard operation is deleteMin.  After throwing away the
minimum element, it repeatedly links trees of equal degree until
no such pairs are left.  The most efficient way to do this is with
an array.  I give two implementations, one using monadic arrays,
the other using accumArray.

In the first implementation, there are three steps.
  1. Allocate an array indexed by degrees.
  2. Insert every tree into this array.  If, when inserting a tree of 
     degree k, there already exists a tree of degree k, link the
     two trees and reinsert the new larger tree.
  3. Transfer the trees into a bag, keeping track of the minimum tree.

>deleteMinFH EmptyFH = error "deleteMinFH EmptyFH"
>deleteMinFH (FH 1 tt f) = EmptyFH
>deleteMinFH (FH n tt f) =
>  let
>    d = log2 (n-1) -- maximum possible degree
>
>    ins a (i, t) =
>        readArray a i >>= \e ->
>        case e of
>          Zero   -> writeArray a i (One t)
>          One t2 -> writeArray a i Zero >>
>                    ins a (i+1, link t t2)

Note that after inserting all the trees, the array contains trees
in the same pattern as the bits of n-1.  Since we know that the
highest order bit of n-1 is one, we know that there is a tree in
the highest slot of the array.

>    getMin a =
>        readArray a d >>= \e ->
>        case e of
>          Zero  -> error "must be One" -- since array is filled as bits of n-1
>          One t -> getMin' a d t EmptyBag 0
>    getMin' a mini mint b i =
>        if i >= d then
>          return ((mini, mint),b)
>        else
>          readArray a i >>= \e ->
>          case e of
>            Zero  -> getMin' a mini mint b (i+1)
>            One t -> if root mint <= root t then
>                       getMin' a mini mint (ConsBag (i, t) b) (i+1)
>                     else
>                       getMin' a i t (ConsBag (mini, mint) b) (i+1)
>            
>  in 
>    runST (newArray (0,d) Zero >>= \a ->
>           applyToAll (ins a) f >>
>           sequence (map (ins a) (getChildren tt)) >>
>           getMin a >>= \ (tt,f) ->
>           return (FH (n-1) tt f))

                         --------------------

The second version of deleteMin uses accumArray to group trees of like
size.  It then performs the linking and all remaining steps purely 
functionally.

>deleteMinFH' EmptyFH = error "deleteMinFH EmptyFH"
>deleteMinFH' (FH 1 tt f) = EmptyFH
>deleteMinFH' (FH n tt f) =
>  let
>    d = log2 (n-1) -- maximum possible degree
>
>    a = accumArray (flip (:)) [] (0,d) (getChildren tt ++ bagToList f)
>
>    doLinks (ts:rest) = startup 0 ts rest
>      where startup i [] [] = []
>            startup i [] (ts:rest) = startup (i+1) ts rest
>            startup i ts [] = combine i ts [] []
>            startup i ts (next:rest) = combine i ts next rest
>
>            combine i [] next rest = startup (i+1) next rest
>            combine i [t] next rest = (i, t) : startup (i+1) next rest
>            combine i (t1:t2:ts) next rest =
>                combine i ts (link t1 t2 : next) rest
>
>    getMin (tt:rest) = foldl chooseMin (tt,EmptyBag) rest
>      where chooseMin (tt1,b) tt2 =
>                if root' tt1 <= root' tt2 then
>                    (tt1,ConsBag tt2 b)
>                else
>                    (tt2,ConsBag tt1 b)
>
>    (new_tt,new_f) = getMin (doLinks (elems a))
>  in
>    FH (n-1) new_tt new_f

                         --------------------

Testing...

>fibToList :: (Ord a) => FibHeap a -> [a]
>fibToList xs = if isEmptyFH xs then []
>               else minFH xs : fibToList (deleteMinFH xs)
>
>fibToList' :: (Ord a) => FibHeap a -> [a]
>fibToList' xs = if isEmptyFH xs then []
>                else minFH xs : fibToList' (deleteMinFH' xs)
>
>makeFH :: (Ord a) => [a] -> FibHeap a
>makeFH xs = foldr insertFH emptyFH xs
>
>fibSort :: (Ord a) => [a] -> [a]
>fibSort = fibToList . makeFH
>
>fibSort' :: (Ord a) => [a] -> [a]
>fibSort' = fibToList' . makeFH
>
>randoms :: Int -> [Int]
>randoms n = take n (iterate (\seed-> (77*seed+1) `rem` 1024) 1967)
>
>test n = fibSort (randoms n) == fibSort' (randoms n)

>--partain
>main = getArgs >>= \ [n] ->
>	putStr (show (test (read n)))