Summary Algorithms and Data Structures
(c) 20040224 Fabian M. Suchanek
http://suchanek.name
This is a summary of the course "Algorithms and Data Structures"
held by Priv.Doz. Dr. Peter Sanders and Dr. Kavitha Telikapalli
in the WS 2003 at Saarland University. Please note that, although
I am also a tutor of this course, this is a student's summary and
not official.
I would like to thank the participants of my tutorial, who helped
finding the bugs in this summary, especially David Steurer,
Gernot Gebhard and Peter Kolloch.
Furthermore, I owe thank to Thomas Schultz, who scanned the largest
part of the summary rigorously for errors.
___ 
Since this text contains some ASCIIart, it should only / \  /
be viewed or printed with a fixed font (eg Courier New). ( ) <
If the font is all right, an "OK" appears here: \___/  \
By reading the following text, you accept that the author does not
accept any responsibility for the correctness or completeness of
this text. If you have any corrections or remarks, please
send me a mail. This is the only way to make the publication of this
summary useful for me, too.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Prerequisites
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Set: An unordered collection of different elements.
S = {s1,...sn}
S = { x  p(x) } is the set of all those x for which p(x) holds
Epsilon of a set S: An element of S.
Notation: eps(S) := x, such that x in S
eps({x}) := x
// This operator has been invented by Hilbert. Since it is
// the only way to convert a set to an element, this summary makes
// use of it.
Size, cardinality of a set: The number of its elements.
Notation: {s1,s2,...sn} = n
Intersection of two sets: The set of those elements which occur in
both of them.
Notation: X /\ Y
Union of two sets: The set of all of their elements.
Notation: X \/ Y
Difference of two sets: The set of those elements which occur in the
first one, but not in the second one.
Notation: X \ Y
Disjoint sets: Sets, the intersection of which is empty.
Picking: Taking an element of a set and deleting it from the set
pickAny(S) := (z:=eps(S), S:=S\{z}, z)
Multiset: An unordered collection of elements. Unlike a set, a multiset
may contain the same element more than once.
Tupel: An ordered collection of elements.
A tuple can be assigned to a set, but not vice versa.
Notation: t=(t[1],t[2],...t[n])
Example: t=(1,'A',abc) is a tupel
Dimension of a tupel: The number of its elements.
Pair: A tupel of dimension 2.
Triple: A tuple of dimension 3.
Subset of a set S: A set, the elements of which are all contained in S.
Notation: S' =< S
Superset of a set S: A set, which contains all elements of S.
Notation: S' >= S
Partition of a set S: A set of subsets S[1],...S[n] of S, such that
each element of S is in exactly one of them.
Cartesian Product of two sets A and B: The set of all possible pairs,
the first element of which is in A and the second element of which
is in B.
Notation: A x B = C
Example: {1,2} x {a,b,c} = { (1,a), (2,a), (1,b), (2,b), (1,c), (2,c) }
nth Power of a set S: The set
S^n := S x S x S x ... S
where 'S' occurs n times.
Relation over two sets A and B: A subset of the cartesian product of A
and B.
Notation:
* (a,b) is an element of the relation R
* r(a,b)
* a r b
* (a,b) is in R
Example: The relation ">" is the set of all pairs of real numbers,
in which the first number is greater than the second.
It holds: ">" = { (1,0), (27.3,12.5), (17,12), ... }
One writes: >(3,2) or 3 > 2 or (3,2) is in ">".
Reflexive relation: A relation R < A x A such that R(a,a) for all a
in A.
Symmetric relation: A relation R < A x A such that
if R(a,b), then R(b,a)
for all a,b in A.
Antisymmetric relation: A relation R < A x A such that
if R(a,b) and R(b,a) then a=b
for all a,b in A.
Asymmetric relation: A relation R < A x A such that
if R(a,b), then not R(b,a)
for all a,b in A.
Complete relation: A relation R < A x A such that
R(a,b) or R(b,a)
for all a,b in A.
Transitve relation: A relation R < A x A such that
if R(a,b) and R(b,c), then R(a,c)
for all a,b,c in A.
Strict relation: A relation R < A x A such that
not R(a,a)
for all a in A.
Linear order on a set A: A relation R < A x A, which is asymmetric,
transitive and complete. Usually, this relation is denoted by '<'.
Example: '<' is a linear order on the set of natural numbers.
Simple order on a set A: A relation R < A x A, which is reflexive,
antisymmetric, transitive and complete. Usually, this relation is
denoted by '=<'.
Example: '=<' is a simple order on the set of natural numbers.
Domain of a relation over A and B: A.
Range of a relation over A and B: B.
Function: A relation f over two sets A and B, such that for every a
in A there is either none or exactly one element b in B with a f b.
Notation for f: f:A>B
Notation for "a f b": f(a)=b.
Notation for the set of all functions f:A>B: B^A
There are B^A = B^A many of these functions f:A>B.
Result of a function f:A>B for an argument x: The element of B
such that f(x)=b.
Real: The set of real numbers.
For a detailed definition and associated operations, cf
"Constraints over the Reals" (cr.txt)
Real+: The set of positive real numbers with 0.
Integer: The set of integer number.
For a detailed definition and associated operations, cf
"Constraints over the Reals" (cr.txt)
Natural: The set of positive integers, including 0.
Interval with the real limits a and b: The set
[a,b] := {x  a =< x =< b}.
Prime: A natural number greater than 1 which can only be divided by
1 and itself.
Order of a function f:Natural>Real: A function g:Natural>Real, such
that there is x0 in Real and c in Real, such that f(x)x0.
Class of an Order f:Natural>Real: The set of functions having f as
an Order.
Notation: O(f)
Notation for "g:Natural>Real is an element of O(f)": g in O(f)
Necessity of a state A for a state B: The fact that B cannot be
without A, B => A.
Sufficiency of a state A for a state B: The fact that, whenever A is,
B is also, A => B. In order to prove the equivalence of two states,
it is sufficient to prove necessity and sufficiency.
iff: if and only if.
Operation, Algorithm: A sequence of instructions, which may have a
result.
Recursive algorithm: An algorithm X, one instruction of which is the
execution of X itself.
Termination of an algorithm: Its property that it cannot be executed
forever.
Partial correctness of an algorithm: Its property that it gives the
right result, if it terminates. Partial correctness does not yet
guarantee that the algorithm terminates.
Total correctness of an algorithm: Its property that it is partially
correct and terminates. In this case, the algorithm does really
what it is intended to do.
Time of an algorithm for an input size n: The order of the
function which calculates the time needed to execute the algorithm,
in dependence upon the input size n. One says "The algorithm is in
time O(g)".
Polynomial in x: An expression of the form
a[0]*x^0 + a[1]*x^1 + ... a[n]*x^n
where x is a real variable and the a[i] are real constants.
// See "Constraints over the Reals" (cr.txt) for more details
Algorithm with a polynomial time: An algorithm, the time of which
can be expressed by a polynomial in the input size.
NPhard problem: A problem, for which no algorithm is known, which
solves the problem in polynomial time.
Space of an algorithm for an input size n: The order of the
function which calculates the number of cells needed to execute
the algorithm, in dependence upon the input size n. One says "The
algorithm is in space O(g)".
RAM, Random Access Machine: A model of a computer with an infinite
amount of cells. A RAM can execute certain instructions.
For a detailed definition cf "Theoretische Informatik" (infod.txt,
German)
Primitive type: A set of finite size.
Boolean: The primitive type {true, false}.
Class, Type: A primitive type or a tuple of types together with a
tuple of names. Usually, there are operations associated to types.
Notation for tuple types: ( : , : , ... )
or
CLASS EXTENDS :
...
PROCEDURE (...) ...
...
The operations can only be called with an instance (s.b.). This
instance is referred to as in the procedure. For convencience,
the string "this." can be left away when accessing components of
. The EXTENDSclause says that the type inherits all components
and operations of the indicated type. If a class inherits an
operation and also defines it, then the newly defined operation is
used for all instances of the extending class.
Data structure, Instance of a type t: If t is a primitive type, then
an element of t, if t is a tuple of types, then a tuple of data
structures of the elements of t.
Notation for the component named c of a datastructure x: x.c
Notation for the operation named f of a datastructure x: x.f
Notation for a datastructure of type t:
(:, :, ...)
or
(, , ...)
Notation for a datastructure of type t constructed from a
datatstructure x of some other type: t(x)
1dimensional array, array, sequence of type t: A type, which is
a tuple of t's.
Notation for data structures: s=
Associated operations:
* s[i] :=
* s.length() :=
* s1 + s2 :=
* s[i..j] :=
nil: A value, which is compatible with all types.
Singly linked list of type t: The type (car:t, cdr:list of type t).
That is: A list data structure is a pair of a data structure and
a list data structure.
CLASS List of T
car : T
cdr : List of T
FUNCTION length() : Integer
RETURN (this==nil) ? 0 : 1+this.cdr.length
PROCEDURE insert(i : Integer, x : T)
IF (i==0) THEN
this:=List(x,this)
ELSE
this.cdr.insert(i1)
PROCEDURE pushFront(x : T)
this.insert(0,x)
PROCEDURE pushBack(x : T)
this.insert(this.length,x : T)
FUNCTION popFront() : T
RETURN this.delete(0)
FUNCTION popBack() : T
RETURN this.delete(this.length1)
FUNCTION last() : T
RETURN (this.cdr==nil) ? this.car : this.cdr.last()
FUNCTION delete(i : Integer) : T
IF i==0 THEN
x:=this.car;
this:=this.cdr;
RETURN x
ELSE RETURN this.cdr.delete(i1)
FUNCTION delete(x : T) : T
IF this==nil RETURN nil
IF this.car==x THEN
RETURN this.delete(0)
RETURN this.cdr.delete(x)
FUNCTION get(i : Integer) : T
IF i==0 THEN RETURN this.car
RETURN this.cdr.get(i1)
PROCEDURE +(l : List of T)
this.cdr:=l
List of type t: A type, which has operations equivalent to
those of a singly linked list of type t.
Doubly linked list of type t: The list
CLASS doublyLinkedList of type T
element : T
prev : doublyLinkedList of type T
succ : doublyLinkedList of type T
with the appropriate operations.
ListStack of type t: The type
CLASS ListStack of T EXTENDS List of T
FUNCTION pop() : T
RETURN this.delete(0)
FUNCTION push(x : T)
RETURN this.insert(0,x)
FUNCTION top() : T
RETURN this.car
FUNCTIOn empty() : Boolean
RETURN this.length==0
Stack of type t: A type, which has operations equivalent to
those of a liststack of type t.
Naive Queue of type t: A list of type t with the following operations
CLASS NaiveQueue of T EXTENDS List of T
this.insert(this.length,x)
FUNCTION dequeue(), q.popFront() : T
RETURN this.delete(0)
FUNCTION head() : T
RETURN this.car
FUNCTION empty() : Boolean
RETURN this.length==0
Queue of type t: A type, which has operations equivalent to
those of a naive queue of type t.
Naive Priority Queue of type t and priority f:t>Real: The type
CLASS NaivePQueue of T with f:t>Real EXTENDS List of T
PROCEDURE insert(x)
IF this==nil  f(this.car)>f(x))
THEN this:=NaivePQueue(x,this)
ELSE this.cdr.insert(x)
FUNCTION deleteMin() : T
RETURN this.delete(0)
PROCEDURE updateKey(x : T), decreaseKey(x : T)
this.delete(x)
this.insert(x)
FUNCTION top() : T
RETURN this.car
Priority Queue of type t and priority f:t>Real: A type, which
has operations equivalent to those of a Naive Priority Queue of type
t and priority f.
Naive heap of type t and priority f:t>Real: A naive priority queue,
which supports the following additional operations:
CLASS NaiveHeap of T with f:t>Real EXTENDS NaivePQueue of T
FUNCTION pred(x : T) : T
IF this==nil THEN RETURN nil
RETURN this.cdr.car==x?this.car:this.cdr.pred(x)
FUNCTION succ(x : T) : T
IF this==nil THEN RETURN nil
RETURN this.car==x?this.cdr.car:this.cdr.succ(x)
PROCEDURE swap(x,y : T) // only neighbors
IF this==nil THEN RETURN
IF this.car==x && this.cdr.car==y THEN
(l.car,l.cdr.car) := (l.cdr.car, l.car)
ELSE this.cdr.swap(x,y)
FUNCTION delete(x : T) : T
IF (this==nil) THEN RETURN nil
IF this.car==x THEN
this:=this.cdr
RETURN x
ELSE RETURN this.cdr.delete(x)
FUNCTION locate(x : T) : T
// Returns x or next smaller element
IF this==nil RETURN nil
IF this.car==x  f(this.cdr.car)>f(x) THEN
RETURN this.car
RETURN this.cdr.locate(x)
Heap of type t and priority f:t>Real: A type, which has operations
equivalent to those of a Naive Heap of type t and priority f.
Reverse of an array X: The array reverse(X) := .
Matrix of type t: A sequence of n sequences of type t. The inner
sequences all have to be of the same length m. The matrix is then
called an n x m matrix. A vector can be seen as n x 1 matrix.
Transposed matrix of a n x m M: The n x m matrix t(M) with
t(M)[j][i] = M[i][j]
Determinant of a 3 x 3 matrix A of type Real: The value
det(A) := a*e*i + b*f*g + c*d*h  c*e*g  b*d*i  a*f*h
if
a b c
A = d e f
g h i
That is: Sum up the products of all \Diagonals, subtract the products
of all /Diagonals.
// For details, see "Constraints over the Reals" (cr.txt)
Position, Index of an element in an array data structure: The number
of its predecessors in the array.
Vector: The type of a tuple of reals. Its operations are defined in the
following.
// For a detailed definition, see "Algebra" (algebra.txt, German)
// or "Analyse de Donnees" (maths.txt, in French)
Sum of two vectors x and y: The vector
x+y := (x[0]+y[0], ..., x[x.length1]+y[x.length1]).
rMultiple of a vector x: The vector vector
r*x := (r*x[0], ..., r*x[x.length1]).
Negative of a vector x: The vector x := (x[0],...x[x.length1]).
Absolute value of a vector: The root of the sum of the squares
of its components.
Notation: x := sqrt(SUM {x[i]^2  i=1..x.length})
Dotproduct of two vectors x and y: The real value
x * y := SUM { x[i]*y[i]  i=1..x.length1 }
The dotproduct measures the "similarity" of the two vectors. It is
0 if the vectors are orthogonal.
// For details, cf "Analyse de Donnees" (maths.txt, in French)
ndimensional array of type t: An array, the elements of which
are (n1)dimensional arrays of type t.
Notation for a[i[0]][i[1]]...[i[n]]: a[i] for an integer vector i
Binary representation of a natural number x: The queue
binrep(x) := (x==0) ? <0> : binrep(x div 2).enqueue(x mod 2)
Natural number of a binary representation b: The natural number
unbin(b) := _unbin(b,1)
_unbin(b,i) := (b==nil) ? 0 :
i*(b.last() mod 2)+unbin( (b.delete(b.length1),b),
i*2 )
Bit: An element of {0,1}.
xor: A function mapping two natural numbers a and b to the natural
number of the binary representation, which is achieved by setting
those bits to 1, for which the binary representations of a and b
differ. It holds:
a xor a = 0
a xor b = c => a = b xor c
a xor b = b xor a
a xor a xor b = b
Random variable: A function, which maps individuals to real numbers.
A random variable represents a characteristic of the individuals.
For each individual, the random variable "knows" the degree of
markedness of this characteristic.
Example: Consider the age of a person. A random variable might
be A, such that A(grand_ma)=99, A(fabian)=23 etc
// For a more detailed definition see "Statistik" (statistik.txt,
// German)
Expectation of random variable X: The sum of the results of X for all
possible individuals, divided by the total number of all possible
individuals. This is a kind of theoretical mean value.
Notation: E(X)=y
Linearity of expectation: The fact that E(X + Y) = E(X) + E(Y).
Probability of a random variable X taking the value y: The number of
all possible individuals i, for which X(i)=y, divided by the total
number of all possible individuals.
Notation: P(X=y)
Probability of X=x if Y=y: The probability that the random variable X
takes the value x, provided that the random variable Y takes the value
y.
Notation: P(X=x  Y=y)
Independant random variables: Two random variables X and Y, such that
P(X=x  Y=y) = P(X=x). That is: Whether Y=y or not, this does not
influence P(X=x).
Probability of X=x and Y=y: The probability that the random variable
Z(i) := (X(i)==x and Y(i)==y) ? 1 : 0
takes the value 1. It holds
P(X=x and Y=y) := P(X=x)*P(X=x  Y=y)
Probability of X=x or Y=y: The probability that the random variable
Z(i) := (X(i)==x or Y(i)==y) ? 1 : 0
takes the value 1. It holds
P(X=x or Y=y) := P(X=x)+P(Y=y)  P(X=x and Y=y)
kth harmonic number: The sum
H(k) := 1/1 + 1/2 + 1/3 + ... + 1/k
It holds that H(k) < ln(n)+1.
Set sum: The function SUM, which maps a set of real numbers to their
sum.
Notation: SUM {x[1],x[2],...x[n]} := x[1]+x[2]+...+x[n]
Lemma, Theorem: A fact.
Sorted sequence of a set A: A sequence S of type A, which contains all
elements of S, such that S[i] Q // Transition function
FUNCTION runOn(S : array of I) : boolean
state := s0
pos := 0
WHILE pos<S DO
state := d(state,S[pos])
pos := pos + 1
RETURN (state in F)
A FSA "reads" a string symbol by symbol. At each symbol read, the
state is changed according to the function d. If the string is over
and the FSA remains in a state marked as "final", then the FSA is
said to "accept" the string. A FSA may be depicted as a graph where
every node stands for a state and the edges show a transition from
one state to another, given a particular input symbol. A FSA may also
be depicted by a transition table, which shows the new state(s) for
a particular current state and a particular input.
// For details, see "Computational Linguistics" (cl.txt) or
// "Theoretische Informatik" (infod.txt, in German)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Graphs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
See also "Graphes et Complexite" (gc.txt, in French)
Graph: The type
CLASS graph
vertices : Set // often noted V
edges : Subset of vertices x vertices // often noted E
A graph is usually seen as consisting of "nodes" (the elements
of V) and "connections" between them (the elements of E). Graphs
are usually drawn as networks, in which the nodes are connected
by arrows. In the following, V is always the set of nodes and
E is always the set of edges. Graphs usually have the following
operations:
* Access to an edge
* Find outgoing edges of a node
* Find an edge
* Insert edge
* Delete edge
* Insert node
* Delete node
Vertex, Node of a graph (V,E): An element of V. The number of nodes
is denoted by n. Usually, the nodes are simply integers 0..n1.
Edge of a graph (V,E): An element of E. The number of edges is
denoted by m. For an edge e=(u,v), one says that e goes from
u to v. Sometimes the edge (u,v) is noted by "u>v". Usually,
the edges are numbered. We define
CLASS edge
startnode : Node
endnode : Node
Reverse edge of an edge (u,v): The edge (v,u).
Going to a node, looking at a node: Shifting attention towards
this node.
Traversing an edge, following an edge: Going from its startnode
to its endnode.
Successor of a node u: A node v, for which there is an edge (u,v).
Outdegree of a node v: The number of edges going out of v.
outdegree(v) := {(v,u) in E}
Indegree of a node v: The number of edges going to v.
indegree(v) := {(u,v) in E}
Selfloop: An edge, the startnode of which is its endnode.
Usually, selfloops are disallowed.
Bidirected graph: A graph, in which for every edge
(u,v), there exists one edge (v,u).
Undirected graph: A bidirected graph, in which the edges (u,v) and
(v,u) are considered equal. This is why the edge is usually
noted by {u,v}. An undirected graph is usually drawn with simple
lines instead of arrows between the nodes.
Degree of a node v: The number of edges involving v.
If the graph is directed, then it holds
degree(v) := indegree(v) + outdegree(v)
In a graph (V,E), the number of nodes with an odd degree is even:
SUM {degree(s)  degree(s) mod 2==1}
+ SUM {degree(s)  degree(s) mod 2==0} = SUM {degree(s)} = 2*E
2*E is even
SUM {degree(s)  degree(s) mod 2==0} is even
=> SUM {degree(s)  degree(s) mod 2==1} is even
=> {degree(s)  degree(s) mod 2==1} is even
Subgraph of a graph (V,E): A graph (V',E') with V'=Real. These
edge costs can be thought of as an amount of money to be paid when
the edge is traversed.
Notation: c(u,v) := c((u,v)) for an edge (u,v)
Weighted graph: A graph with edge weights. We define
CLASS WeightedGraph EXTENDS Graph
c:FUNCTION edges>Real
Path in a graph (V,E): A sequence of nodes , such that
v[i] in V and (v[i],v[i+1]) in E. A path can also be given as the
sequence of the respective edges: <(v[0],v[1]), (v[1],v[2]), ...
(v[k1],v[k])>. Usually, a path is noted by v[0]>v[1]>...v[k] or
as v[0]~~>v[k].
Start node of a path: Its first element.
End node of a path: Its last element.
Node reachable from a node s: A node t, for which there is a path from
s to t.
Length of a path in a graph: The number of its nodes minus 1.
Notation: length() := k
Weight, Cost of a path p in a weighted graph: The sum of the
costs of the edges of p.
Notation: c() := SUM {c(v[i],v[i+1])  i=0..k1}
Simple path: A path, the nodes of which are all distinct.
Walk: A path, which is not simple.
Bottleneck weight of a path p in a weighted graph: The minimum of
the weights of the edges of p.
Cycle: A path, the start node of which is at the same time its end
node.
Simple cycle: A cycle, the nodes of which are all different
except for the start and end node.
Hamiltonian cycle: A simple cycle which contains all nodes.
Finding a Hamiltonian cycle is an NPcomplete problem.
Euler cycle: A cycle which contains all edges once.
Connected graph: A graph, in which there is for two nodes always a
path from one to the other, if one adds to each edge (u,v) its
reverse edge (v,u). This definition simply says that the graph does
not consist of separated components.
Strongly connected graph: A connected graph, in which there is a path
from every node to every other node.
Complete graph: A graph (V,V x V). That is: Every node is connected to
every other node.
Eulerian graph: An undirected graph which contains an Euler cycle. A
graph is eulerian, iff
* it is connected and
* every node has even degree
Necessity:
* If the graph was not connected, there could not be a cycle touching
all its edges. Therefore, the graph must be connected.
* If the Euler cycle enters into a node, it must also quit this node
again. Since every edge is only to be visited once, each passing
of a node requires 2 separate edges.
Sufficiency:
By the following algorithm, an Euler cycle can be found if the above
conditions are met:
FUNCTION eulerCycle_(VAR (V,E):Graph,s: V): List of E
IF s==nil THEN s := eps(V)
result : List of E := <>
u : V := s
REPEAT
next : V := eps({v  (u,v) i E})
List.pushBack((u,next))
E := E \ {(u,next)}
u := next
UNTIL u==s
WHILE Ex i: degree(result.get(i).endNode) mod 2==1 DO
s := result.get(i)
small : List of E := eulerCycle((V,E),s.endNode)
small.cdr := s.cdr;
s.cdr := small
RETURN result
FUNCTION eulerCycle((V,E):Graph): List of E
RETURN eulerCycle_((V,E):Graph,eps(V))
Start at any node. Follow an edge to another node and delete this
edge. Continue this way until you get to a node with no more edges.
Since every node has an even degree, each node can be entered and
left. If one gets to a node with no exit, then this must be the
start node, because it is the only one with an uneven degree.
Start at another node and do the same. Continue until no more edges
are left. This way, many cycles are created. Every one of them has
at least one node in common with another cycle. If this were not
the case, the graph would not be connected. These common nodes
can be used to glue together the cycles, creating one large
cycle  the Euler cycle. The algorithm runs in time O(E).
acyclic graph: A graph which does not contain any cycles.
Undirected tree: An undirected graph, in which one of the following
properties holds:
1. There is exactly one path from one node to another
2. The graph is connected and has exactly n1 edges
3. The graph is connected and acyclic.
These properties are equivalent:
* 1 => 2
If there is a path from one node to another, then the graph is
connected. If there is a path from one node to another, then
the graph cannot have less than n1 edges, since every node needs
one edge. If the graph had more than n1 edges, it would have
a superfluous edge, which could be used to construct an alternative
path.
* 2 => 3
If the graph is connected, then it is connected. The graph has
n1 edges. Suppose it contains a cycle. Then we can delete one
edge without destroying the connectedness. Then the graph contains
n2 edges, which can never connect all nodes. This contradiction
implies that the graph cannot contain a cycle.
* 3 => 1
Since the graph is connected, there exists one path from one node
to another. Suppose there were two paths from one node u to another
node v. Then u~~>v~~>u would be a cycle. Since 3 says that the
graph does not contain any cycle, each path must be unique.
Forest: An undirected acyclic graph. A forest can be seen as consisting
of multiple trees (which is quite a reasonable definition). Such a
graph can at most have n1 edges. If it had more, it would contain a
cycle, see definition of "undirected tree".
Directed tree: A directed graph, which has a node r, such that there is
exactly one path from r to every other node. This entails
* Connectedness
If the graph was not connected, there could not be this path from
r to a node in another component.
* No cycles
If the graph allows a cycle around a node u, then there are
infinitley many paths from r to u, namely r~~>u, r~~>u~~cycle~~>u,
etc.
* At most one incoming edge
Every additional incoming edge would allow another path from r.
If all edges are reversed in such a graph, then this is also called a
tree.
Root of a directed tree: The node, from which there is a path to every
other node.
Parent of a node u in a directed tree: The node for which there is
an edge to u. The parent of the root is the root.
Child of a node u in a directed tree: A node v for which there is an
edge (u,v). A node may have multiple children.
Leaf of a directed tree: A node which has no children.
Interior node of a directed tree: A node which has children.
Ordered directed tree: A directed tree, in which each node maintains
an ordering of its children.
DAG, directed acyclic graph: A graph without cycles. This type of
graphs is for instance used in multiple inheritance "trees".
A graph can be tested for being a DAG by the following algorithm:
FUNCTION isDAG((V,E) : Graph) : Boolean
WHILE Ex v in V: outdegree(v)==0 DO
V := V \ {v}
E := E \ {(u,v}  u in V} \ {(v,u}  u in V}
RETURN (V=={})
That is: While there is still a node without outgoing edges, delete
this node. Iff no node remains, then the graph is a DAG. This is
because if the graph contains a cycle, then this cycle is supported
by outgoing edges. Each node protects its predecessor from being
deleted.
// For an application of DAG, see "Representing Ontological Structures
// in CLOS, Java and FAST" (ba.doc)
Graph representation: A type, which provides the operations of graph
in an efficient way on a computer.
// Since the operations can mostly be implemented trivially,
// this summary only gives their runtimes
Edge sequence representation: The graph representation
CLASS edgeSequence
edges : List of Edge
nodes : List of Node
This class uses a list of all edges <(u1,u2),(u5,u2),
...> and a list of all nodes of the graph. (The node array is necessary
to store nodes which have no edges.) Operations with their time:
* Access to an edge: O(m)
* Find outgoing edges of a node: O(m)
* Find an edge: O(m)
* Insert edge: O(1)
* Delete edge: O(m) (first find the edge)
* Insert node: O(1)
* Delete node: O(n)
Additional information for edges or nodes can be stored in partner
arrays of the two arrays.
Space: O(m+n)
Adjacency array: The graph representation
CLASS adjArray
edgeArray : array of Node
nodeArray : array of Integer
This class uses two arrays: The "edge array" stores, for each node u,
the end nodes of all edges going out of u. The "node array" stores,
for each node u, the index of the cell in the edge array, which
contains the first outgoing edge of u. The last element of the node
array stores the length of the edge array.
Operations with their time:
* Find outgoing edges of a node: O(1)
* Find an edge: O(m)
* Insert edge: O(m+n), needs to create a new edge array
* Delete edge: O(m), a second node array can be used, which stores
the index of the last edge of a node.
* Insert node: O(n), needs to create a new array
* Delete node: O(n), the node array needs to be rearranged
Additional information for edges and nodes can be stored in partner
arrays of the two arrays.
Space: O(m+n)
Adjacency list: The graph representation
CLASS adjList
list : array of List of Node
Here, each array element with index i stores a list of the
successors of node i. Operations with their time:
* Find outgoing edges of a node: O(1)
* Find an edge: O(m)
* Insert edge: O(1)
* Delete edge: O(m) (first find the edge)
* Insert node: O(n), needs to create a new array
* Delete node: O(n), the array needs to be rearranged
Space: O(m+n), but three times as much as adjacency arrays.
Adjacency matrix: The graph representation
CLASS adjMatrix
matrix : array of array of Bit
The matrix contains a 1 at position i,j iff there is an edge from node
i to node j. Operations with their time:
* Access to an edge: O(1)
* Find outgoing edges of a node: O(n)
* Find an edge: O(1)
* Insert edge: O(1)
* Delete edge: O(1)
* Insert node: O(n^2), needs to create a new array
* Delete node: O(n^2), needs to create a new array
Space: O(n^2)
If an adjacency matrix is raised to the kth power, the position
i,j indicates the number of paths of length k between nodes i and
j. For a definition of kth power and a proof of the path property,
see "Graphes et Complexite" (gc.txt, French).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Graph traversal
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Graph traversal: The process of systematically going through all edges
and nodes of a graph. This process starts in a certain start node.
Shortest path from node u to node v: The path with the start node u
and the end node v, which has the smallest length of all possible
such paths.
Shortest subpath lemma: The fact that if is the
shortest path from u to v, then is the shortest path
from u to v'. Proof: Suppose there was a shorter path from u to
v'. Then this path could be used to go from u via v' to v. This
path would be shorter than the original one, leading to a
contradiction.
BFS, Breadthfirst search: A graph traversal, which first looks at the
start node s, then at all nodes reachable from s by one edge, then
all nodes reachable from s by two edges etc.
Algorithm:
PROCEDURE bfs((V,E) : Graph, s : Node)
q : Queue of Node :=
marked : Array[V] of Boolean :=
marked[s] := true
WHILE q != nil DO
u := q.popFront()
FOREACH (u,v) in E DO
IF marked[v]==false THEN
q.pushBack(v)
marked[v]:=true
The idea of this algorithm is to maintain a queue, which
contains the nodes to be explored next. Since the children of
an explored node are enqueued at the end of the queue, it is
guaranteed that the nodes of one level are finished before
starting with the next level. This algorithm only works for
nodes which are reachable from s.
The algorithm runs in time O(m+n), because every edge and
every node is explored exactly once.
DFSAction: The following type:
CLASS DFSAction
PROCEDURE traverse((u,v) : Edge)
// To be executed if the edge (u,v) is traversed
PROCEDURE root(s : Node)
// To be executed if the node s demarks a new start node
PROCEDURE backtrack((u,v) : Edge)
// To be executed if the edge (u,v) is traversed in the
// opposite sense when there is nothing interesting at node v
Such a DFSAction can be plugged into a graph traversal
algorithm. The algorithm will execute the operations of the
DFSAction in the order it traverses the graph.
// This is the objectoriented interpretation of the blue lines
// in the algorithms in the slides
DFS, Depthfirst search: A graph traversal, which first looks at the
start node and then does depthfirst searches for all children
of the start node.
Algorithm:
PROCEDURE dfs((V,E) : Graph, s : Node, a : DFSAction)
S : Stack of Node := nil
markedNode : Array[V] of Boolean :=
markedEdge : Array[E] of Boolean :=
FOREACH s in V DO
IF markedNode[s]==false THEN
markedNode[s] := true
S.push(s)
a.root(s)
WHILE S != nil DO
v := S.top()
IF Ex e=(v,w) in E : markedEdge[e]==false THEN
a.traverse(v,w)
markedEdge[e]:=true
IF markedNode[w]==false THEN
markedNode[w]:=true
S.push(w)
ELSE
w:=S.pop()
a.backtrack(S.top(),w)
This algorithm maintains a stack, the top of which is the node
currently being explored. The first edge of the current node
is taken and followed. This continues until a node is found,
which has already been visited (visted nodes are marked). Then
the algorithm "backtracks" to the last node in the stack, which
still has unmarked edges. Since this procedure only finds nodes
reachable from the start node, all nodes are given the chance of
being the start node ("FOREACH s in V..."). If a node is found
which has not yet been visited by a previous round, it is chosen
as a new start node. Note that S.top() has to be nil if the stack
is empty, because backtrack will also be called for the start node.
The algorithm runs in time O(m+n), because every edge and
every node is explored exactly once. If this algorithm is to be
run on an undirected graph (represented as a bidirected graph),
then every marking of an edge also has to mark the reverse edge.
The marking of the visited edges can be avoided if adjacency arrays
are used: For each node, one stores not only the index of its edges
in the edge array, but also the index of the last edge visited. One
visits the edges one after the other and increases the
lastvisitededgepointer each time. Thus, each node knows which edge
to visit next without the need to mark the edges.
New node of a graph traversal: A node which has not yet been visited.
In DFS, these nodes are unmarked.
Open node of a graph traversal: A node which has been visited and will
be checked out again for successor nodes. In DFS, these nodes are
marked and kept on the stack.
Finished node, Closed node of a graph traversal: A node which has been
visited and will not be checked again. In DFS, these nodes are marked
and are no longer on the stack.
DFStree of a graph (V,E): The graph (V,E'), where E' is the subset of
those edges over which DFS backtracks. E' can be obtained by plugging
the following DFSAction into the DFSAlgorithm:
CLASS DFSTree EXTENDS DFSAction
parent : Array[n] of Node
PROCEDURE root(s : Node)
parent[s] := s
PROCEDURE backtrack((v,w) : Edge)
parent[w] := v
When an edge (u,v) is backtracked, then the node u stores a
parentpointer to the node v. The set E' can now be derived from
the parentpointers:
E' := {(u,v)  parent[v]==u}
DFSNumbering: A DFS, which assigns to each node a number
indicating the time when it is explored. DFSNumbering can be done by
plugging the following DFSAction into DFS:
CLASS DFSNumbering EXTENDS DFSAction
dfsNum : Array[n] of Integer := <1,...>
d : Integer
PROCEDURE root(s : Node)
dfsNum[s] := 0
d := 1
PROCEDURE traverse((u,v) : Edge)
IF dfsNum[v]==1 THEN
dfsNum[v] := d
d := d + 1
Whenever an edge (u,v) is traversed and has not yet been seen, then
the DFSNumber d is assigned to v and d is increased.
Tree edge of a DFSrun: An edge leading to a new node. That is: For the
edge (u,v), v is unmarked. Tree edges make up a spanning tree.
Spanning tree of a connected graph (V,E): A tree (V,E'), where E'=dfsNum[v]. Back edges lead from a node to its
(grand)parent in the spanning forest.
Forward edge of a DFSrun: An edge leading to a closed node with a
greater dfsNum. That is: For the edge (u,v), v is marked and no longer
on the stack and dfsNum[u]dfsNum[v]. cross edges connect
two nodes, none of which is the (grand)parent of the other in the
spanning forest. A DFS in an undirected graph cannot encounter any
cross edges: Since v is marked, it must have been visited. Since
dfsNum[u]>dfsNum[v], u must have been unvisited at this time. Now
all edges of v will be explored, in particular also {v,u}, which
exists in the undirected graph as well. Thus, {v,u} becomes a
tree edge  and cannot be a cross edge.
Spanning forest by network problem: The problem of, given a set V of
n nodes, O(n) space in memory and m edges one after the other,
calculating a spanning tree of the resulting graph. The problem is
that the edges cannot just all be read and inserted into the graph,
because one only has O(n) space. The problem can be solved as
follows: One reads a set E[1] of 2*n edges. Then, one calculates
a spanning forest of (V,E[1]) by DFS. This takes time O(n+m), i.e time
O(n + n)=O(n) here. It results in the tree edges F[1], where F[1] <
V < E[1]. Now, one reads the next set of edges E[2], and computes
again a spanning forest of (V,F[1] \/ E[2]). This is repeated until
all edges have been read. If there are m edges, then this means
repeating an O(n) process m/n times. The resulting total running time
is O(n+m).
Topological Sorting of a graph: The construction of a list
of all nodes, such that any edge of the graph goes from a prior
node to a later node. To sort a graph topologically, it must not
contain any cycles (it must be a DAG).
A topological sorting can be achieved by plugging the following
DFSAction into the DFS:
CLASS TopSort EXTENDS DFSAction
topologicalOrder : List of Node := nil
PROCEDURE traverse((u,v) : Edge)
IF (u,v) is a back edge THEN THROW exception("Cycle found")
PROCEDURE backtrack((u,v) : Edge)
topologicalOrder.pushFront(v)
Whenever an edge is backtracked, its end node is pushed on the
stack. Thus, only closed nodes are put on the stack. The stack
is in topological order, because whenever a node u is pushed on
it, there is never an edge from another node v on the stack to u.
Proof: Assume there was such an edge
* If it was a tree edge, then v would not be closed and would
not be on the stack
* If it was a back edge, then the algorithm would have signalled
a cycle
* If it was a forward edge, then v would not be closed and would
not be on the stack
* If it was a cross edge, then u would already be closed
Connected component of a graph G: A connected, vertex induced subgraph
of G. That is: A subset of the nodes of G, together with all their
edges among each other, which is connected. Intuitively, a connected
component is a number of nodes, which are all somehow connected by
edges of any directions.
SCC, Strongly connected component of a graph G: A maximal, strongly
connected, vertex induced subgraph of G. Thus, a SCC is a subset of
the nodes of G, together with all their edges. Within this subgraph,
there is a path from every node to every other node.
The SCC can be determined by plugging the following DFSAction
into DFSNumbering:
CLASS SCC EXTENDS DFSNumbering
openNodes : Stack of Node := nil
openComponents : Stack of Node := nil
// These are nodes, each of which represents a component
// A component is on this stack, if there may still be new
// elements for it
component : Array[n] of Node :=
// This array stores for each node one component representative
PROCEDURE root(s : Node)
super.root(s)
openNodes.push(s)
openComponents.push(s)
PROCEDURE traverse( (u,v) : Edge)
super.traverse((u,v))
IF (u,v) is a tree edge THEN
// Every new node is seen as a component per se
openNodes.push(v)
openComponents.push(v)
ELSE IF v in openNodes THEN
// v is already open, we have a back edge and a cycle.
// Mark all nodes on the cycle as belonging to the
// same component
WHILE dfsNum[v] < dfsNum[openComponents.top()] DO
openComponents.pop()
PROCEDURE backtrack( (u,v) : Edge)
IF v == openComponents.top() THEN
// Backtracking over a component representative
// Assign this representative to all open nodes
openComponents.pop()
REPEAT
u:=openNodes.pop()
component[u] := v
UNTIL u == v
This algorithm works due to the following SCCCycleLemma.
SCCCycleLemma: The fact that two nodes are in the same SCC iff they
are on a cycle.
Necessity:
If two nodes u and v are in the same SCC, then there is a path from
every node in the component to every other node in the component.
In particular, there is a path from u to v and a path from v to u,
which form a cycle together.
Sufficiency:
If two nodes lie on a cycle, then there is a path from every one
of these two to the other one. As a result, they form a strongly
connected component.
Topologically sorted SCCalgorithm: A modified SCCalgorithm, which
uses numbers as component representatives. The numbers are assigned
in such a way that if each SCC is collapsed to one node, the numbers
indicate a topological sorted sequence of the nodes. The backtrack
procedure of the SCCalgorithm needs to be replaced as follows:
CLASS TopSortSCC EXTENDS SCC
component_counter := 0
PROCEDURE backtrack( (u,v) : Edge)
super.backtrack((u,v))
IF v == openComponents.top() THEN
openComponents.pop()
REPEAT
u:=openNodes.pop()
component[u] := component_counter
UNTIL u == v
component_counter := component_counter + 1
This procedure assigns a component number whenever a component is
closed. Since the component numbers are assigned in the wrong order,
they need to be reversed after the DFS finished:
PROCEDURE renumber()
FOREACH v in V DO
component[v] := component_counter  component[v]
Twoedgeconnected nodes in an undirected graph: Two nodes, which
lie on a cycle.
Twoedgeconnected component of an undirected graph: A set of twoedge
connected nodes. This means: If one cuts one edge of a twoedge
connected component, all of its nodes still remain connected.
Since SCCs in undirected graphs are twoedgeconnected components, the
SCC algorithm can be used to find maximal twoedgeconnected
components.
Biconnected nodes in an undirected graph: Two nodes, which lie on a
simple cycle. Biconnectedness implies twoedgeconnectedness.
Biconnected component of an undirected graph: A set of edges of the
graph, which lie on a simple cycle. The nodes of any biconnected
component are biconnected and also twoedgeconnected. That is:
If one takes out one of the nodes of a component, the others
still remain connected. Biconnected components can be found by
plugging the following DFSAction into the DFSNumbering algorithm:
CLASS Biconnected EXTENDS DFSNumbering
openComponents : Stack of Edge := nil
openEdges : Stack of Edge := nil
components : Array[m] of Edge :=
PROCEDURE traverse((u,v) : Edge)
super.traverse((u,v))
openEdges.push({u,v})
IF {u,v} is a tree edge THEN
openComponents.push({u,v})
ELSE IF {u,v} is a back edge THEN
WHILE dfsNum[v] < dfsNum[openComponents.top()[1]] AND
dfsNum[v] < dfsNum[openComponents.top()[2]] DO
openComponents.pop()
PROCEDURE backtrack((u,v) : Edge)
super.backtrack((u,v))
IF {u,v} == openComponent.top() THEN
openComponents.pop()
REPEAT
e := openEdges.pop()
component[number(e)] := {u,v}
UNTIL e == {u,v}
This algorithm works very similar to the SCCalgorithm: Both depend on
cycles: Each tree edge creates a new component. Whenever a backedge (a
cycle) is found, the involved items are merged to one component.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Shortest paths
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SSSP Problem, Single Source Shortest Path Problem: The problem of,
given a start node in a graph, finding all shortest paths from this
start node to all other nodes.
SSSPTreeLemma: The fact that, if there is a graph and a start node s,
such that all nodes are reachable from s and such that the graph does
not contain any negative cycles, then there is a spanning tree, all
paths of which are shortest paths.
Proof: Suppose we mark all edges in the graph, which contribute to
a shortest path from s to any other node. This shortest path
always exists for each node because each node can be reached from s.
So these edges will reach all of the nodes of the graph. Now suppose
these edges contain a cycle s ~~> u ~~> v > u. Then v > u cannot
contribute to any shortest path because u can be reached from s by
the much shorter path s~~>u.
Hence, the solution to a SSSP problem can be given by a spanning tree.
This spanning tree can be given by determining for each node a parent
node.
Negative cycle in a weighted graph: A cycle, in which the sum of the
edge costs is negative.
Shortest path length for two nodes u and v: The value
* my(u,v) := +oo if there is no path from u to v
* my(u,v) := oo if there is a path from u to v which contains
a negative cycle
* my(u,v) := min({c(p) : p is a simple path from u to v})
else
It makes sense to assign my(u,v)=oo for negative cycles, because
if u~~>w~~>w~~>v is a path with the negative cycle w~~>w, then
this cycle can be passed infinitely often, lowering the path costs to
oo. The value my(u,v) will be used to argue about algorithms for the
SSSP problem, but it will not be used in the algorithms themselves,
because the value of my(u,v) is usually unknown.
Minimal distance of a node v: The length of the shortest path from a
given start node s to v, i.e. my(s,v).
Generic SSSP algorithm: The following blueprint for an algorithm which
solves a SSSP problem:
FUNCTION genericSSSP((V,E,c) : WeightedGraph, s : Node) :
(Array of Real,Array of Node)
d : Array[V] of Real := <+oo,...> // d[i] = Tentative distance
from Node s to Node i,
this distance is always
greater than the length
of the shortest path from
s to i
d[s] := 0
parent : Array[V] of Node := // Stores for each node a
pointer to the parent
WHILE Ex v in V : d(v) != my(s,v) DO // This condition is a
theoretic one, it cannot
yet be checked
Choose e in E
relax(e)
RETURN (d, parent)
FUNCTION relax((u,v) : Edge) : boolean
IF d(u)+c(u,v) < d(v) THEN
// We can reach u with cost d(u), we can pass
// (u,v) with cost c(u,v), hence v can be reached with
// cost d(u)+c(u,v). Store u as a new parent of v
parent[v] := u
d[v] := d(u)+c(u,v)
RETURN true
RETURN false
This algorithm uses an array d of tentative distances. d[i]=x means
that a path of length x from s to node i has been found. The algorithm
continuously "relaxes edges" until every node has its optimal
distance. Relaxing an edge means: Check whether the edge is useful for
reaching the end node of the edge. This is the case if the distance
of the start node plus the cost of the edge is smaller than the
tentative distance of the end node. Then, the start node is saved as
a predecessor of the end node and d[i] is updated.
RelaxationLemma: The fact that in the genericSSSP procedure, the
tentative distance d[i] of a node i is its minimal distance, if the
edges of the shortest path from the start node to node i are
relaxed in the order given by the path. This lemma is quite
intuitive: We pick the shortest path from the start node to node i.
Then we follow this path and sum up the costs of the edges. This
is what the algorithm does, it always adds the cost of edge (u,v)
to the edge sum calculated so far in d[u] and stores the result in
d[v]. Then d[i] will contain the sum of all the edges of the
shortest path, and hence nothing else than the minimal distance of i.
ParentCycleLemma: The fact that, if a node v lies on a cycle of the
parentpointers of the genericSSSP procedure, then the minimal
distance of v is oo. Proof: Suppose the cycle of parentpointers
is s<~~u<~~v<~~u. Then these parent pointers have been installed
because c(s~~>u~~>v~~>u) < c(s~~>u), i.e. u can be reached with
less costs if one travels from s to u, from u to v and from v to u.
Now it is
c(s~~>u~~>v~~>u) < c(s~~>u)
<=> c(s~~>u) + c(u~~>v~~>u) < c(s~~>u)
<=> c(u~~>v~~>u) < 0
So u~~>v~~>u is a negative cycle and my(s,v)=oo.
BellmanFordAlgorithm: The genericSSSP procedure, where the main
loop has been replaced by:
FOR i := 1 to V1 DO
FOREACH e in E DO relax(e)
// Check for edges which can still be relaxed
relaxable : Set of Node := {}
FOREACH (u,v) in E DO
IF relax(u,v) THEN
d[v]:=oo
relaxable := relaxable \/ {v}
// Spread the distance oo to all nodes reachable from
WHILE relaxable != {} DO
u := pickAny(relaxable)
FOREACH (u,v) in E with d[v]!=oo DO
d[v]:=oo
relaxable := relaxable \/ {v}
This algorithm relaxes V1 times all edges of the graph. This
suffices to find the shortest path to each node: The shortest path to
any node i has at most V1 elements. In the first loop execution,
its first edge will be relaxed (among others). In the second loop
execution, the first edge will not change any more, but the second
edge will be relaxed (among others). By help of the relaxation lemma,
the shortest path is found.
After having done the relaxation, those edges which could still be
relaxed once more deserve attention: They lie on a negative cycle.
This algorithm works for arbitrary edge weights. It could be
improved by relaxing only those edges, the start node of which
received a different d[i]. It runs in time O(V*E).
SometoSomeProblem: The problem of, given a weighted graph (V,E,c)
and two subsets S < V and T < V, finding a shortest path with the
property that it starts in a node in S and ends in a node in T. This
problem can be solved by adding an artificial source node s, which
is connected by a 0weigh edge to all nodes in S. Furthermore, an
artificial node t is added and all nodes in T are connected to it.
Then, BellmanFord can be run on ((V,E,c),s). The shortest path from
s to t (as found by BellmannFord) is a path with the desired property.
Arbitrageproblem: The problem of, given n currencies with their
exchangerates, finding a tour of currencies, such that exchanging an
amount of money to the first currency, then the result to the second
etc. and finally back to the initial currency will result in an
increase of money. This problem can be interpreted as a graph problem:
Each currency is a node and each possibility of exchange between two
currencies is an edge. The edges are weighted by the exchange rate.
The question is: Does there exist a cycle in the graph, such that the
product of the edge weights is greater than 1? Since the known graph
algorithms are designed to minimize the sum of the edge weights
instead of maximizing their product, we use the following fact:
x[1] * x[2] * ... * x[k] > 1
<=> 1/x[1]*1/x[2]* ... *1/x[k] < 1
<=> log(1/x[1]*1/x[2]* ... *1/x[k]) < log(1)
<=> log(1/x[1]) + ... log(1/x[k]) < 0
If we replace the edge weights by log(1/x[i]), then the problem
reduces to finding a negative cycle. BellmanFord can be used to
identify such negative cycles.
Solution to the SSSP problem on DAGs: The genericSSSP procedure,
where the main loop has been replaced by:
topologicalOrder := topologicalSort(V,E,s)
FOREACH v in topologicalOrder DO
FOREACH (v,w) DO relax(v,w)
If the nodes are in topological order, it follows from the
relaxation lemma, that all shortest paths are found. The runtime is
O(V+E).
Dijkstra's Algorithm: The SSSP algorithm, where the main loop has
been replaced by:
Q : PriorityQueue of Node with Priority d
Q.insert(s)
WHILE Q != nil DO
u := Q.deleteMin()
FOREACH (u,v) in E DO
IF relax(u,v)==true THEN
IF v in Q THEN Q.decreaseKey(v)
ELSE Q.insert(v)
This algorithm is a kind of greedy BFS: It applies the relaxfunction
first to those nodes, which promise to be part of a short path, i.e.
those with a small d[i]. Relax candidates wait in a priority queue, in
which those with a small d[i] are in front. The smallest node is taken
and its outgoing edges are relaxed. In doing so, the end nodes of these
relaxed edges also become relax candidates. They are put into the queue
according to their temptative distance. If a node becomes a relax
candidate, when it is already a relax candidate, its position
in the queue is updated according to the new d[i] (Q.decreaseKey).
The first node of the queue does already have its minimal distance
d[i]==my(s,i). Proof: Suppose there was a shorter path from s to i.
Then some nodes on this path will already have been explored by the
algorithm. Be j the last node on this path explored by the algorithm.
Since this path is shorter than the one found by the algorithm,
d[j] v[1] 0> v[2] 0> v[3] ...
...where edges (v[i],v[j]) have to be added for all j>i. These
edges have weight c(v[i],v[j])=ji. If Dijkstra starts exploring this
graph from v[0], it will explore all edges (v[0],v[i]) and add all
nodes v[i] to the queue, each with an temptative distance of
d[v[i]]=i. Then it follows the smallest edge to v[1] and is forced to
correct all temptative distances d[v[i]] by 1 etc.
Limousineproblem: The problem of, given a map of a city with only
eastwest or northsouth streets, finding a cheap route from one
street intersection to another, where the cost of a route is defined
as the the length of the route plus the number of turns times a
constant C. That is: Find a short path from one point to another,
avoiding turns. The problem can be interpreted as a graph problem,
with the street intersections being the nodes and the street parts
the edges. The edge weights are the lengths of the street parts.
In order to weight the turns, we split up the graph into two
parts: The graph of all eastweststreets and the graph of all
northsouthstreets. For each intersection of a northsouthstreet
with a eastweststreet, we add an edge of weight C from the one
graph to the other. Now, shortest paths can be calculated e.g. by
Dijkstra and their cost will automatically include the cost for
turns.
SSSPinSCCproblem: The problem of solving the SSSP in a graph (V,E)
with SCCs of at most a constant size C. This problem can be solved
in time O(V+E): First, all SCCs need to be calculated
(O(V+E)). Next, the shortest paths from all nodes within the SCC
to all other nodes within the SCC are calculated. Due to the constant
number of nodes in the SCC, this takes constant time. Then, the SCCs
can be "unwringled": We make a copy A' of all nodes A in a SCC, and
delete all edges in the SCC. We only keep the edges entering the
SCC (connecting them to their node in A) and those leaving the SCC
(connecting them from their node in A'). Now, we add edges between A
and A', which mirror the length of the shortest path from each node in
A to each node in A'. Thereby, the whole graph becomes acyclic. Now,
the SSSPproblem can be solved in time O(V+E).
Leftrightminimum of a sequence s of type Real: A value, which is
smaller than all values with a lower index.
Example: s=<5,3,8,4,2,7>
^ ^ < leftrightminima
Probability of a value of being a leftrightminimum: The value
P(x[i]==min({x[1],x[2],...x[i]})) = 1/i
where x is the sequence and i is the position of the value,
counting 1based.
Expectation of the number of leftrightminima: The value
P(s[1] is LRmin or s[2] is LRmin or ...)
= 1/1 + 1/2 + 1/3 + ... 1/s.length
= H(s.length) < ln(s.length)+1
where s is the sequence of type Real.
Expectation of the number of decreaseKey operations in Dijkstra's
Algorithm on a graph (V,E): A value smaller than V*ln(E/V).
Proof: Consider a node v with all its entering edges e[i]. We
give numbers to the edges in the order in which Dijkstra explores
them: e[1] is the edge explored first etc. Let d[i] be the
tentative distance of the start node of e[i]. Then a
decreaseKeyoperation takes place for every leftrightminimum
of the sequence
The number of leftright minima is H(k). The first node is always a
leftrightminimum, but causes it an insertoperation instead of a
decreaseKeyoperation. hnce we have to count H(k)1 decreaseKey
operations per node.The expectation for the number k of edges of
a node is E/V. Since we have V nodes, the expectation of
the number of decreaseKeyoperations is
V*(H(E/V)1) < V*ln(E/V)
Discrete function: A function, the range of which is finite.
Discrete priority queue: A priority queue, the priority function of
which is discrete. Usually, the priority function maps to a finite
subset of the natural numbers.
Monotone priority queue: A priority queue, in which there is never an
element inserted before the element, which has been deleted most
recently. That is: "Minima never decrease", the first element of
the queue will always have a lower priority value than any element
which will take up this first place later. The priority queue of
Dijkstra's algorithm is monotone, because the priority values d[i]
increase with the distance from the start node.
"Range C restricted priority queue": A priority queue, where the
difference of the smallest and the largest priority value never
exceeds the value C. The priority queue of Dijkstra is restricted
to a range of the maximum edge weight: Whenever the first node i
is taken out of the queue in order to relax its edges, the new nodes
resulting from this relaxation will have a priority value of
d[i]+c(e), which is smaller than d[i]+max({c(e) : e in E}). Thus,
the maximum difference of the smallest and the greatest node is
the maximum edge weight.
// I don't know if this has a name, but I feel we need it for
// the bucket queue.
Bucket queue of type t with priority f: A discrete, monotone and
range C restricted priority queue of type t and priority f, which
is defined by
CLASS BucketQueue of t
B : array[C+1] of List of t :=
minIndex : Integer := 0
PROCEDURE insert(x : t)
B[f(x) mod (C+1)].pushFront(x)
FUNCTION deleteMin() : t
minIndexBound := minIndex + C + 1
WHILE B[minIndex mod (C+1)]==nil DO
IF minIndex > minIndexBound THEN RETURN nil
minIndex := minIndex + 1
RETURN B[minIndex mod (C+1)].popFront()
PROCEDURE decreaseKey(x : t)
insert(x)
The "" can be done efficiently, if each element stores a
reference to its list entry in the bucket queue. The bucket queue
keeps C+1 "buckets" (lists) of elements. minIndex points to the
bucket, which contains the element with the lowest priority value.
The bucket right to B[minIndex] contains the elements, the priority
value of which is the minimum priority value + 1 etc. If the end
of the array is reached, the sequence wraps to the beginning, i.e.
an element x can always be found at position f(x) mod (C+1). If an
element is inserted, it is added to the list of its corresponding
bucket. If the minimum is to be deleted, minIndex seeks to the right
until a new nonempty bucket has been found, which contains the
actual lowest priority value element. This element is removed and
returned. Since the queue is monotone, new elements will always be
inserted to the right of minIndex (modulo C+1). Since the queue is
range restricted, the greatest possible difference between the lowest
priority value and the largest priority value is C. Hence, a new
element with a very large priority value will never overrun the
invisible limit of the minIndexbucket, since its distance to
minIndex is C and the queue has a length of C+1. If m is the number of
elements in the bucket queue, then it needs a space of O(m+C) and a
time of O(C) for insert and decreaseKey. If used in Dijkstra, the
time for insert and decreaseKey can be reduced further as follows: Be
m the weight of the lightest edge
m := min({c(e)  e in E).
Determine the ratio of the largest edge weight to the smallest edge
weight
r := max({c(e)  e in E}) / m
Now, we use r+1 buckets of "width" m, that is the first bucket takes
nodes v with 0 < d[v] < m, the second bucket takes nodes v with
m =< d[v] < 2*m etc. Then the insertprocedure looks as follows:
PROCEDURE insert(x : Node)
B[d[x]/m mod (r+1)].pushFront(x)
...and the other operations are adjusted appropriatly. This algorithm
works because the d[i] have at least a distance of m among each other,
since they are sums of edge weights and the minimum edge weight is m.
msd, most significant distinguishing bit of two integers x and y: The
0based index of the leftmost bit of the binary representation of x
and y, for which x and y are different. Example:
msd(53,59):
53 = 00110101
59 = 00111011
^ msd, third bit from right
The msd can be calculated as follows:
msd(x,y) := (x==y) ? 1 : round_down(log2(x XOR y))
If the numbers are in a constant range, then log2(x) can be calculated
in O(1). The msd can be seen as a logarithmic distance.
Radix heap of type t with priority f: A discrete, monotone and
range C restricted priority queue of type t and priority f, which
is defined by
CLASS RadixHeap of t and priority f
K : integer := 1 + round_down(log(C))
B : array[1,...K] of list of t
minVal : integer := 0
PROCEDURE insert(x : t)
B[min(K,msd(f(x),minVal))].pushFront(x)
PROCEDURE decreaseKey(x : t)
insert(x)
FUNCTION deleteMin() : t
i := min({j : j in {1,...K} and b[j] != nil })
IF i==1 THEN RETURN B[1].popFront()
minVal := min({ f(x) : x in B[i]})
minElement := eps({x f(x)==minVal})
B[i].delete(minElement)
FOREACH v in B[i] DO
decreaseKey(v)
RETURN minElement
The radix heap works similar to a bucket queue, but the buckets are
"wider" for areas far away from the minimum. Thus, two elements x1 and
x2 may fall into two separate buckets, if their priority is close to
the element with the lowest priority value. But they may fall into the
same bucket, if their priority is far away from the minimum priority.
A new element is inserted in function of the logarithmic distance of
its priority value to the lowest priority value. A priority value is
changed (decreaseKey) by deleting the element and reinserting it.
Deleting it may be implemented efficiently if the element stores a
reference to its position. Bucket 1 stores elements, the priorities
of which are equal to the lowest priority value (see definition of
msd). If the minimum is to be deleted from the queue, the first
nonempty list is seeked and its minimum is returned. Then, the
queue has to be rearranged: Each element of the minimum bucket is
deleted and reinserted. Thereby, the old minimum value (the one just
removed) is taken as a measure for msd. All elements to the right
of bucket B[i] are still in the right order. Proof: Be oldMin the
old minimal priority value, according to which the elements are spread
over the buckets. Be x the element in bucket B[i] to be removed. Be y
an element in the bucket B[i+k]. Then we have to prove, that y still
remains in B[i+k], even if x is taken as a measure for msd:
i+k = i+k
log2(yoldmin) = log2(xoldmin) + k
log2(yoldmin)log2(xoldmin) = k
log2((yoldmin)/(xoldmin)) = k
(yoldmin)/(xoldmin) = 2^k
yoldmin = 2^k*(xoldmin)
y = 2^k*(xoldmin) + oldmin
Now the new position of y is
log2(yx)
= log2(2^k*(xoldmin)+oldminx)
= log2((2^k1)*(xoldmin))
= log2(xoldmin) + log2(2^k1)
~= i + k
High Part of an integer x with a bitbreadth K: The value
highPart(x,K) := x div 2^(K/2)
This is the value, which can be obtained by cutting away the right
half of the binary representation of x.
Low Part of an integer x with a bitbreadth K: The value
lowPart(x,K) := x mod 2^(K/2)
This is the value, which can be obtained by cutting away the left
half of the binary representation of x.
Van Emde Boas Search Tree, VEBTree of bitbreadth K: A type, which
allows the storing and retrieval of natural numbers. The log2 of
the numbers may not exceed K. K must be a power of 2. The VEBTree
is either
* the value nil OR // VEBtree is empty
* an integer value OR // VEBtree just stores 1 value
CLASS VEB
minM : integer // stores minimum value of values contained
maxM : integer // stores maximum value
highBits : VEBTree of bitbreadth K/2
lowBits : Array[2^(K/2)] of VEBTree of bitbreadth K/2
FUNCTION minVal() : integer
// Accessor function for minimal value
IF this==nil THEN RETURN +oo
IF this== e THEN RETURN e
RETURN minM
FUNCTION maxVal() : integer
// Accessor function for maximal value
IF this==nil THEN RETURN oo
IF this== e THEN RETURN e
RETURN maxM
PROCEDURE insert(x : integer)
IF this==nil THEN
// make trivial VEBtree of 1 element
this:=x
RETURN
IF this== e THEN
// convert trivial VEBtree to complex one
minM := e
maxM := e
highBits.insert(highPart(e,K))
lowBits[highPart(e,K)].insert(lowPart(e,K))
// Insert x
IF lowBits[highPart(x,K)] is empty THEN
highBits.insert(highPart(x,K))
lowBits[highPart(x,K)].insert(lowPart(x,K))
IF xmaxM THEN maxM := x
PROCEDURE remove(x : integer)
IF this==x THEN
this:=nil
RETURN
subtree := lowBits[highPart(x,K)]
subtree.remove(lowPart(x,K))
IF subtree is empty THEN
highBits.remove(highPart(x,K))
// Did we destroy our own maximum or minimum?
IF x==minM THEN minM := lowBits[highBits.minVal()].minVal()
IF x==maxM THEN maxM := lowBits[highBits.maxVal()].maxVal()
// Check out whether we are trivial
IF highBits== x1 AND
lowBits[x1]== x2 THEN
this := x1*2^K+x2
PROCEDURE locate(x : integer) : {0,...2^K1}
// Returns x if x is in the tree
// Otherwise, return next greater element which is in the tree
// If there is none, return +oo
IF this=={} THEN RETURN +oo
IF this== x THEN RETURN x
IF x>maxVal() THEN RETURN +oo
IF x
F : Set of Node := {} // This set will keep the nodes, which already
have their final distance
Q.insert(s)
WHILE Q != nil DO
WHILE F != {} DO
u := pickAny(F)
FOREACH (u,v) in E DO
IF relax(u,v)==true THEN
IF v in Q THEN Q.decreaseKey(v)
ELSE Q.insert(v)
F := F \/ deleteSettled(Q)
FUNCTION deleteSettled(Q : RadixHeap of type t with priority f)
: Set of t
// Deletes and returns all elements of Q having a minimal
// priority value
F' : Set of t := {}
i := min({j : j in {1,...K} and Q.B[j] != nil })
IF i==1 THEN
F' := { x : x in Q.B[1] }
Q.B[1] := nil
RETURN F'
minVal := min({ f(x) : x in B[i]})
// Run through all nodes in B[i], select the one with d[v]=minVal
// as well as all those nodes which also already have their final
// distance
FOREACH v in Q.B[i] DO
IF f(v) < minVal + minIn(v) THEN
F' := F' \/ {v}
ELSE
decreaseKey(v)
RETURN F'
This variant of Dijkstra always reloads a couple of nodes from the
queue, which already have their final distance. Then the algorithm
works on these interesting nodes (kept in F) and afterwards refills
F. The algorithm knows that the minimal node v in the queue already
has its final distance. If there is another node w after v in the
queue, then the algorithm looks at the ingoing edges of w: If the
smallest of them comes from a node u, which is before v (i.e.
d[u]=d[w]smallestEdge < d[v]), then this node u already had its
final distance, hence also w must have its final distance. Thus,
w can be taken into F. By this "caching", the algorithm avoids
unnecessary queuereorderingoperations and runs in time O(m+n).
Proof: We have to count the expected number of moves of one node
in the queue. Be C the maximum edge weight, be K the length of the
bucket queue, i.e. K:=1+log(C). Any node v will only move down
to minIn(v), because then it is caught and put to F. In the radix
heap, this means that the buckets the node will pass on its way from
bucket K+1 to its exit position is K + 1  log(minIn(v)). Assume that
log(c(e)) is always rounded down:
E(#moves) =< SUM {E(K + 1  log(minIn(v)))  v in V }
=< n + SUM {E(K  log(c(e)))  e in Edges}
Let's look at one edge e. The term E(K  log(c(e))) is
P(log(c(e))==K1)*1
+ P(log(c(e))==K2)*2
...
+ P(log(c(e))==0)*K
= P(c(e) in [2^(K1),2^K1] )*1
+ P(c(e) in [2^(K2),2^(K1)1] )*2
...
+ P(c(e) in [0,1] )*(K1)
= 1/2 * 1
+ 1/4 * 2
...
+ 1/2^K * K
=< 2
Thus, the big term reduces to
=< n + m*2
// This proof has been invented by Thomas Schultz, thanx!
APSP Problem, All Pairs Shortest Paths Problem: The problem of, given
a graph, finding the shortest path from each node to each other node.
This problem could be solved by solving the SSSP problem for every
node. In the case of negative edgeweights, BellmanFord is required.
If one would like to use the faster Dijkstra, one has to change
the negative edges first. This change may not destroy eventual
shortest paths. This is done by nodepotentials.
Nodepotential in a weighted graph (V,E,c): A function pot:V>Real, such
that for all pairs of nodes u,v
pot(u)+c(u,v)pot(v) >= 0
This corresponds to "lifting up" all edge weights in order to achive
nonnegative edge weights.
Reduced cost in a weighed graph (V,E,c) with a nodepotential pot:
The function tau:E>Real, such that
tau(u,v) = pot(u)+c(u,v)pot(v)
This is a new cost function, which follows the "lifting up" of the
nodepotential.
Reduced cost of a path p in a weightd graph with nodepotntial: The
sum of the reduced costs of all edges of p, noted tau(p).
Reduced Cost Lemma: The fact that shortest paths in a weighted graph
(V,E,c) are the same as in the weighted graph (V,E,tau), where
tau is the corresponding reduced cost function. Proof:
tau(p) = tau(p[1],p[2]) + tau(p[2],p[3]) + ...
= pot(p[1])+c(p[1],p[2])pot([p2]) +
pot(p[2])+c(p[2],p[3])pot([p3]) +
...
= pot(p[1])+ c(p) pot(p[p.length1])
Since for a shortest path between two nodes, pot(p[1]) and
pot(p[p.length1]) are constant, shortest paths with respect to tau
are exactly the shortest paths with respect to c.
Reduced Cost Node Potential Lemma: The following fact: If a
weighted graph (V,E,c) does not contain any negative cycles and
if s is a node, from which all other nodes can be reached and
the node potential is chosen
pot(v) := my(s,v),
then tau(e)>0 for all edges e. Proof:
For any edge (u,v), it holds that
my(s,u)+c(u,v) >= my(s,v)
<=> my(s,u)+c(u,v)my(s,v) >= 0
<=> pot(u) +c(u,v)pot(v) >= 0
<=> tau(u,v) >= 0
APSP finder: The following algorithm, which solves the APSP problem:
FUNCTION APSPfinder((V,E,c) : WeightedGraph) : Array of Array of
Node
// Returns a 2dimensional array, in which at position i,j
// there is the node to go to if one wants to go from
// node i to node j
Result : Array[V][V] of Node
// Introduce a new node s
V := V \/ {s}
// Introduce edges from s to all nodes
E := E \/ {(s,v)  v in V\{s}}
// These edges have weight 0
c := c \/ {((s,v),0)  (s,v) in E}
// Run BellyFord to get minimal distances from s
// to all nodes
d := BellmanFord((V,E,c),s).d
// Introduce the node potential to be the minimal
// distance to s
pot := {(v,d[v])  v in V}
// Define tau as usual
tau := {((u,v),pot(u)+c(u,v)pot(v))  (u,v) in E}
// Run Dijkstra for each node
FOREACH v in V DO
Result[v] := Dijkstra((V,E,tau),v).parent
RETURN Result
This procedure first finds some node potential function. Then
it defines a reduced cost function based on this node potential.
Now, Dijkstra can run, because no negative edges are left.
Dijkstra runs for every node, thus producing a total running time
of O(n*m + n*(m+n*log(n))) = O(n^2*log(n)).
Heuristic in a weighted graph (V,E,c) for a target node t: A function
f:V>Real, such that f(v)==0, but this is not mentioned in the slides
// (?)
A*: The following algorithm, which finds a shortest path from one
node to another in a weighted graph:
FUNCTION A*((V,E,c) : WeightedGraph, s : Node, t: Node,
f : Heuristic for (V,E,c) and t) : Sequence of Node
d : Array[V] of Real := <+oo,...>
d[s] := 0
parent : Array[V] of Node :=
Q : PriorityQueue of Node with Priority d(v)+f(v) := <>
Q.insert(s)
WHILE !Q.empty() DO
u := Q.deleteMin()
FOREACH (u,v) in E DO
IF relax(u,v)==true THEN
IF v in Q THEN Q.decreaseKey(v)
ELSE Q.insert(v)
RETURN
FUNCTION relax((u,v) : Edge) : boolean
IF d(u)+c(u,v) < d(v) THEN
parent[v] := u
d[v] := d(u)+c(u,v)
RETURN true
RETURN false
This algorithm picks those nodes first, which not only have a close
distance to s, but also promise to have a short distance to t. This
algorithm is equivalent to Dijkstra with a node potential of
pot(v) = f(v). It will only work if tau(e)>=0. This is the case
if
tau(u,v) >= 0
<=> pot(u)+c(u,v)pot(v) >= 0
<=> f(u,t)+c(u,v)+f(v) >= 0
<=> c(u,v)+f(v) >= f(u)
This is a kind of triangleinequality.
// For another A* approach, see "Intelligence artificielle"
// (ia.txt, French)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Maximum Flows
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Capacity graph: The type
CLASS capacityGraph EXTENDS Graph
capacity : Node > Real+
sourceNode : Node
targetNode : Node
This is a weighted connected graph, the cost function of which is
called "capacity". The capacity must be positive. Such a graph can
be thought of as a network of tubes, where the diameter of the tubes
is given by the capacity of the edges.
Common Notation: (V,E,cap,s,t)
Flow in a capacity graph (V,E,cap,s,t): A function f:E>Real, such that
* f(e) >= 0 for all e in E
* f(e) =< cap(e) for all e in E
* SUM {f(u,v)  (u,v) in E } = SUM {f(v,w)  (v,w) in E}
for all v in V\{s,t}
The flow can be thought of as water which flows through the tubes.
cap(e) gives the maximum number of liters per second, whereas f(e)
gives the actual number of liters per second. The last condition
states that the amount of water flowing into a node must also
leave it again. This does not apply to s and t.
Flow graph: A capacity graph with a flow. A flow graph is usually
drawn like a normal graph with the notion "cap(e)/f(e)" with each
edge.
Excess of a node v in a flow graph (V,E,cap,s,t,f): The difference of
the flow of its incoming egdes and the flow of its outgoing edges.
excess(v) = SUM {f(u,v)  (u,v) in E }  SUM {f(v,w)  (v,w) in E}
By definition, the excess of all nodes except for s and t is 0. The
excess can be thought of the amount of water running out of a node.
("Ueberlauf").
Value of a flow in a capacity graph (V,E,cap,s,t): The excess of t.
Notation: val(f) = excess(t)
Maximal flow in a capacity graph: A flow, such that there is no flow
with a greater value. The goal of this section is finding a maximal
flow.
stExcessLemma: The fact that the excess of the source node in a
capacity graph is the negative excess of the target node.
excess(s) = excess(t) in (V,E,cap,s,t,f)
Proof:
excess(s) + excess(t) = SUM {excess(v)  v in V }
because excess(v)=0 for every v with v!=s, v!=t.
SUM {excess(v)  v in V } = 0
because every f(e) appears twice in the SUM: Once as an input
to a node and then (with inverse sign) as its output. The only
values which do not go away are excess(s) and excess(t)
=> excess(s) + excess(t) = 0
=> excess(s) = excess(t)
Saturated edge in a flow graph (V,E,cap,s,t,f): An edge the capacity
of which is its flow. f(e) = cap(e)
Cut in a graph (V,E): A subset of the nodes V.
Edge leaving a cut S: An edge, the start node of which is in S and the
end node of which is outside S.
Edge entering a cut S: An edge, the end node of which is in S and the
start node of which is outside S.
Edge of a cut S: The set of its entering and leaving edges. A cut can
also be given by its edges, if the graph is connected.
stcut in a capacity graph (V,E,cap,s,t): A cut, which contains s and
does not contain t.
Capacity of a stcut S: The sum of the capacities of the edges
leaving S.
Notation: cap(S) = SUM {cap(u,v)  u in S, v in V\S}
with the capacity graph being (V,E,cap,s,t)
Saturated cut in a flow graph (V,E,cap,s,t,f): An stcut, such that
* all leaving edges are saturated, i.e. f(e)=cap(e)
* all entering edges have flow 0, i.e. f(e)=0
Cutcapacitylemma: The fact that the capacity of a cut is always
greater or equal to the value of a flow.
Proof:
val(f) = excess(t) = excess(s) = SUM {excess(u)  u in S }
= SUM {f(e)  e leaves S}  SUM {f(e)  e enters S}
= SUM {f(e)  e leaves S}  0
=< SUM {cap(e)  e leaves S}
= cap(S)
Saturatedcutlemma: The fact that if there exists a saturated cut in a
flow graph, then the value of the flow is the capacity of the cut.
Proof: See proof of Cutcapacitylemma: For a saturated cut, the "=<"
is a "=".
Saturatedcutmaximalflowlemma: The fact that if there exists a
saturated cut in a flow graph, then the flow is maximal. Proof:
Assume there was a flow with a greater value. Then this value would
contradict the cutcapacity lemma.
If a flow is maximal, then there exists a saturated cut, but there
may also be unsaturated cuts.
Residual network of a flow graph (V,E,cap,s,t,f): The capacity graph
(V,E',cap',s,t), where E' is
{e  e is an unsaturated edge in E } \/
{(v,u)  (u,v) is an edge with flow in E}
Thus, E' takes up all unsaturated edges of E and in addition takes up
the reverse edges of all flow carrying edges. For edges of the first
kind, one sets:
cap'(e) = cap(e)  f(e)
For edges of the second kind, one sets
cap'(v,u) = f(u,v)
The residual network displays where there can be additional flow in
the flow graph. This can happen in two kinds of edges: If an edge is
not saturated, flow can be added. If an edge already carries flow,
then this flow can be undone. This is mirrored by the edges of the
second kind.
Algorithm:
FUNCTION residualNetwork((V,E,cap,s,t,f) : FlowGraph) :
CapacityGraph
result : CapacityGraph := (V,E'={},cap'={},s,t)
FOREACH (u,v) in E DO
IF cap(u,v)>f(u,v) THEN
E' := E' \/ {(u,v)}
cap' := cap' \/ {((u,v),cap(u,v)f(u,v))}
IF f(u,v)>0 THEN
E' := E' \/ {(v,u)}
cap' := cap' \/ {((v,u),f(u,v))}
RETURN result
Minimal cut in a capacity graph (V,E,cap,s,t): A stcut with a minimal
capacity. That is: A subset of V, which
* contains s
* does not contain t
* has a minimal sum of the capacities of the leaving edges
Residualpathlemma: The fact that if there is a path from s to t in a
residual network of the flow graph (V,E,cap',s,t,f), then the flow
cannot be maximal. Proof: The residual network reflects flow that can
be added. If one adds flow along the path, then the overall flow
increases. Vice versa, if there is no such path, then no flow can be
added and the flow is maximal.
MaxFlowMinCutTheorem: The fact that if there is no path from s to
t in the residual network of a flow graph (V,E,cap,s,t,f), then the
set of all nodes reachable from s in the residual network constitutes
a minimal cut. Proof: Be S the stcut, which contains all nodes
which are reachable from s in the residual network Gf. We know that
all edges e leaving S in the original graph must have maximal
flow f(e)=cap(e), because these edges do not appear in Gf.
Furthermore, we know that all edges entering S in the original
graph must have flow 0, because their reverse edge does not appear
in Gf. As a result, S is saturated. It follows by the saturatedcut
maximalflowlemma that val(f)=cap(S). Hence, there cannot be a
stcut with a smaller capacity that cap(S), because val(f) has to
go through it. It follows that cap(S) is minimal, S is a minimal cut.
Capacity of a path in a capacity graph: The minimum of the capacities
of the edges in the path.
Augmenting path in a flow graph G: A path from the source node to the
target node in the residual network of G. This path shows that the
total flow can still be increased by increasing the flow along this
path.
MaxResCapPathLemma, Maximal Augmenting Path Lemma: The fact that
an augmenting path with a maximal capacity can be found in time
O(m*log(m)) in a residual network, where m is the number of edges.
Proof by algorithm:
FUNCTION findMaxAugPath((V,E,cap,s,t) : CapacityGraph) :
List of Node
m := E
// Sort edges by their capacity, big capacity first
:= sortBy(E,cap,DECREASING)
// Do a binary search to find the minimum set of edges,
// which still allows a path from s to t
lasti:= m
i := m/2
REPEAT
// Find a path with edges e[1],...e[i]
E' := {e[0], ...e[i]}
p := findPath((V,E'),s,t)
// If there exists a path, try smaller i
// else try greater i
tmp:=i
IF p==nil THEN
i:=i+abs(lastii)/2
ELSE
i:=iabs(lastii)/2
lasti:=tmp
UNTIL abs(lastii)=<1 && p!=nil  i>=m
RETURN p
This algorithm tries to build up a path of the most powerful
edges. To do so, it sorts these edges by their capacity and then
does a binary search to find out how many of these edges are needed
for a maximal augmenting path. The time is in O(n+m)*O(log(m)) =
O(m*log(m)+n*log(m)) = O(m*log(m)), since m>=n in a connected
capacity graph.
AugmentationsLemma: The fact that in an integer capacity graph with m
edges and a maximum flow value of maxVal, O(m+m*log(maxVal/m))
augmentations of an initial 0flow suffice to find the maximal flow.
Proof: (?)
FordFulkerson Algorithm: The following algorithm, which finds a
maximal flow in a capacity graph:
FUNCTION fordFulkerson((V,E,cap,s,t) : CapacityGraph) :
FUNCTION E>Real
n := V
m := E
f : FUNCTION E>Real := {(e,0)  e in E}
REPEAT
(V,E',cap',s,t) := residualNetwork(V,E,cap,s,t,f)
p := findPath((V,E),s,t)
IF p==nil THEN
// If t is not reachable from s in Gf, we have a
// maximal flow
RETURN f
// Find out how much we can increase
incVal := min({cap'(e)  e in p})
// Increase flow along p
FOREACH e in p DO
f(e) := f(e)+incVal
UNTIL 42==0
This algorithm bases on the residual network: If there is a path from
s to t in it, then there is a possibility to push flow from s to t. The
algorithm finds such paths and increases the flow along them. With
integer capacities, the algorithm has a running time of
O((m+m*log(maxVal/m))*m*log(m)), where maxVal is the value of the
maximum flow. Reason: According to the Maximal Augmenting Path Lemma,
a maximal augmenting path can be determined in time O(m*log(m)).
This procedure is necessary for every flow augmentation. According to
the AugmentationsLemma, O(m+m*log(maxVal/m)) of these augmentations
suffice. Thus, the claimed running time results. With noninteger
capacities, the algorithm may run forever.
Example for this (?)
BFSdepth of a node in a graph (V,E) relative to a start node s: A
natural number, which identifies the time when the node was found by a
BFS starting from s. Thus, s gets the number 0 and each "sphere"
discovered by BFS gets a new number. Algorithm:
FUNCTION BFSDepth((V,E) : Graph, s : Node) : Array of Integer
// Essentially do a BFS
q : Queue of Node :=
marked : Array[V] of Boolean :=
marked[s]:=true
depth : Array[V] of Integer := <0,...>
WHILE q != nil DO
u := q.popFront()
marked[u] := true
FOREACH (u,v) in E DO
IF marked[v]==false THEN
q.pushBack(v)
marked[v]:=true
// The depth of a node is the depth of its predecessor + 1
depth[v]:=depth[u]+1
RETURN depth
Reverse BFS Depth: A BFSDepthAlgorithm carried out on a graph, the
edges of which have been reversed.
Algorithm:
FUNCTION reverseBFSDepth((V,E) : Graph, s : Node) : Array of Integer
RETURN BFSDepth((V,{(v,u)  (u,v) in E),s)
Layered subgraph of a graph (V,E) relative to a start node s: A
graph (V,E'), where E' contains only those edges, which connect
a node of BFSdepth i with a node of BFSdepth i+1.
Algorithm:
FUNCTION layeredSubgraph((V,E) : Graph, s : Node) : Graph
depth := BFSDepth((V,E),s)
E' := {}
FOREACH (u,v) in E DO
IF depth[u]+1==depth[v] THEN E' := E' \/ {(u,v)}
RETURN (V,E')
This algorithm divides the graph into "spheres", which indicate the
distance of the nodes to the source node. Any edge, which goes within
one sphere or which goes backwards is cancelled. This algorithm does
not destroy shortest paths from s to any target node.
Blocking Flow in a capacity graph: A flow, which uses the capacities
in such a way, that no other flow may go from the source node to the
target node. That is: For every path from the source node to the
target node, there is a saturated edge in it.
Algorithm:
FUNCTION blockingFlow((V,E,cap,s,t) : CapacityGraph) :
FUNCTION E > Real
p : List of Nodes :=
f : FUNCTION E>Real := {(e,0)  e in E}
REPEAT
IF p.last()==t THEN
// We found a path from s to t, breakthrough
// Calculate maximal increase
incVal := min({cap(e)f(e)  e in p})
// Increase flow along path
FOREACH e in p DO
f(e) := f(e)+incVal
IF f(e)==cap(e) THEN E := E\{e}
// Reset path
p :=
ELSE
// Extend path if possible
IF Ex (p.last(),w) in E THEN p.pushBack(w)
ELSE
// If there is no more edge, we're done
IF p.last()==s THEN RETURN f
ELSE
// Retreat (Backtrack), delete the edge from
// the path and the set of edges
E := E\{ (p.get(p.length()2),p.last()) }
p.delete(p.last())
UNTIL 42==0
This algorithm works similar to a DFS. It tries to find a path from
s to t. Then it sends a flow through it and deletes the saturated
edges. Next, it finds another possible path. In doing so, it simply
destroys useless or saturated edges. As a result, no path from s to
t with an unsaturated edge is left. The running time results from
the number of operations:
* Number of breakthroughs =< m, because each breakthrough saturates
one edge
* Number of retreats =< m, since one edge is removed
* Number of extend operations =< Number of retreats + n*Number of
breakthroughs, because each extend
* is later undone by a retreat OR
* is part of a breakthroughpath
For the first case, only as many extends as retreats can appear.
In the second case, a breakthroughpath may involve at most n
extends. Thus, there may be all in all n times as many extends as
breakthroughs. Since each extend falls either in the first case or
the second, the claimed upper bound results.
Thus, the running time is in O(m+m+m+n*m)=O(n*m). If the capacities
of the edges are only 1 or 0, then the algorithm even runs in O(m+n),
since every breakthrough saturates all edges on the path.
Note: A blocking flow is not necessarily a maximal flow.
Dinic's Algorithm, Dinitz' Algorithm: The following algorithm, which
finds a maximal flow in a flow graph
FUNCTION dinic((V,E,cap,s,t) : CapacityGraph) : FUNCTION E>Real
f : FUNCTION E>Real := {(e,0)  e in E}
REPEAT
Gf := residualNetwork(V,E,cap,s,t,f)
Lf := layeredSubgraph(Gf,s)
// In the slides, a reverse BFS is done! (why(?))
fb := blockingFlow(Lf)
f := f + fb
UNTIL fb is 0function
RETURN f
This algorithm sees a maximal flow as a sum of blocking flows. The
blocking flows are calculated in the residual network, because this
graph reflects the edges which still allow flow. In calculating
the blocking flows, the algorithm does not care for the edges within
one "sphere" in the sense of the layered subgraph. These edges do
not contribute efficiently to a flow. As soon as the direct edges
are saturated, the residual network changes and the formerly
discarded edges come back into play. The loop is executed at most
n times (proof to follow(?)), and each loop requires the calculation
of the residual network (O(m)) and a blocking flow (O(n*m)),
resulting in a total time of O(m*n^2). If the capacities are either
0 or 1, then the algorithm runs in O((m+n)*sqrt(m)) (why(?)).
Bipartite graph: A graph (V,E), in which there is a subset V' of V,
such that there is no edge within V' and no edge within V\V'. This
graph reflects the idea of two different kinds of nodes: V' and
V\V'. Edges always connect nodes of one kind to nodes of the other
kind. Notation: (V' \/ V\V',E)
IntegerFlowTheorem: The fact that in a flow graph with integer
capacities, there is a maximal flow with integer flow values.
Proof: Use FordFulkerson to construct the maximal flow and
increase the flow by 1 in each iteration.
Matching in an undirected graph (V,E): A subset M of E, such that
(V,M) has maximum degree 1. This corresponds to a pairwise mapping of
nodes to other nodes.
Observations concerning flows:
* A maximal flow can contain cycles. Consider e.g.
,1000/1000> a
t <1/1 s <1000/1000'
* There always exists a cyclefree maximum flow, because the flow in
each cycle can be reduced by the minimum flow of its edges without
disturbing the maximal flow.
* If all edges in a capacity graph have unique capacities, this does
not mean that the maximal flow is unique. Consider e.g.
s 1/1> a 2/1> b 4/1> t
'3/0>
* If one replaces each edge by an undirected edge in a capacity graph,
then the maximal flow may change. Consider e.g.
s <1/0 t
* If all edge capacities in a capacity graph are multiplied by a
positive constant c, then the minimal cut remains unchanged. Proof:
Be cap'(e)=c*cap(e) the new capacity function. Consider the minimal
cut S and any other cut S'. It holds
cap(S) < cap(S')
c * cap(S) < cap(S') * c
c * SUM {cap(e)  e leaves S} < SUM {cap(e)  e leaves S'} * c
SUM {c*cap(e)  e leaves S} < SUM {c*cap(e)  e leaves S'}
SUM {cap'(e)  e leaves S} < SUM {cap'(e)  e leaves S'}
cap'(S) < cap'(S')
Thus, even in the changed graph, S has a smaller capacity than S'
and hence S is still the minimal cut.
* If all edge capacities in a capacity graph are increased by a
constant c, then the minimal cut may change. Consider e.g. a minimal
cut S with leaving edges e1 and e2, each with capacity
cap(e1)=cap(e2)=1. Now consider a nonminimal cut S' with the single
leaving edge e3 with cap(e3)=3. If we now add c=1000 to each
capacity, the new capacity of S is
cap(S) = cap(e1)+1000 + cap(e2)+1000 = 2002
whereas the capacity of the formerly nonminimal cut S' becomes
cap(S') = cap(e3)+1000 = 1003
* If all edge capacities in a capacity graph are increased by a
constant c, then the value of the maximal flow does not necessarily
change by a multiple of c. Consider e.g.
s 3/2> a 1/1> b 1/1> t
'1/1>
The flow is 2 here. If the capacities are increased by c=2, then
the two edges entering t triple their capacity to 6 all in all, but
the edge leaving s will only allow a flow of 5. Thus the flow will
increase from 2 to 5, which is not a multiple of c.
Spamming problem: See your mailbox.
In terms of graph theory: The problem of, given a network with one
spamming node q and a set S of peaceful users, determining the
connections where a spam filter needs to be installed, such that all
users are protected and a minimal amount of money is spent. Installing
a spam filter between the network nodes u and v costs c(u,v). The
communication network can be seen as a capacitygraph (V,E,cap) with q
in V, and S < V. We create an artificial target node t, which we
connect to all nodes of S with a capacity of +oo. The goal is to find
a minimal qtcut: If we block all edges leaving this cut, we get a
total isolation of q with a minimal effort. The minimal cut can be
found by finding a maximum flow, creating the residual network and
determining the nodes reachable from q. This cut does not contain any
node of S: The edges from S to t will certainly also appear in the
residual network, due to their high capacity. Thus, if any node in S
was reachable from q, then t would also be reachable from q and the
calculated flow cannot be the maximal flow.
Most vital arc, MVA of a capacity graph: An edge, the deletion of which
causes the largest decrease in the value of the maximal flow.
Some observations:
* A MVA does not need be an edge with a maximal capacity.Consider e.g.
a <100000/0 s 1/1> t
* A MVA does not need be the one with a maximal flow. Consider e.g.
,1000/1000> a
t <1/1 s <1000/1000'
* A MVA does not need to be leaving a minimal cut. Consider e.g.
s 1/1> b 10000/2> t
'1/1> a 1/1>
* A graph may contain several MVAs. Consider e.g.
s 1/1> a 1/1> t
Bipartite Matching Problem: The problem of, given a bipartite graph
(A \/ B,E), finding a matching. This corresponds to the idea of
creating couples (a,b) of nodes of A and nodes of B, such that a and b
are connected. The bipartite matching problem can be reduced to the
maximum flow problem: One adds a source node s and connects it
to all nodes of A. Then one adds a target node t and connects all
nodes of B to it. All edges are given a capacity of 1. Then the
problem is to find a maximal flow, in which all edge flows are
1 or 0: The 1edges are to keep and the 0edges are to throw away.
According to the IntegerFlowTheorem, such a flow does exist.
Dinner problem: The problem of, given n families with a[i] members and
m tables with b[j] chairs, finding a seating order such that no two
families sit together at one table. This problem can be solved
similarly to the Bipartite Matching Problem: One considers each family
a node and each table a node. One adds a source node s and edges with
capacity a[i] to every family i. One adds a target node and egdes with
capacity b[i] from every table to t. Then one adds edges with capacity
1 from every family to every table. Now, finding a maximal flow with
integer capacities in this graph will highlight some edges between the
families and the tables, indicating where the family members have to
go.
Preflow in a capacity graph (V,E,cap,s,t): A function f:E>R, such
that
* f(e) >= 0 for all e in E
* f(e) =< cap(e) for all e in E
* SUM {f(u,v)  (u,v) in E } >= SUM {f(v,w)  (v,w) in E}
for all v in V\{s}
The constraint of no excess is relaxed: A node may have positive
excess. Thus, a preflow can be seen as a trial of a flow. As a
preflow is very similar to a flow, a preflow is often treated like
a flow, especially concerning the existence of the residual network
and the notions of saturation and excess.
Preflow graph: A capacity graph with a preflow. Such a graph can be
imagined as a system of tubes with balloons in between. If there is
water preflow in a tube, the balloon at the end of this tube inflates.
It is under pressure and needs to be released.
Pushoperation in a preflow graph: The operation of increasing the
preflow across a given edge. As a precondition, the edge capacity must
allow this increase and the start node of the edge must have an excess
greater or equal to the value of the added preflow. The following
algorithm does this and adjusts the residual network:
FUNCTION push((u,v) : Edge, d : Real,
(V,E,cap,s,t, var f) : PreflowGraph,
(V,var E',s,t,var cap') : ResidualNetwork)
ASSERT d>0
ASSERT cap(u,v)f(u,v)>=d
ASSERT excess(u)>d
f(u,v) := f(u,v) + d
// update excess, if it is stored as an array
// excess(u):=excess(u)d
// update residual network (V,E',cap',s,t)
cap'(u,v):=cap'(u,v)d
// Eventually delete (u,v)
IF cap'(u,v)==0 THEN E':=E'\{(u,v)}
// Eventually add (v,u)
IF not (v,u) in E' THEN E':=E'\/{(v,u)}
// Set value of residual capacity
cap'(v,u) := f(u,v)
Saturating pushoperation: A pushoperation, which saturates an edge.
Levelfunction in a residual network (V,E,cap,s,t): A function
d:V>Natural, such that
* d(s)=V
* d(t)=0
* in the residual network of the flow graph, there is no edge (u,v)
such that d(u)=d(v)+1. Steep edges are not permitted by the
definition of d.
Downward edge in a residual network with a level function d: An edge
(u,v), where d(u)=d(v)+1.
Horizontal edge in a residual network with a level function d: An
edge (u,v), where d(u)=d(v).
Upward egde in a residual network with a level function d: An
edge (u,v), where d(u)=d(v)1.
Active node in a preflow graph: A graph with positive excess.
Eligible edge in a residual network with a level function d: An edge,
the start node of which has a higher level than the end node. Since
step edges are forbidden, this corresponds to a downward edge.
Relabeling a node v in a residual network with a level function d:
Increasing its value d(v). This corresponds to lifting the node in
the landscape of the graph.
Generic PreflowPushAlgorithm, PPAlgorithm: The following algorithm,
which finds a maximal flow in a capacity graph:
FUNCTION genericPreflowPush((V,E,cap,s,t) : CapacityGraph) :
FUNCTION E > Real
// Preflow
f : FUNCTION E>Real := {(e,0)  e in E}
n := V
m := E
// Node levels
d : FUNCTION V > N := {(v,0)  v in E}
d(t) := 0
d(s) := n
// Calculate residual network
(V,E',cap',s,t) := residualNetwork(V,E,cap,s,t,f)
// Saturate each source edge by a preflow
FOREACH (s,u) in E DO
push((s,u),cap(s,u) ,(V,E,cap,s,t,f),(V,E',cap',s,t))
// Select active node
WHILE Ex v in V\{s,t} : excess(v)>0 DO
// Select eligible edge
IF Ex (u,v) in E' : d(u)>d(v) THEN
// What flow can we push along it?
d := min(excess(u),cap'(u,v))
push(e,d ,(V,E,cap,s,t,f),(V,E',cap',s,t))
ELSE
// We have excess, but no downward edge > relabel u
d(u) := d(u) + 1
RETURN f
This algorithm first puts as much preflow as possible on each edge,
which leaves s. Then it selects a node with positive excess and
distributes the excess among the outgoing downward edges. If there
are no such downward edges, the node is risen, until downward edges
appear. Finally, the superfluous flow will flow back to s. The
algorithm leaves much space for heuristics: An active node has to
be selected in some reasonable manner and eligible edges have to be
selected.
Breaking the ties in a PPAlgorithm: Selecting the active nodes and
the eligible edges.
arbitrary preflowpushalgorithm, arbitrary PPAlgorithm: The
PPAlgorithm, where the active nodes and the eligible edges are
selected randomly.
SteepEdgesLemma: The fact that the PPAlgorithm does not produce
any steep edges. Proof:
If there is a node on a lower level than a node v, then v can push
its flow into this lower node until the edge is saturated.
Then, the edge disappears in the residual network.
ActiveNodePathLemma: The fact that from every active node in a
preflow graph, there is a path in the residual network to the source
node. This lemma guarantees that we can always get rid of superfluous
preflow by pushing the flow towards the source node. Proof:
Be G=(V,E,cap,s,t,f) the preflow graph.
Be Gf=(V,E',cap',s,t) the residual network of G.
Be v the node to be analyzed.
Be S the set of nodes reachable in Gf from v
S := {w  Ex path in Gf from v to w}
Be T the set of the remaining nodes
T := V\S
Then the total excess within S can be calculated as the difference of
the inflow and the outflow
SUM {excess(v)  v in S} = SUM {f(e)  e enters S} 
SUM {f(e)  e leaves S}
In Gf, there is no outgoing edge of S
~ (Ex (u,v) in E' : u in S && ~(v in S))
Hence, the flow coming into S must be 0:
All e in E : e enters S => f(e)=0
Thus, the sum of the excesses in S must be negative
SUM {excess(v)  v in S} = 0  SUM {f(e)  e leaves S} < 0
Since s is the only node, which is allowed to have negative excess,
s must be in S.
MaxLevelLemma: The fact that during the PPAlgorithm, no node attains
a level above or equal to 2*n, where n is the number of nodes. This
lemma works as an upper bound for the overall height of the graph
landscape. Proof:
Look at a node v.
By the ActiveNodePathLemma, there is a simple path from v to the
source node s in the residual network.
Such a path has at most n1 nodes.
By the SteepedgesLemma, the difference in height of two
neighboring nodes is at most 1.
Thus, the height covered by the path is at most n1.
The height of the source node s is always d(s)=n.
Since the sum of both is d(s)+(n1)*1 = n+n1 = 2*n1,
d(v) cannot overrun 2*n.
PPPartialCorrectnessLemma: The fact that, if the PPAlgorithm
terminates, then its result is a maximal flow. This is a claim of
partial correctness. Proof:
First of all, the result of the algorithm is a flow, because the
algorithm only terminates if the excess of all normal nodes is 0.
In order to show that it is maximal, it suffices to show that
there is no path from the source node to the target node in
the residual network (by the ResidualPathLemma).
Since d(s)=n and d(t)=0 and since the length of the path is
at most n1, there would have to be steep edges in this path in the
residual network.
This contradicts the SteepEdgesLemma.
MaxRelabelsLemma: The fact that the PPAlgorithm executes at most
2*n^2 relabel operations. Proof: There are n nodes and each of them
is relabled at most 2*n times (by the MaxLevelLemma). Thus, the
total number of relabeling operations can be 2*n^2 at most.
SearchcostLemma: The fact that total time for the search for eligible
edges in the PPalgorithm can be reduced to O(m*n). Proof: If each
active node maintains a pointer to its edges, the eligible edges can
be selected one after the other, avoiding a search. If the node is
relabeled, its pointer needs to be reset. No edge e=(u,v) prior to
the current one may ever become eligible again: This could only
happen if flow was sent back to the node via the reverse edge (v,u).
Since flow only goes downward, (v,u) would be downward, (u,v) would be
upward and hence not eligible. Since a node is at most relabeled 2*n
times, the total time needed for searching eligible edges amounts to
SUM {2*n*degree(v)  v in V} = 4*n*m in O(n*m).
MaxSatPushLemma: The fact that the PPAlgorithm executes at most
n*m saturated pushes. Proof:
A saturating push removes an edge (u,v) from the residual network.
In order to execute a second saturating push, the edge has to give
away flow first.
This can be done by pushing flow across (v,u).
In order to push flow across (v,u), v has to be lifted above u.
By the MaxRelablesLemma, this vice versa overrunning can only
happen n times.
Since there are m edges, there can be at most n*m saturating pushes.
MaxNonsatPushLemma: The fact that the PPAlgorithm executes at
most (1+m)*2*n^2 nonsaturating pushes. Proof:
Use a "potential function" phi := SUM {d(v)  v is active}
Initially, phi is 0.
A relabel operation increases phi by 1.
By the MaxRelablesLamma, there are less than 2*n^2 relabels.
A saturating push increases phi by at most 2*n, because one new
node becomes active and its level is at most 2*n.
By the MaxSatPushLemma, there are at most n*m saturating pushes.
Hence, the total increase of phi is at most 2*n^2*1 + 2*n*n*m
A nonsaturating push decreases phi by at least 1, since one node
pushes its flow away and becomes inactive.
phi will never be negative, because it is the sum of positive values.
Hence, the maximal number of nonsaturating pushes is (1+m)*2*n^2.
Since the nonsaturating pushes make up the greatest time factor,
the PPAlgorithm runs in time O(n^2*m).
PPTerminationLemma: The fact that the PPAlgorithm terminates.
Proof: In each loop, the algorithm either does a push or a relabel
operation. The relabeloperations are bound by the MaxrelabelsLemma and
the pushoperations are bound by the MaxSatPushLemma and the
MaxNonsatPushLemma. Hence, the algorithm terminates.
PPTotalCorrectnessLemma: The fact that the PPAlgorithm is totally
correct. This follows from the PPPartialCorrectnessLemma and the
PPTerminationLemma.
HighestLevelPreflowPushAlgorithm, HighLevelPPAlgorithm: A
PPAlgorithm, which selects always the active nodes at the highest
level. To do so, it maintains a priority queue with priority d, which
hands out nodes with a high level first.
Number of nodes dominated by a node v in PPalgorithm: The number
of nodes, which have a level equal to or less than d(v).
dominatedNodes(v) := {u  d(u) =< d(v), u in V}
Tuning parameter of a PPalgorithm runtimeproof: A real value named K.
This value will be chosen to be K := sqrt(m) in the
HighLevelPPTheorem.
Scaled Number of nodes dominated by a node v in PPalgorithm: The
number of dominated nodes, divided by K:
d'(v) := dominatedNodes(v)/K
d* of a PPalgorithm: The highest level.
d* := max({d(v)  v is active})
phase of a PPalgorithm: The set of all pushoperations, which are
carried out between two changes of d*.
Expensive phase of a PPalgorithm: A phase, which involves more than
K nonsaturating pushes.
Cheap phase of a PPalgorithm: A nonexpensive phase. K serves as a
parameter for distinguishing expensive and cheap phases.
MaxPhasesLemma: The fact that a PPalgorithm has at most 4*n^2
phases. Proof: Every change of d* starts a new phase. d* can only be
increased by Relabeloperations. According to the MaxRelabelsLemma,
there are at most 2*n^2 relabeloperations. Thus, d* is only
increased 2*n^2 times  and hence cannot be decreased more than 2*n^2
times. As a result, there is a total upper bound of 4*n^2 changes.
This also bounds the number of nonsaturating pushes in cheap phases,
because there are at most 4*n^2 phases and every cheap phase has at
most K nonsaturating pushes.
Domination potential function of a PPalgorithm: The function
phi := SUM {d'(v)  v is active}
PhiBoundLemma: The fact that the value of the
Domination potential function phi of a PPalgorithm is always >= 0
and initially =< n^2/K. Proof: Since phi is a sum of positive
values, it is greater than or equal to 0. Since initially, the source
node dominates over all nodes (n) and all the others (n1) dominate
over the others (n2), we get a sum of dominated nodes of
n+(n1)*(n2) =< n^2. Since the domination potential function is this
value, divided by K, the claimed bound results.
PhiIncreaseLemma: The fact that each relabeloperation and each
saturating push increases the domination potential function phi of a
PPalgorithm by at most n/K. Proof: The relabel operation changes
phi by increasing the d'(v) of one node. A saturating push also
changes phi by the d'(v) of one node, because the flow receiving node
may become active and thus change d'. In general, it holds for any
node v that
d'(v) := dominatedNodes(v)/K =< n/K
Thus, the above operations can change phi only by at most n/K.
PhiDecreaseLemma: The fact that the domination potential function
phi of a PPalgorithm is not increased by a nonsaturating push
operation. Proof: If a push across the edge (u,v) is nonsaturating,
then u got rid of its excess and does no longer belong to the set of
active nodes. By contrast, v may become active. But since v is below
u, its d'(v) will be smaller than d'(u). Hence, the potential will
not increase, but rather decrease.
// It seems to be the wrong way round in the slides.
ExpensivePhaseLemma: The fact that an expensive phase of a
HighLevelPPalgorithm decreases the domination potential by at
least the number of nonsaturating pushes, which is >K. Proof:
In the HighLevelPPAlgorithm, the nodes at level
d* are chosen first. Since d* does not change during a phase, such
a phase will empty the level d*. By the phidecreaseLemma, each
nonsaturating push decreases the potential by at least 1. Since
by definition, an expensive phase compromises at least K
nonsaturating pushes, phi decreases by at least K.
HighLevelPPTheorem: The fact that the HighLevelPPAlgorithm runs
in time O(n^2*sqrt(m)). Proof:
By the PhiBoundLemma, phi is initially smaller than n^2/K.
By the MaxRelabelsLemma, there are at most 2*n^2 relabeloperations.
By the MaxSatPushesLemma, there are at most n*m saturating pushes.
By the PhiIncreaseLemma, a relabeloperation or a saturating push
increases phi by at most n/K.
Hence, phi is always below (2*n^2+n*m)*n/K + n^2/K =
(2*n^3+n^2+m*n^2)/K
Thus, phi can only be decreased (2*n^3+n^2+m*n^2)/K times.
By the ExpensivePhaseLemma, each expensive phase decreases phi by
the number of nonsaturating pushes.
Thus, there can be only (2*n^3+n^2+m*n^2)/K nonsaturating pushes in
all expensive phases.
By the MaxPhasesLemma, there are at most 4*n^2*K nonsaturating
pushes in all cheap phases together.
This results in a total sum of nonsaturating pushes of
(2*n^3+n^2+m*n^2)/K + 4*n^2*K
This value can be minimized by chosing K = sqrt(m), resulting in a
total number of nonsaturating pushes of O(n^2*sqrt(m))
Since the nonsaturating pushes make up the greatest time factor of
the algorithm, the whole algorithm is in time O(n^2*sqrt(m)).
Heuristic improvement of the PPAlgorithm: A modification of the
PPAlgorithm, which aims at a better performance.
Aggressive Relabeling, Local Relabeling: A heuristic improvement of
the PPAlgorithm, which relabels a node v by
d(v) := 1 + min({d(w)  (v,w) in E})
Instead of increasing d(v) one by one until d(v) overtakes another
d(w), the value is increased immediately in such a way that it
overtakes the smallest neighbor. This works like a sequence of
normal relabelings and saves unneccessary steps.
Global Relabeling: A heuristic improvement of the PPAlgorithm,
which sets d to
d := reverseBFSDepth(Gf,t)
where Gf is the residual network of the graph and t is the target
node. This happens at the beginning and then periodically every O(m)
phases. Thus, the nodes are given their "optimal" level, which
encourages the flow towards the target node.
BackBeam: A heuristic improvement of the PPAlgorithm, which treats
nodes with a level greater than n apart. Nodes of this kind cannot be
connected to the target node in the residual network, as this would
entail steep edges. They may only send their flow back towards the
source node.
Gap heuristics: A heuristic improvement of the PPAlgorithm, which
bases on the observation that if there is an empty level, then the
nodes above this level will send their flow back towards the source.
If there is a level with no nodes (a "gap"), then there cannot be any
edge in Gf from the upper level the lower level (no steep edges).
Hence, the only paths from nodes in the upper level go to the source
node (at level n). Now, the nodes which still have excess will
struggle to rise their level until they can send back the flow to
the source. This process can be shortened by the following rule
IF Ex i : ~(Ex v in V: d(v)=i) THEN
FOREACH v in V: d(v)>i DO
d(v):=n
Thus, a number of previsible steps is avoided.
CostCapacityGraph: The type
CLASS costCapacityGraph EXTENDS Graph
cost : Edge > Real
capacity : Node > Real+
supply: Edges>Real,
such that SUM {supply(v)  v in V} = 0
This graph can be seen as a transportation network, in which the nodes
supply a good (supply(v)>0) or demand that good (supply(v)<0).
Meanwhile, the edges have a maximum amount of that good, which they
can carry (capacity). They also have costs for every unit of the good
crossing the edge. The goal is to find the least expensive flow of
goods from the suppliers to the demanders.
Supply node in a cost capacity graph: A node with a positive supply
value.
Demand node in a cost capacity graph: A node with negative supply
value.
Feasible flow, flow in a cost capacity graph (V,E,c,cap,supply): A
function f:E>Real, such that
* f(e) >= 0 for all edges e
* f(e) =< cap(e) for all edges e
* excess(v) = supply(v) for all nodes v
All notions relative to flows in capacity graphs are taken over,
particularly "excess" and "saturation". A flow in a cost capacity
graph can be found by finding maximal flows: We add a source node s
and connect s to all nodes v with positive supply with a capacity of
supply(v). Then we add a target node t and connect all nodes v with a
negative supply to t with a capacity of supply(v). A maximal flow in
this capacity graph outlines a feasible flow.
Cost of a flow f in cost capacity graph (V,E,c,cap,supply): The sum of
the edge costs, weighted by the flow:
c(f) := SUM {c(e)*f(e)  e in E}
MinCostFlow in cost capacity graph: The flow with the least cost.
Residual Cost network of a cost capacity graph (V,E,c,cap,supply) with
a flow f: The cost capacity graph (V,E',c',cap',supply) with
* E' =
{e  e is an unsaturated edge in E } \/
{(v,u)  (u,v) is an edge with flow in E}
* cap'(u,v) = cap(u,v)f(u,v) if (u,v) unsaturated in E
cap'(v,u) = f(u,v) if (u,v) has flow in E
* c'(u,v) = c(u,v) if (u,v) unsaturated in E
c'(v,u) = c(u,v) if (u,v) has flow in E
This corresponds to a residual network of a capacity graph, to which
the approproiate edge costs have been added. If a flow can be undone,
then the edge costs in the residual cost network are negative to
indicate a potential decrease in cost. Algorithm:
FUNCTION residualCostNetwork((V,E,c,cap,supply) : CostCapacityGraph,
f : FUNCTION E>Real) : CostCapacityGraph
result : CostCapacityGraph := (V,E'={},cap'={},c'={},supply)
FOREACH (u,v) in E DO
IF cap(u,v)>f(u,v) THEN
E' := E' \/ {(u,v)}
cap' := cap' \/ {((u,v),cap(u,v)f(u,v))}
c' := c' \/ {((u,v),c(u,v))}
IF f(u,v)>0 THEN
E' := E' \/ {(v,u)}
cap' := cap' \/ {((v,u),f(u,v))}
c' := c' \/ {((v,u),c(u,v))}
RETURN result
MinCostCycleLemma: The fact that the existence of a MinCostFlow in
a cost capacity graph is equivalent to the absence of negative cycles in
the residual cost network. Proof: (?).
CycleCancelingAlgorithm: The following algorithm, which finds a
MinCostFlow in a cost capacity graph:
FUNCTION cycleCancel((V,E,c,cap,supply) : CostCapacityGraph) :
FUNCTION E>Real
f := feasibleFlow(V,E,c,cap,supply)
Gf := residualCostNetwork((V,E,c,cap,supply),f)
WHILE Ex cycle C in Gf: c(C)<0 DO
minVal := min({cap'(e)  e in C})
FOREACH e in C DO
f(e):=f(e)+minVal
Gf := residualCostNetwork((V,E,c,cap,supply),f)
RETURN f
This algorithm bases on the MinCostCycleLemma: Any feasible flow is
taken and then all negative cycles in it are removed. In doing so, the
flow still remains a valid feasible flow.
Transportation problem: The problem of, given a cost capacity graph
with infinitely large capacities, finding a MinCostFlow. This
results in a focus on the costs, ignoring the edge capacities.
MinimumCostBipartiteMatchingProblem: The transportation problem in
a bipartite cost capacity graph (A \/ B,E,cap,c,supply), where
* supply(a)=1 for all a in A
* supply(b)=1 for all b in B
* cap(e)=+oo for all e in E
For this problem, the desired mincostflowvalues may only be integer,
i.e. f(e)=0 or f(e)=1 for all edges. This mirrors the idea of
suppliers in A, which want to get their goods to the demanders in B
with a minimal transportation cost.
Weight of a matching in an undirected cost graph (V,E,c): The sum of
the weights of the edges.
c(M) := SUM {c(e)  e in M} for a matching M
MaximumWeightMatching: The problem of, given a weighted graph
(V,E,c), finding a matching with a maximum weight.
ApproximateWeightedMatchingAlgorithm: The following algorithm, which
finds a matching in a weighted graph, the weight of which is at least
half of the weight of a maximumweightmatching:
FUNCTION appWeightedMatching((V,E,c) : UndirectedCostGraph) :
Subset of E
M' : List of Edge := <>
WHILE E' != {} DO
// Pick any node
v:= eps(V)
// Build up a heavy path greedily
WHILE degree(v)>0 DO
// Find the edge leaving v with a maximum weight
(v,w) := eps({(v,u)  c(v,u)=max({c(v,u)(v,u) in E})})
M'.pushBack(v,w)
V := V\{v}
v := w
// Now we have a set of edges, try building a matching
// of all those at even positions
weight1 := SUM {c(M'.get(i*2))  i=0..M'.length/2}
// Try building a matching with all odd edges
weight2 := SUM {c(M'.get(i*2+1))  i=0..M'.length/2}
// Return the matching with the larger weight
IF weight1>weight2 THEN
RETURN {M'.get(i*2)  i=0..M'.length/2}
ELSE
RETURN {M'.get(i*2+1)  i=0..M'.length/2}
This algorithm greedily collects edges to paths. It always choses the
heaviest edge to go on. Having taken a node, it deletes the node in
order to avoid using it again. If it cannot go on, it restarts at
another node. Then it tries the two possible matchings contained
in the list:
a > b ==> c > d ==> e > f ...
The two matchings are the set of ">" edges and the set of "==>"
edges. The algorithm takes the matching, which is more heavy. Let's
look at a maximal weight matching M*. Consider any edge e* in M*.
Either e* is in M', then c(e*) is in c(M*) as well as in c(M'). Now
assume e* is not in M'. Let v be the node of e*, which is deleted
first by the algorithm (either the start node or the end node of e*).
If the algorithm arrives at v, the edge e* is still alive. Now the
algorithm will pick the heaviest edge leaving v. As, according to the
assumption, e* is not picked, another edge e' is picked with
c(e')>c(e*). Now M* cannot contain e', because M* already contains
e* and each node has degree 1 in M*. Hence, M' contains a heavier
edge than M*. If this argument is repeated for all edges, it results
that the sum of the weights of the edges contained in M' is more than
the weight of M*. Hence, the weight of the "more heavy half" of M'
is greater than c(M*)/2.
MulticastingNetwork: A weighted directed acyclic graph (V,E,cap)
with a source node s in V and a target node set T < V. All nodes
except for the source node have the same outdegree. (?)
Rate of a multicastingnetwork: The capacity of the smallest of the
minimal stcuts, where t runs through all the target nodes.
h = min{cap(S)  S is a minimal stcut with t from T}
Information flow in a multicastingnetwork (V,E,cap,S,T): A function
f:E>{0,1}, which maps each edge to a bit, such that
f(u,v) = f(v,w) for all u,v,w in V\T\{s}
If f(u,v)=b, one says that node u sends bit b to node v.
MulticastingProblem: The problem of, given integer numbers r and t,
creating a multicastingnetwork with rate r, such that the source
node can send a maximal number of bits in the same time, which are
all received by the t target nodes. The problem is that each node only
has r leaving edges and that two nodes may only be connected by one
edge. Thus, intermediate nodes have to be created to multiply the
information.
Coding in a multicastingnetwork: The idea that the source node does
not send only the bits it should send, but applies some function
to them. The target nodes undo this function and get the real bits
again. With coding, h bits can be sent at a time, if h is the rate
of the network. ___________
Example: __/  Here, each target node
_____a_______ v1 _____ t1  receives either bit a and
/ /  b, or one of them and the
s a xor b v2 === t2  value a xor b. By one bit
\_____b_______ /  and the value xor b, the
'v3 === t3  other bit can be
_______ calculated.
Linear Coding: Coding in a multicastingnetwork, which involves
vectors associated to each edge (?).
The essential thing seems to be that the coding can be done in
polynomial time. (?)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Strings
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Sorting algorithm: An algorithm, which produces the sorted sequence
of a given set.
Rank of an element x in a set S: The index of x in the sorted
sequence of S.
Bucket sort: The following sorting algorithm for a set of natural
numbers in a fixed range:
FUNCTION bucketSort(S : set of Integer) : array of Integer
// Buckets count occurences of numbers
bucket : array[max(S)] of Integer :=<0,0,...0>
FOREACH i in S DO
bucket[i] := bucket[i] + 1
// Collect buckets to form the ordered sequence
result : array[S] of Integer
resultIndex := 0
FOR bucketIndex := 0 TO bucket.length()1 DO
WHILE bucket[bucketIndex] > 0 DO
result[resultIndex] := bucketIndex
bucket[bucketIndex] := bucket[bucketIndex]  1
RETURN result
This algorithm has one bucket for each number which may appear in S.
It the goes through all the numbers of the set and throws them into
the corresponding buckets (by increasing a counter). Then, it goes
through all the buckets and constructs the sorted sequence. This
algorithm takes time O(S+max(S)).
Char, Character: A symbol. Usually, a character is encoded by a
natural number.
Alphabet: A finite set of characters.
Ordered alphabet: An alphabet with a linear order on it. By the help
of equality, such an alphabet also posesses a simple order.
Constant alphabet: An ordered alphabet of a fixed and finite size.
Integer alphabet: The alphabet {1,...K} for some natural number K.
Any constant alphabet can be represented by an integer alphabet.
When an algorithm uses an integer alphabet, one describes its
runtime also in dependence of K.
// The paper on the Skewalgorithm defines an integer alphabet
// slightly different
Sentinel character of an ordered alphabet A: A special character,
which is smaller than all characters of A. Thus, the character cannot
be in A. "Smaller than" is to be read in the sense of the linear
order of A. In an integer alphabet, the sentinel character is usually 0.
// This corresponds to the final zero of CStrings and simplifies
// sorting
String over an ordered alphabet A: A sequence of elements of A, plus
the sentinel character of A. We will assume that A be an integer
alphabet with the sentinel character 0. The sentinel character is not
noted.
Notation: S := S.length()
S =
S[i] := for i>=S
S.alphabet := A
"abc" := <'a','b','c'>
S1 + S2 := + S2 // "concatenation"
Substring from i to j of a string S: The subsequence of S, starting
at position i and extending to position j1, plus the sentinel character.
Notation: S[i:j) :=
S[i:j).length() = ji
// Like for intervals, the ")" is probably supposed to indicate
// that the character at position j does not belong to the substring.
Prefix of size i of a string S: The substring from 0 to i of S.
"iSuffix of a string S", Suffix at position i in S: The substring
from i to S.length() of S.
Notation: S[i..]
Lexicographic order: The following linear order on the set of Strings:
S1 < S2 :<=> S1[0] < S1[0] 
(S1[0]==S2[0] &&
S1[1..] < S2[1..]) )
Pattern: A (short) string.
Occurence of a pattern P in a string S: A substring of S, which equals
P.
Distinguishing prefix of a string S in a set of strings: The shortest
prefix of S, which is not a prefix of another string in the set. If
for a string S, there is no such prefix, then the distinguishing prefix
of S is S itself.
String sorting problem: The problem of, given a set of strings, producing
a sorted sequence of them. A string sorting algorithm only needs to look
at the distinguishing prefixes of the strings.
Standard sorting algrithm: A sorting algorithm, which can sort any
ordered set. For string sorting, it is possible to perform better than
Standard sorting algorithms, because one can exploit the internal
structure of strings and the recusrive definition of their order.
Standard sorting algorithms do O(n*log(n)) element comparisons, where
n is the number of input elements. Depending on the size of the strings,
each comparison may again contribute a significant time factor.
MultikeyQuicksort: The following string sorting algorithm:
FUNCTION multiKeyQuickSort(S : Set of String, p : Integer := 0)
: array of String
ASSERT All strings of S have the same prefix of length p
// Trivial cases
IF S==0 THEN RETURN <>
IF S==1 THEN RETURN
// Choose Pivotelement
pivot := eps(S)
S1 := {s in S  s[p] < pivot[p]}
S2 := {s in S  s[p] == pivot[p]}
S3 := {s in S  s[p] > pivot[p]}
RETURN multiKeyQuickSort(S1,p) +
multiKeyQuickSort(S2,p+1) +
multiKeyQuickSort(S3,p)
This algorithm get a set of strings S, which already share the first p
characters. It chooses one arbitrary string of S, called the "Pivot
element". Then it partitions S into 3 subsets by looking at the pth
character: Depending on whether this character is smaller than, equal
to or greater than the pth character of the pivot element, the
string goes to S1, S2 or S3. Due to the nature of the lexicographical
order, this means that all strings in S1 are smaller than the pivot
element, all strings in S2 are greater than the pivot element and
nothing can be said about the strings of S2. Then, each group is
sorted recursively, where p advances only for the group S2. There
are two cases: Let's look first at one string s, which survives each
recursive call in the set S1. Assuming a perfect choice of the
Pivotelement, this set is halved each time. Thus, recursion depth
cannot exceed log(n), where n is the number of strings. Since about
O(n) strings face the same fate as s, and each of them requires the
comparisons, runtime is at least O(n*log(n)). In the other case, the
strings share long common prefixes. Then, all strings go to S2. The
algorithm will only stopped after having looked at all characters of
the distinguishing prefixes of the strings. Let D denote the sum of
the lengths of the distinguishing prefixes, then the runtime is
O(D+n*log(n)). It can be shown that for a random choice of the
Pivotelement, the same runtime order results.
LSDfirst Radixsort, Least significant digit first Radix sort: The
following string sorting algorithm:
FUNCTION LSDradixSort(S: Set of String) : array of String
ASSERT eps(S).alphabet is an integer alphabet
alphsize := eps(S).alphabet
buckets : array[2][alphsize] of List of String <<<>,<>,..>..>
FOREACH s in S DO
bucket[0][0].pushBack(s)
switch := 1
FOR len := max({s  s in S}) DOWNTO 0 DO
switch := (switch+1) mod 2
FOR i := 0 TO alphsize DO
WHILE bucket[switch][i] != <> DO
s := bucket[switch][i].popFront()
bucket[(switch+1)mod 2][s[len]].pushBack(s)
result : array[S] of String
j := 0
FOR i := 0 TO alphsize DO
WHILE bucket[switch][i] != <> DO
result[j] := bucket[switch][i].popFront()
j := j + 1
RETURN result
This procedure does a bucket sort for each position ("digit") of the
strings. It starts with the rightmost (the "least significant"), i.e.
the last character of the longest string (the others automatically
have the 0 there). It puts the strings into the buckets. Then, it
goes through all buckets (in the right order) and puts the strings
into new buckets according to their butlast character. This continues
until all positions have been sorted. The running time in O(N), if N
is the total number of characters. This algorithm does not make use of
the fact that the distinguishing prefix (which may be much shorter
than the string) suffices for sorting.
MSDRadixSort, most significat digit radix sort: The following string
sorting algorithm:
FUNCTION MSDRadixSort(S : List of String, p : Integer := 0)
: List of String
ASSERT All strings share the first p characters
ASSERT eps(S).alphabet is an integer alphabet
alphsize := S.alphabet
IF S.length(),<>,...>
FOREACH s in S DO
bucket[s[p]].pushFront(s)
FOR i:=0 TO alphsize1 DO
bucket[i] := MSDRadixSort(bucket[i],p+1)
RETURN bucket[0] + bucket[1] + ... + bucket[alphsize1]
This algorithm does repeated bucketsorts like LSDRadixSort, but it
starts at the first character: It puts the strings into buckets
according to their first character and then sorts the buckets
recursively. The parameter p indicates the next position to analyze.
If the size of the set to sort is smaller than the size of the
alphabet, the recursion breaks and multikeyQuicksort is used. This
happens in order to avoid that initializing the buckets to <> takes
longer than the actual sorting process. Taking this into
consideration, the runtime is dominated by the FOREACHloop. The
total number of pushFront()operations for all strings and through all
recursions is D, if D is the sum of the lengths of the distinguishing
prefixes. This is because after having bypassed the distinguishing
prefixes, all strings are in oneelementsets. Thus, running time is
at least O(D). In addition, a break of the recursion causes a call to
multikeyQuicksort, which has a runtime of O(n*log(n)), where n is the
number of strings. MultikeyQuicksort is only called on alphsize
strings and thus, the overall runtime is O(D+n*log(alphsize)).
Suffix Sequence of a string S: The sorted sequence of all suffixes of
S.
Suffix sorting problem: The problem of, given a string S, creating
its suffix sequence.
Example: S="banana"
Suffixes ~~~> Sorted suffixes
Pos Pos Rank
0 banana 5 a 0
1 anana 3 ana 1
2 nana 1 anana 2
3 ana 0 banana 3
4 na 4 na 4
5 a 2 nana 5
Suffix array of a String S: An array A of integers, where A[i] is
the starting position of the suffix with rank i in the sorted
sequence.
Notation: SuffixRankToPos[i]=j
means: The starting position of the ith suffix is j
Example: S = "banana"
SuffixRankToPos = <5,3,1,0,4,2>
^ The second suffix in the sorted
sequence is the 3suffix
S[3..] = "ana"
Inverse suffix array of a string S: An array A of integers, where A[i]
is the rank of suffix S[i..].
That is: The array can be read as an "annotation" for the string S,
which notes for each suffix its rank in the sorted sequence of all
suffixes
Notation: SuffixPosToRank[i]=j
means: The rank of the isuffix is j
Example: S="banana"
SuffixPosToRank = 325140
^ Suffix "ana" is has rank 1 in the sorted
sequence of all suffixes (cf. ex. above)
The suffix array can be translated to the inverse suffix array and
vice versa by doing
PROCDURE swap(A : array of Integer) : array of Integer
B : array[A.length] of Integer
FOR i:=0 TO A.length1 DO B[A[i]] := i
RETURN B
// I find the notions "suffix array" and "inverse suffix array" hard
// to distinguish. Therefore, I chose different names.
// SuffixRankToPos means: I put the rank as an index, so I get the
// starting position as a result. (and vice versa)
Lexicographic names for a set S: A function f:S>A so that
S[i] f(S[i])
triples2 : array[S/3] of String :=
// Put these triples together, sort them
triples1&2Sorted : array of String
:= LSDradixSort({s  s in triples1 or triples2})
// Create lexicographic names for triples1&2Sorted
lex : FUNCTION triples1 \/ triples2 > Integer := {}
lex(triples1&2Sorted[0]):=0
FOR i:=1 TO triples1&2Sorted.length1 DO
IF triples1&2Sorted[i]!=triples1&2Sorted[i1] THEN
// Different from previous => new lexicographic name
lex(triples1&2Sorted[i]) := lex(triples1&2Sorted[i1])+1
ELSE
// Same as previous => same lexicographic name
lex(triples1&2Sorted[i]) := lex(triples1&2Sorted[i1])
// Translate triples1 and triples2 to their lexicographic
// counterparts
triples1Lex : array[S/3] of Integer
triples2Lex : array[S/3] of Integer
FOR i:=0 TO S/3 DO
triples1Lex[i] := lex(triples1[i])
triples2Lex[i] := lex(triples2[i])
// Concatenate triples1Lex + triples2Lex
triples1&2Lex : array of Integer := triples1Lex+triples2Lex
// Create their suffix sequence
triples1&2LexSuffixRankToPos : array of Integer
// Do we have all distinct elements?
IF lex(triples1&2Sorted[triples1&2Sorted.length1] ==
triples1&2Sorted.length1 THEN
// Their lexical number is their SuffixPosToRank,
// so do a swap swap to get the SuffixRankToPos
triples1&2LexSuffixRankToPos := swap(triples1&2Lex)
ELSE
// Use recursion
triples1&2LexSuffixRankToPos := skew(triples1&2Lex)
// Convert this back to the original suffixes
suffixes1&2Sorted : array[S*2/3] of (Suf:String,Pos:integer)
FOR i := 0 TO S*2/31 DO
pos := triples1&2LexSuffixRankToPos[i]
IF pos < triples1Lex.length THEN pos := pos*3+1
ELSE pos := (postriples1Lex.length)*3+2
suffixes1&2Sorted[i] := (S[pos..],pos)
// Now take care of the suffixes starting at positions i,
// where i mod 3 = 0: Get them and sort them. The suffix of
// each of these suffixes has already been sorted in the previous
// rounds, so we only need to look at the first character and
// the position of the following suffix in the sorted sequence
// (This is neglected here)
suffix0 : Set of (String,integer) := {(S[i..],i)i mod 3 == 0}
suffix0Sorted : array of (Suf:String,Pos:integer) := sort(suffix0)
// Now merge suffixes1&2Sorted and suffixes0Sorted
// This merging can be done efficiently by taking into account
// the subsuffixes, which are already sorted
// (This is ignored here)
result : array[S] of integer
i := 0
index0 := 0
index1&2 := 0
WHILE i < S DO
IF suffix0Sorted[index0]
SuffixRankToPos := skew(S)
L := 0
FOR i:=0 TO S1 DO
// Look at the suffix starting at i
// Does it have any predecessor?
IF SuffixPosToRank[i]>1 THEN
// Find its predecessor in the sorted sequence
pred := SuffixRankToPos[SuffixPosToRank[i]1]
// Count shared characters, starting at L (!)
WHILE S[i+L]==S[pred+L] DO
L := L + 1
lcp[SuffixPosToRank[i]] := L
L := L  1
IF L < 0 THEN L := 0
RETURN lcp
The main loop is executed S times. In each execution, L may be
increased by the inner WHILEloop and is then decreased by 1. Assume
the first WHILElooprun increases L to n. As L cannot bypass n, each
following WHILElooprun can at least increment L by one. Thus, we
still have a total time of S+S, i.e. O(S). If another WHILE
looprun wants to increase L more than once, say k times, then the
first WHILElooprun has to give up exactly k incrementations. But all
in all, the total number of incrementations will not bypass n. Thus,
the algorithm needs time O(S). L only needs to be decreased by 1
each time due to the following LCPContinuitylemma.
LCPContinuitylemma: The fact that, given a string S and a position i,
the isuffix has more shared characters with its predecessor in the
sorted sequence than the (i1)suffix minus 2.
lcp[SuffixPosToRank[i]] > lcp[SuffixPosToRank[i1]]  2
Example: S="banana", i=3
lcp[SuffixPosToRank[3]] > lcp[SuffixPosToRank[2]]  2
lcp[4] > lcp[5]  2
1 > 2  2
Picture:
SuffixSeq Rank lcp
A 0 0
Ana 1 = SuffixPosToRank[3] 1 ,
ANAna 2 3 
banana 3 0  1 > 2  2
na 4 0 
NAna 5 = SuffixPosToRank[2] 2 '
"ana" has more shared characters with its predecessor "a" than
the suffix "nana" has with its predecessor "na", minus 2.
Proof: If the (i1)suffix shares less than 1 character with its
predecessor, then the lemma follows immediatly, because the lefthand
side of the inequality is positive, while the righthand side is
negative. Now let's suppose the (i1)suffix shares L>1 characters
with its predecessor, i.e L=lcp[SuffixPosToRank[i1]]>1. In the
example, i1 is 2 and the 2suffix "nana" shares 2 characters with its
predecessor "na". Be j the position of this predecessor in the string,
i.e. j=SuffixRankToPos[SuffixPosToRank[i1]1], which is 4 in this
case. Since the jsuffix and the (i1) suffix share L characters, the
(j+1)suffix and the isuffix must share L1 characters ("a" and "ana"
share 21=1 characters). Since the jsuffix preceded the (i1)suffix
in the sequence and they share L characters, the (j+1)suffix must
also precede the isuffix. In the example, the suffix "a" precedes the
suffix "ana", although there may be other suffixes in between. The
goal is to show that the isuffix shares more than L2 characters with
its immediate predecessor. We already showed that the isuffix shares
L1 characters with the (j+1)suffix, which is a (possibly far)
predecessor. But since the suffixes are sorted, all suffixes between
the (j+1)suffix and the isuffix also have to share at least those
L1 characters. In particular, the direct predecessor of the isuffix
must share at least L1 characters with the isuffix, which is the
claim.
// I owe thank to David Steurer, who helped me understanding this
// lemma.
Labelled tree: The type
CLASS LabelledTree EXTENDS Graph
ASSERT the graph is a tree
root : Node
label:Edges>X
d:Node>Y
This is a directed tree with a function label, which equips each edge
with a name or symbol and a function d, which equips each node with a
name or symbol.
Suffixtree of a string S: A labelled tree, so that
* its edges are labelled with substrings of S
* there is one leaf for each suffix of S
* the concatenation of the edge labels on a path from the
root to a leaf gives the suffix of the leaf
* labels of edges leading to siblings have a distinct first
character
* the label of a node is the length of the concatenations of the
labels of the edges on the path from the root to the node.
The suffixtree can be constructed by hand by the following (simple)
algorithm:
* Do a suffix sort of S, taking into account the sentinel character.
* Introduce a root node and add all suffixes of S as children
* While there is still a suffix x!=<0> which is a prefix of sibling
suffixes of x
* chop off the trailing 0
* add 0 as a child of x
* make all sibling suffixes with prefix x children of x
* chop off x of all these new children
* Label the edges with the node they lead to and compute d
Example: S="banana"
Suffix sort Make tree
root
0 banana0 ~~~~~> 6 0 ~~~~~> > 0
1 anana0 5 a0 > a > 0
2 nana0 3 ana0  '> na > 0
3 ana0 1 anana0  '> na0
4 na0 0 banana0 > banana0
5 a0 4 na0 > na > 0
6 0 2 nana0 '> na0
The suffix tree can also be constructed by the following (complicated)
algorithm:
FUNCTION suffixTree(S : String) : LabelledTree
SuffixRankToPos := skew(S)
lcp := lcp(S)
T : LabelledTree := (V={root},E={},label={}, d={})
d(root) := 0
v := root
FOR i := 0 TO S1 DO
// Insert the suffix with rank i into the tree
WHILE d(parent(v)) > lcp[i] DO v := parent(v)
IF d(v) > lcp[i] THEN
// Split the edge
v' := newNode(V)
d(v') := lcp[i]
e1 := newEdge(E,parent(v),v')
label(e1) := label(parent(v),v)[0..lcp[i]]
e2 := newEdge(E,v',v)
label(e2) := label(parent(v),v)[lcp[i]+1..]
E := E \{(parent(v),v)}
v := v'
// Now d(v) == lcp[i]
v'' := newNode(V)
d(v'') := SSuffixRankToPos[i]
e := newEdge(E,v,v'')
label(e) := S[SuffixRankToPos[i]+lcp[i]..]
RETURN T
This algorithm inserts the suffixes according to their sorted sequence
into the tree. Therefore, it can build the edges from the root to the
leaves. Only if a suffix is not a prolongation of its predecessor, the
algorithm needs to go up in the tree. If it finds a suitable node
(lcp[i]==d[i]), it simply adds the new suffix. If not, it breaks the
edge in two: The first half gets the first part of the label, the
second gets the second part. In between, a new node is added. To this
new node, the suffix can then be added as usual. This algorithm works
similar to DFS. The main loop loops n:=S times. Since each edge is
created once and backtracked at most once, we get an additional time
consumption of O(2*n), which does not influence the overall running
time of O(n).
Multikey Suffix Quicksort: The Multikey Quicksort algoroithm, applied
to the suffixes of a string. For the string S with n=S and sigma=
S.alphabet, the following holds:
* It is not necessary to make copies of all suffixes, because the
algorithm can as well work on pointers to the suffixes in S.
* If sigma>n, then one can design S such that the sum D of
the lengths of the distinguishing prefixes is in O(n): One just
creates S such that it does not contain equal characters. As a
result, the algorithm does not need any recursive call.
* For two nonoverlapping substrings S[i:i+k) and S[j:j+k), the
probability that they are equal is sigma^k. Proof: For
two randomly chosen characters, the probability of being equal is
1/sigma. If the strings are to be equal, their first character has
to be equal and their second character has to be equal etc. Since
the choices of characters are independent, it is
P(S[i+1]=S[j+1]  S[i] = S[j] ) = P(S[i+1]=S[j+1])
Thus, the overall probability amounts to the product of all
single probabilities: 1/sigma*1/sigma*...1/sigma = sigma^k.
* Even if two strings overlap, the probability of being equal is
still sigma^k, because the strings are independently chosen.
* The probability of two equal substrings of length k in S is bounded
by n^2*sigma^k. The substring starting at position 0 could be
equal to (nk1) other substrings. The substring at position 1
could be equal to (nk2) substrings etc.. All in all,
(nk +1)(nk)/2 < n^2 pairs could be equal. For each of them, the
probability is sigma^k. Thus, the probability of any two equal
substrings is less than n^2*sigma^k.
* The probability that there are any two suffixes with a common prefix
of length 3*log(n)/log(sigma) or longer is at most 1/n. The bound
n^2*sigma^k for the probability of equal substrings of length k
also applies to substrings of length k or longer, because every
equal substrings of length >k also have to match within the first
k characters. Thus, we can set k=3*log(n)/log(sigma) and obtain
n^2*sigma^(3*log(n)/log(sigma)) = n^2*sigma^(log(n^3)/log(sigma))
= n^2*n^3 = 1/n
* The average case for a multikey suffix quicksort run is bounded by
O(log(n)*n). Two cases can be distinguished:
* There exist at least two suffixes, which share a prefix of length
at least 3*log(n)/log(sigma). The probability of this case is
bounded by 1/n. These two strings will force the algorithm to keep
them in the same bucket until their equal parts have been
bypassed. The sum D of the lengths of all equal parts is bounded
by n^2. Thus, we get the running time of O(D+n*log(n)) =
O(n^2+n*log(n)).
* If there do not exist any two substrings with an equal prefix of
3*log(n)/log(sigma), then D is bounded by n*3*log(n)/log(sigma).
The running time is thus O(D+n*log(n)) = O(n*3*log(n)/log(sigma)
+ n*log(n)) = O(n*log(n)).
The average case running time is the sum of the special case running
times, weighted with their probabilities:
O( (n^2+n*log(n)) * 1/n + n*log(n) * (11/n) )
= O(n*log(n))
The prefix length 3*log(n)/log(sigma) has been chosen arbitrarily:
If one choses a smaller value, the first case becomes less probable,
but the second case becomes worse. Vice versa, if one choses a
greater value, then the second case becomes better, but the first
case becomes more probable.
BurrowsWheelerTransform of a string S: The string BW(S), which is
obtained as follows: Create a matrix of size (S+1) x (S+1) and
write S with its sentinel character in the first line. Then rotate
this line one to the left with the first character wrapping to the
last position and write this into the second line. Continue until the
matrix is filled up. Now sort the lines alphabetically. The rightmost
column of the matrix, read topdown, is the BurrowsWheeler
Transformation of S. Since the sentinel character is smaller than all
other characters, the sorting of the matrix lines boils down to the
sorting of the suffixes of S. Hence, BW(S) can also be calculated as
the sequence of the characters, which precede the suffixes of the
sorted suffix sequence. The algorithm would be:
FUNCTION BW(S : String) : array of S.alphabet
SA := skew(S)
result : array[S+1] of S.alphabet
result[0] := S[S1]
FOR i:= 0 TO S1 DO
result[i+1]:=S[ (SA[i]1+S) mod S ]
RETURN result
Matching: The process of finding the position(s) in a string, where a
pattern occurs as a substring.
BruteForceMatching: The following (simple) algorithm for finding
a pattern P in a string S:
FUNCTION bruteForceMatching(P,S: String) : Set of Integer
result : Set of Integer := {}
FOR j := 0 TO SP DO
IF S[j..j+P] == P THEN result := result \/ {j}
RETURN result
This algorithm just checks out every position of S. If P matches with
all of its characters, then the position is added to the resultset.
The algorithm runs in time O(P*S).
SuffixTreeMatchingAlgorithm: The following algorithm, which finds a
pattern P in a string S by using the suffix tree of S:
FUNCTION suffixTreeMatch(P,S : String,(V,E,root,label,d) : SuffixTree)
: Set of Integer
v := root
i := 0 // Index of first character, that has not yet
been matched
WHILE i < P DO
IF Ex (v,w) in E : label(v,w) == P[i..i+label(v,w)1] THEN
i := i + label(v,w)
v := w
ELSE
RETURN {}
RETURN { S  d(v')  v' is a (grand)child of parent(v)
and v' is a leaf}
This algorithm goes down in the suffix tree with edges matching the
pattern. When the pattern has been eaten, all child suffixes of the
current node are occurences. It has a runtime of O(P+n), where n is
the number of occurences of P in S. The convenient thing about this
algorithm is that once the suffix tree of S has been calculated,
different patterns can be searched very efficiently.
FSAMatchingAlgorithm: The following algorithm, which uses a FSA to
find a pattern P in a string S:
FUNCTION FSAMatching(P,S : String) : Set of Integer
result : Set of Integer := {}
myFSA : FSA := (Q : {0,...P}, // States == positions in P
I : S.alphabet, // Input will be S
s0 : 0, // Start state = 0
F : {P}, // End state = length of P
d) // Transition function
myFSA.d(q,c) := max({ j  j=0...P,
P[0..j) == P[qj:i) + })
// with max({}) := 0
result : Set of Integer := {}
state := s0
pos := 0
WHILE pos<S DO
state := myFSA.d(state,S[pos])
pos := pos + 1
IF state in myFSA.F THEN result := result \/ {posP}
RETURN result
This algorithm first constructs a FSA for the pattern. The FSA has a
state for each position of P. It is fed with the string S. With each
character, it remains in the start state. When the first character of
P occurs, the FSA advances one state. If the second character occurs,
the FSA advances one state and so on, until the final state has been
reached, meaning that the pattern occured. If the FSA encounters a
character c which disturbs the pattern, it returns to a previous
state q'. This state is chosen such that it allows continuing a
new start of the pattern, which may lie in the pattern itself. If the
current state is q, then the new state is the length of the longest
prefix of P, which is a suffix of the part of the pattern read so
far (P[0..q]) plus the character just read (c).
Example: P = "bonbons"
S = "bonbonbons"
^ The FSA runs till here (q=6) until it
encounters a mismatch (c='b').
Now, the FSA will not restart matching from P[0], but it will
continue with matching the second 'b' of the pattern with the third
'b' of the string: It knows that the last three characters of the
matched part "bonbon" may be the starting of another match of
"bonbons". Here, the FSA would be
0 b> 1 o> 2 n> 3 b> 4 o> 5 n> 6 s> 7
/ \ ^b'
'b'
with d(q,c) = 0 for all transitions not mentioned and F={7}
This FSA needs the transition function to be stored, occupying a space
of O(S.alphabet*P). To calculate d, a time of
O(S.alphabet*P^3) is needed, if this is done naively. The running
time of the loop is O(S). Although the FSA can the be used to search
in different strings for the same P, this algorithm is rather
inefficient.
MorrisPrattFailureFunction of a pattern P: A function
F:{1,..P}>{0,..P},
which returns for each position i=0..P the length of the longest
prefix of S, which is a suffix of P[0:i).
F(i) := max({j : j=0..i1, P[0:j) == P[ij:i) })
Here, it is max({}):=0. Often, F is stored as a sequence
<0,F(1),F(2),...>.
Example: P="bonbon"
F= _000123
^ The 2 preceding characters of this position
("bo") match with the pattern start.
The MorrisPrattFailureFunction can be calculated as follows:
FUNCTION MorrisPrattFF(P : String) : array of Integer
F : array[P+1] of integer
posInP := 1
prefixLen := 0
F[1] := 0
// Run over P
WHILE posInP < P DO
// Run over a matching with P[0..prefixLen]
WHILE posInP < P && P[posInP] == P[prefixLen] DO
posInP := posInP + 1
prefixLen := prefixLen + 1
F[posInP] := prefixLen
// No matching found, try next char
IF prefixLen == 0 THEN posInP := posInP + 1
// Partial match found, try to use this as a new beginning
ELSE prefixLen := F[prefixLen]
RETURN F
This algorithm does a matching of the pattern in the pattern itself.
It runs over the whole pattern and with each character, tries to match
with the start of P. While this succeeds, F is filled appropriately.
As soon as the complete match fails, prefixLen jumps back to the
longest prefix, which has been found  and this is given by F. During
this algorithm, prefixLen is always smaller than posInP, which allows
the assignment prefixLen := F[prefixLen]. Since posInP only increases,
the running time is in O(P).
MorrisPrattMatchingAlgorithm: The following algorithm, which finds a
pattern P in a string S by the MorrisPrattFailureFunction of P:
FUNCTION MorrisPrattMatching(S,P : String, F : MorrisPrattFF(P))
: Set of Integer
result : Set of Integer := {}
posInS := 0
matchedChars := 0
// Run over the whole string S
WHILE posInS + (PmatchedChars) < S DO
INVARIANT A: P[0:matchedChars) == S[posInSmatchedChars:posInS)
INVARIANT B: occurences before posInSmatchedChars are in
// Run over a sequence of matching characters
WHILE matchedChars < P && P[matchedChars]==S[posInS] DO
matchedChars := matchedChars + 1
posInS := posInS + 1
// All chars matched? > husch husch ins Koerbchen
IF matchedChars == P THEN
result := result \/ {posInS  P}
// No chars matched at all? > try next char in S
IF matchedChars == 0 THEN
posInS := posInS + 1
ELSE
matchedChars := F[matchedChars] // Seek back within pattern
RETURN result
This algorithm goes through the string once. With each character, it
tries to match the beginning of P. When the start of P has been found,
the algorithm tries to match the whole string. While doing so,
counts the number of characters which already match. If
the complete match fails, it gives up the complete match. But F tells
it, how many characters it can nevertheless counted as found, because
while scanning the incomplete pattern, a new pattern may have started
in between. Invariant A says that counts the number of
chars, which match from the beginning of P with the suffix of the
currently read prefix of S. This is maintained trough all statements of
the loop. Invariant B says that no occurences are left out. By
invariant A, an occurence must be signalled by ==P.
This case is cared for by the IFstatement. This algorithm runs in
time O(S), because the difference of posInS  matchedChars increases
with each main loop execution. Having precomputed F in time O(P), F
can be used for a search in different strings. Now it becomes obvious,
why F[i] considers the substring preceding i instead of the substring
including position i: If, during the comparison, a character at
position i mismatches, then one knows that the preceding substring was
all right. If F[i] is used now, it must refer to this preceding
substring.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Hashing
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Key for a type t: A function key:t>tkey, which maps each instance of
the type to an element of a keyset tkey, such that no two instances
get the same key.
ListDictionary of type t: The type
CLASS ListDictionary of T EXTENDS List of T
key : FUNCTION T > tkey
M : Integer := 0 // Size
PROCEDURE insert(e : T)
this.pushFront(e)
M := M + 1
PROCEDURE remove(k : tkey)
IF this==nil THEN RETURN
IF key(this.car)==k THEN
this:=this.cdr
M := M  1
ELSE this.cdr.remove(k)
FUNCTION find(k : tkey) : T
IF this==nil RETURN nil
IF key(this.car)==k THEN RETURN this.car
RETURN this.cdr.find(k)
Dictionary of type t: A type, which has operations equivalent to
those of a listdictionary of type t.
Arraydictionary of type t: An array of type t, which has the
operations of a dictionary.
HashFunction for an arraydictionary of type t: A function
h:tkey>Natural, which maps an element of t given by its key to
a position in the array.
Collision conflict in an arraydictionary: The phenomenon that two
elements are mapped to the same position by the hashfunction.
Closed HashDictionary, HashDictionary with Chaining of type t: An
array dictionary of type t, which maintains a list for each array
position. If two elements are mapped to the same position, they are
simply stored both in the corresponding list. The definition for a
HashDictionary of type t with key key and hashfunction h is
the following:
CLASS ChainingHash of type T
m : Integer // Number of possible values of
key : FUNCTION T > tkey
h : FUNCTION tkey > {0,...m1}
d : array[this.m] of List of T
M : Integer := 0 // Number of elements
PROCEDURE insert(e : T)
d[h(key(e))].pushFront(e)
M := M + 1
PROCEDURE remove(k : tkey)
d[h(key(e))].remove(k) // Using the ListDictionary remove
M := M  1
FUNCTION d.find(k : tkey) : T
RETURN d[h(key(e))].find(k) // Using the ListDictionary find
Inserting needs constant time, whereas removing and deleting may in
the worst case take a time proportional to the number of elements in
the list. But the expected time of removing and finding is O(1), if
the size of the array d.size and the number of elements d.M in it
are about the same. Proof: The expectation of the length of a list
has to be determined. Be X:= d[h(k)].length for any key k, then
X={e in d.M : h(key(e)) == h(k) }
We define the random variable
X(e) := h(key(e))==h(k)?1:0
Then it holds
E(X) = E(SUM {X(e)  e in d.M}) = SUM {E(X(e))  e in d.M}
= d.M*P(Xe=1) = d.M*1/d.size = d.M/d.size
If the array has approximately the size of the number of elements,
then the expected list length reduces to 1.
// Warning: The notions "closed hashing" and "open hashing" are
// sometimes used vice versa! Therefore, this summary uses the
// nonambiguous terms.
Universal Hashfunction family for an arraydictionary d: A set
of hashfunctions, where the probability of two elements being
mapped to the same position is less than or equal to 1/d.size:
P(h(key(e1))==h(key(e2))) =< 1/d.size
"Vectorhashfunction" for a dictionary d: The hashfunction
h[a](x[1],...x[k]) = SUM {x[i]*a[i]  i=0..k} mod d.size
where
* the keyfunction maps each element to a vector of
{0,...d.size1}^k. This is why the input to h is (x[1],...x[k])
* a=(a[1],...,a[k]) is a vector of {0,...d.size1}^k, which
characterizes h
* d.size must be a prime number
The set of all these hashfunctions is universal. Proof: We consider
two keyvectors x=(x[1],...x[k]) and y=(y[1],...y[k]), which differ
in at least one position, say j: x[j]!=y[j]. There are
{0,...d.size1}^k different a and hence {0,...d.size1}^k different
h[a] in this family. We want to count those, which map x and y to
the same position:
h[a](x) = h[a](y)
SUM {x[i]*a[i]} mod d.size = SUM {y[i]*a[i]} mod d.size
0 mod d.size = SUM {a[i]*(y[i]x[i])} mod d.size
a[j]*(y[j]x[j]) mod d.size = SUM {a[i]*(y[i]x[i])  i=0..k, i!=j}
mod d.size
a[j] mod d.size = SUM {a[i]*(y[i]x[i])  i=0..k, i!=j}
/ (y[j]x[j]) mod d.size
Since a[j] {0,...2^201}
h[t](x) := (x mod 2^20) xor t(x div 2^20)
where t is any function
t: {0,...2^121} > {0,...2^201}
Picture: We get a binary input x of 12 higher and 20 lower bits:
xxx12xxxx xxxxxxx20xxxxxxxxxxx __ < A := (x mod 2^20)
  xor
' t > ttttttt20ttttttttttt __ < B := t(x div 2^12)

rrrrrrr20rrrrrrrrrrr <' < A xor B
The set of all these hashfunctions is universal. Proof: We have to
show that for two different given keys x and y and a random t, the
probability of x and y being mapped to the same cell is less than or
equal to 1/d.size, where d.size=2^20:
All x: All y: Random t: x!=y => P(h[t](x)==h[t](y))
There are two cases:
* If the high parts of x and y are equal, then t(x div 2^20) and
t(y div 2^20) will evaluate to the same result. Since the low parts
of x and y are different, this difference will be passed directly to
the result. Thus, the probability of x and y being mapped to the
same cell is 0 in this case.
* If the high parts are different, we need to find those t, which map
x and y to the same cell, i.e.
h[t](x) = h[t](y)
<=> (x mod 2^20) xor t(x div 2^20) = (y mod 2^20) xor t(y div 2^20)
We add an (y mod 2^20) left and right
<=> (y mod 2^20) xor (x mod 2^20) xor t(x div 2^20) =
(y mod 2^20) xor (y mod 2^20) xor t(y div 2^20)
<=> (y mod 2^20) xor (x mod 2^20) xor t(x div 2^20) = t(y div 2^20)
x and y are fixed. Hence, the left part of the equation is also
fixed. Now x and y are mapped to the same cell, iff t(y div 2^20)
takes exactly this fixed value. Since t is chosen randomly and
t(y div 2^20) has the choice among 2^20 values, the probability is
1/2^20, which is 1/d.size.
// Why is the left side of the equation fixed, as it contains t (?)
Thus, the overall probability of the two keys being mapped to the
same cell is less than 1/d.size, the family is universal.
Open HashDictionary, HashDictionary with Probing: An array dictionary
of type t, which seeks for a new free place in the array if an element
is mapped to an occupied cell.
HashDictionary with Linear Probing: A hashdictionary with probing,
which uses the next free array cell whenever a new entry collides with
a present entry. If the last array cell is bypassed, the search
continues at the beginning.
HashDictionary with Bounded Linear Probing: A hashdictionary with
linear probing, which maintains additional space at the end of the
array. Such a dictionary type t is defined as
CLASS BoundedLinearProbingHash of T
m : Integer // Number of possible values of
key : FUNCTION T > tkey
h : FUNCTION tkey > {0,...m1}
key : T
size : Integer // Real size
d : array[this.size] of T
m' : Integer := size  d.m
PROCEDURE insert(e : t)
// The element should go here
pos:=h(key(e))
// Seek for next free place
WHILE this.d[pos]!=nil DO pos := pos+1
IF pos==this.size1 THEN THROW Exception("Too many entries")
// Put it in the free position
this.d[pos]:=e
FUNCTION find(k :tkey) : t
// Should be here:
pos:=h(k)
// Seek pos and subsequent places
WHILE this.d[pos]!=nil DO
IF key(d[pos])==k THEN RETURN this.d[pos]
pos := pos + 1
RETURN nil
PROCEDURE remove(k : tkey)
// First seek the element
pos:=h(k)
IF this.d[pos]==nil RETURN
WHILE key(d[pos])!=k DO
pos := pos + 1
IF this.d[pos]==nil THEN RETURN
// We plan a hole for this.d[pos].
// Seek next element which could fill it
pos' := pos + 1
WHILE this.d[pos']!=nil DO
// Would this.d[pos'] have liked to go to this.d[pos]?
IF h(key(d[pos'])) =< pos THEN
// Put this.d[pos'] where it wants to be
this.d[pos]:=d[pos']
pos := pos'
pos' := pos' + 1
this.d[pos]:=nil
Thus, there is never a hole between the actual position of an element
and the wouldliketobeposition of an element:
All j: j > h(key(d[i])) && j < i => d[j]!=nil
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Minimum Spanning trees
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Weight of a spanning forest of a weighted graph: The sum of the
weights of its edges.
WUC Graph: Weighted undirected and connected graph.
Minimum spanning tree, MST of a WUC graph (V,E,c): A spanning tree,
which has a minimal weight among all spanning trees of (V,E,c). c
must be positive.
Minimum spanning forest of a weighted undirected graph (V,E,c): A
spanning forest, which has a minimal weight among all spanning
forests of (V,E,c).
CutProperty: The fact that the lightest edge of a cut in a WUC graph
can be used in a minimal spanning tree of that graph. Proof: Be (V,E,c)
a WUC graph, be S a cut, i.e. a subset of the nodes.
Be e* the minimal edge. Be R the connected component which is
connected by e* to S. Either e* is in the MST, then the CutProperty
holds anyway. If e* is not in the MST, then any other edge e between
R and S has to be in the MST. But since S and R are both connected
components, both e and e* will do the job of connecting S and R.
Hence, e can be replaced by e* in the MST and the CutProperty holds.
CycleProperty: The fact that the heaviest edge on a cycle in a
WUC graph is not needed in a MST. Proof: Suppose a MST uses the
heaviest edge e* of a cycle. Then any other edge e of this cycle must
miss in the MST. If e* is replaced by e in the MST, the result is
still a MST.
JarnikPrimAlgorithm, PrimAlgorithm: The following algorithm, which
constructs a MST of a WUC graph:
FUNCTION prim((V,E,c) : WUCGraph) : Set of Edge
n := V
// Edge to the predecessor of each node in the MST
pred : array[n] of Edge
// Distance of each node to the MST
dist : array[n] of Integer := <+oo,...>
// Nodes to be examined
Q : Priority Queue of Node with priority dist := <>
// Start with any node
s0 := eps(V)
Q.insert(s0)
WHILE Q != <> DO
s := Q.deleteMin()
dist[s] := 0
FOREACH (s,v) in E DO
IF c(s,v) < dist[v] THEN
dist[v] := c(s,v)
pred[v] := (s,v)
IF v in Q THEN Q.decreaseKey(v) ELSE Q.insert(v)
RETURN { pred[v]  v in V ,v != s0 }
This algorithm works similarly to Dijkstra: It starts with a partial
MST T, which initially only consists of an arbitrary start node s0.
Then, it greedily adds the smalles edge leaving T to T. Thereby, T
grows until all nodes have ben eaten. Similarly to Dijkstra, the
algorithm uses time O(n+m) plus the time of the priority queue, which
is O(m+n*log(n)) all in all with a Fibonaccy Priority Queue (not
treated here) or O((m+n)*log(C)) with a RadixHeap.
UnionFindType: A type which can memorize a partition of a set of
Integers M={1,...n}, where n is a given natural number. A partition
is a set of subsets of M. In each subset, on element is chosen as a
"leader", which represents the subset. For each element, one
establishes parentpointers, lead to the leader of the current subset.
The type is
CLASS UnionFind
parent : array[n] of Integer := <1,2,3,...n>
age : array[n] of Integer := <0,0,...>
FUNCTION find(i : {1,...n}) : {1,...n}
// Returns the leader of the subset, which i is in
IF parent[i]==i THEN RETURN i
i' := find(parent[i])
// Now as we found the leader, store it directly as a shortcut
// in parent, avoiding the search for a future query
// ("Path compression")
parent[i] := i'
RETURN i'
PROCEDURE _link(i, j : {1,...n})
// Links two leaders to a new subset
ASSERT i and j are leaders of different subsets
// Have the "older" leader win
IF age[i] follows for a given element the
parentpointer until it finds the subsetleader (the one which points
to itself). In order to save this search for future calls, the parent
pointer of the element is directly bent to the leader. The procedure
connects two subsets to a new subset. Thereby, it just bends
the pointer of the younger leader. The older leader then becomes even
older. The age of a leader measures the size of its subset. By always
chosing the older leader, it is guaranteed that those elements, which
are in the greater (older) subset do not need an additional step in
, whereas the elements in the smaller set do. But they are less,
which improves the overall runtime. The runtime for all operations
is nearly O(1).
// The slides do not provide a proof for this runtime
UnionFindDatastructure: An instance of the UnionFindType.
KruskalAlgorithm: The following algorithm, which constructs a MST of
a WUC graph:
FUNCTION kruskal((V,E,c) : WUCGraph) : Set of Edges
n := V; m := E
T := UnionFind(n)
result : Set of Edge := {}
Esorted : array[E] of Edge := sortBy(E,c,INCREASING)
FOR i:=0 TO E1 DO
(u,v) := Esorted[i]
IF T.find(u) != T.find(v) THEN
result := result \/ (u,v)
T.union(u,v)
RETURN result
This algorithm grows the MST by joining subtrees. Each subtree is seen
as a subset of V. These subsets are maintained in the UnionFind
Datastructure. First, every node is a tree of its own. The algorithm
selects the lightest edge available and connects the two trivial
subtrees to a new tree. It goes on like this, always chosing the next
edge according to its weight: If it connects two subtrees, it is
absorbed into the tree. If it is a connection within a tree, the
cycleproperty allows discarding the egde. Since and
run in constant time, the dominating time factor is the sorting of
the edges, which makes the overall runtime O(m*log(m)).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Geometry in the Plane
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This is a section on twodimensional objects.
// All points and vectors are twodimensional.
// All sets of points defined here also form a type with the
// appropriate components (e.g. Circle =(radius: Real, center: Point)).
// All of these types have an operation "contains(p :Point)",
// which determines whether p is part of the point set.
// All real values or properties defined here in dependence upon
// some other objects implicitely define a function
// (e.g. "right turn: A sequence of Points x,y,z, such that..."
// defines the function rightturn(x,y,z : Point) : Boolean).
// These functions can be referred to in colloquial English
// (e.g. "x, y and z form a right turn" instead of
// rightturn(x,y,z))
Circular array, circular sequence: A sequence, the last element of
which is its first element and all other element of which are
distinct. For convenience, the last element is left away and one
defines
(x[1],...x[n]) := (x[2],...x[n],x[1])
x[i] := x[(n+i) mod n]
Point: The type (x:Real,y:Real). This type inherits all operations
of vectors.
Coordinate of a point: One of its elements.
Distance of two points p and q: The value pq.
Plane: The set of all twodimensional points. One thinks of the plane
as a blackboard, such that the notions "left" (small xcoordinate),
"right" (large xcoordinate), "down" (small ycoordinate) and
"up" (great ycoordinate) get meaning.
Area: A subset of the plane.
Epsilonenvironment of a point p: The set of all points having a distance
less than epsilon to p: { p'  pp'0: Ex p': (p' not in area && pp'0} is a subset of the
area.
Line segment between two points p and q: The set
{ p+r*(qp)  0 =< r =< 1}
= { r*p+(1r)*q  0 =< r =< 1}
Notation: pq with a bar above them, here: pq
It is pq = qp
Arrow: A pair of points.
Notation: pq with an arrow above them, here: p>q
°: An abbreviation for "*pi/180"
Angle of three points (p,q,r): The value
angle(p,q,r) := (atn((r.yq.y)/(r.xq.x))  atn((p.yq.y)/(p.xq.x))
+ 2*pi) mod (2*pi)
Example:
pq 90°

r
pqr 180°
r

pq 270°
// This definition is quite arbitrary
Turn matrix of three points p,q,r: The matrix
1 p.x p.y
turn(p,q,r) := 1 q.x q.y
1 r.x r.y
Left turn: A sequence of three points (p,q,r), such that
angle(p,q,r)>180°. In this case, det(turn(p,q,r))>0.
Right turn: A sequence of three points (p,q,r), such that
angle(p,q,r)<180°. In this case, det(turn(p,q,r))<0.
Collinear points: A sequence of three points (p,q,r), such that
angle(p,q,r)==180° or angle(p,q,r)==0°. In this case,
det(turn(p,q,r))==0.
Point left of an arrow (p,q): A point r, such that angle(p,q,r)>180°.
Point right of an arrow (p,q): A point r, such that angle(p,q,r)<180°.
Intersection of two line segments: A point which belongs to both
line segments. Given two line segments
{ a+r*(ba)  0 =< r =< 1} and { c+r*(dc)  0 =< r =< 1},
their intersection is the point a+r'*(ba), such that 0=b and a and b
lie on opposite sides of c>d.
Line through the points p and q: The set
{ p+r*(qp)  r in Real }
= { r*p+(1r)*q  r in Real }
Line through the point p in the direction of vector d: The set
{ p+r*d  r in Real }
This is the line through the points p and p+d.
Line with the real parameters a, b, c: The set
{ (x,y)  a*x + b*y = c }
This is the line through (0,c/b) and (1,(ca)/b) and the line
through (0,c/b) in the direction (1, a/b ).
Halfplane with the real parameters a,b,c: The set
{ (x,y)  a*x + b*y =< c }
This is the set of points lying to one side of the line
given by these parameters.
Direction of a halfplane (a,b,c): The vector (a, b). It points to
the direction where the halfplane is unbounded.
Sequence of points in clockwise order: A sequence of points, each of
which forms a right turn with its two successors.
Sequence of points in counterclockwise order: A sequence of points,
each of which forms a left turn with its two successors.
Polygon vertices: A circular sequence P of points, such that no two
line segments p[i]p[i+1] p[j]p[j+1] intersect (except for equal
points) and such that no point is collinear with its two sucessors.
Edge of polygon vertices: A line segment of a vertex and its successor.
Polygon of polygon vertices P: The set
{ q   (+oo,+oo)q /\ E  mod 2 == 1 },
where E is the set of all edge points of P. That is: If I take a point
exterior to the boundary and draw a line to q, then this line crosses
the boundary an odd number of times. A polygon is mostly given by its
vertices. We define the type Polygon as being a circular array.
Circle with radius r around a center point p: A set of points, all of
which have the distance r to p: { q  qp==r }
A circle is mostly given by its radius and center.
Disc of a circle C: The set
{ q   (+oo,+oo)q /\ C  mod 2 == 1 },
that is the "inner area" of the circle. A disc is mostly given by its
circle.
Interior of a circle: Its disc.
Convex area: An area S, such that for any two points p and q in S,
the line segment pq is a subset of S. This corresponds to an area
without "inward corners". The intersection of convex areas is again
convex.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Convex Hull
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Site: A marked point. This notion is used to distinguish certain
given points from the other infinitly many points of the plane.
Convex hull of a set of sites S: The smallest convex area containing
S. If one thinks of the sites as nails in a board, then the convex
hull corresponds to a rubber band (Gummiband) embracing all of them.
If S is finite, then the convex hull is a convex polygon, which is
mostly given by its vertices in clockwise order.
Notation: CH(S)
Reduction of the sorting problem to the convex hull problem: The
following algorithm, which solves a sorting problem if a convex hull
problem can be solved:
FUNCTION sortByCH(A : array of Real) : array of Real
points : array[A.length] of Point
FOR i:=0 TO A.length1 DO
points[i] := (A[i],A[i]^2)
ch : circular array[A.length] of Point := convexHull(points)
minx := ch[0].x
FOR i:=0 TO A.length1 DO
IF ch[i].xq THEN
counter := counter + 1
IF counter==P2 THEN
edges := edges \/ (p,q)
result := array[edges.length] of Point
result[0]:= eps(edges).startNode
FOR i:=1 TO edges1 DO
result[i] := eps({ q  (result[i1],q) in edges})
RETURN result
This algorithm considers each pair of sites and checks whether all
other sites lie right of the line egment defined by the pair. If so,
it considers the pair an edge of the convex hull. Last, the set of
edges is transformed to the sequence of vertices. This algorithm runs
in time O(P^3).
Incremental algorithm: An algorithm, which solves aproblem first for
one constraint and then adds the other constraints one after the
other.
Upper Hull of a set of sites: The subsequence of the convex hull
polygon vertices, in which each site has a larger xcoordinate
that its predecessor.
Graham's Scan: The following incremental algorithm, which computes
the convex hull of a set of sites P:
FUNCTION halfHull(P : array of Point, X : {left, right})
: array of Point
result : List of Point := <>
result.push(P[0])
result.push(P[1])
FOR i:=2 TO P.length1 DO
WHILE result.get(result.length2),
result.get(result.length1),
P[i]
do not form a X turn DO
result.popBack()
result.pushBack(P[i])
RETURN Array(result)
FUNCTION graham(P : Set of Point) : Polygon
P' : array[P] of Point := sortBy(P,x,INCREASING)
upperHull := halfHull(P', right)
lowerHull := halfHull(P', left)
RETURN Polygon(upperHull + reverse(lowerHull))
This algorithm first sorts the sites from left to right, i.e. by
their xcoordinate. It considers the boundary of the convex hull
polygon of P to consist of two parts: The "upper" hull (only right
turns from left to right) and the "lower hull" (only left turns from
left to right). Both are calculated by adding one point after the
other to the hull half. If the new site does not form the correct
turn with its two predecessors, then the predecessor is kicked out.
Then the two hulls are connected. This algorithm runs in time
O(n*log(n)).
b_______c e Adding e will remove f of the
/ \ / convex hull
a/ \/
f
Jarvi's March: The following algorithm, which computes
the convex hull of a set of sites P:
FUNCTION jarvi(P : Set of Point) : Polygon
result : List of Point := <>
p0 : Point := (0,+oo)
// Pick the lowest point
p1 : Point := lowestPoint(P)
P := P \ lowestPoint(P)
// Loop until we are back at the beginning
REPEAT
// Find p2, such that angle(p0,p1,p2) is maximal
p2 := eps(P)
FOREACH p in P DO
IF angle(p0,p1,p) > angle(p0,p1,p2) THEN
p2 := p
P := P \ {p2}
(p0,p1) := (p1,p2)
UNTIL p1==result.get(0)
RETURN Polygon(result)
This algorithm first picks the lowest site p1 of the set. Then it
checks which other site makes up a maximal angle with p1(0,+oo).
This site is chosen as a next site in the vertex sequence. Then
it checks again which remaining site makes up a maximal angle
with the newly added site and the predecessor of it. This continues
until the start site point has been met. If n is the number of sites
and h is the number of convexhullvertices, then this algorithm
runs in time O(n*h).
Tangent of a point p to a connected area A: A line through p, such
that the line has exactly one point q in common with A and the
whole area lies to the right (or left) of p>q. The left tangent
of a convex polygon P (i.e. where the polygon will lie to the right
of it) can be calculated in time O(log(P.length)) as follows:
FUNCTION leftTangentTouchPoint(P : Polygon, p1:Point) : Integer
ASSERT p1 lies to the left of the polygon, P convex
// phi(i) is the angle which P[i] forms with p1 and a
// point in the negative yinfinity
phi : FUNCTION Integer > Real,
phi(i) := angle(P[i],p,(0,oo))
// We need to find the maximum of phi
// Define some useful functions on phi
phiGrows : FUNCTION Integer > Boolean,
phiGrows(i) := phi(i1) < phi(i) < phi(i+1)
phiShrinks : FUNCTION Integer > Boolean,
phiShrinks(i) := phi(i1) > phi(i) > phi(i+1)
isMax : FUNCTION Integer > Boolean,
isMax(i) := phi(i1) < phi(i) > phi(i+1)
isMin : FUNCTION Integer > Boolean,
isMin(i) := phi(i1) > phi(i) < phi(i+1)
// When does P[i]...P[j] contain the maximum?
containsMax : FUNCTION Integer x Integer > Boolean
containsMax(i,j) :=
phiGrows(i) && phiShrinks(j) 
phiGrows(i) && phiGrows(j) && phi(i)>phi(j) 
phiShrinks(i) && phiShrinks(j) && phi(i) i (circular)
We look at two indices i and j. By the help of phi(i) and phi(j),
we can tell whether the maximum lies between i and j: If phi(i)
decreases and phi(j) increases, then the maximum must be in between.
The maximum is also in between, if both phi(i) and phi(j) increase,
but phi(i) > phi(j). A number of other cases applies similarly.
Now we do a binary search with i and j: If the maximum is in between,
we look at the index k:=(i+j)/2. If the maximum is not in between,
we swap i and j. Since the difference ij is always divided by two,
the total running time is O(P.length).
// This is David Steurer's solution. Thanx!
Chan's algorithm: The following algorithm, which computes
the convex hull of a set of sites P:
FUNCTION partialHull(P : Set of Point, m : Integer) : Polygon
r := P / m
// Partition P into r sets of at most m elements
PP : array[r] of Set of Point := <{},...{}>
FOR i:=0 TO r1 DO
FOR j:=0 TO m DO
IF P=={} BREAK
PP[i] := PP[i] \/ {pickAny(P)}
// Compute the convex hulls of each subset
PPhulls : array[r][m] of Point
:=
// Do Jarvi's March on polygons
p0 := (0,+oo)
p1 := lowestPoint(P)
p2 := p0
result : List of Point :=
FOR i:=1 TO m DO
// Find p2, such that angle(p0,p1,p2) is maximal
FOR j:=1 TO r DO
left := PPhulls[j][leftTangentTouchPoint(p1,PPhulls[j])]
right :=PPhulls[j][rightTangentTouchPoint(p1,PPhulls[j])]
IF angle(p0,p1,left) > angle(p0,p1,p2) THEN
p2 := left
IF angle(p0,p1,right) > angle(p0,p1,p2) THEN
p2 := right
IF p2==lowestPoint(P) THEN RETURN Polygon(result)
result.pushBack(p2)
(p0,p1) := (p1,p2)
RETURN nil
FUNCTION chan(P : Set of Point) : Polygon
t := 1
REPEAT
m := min(2^(2^t), P)
result := partialHull(P,m)
UNTIL result != nil
RETURN result
The partialHull function partitions the set of sites into subsets
of size m and computes the convex hull of each of them. By Grahams
scan, this takes time n/m * O(m*log(m)) = O(n*log(m)). Then, it runs
Jarvi's March on these polygons. Given that leftTangentPoint runs in
time O(log(m)), this takes time O(h*n/m*log(m)), where h is the number
of vertices of the final convex hull. If m=h, then the algorithm
runs in time O(n*log(h)). Since h is unknown, it is found by repeated
squaring in the chanfunction: First m=4 is tried, then m=16, m=256,
m=65536 etc. The tth iteration takes time O(n*log(2^(2^t))) =
O(n*2^t). If m>h, then the partialHullfunction will find the hull, so
the iteration is executed at most log(log(t)) times. The overall time
amounts to SUM {n*2^t  t=1...log(log(h))} in O(n*log(h)).
_
_ : Compute the small polygons by Graham, then
/:\ _ _ compute an allenclosing polygon by Jarvi
\:/ /:\
\:/
Polygonmerging: The following algorithm, which merges two disjoint
convex polygons to a new one:
FUNCTION polygonMerge(P1, P2 : Polygon) : Polygon
ASSERT P1 and P2 convex and disjoint, P1 left of P2
// Find upper tangent
i1 : Integer := rightMostPointIndex(P1)
i2 : Integer := leftMostPointIndex(P2)
WHILE not (P2 is right of P1[i1]>P2[i2] &&
P1 is right of P1[i1]>P2[i2]) DO
WHILE not P2 is right of P1[i1]>P2[i2] DO
i2 := i2 + 1
WHILE not P1 is right of P1[i1]>P2[i2] DO
i1 := i1  1
// Find lower tangent
j1 : Integer := rightMostPointIndex(P1)
j2 : Integer := leftMostPointIndex(P2)
WHILE not (P2 is left of P1[j1]>P2[j2] &&
P1 is left of P1[j1]>P2[j2]) DO
WHILE not P2 is left of P1[j1]>P2[j2] DO
j2 := j2 + 1
WHILE not P1 is left of P1[j1]>P2[j2] DO
j1 := j1  1
// Merge the halves and tangents to the result
result : List of Point
REPEAT
result.pushBack(P[j1])
j1 := j1 + 1
UNTIL j1==i1
REPEAT
result.pushBack(P[i2])
i2 := i2 + 1
UNTIL i2==j2
RETURN Polygon(result)
This algorithm assumes P1 to be left of P2. We find the rightmost
vertex i1 of P1 and the leftmost vertex i2 of P2. They so to speak face
each other. Then we move i1 and i2 up alternatedly, until we found
the upper tangent. The condition "P1 lies right of P[i1]>P[i2]" can
be checked by looking if both P1[i11] and P1[i1+1] lie to the right
of P[i1]>P[i2]. We do the same to find the lower tangent and finally
merge the left half of P1, the upper tangent, the right half of P2
and the lower tangent. The runtime is in O(n), where n is the total
number of vertices.
_______c z________ The run will be ax, ay, az,
/ \b y/ / bz, cz, done
   /
\ /a x\_____/
\_____/
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Line segment intersection
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Sweep line: An arrow (i,oo)>(i,+oo), which is moving across
the plane with increasing i. Sweep line techniques are used for
solving problems in the plane. Thereby, the problem has already been
solved to the left of the sweep line.
Event point of a sweep line: A site, which will be part of the sweep
line at some time. The sweep line always jumps from one event point
to the next.
Dead area of a sweep line: An area left of it.
Active area of a sweep line: An area which has points in common with
the sweep line.
Dormant area of a sweep line: An area to the right of the sweep line.
Line intersection sweeping: The following algorithm for calculating
the set of all intersection points of a set of line segments:
FUNCTION lineIntersect(L : Set of LineSegment) : Set of Point
result : Set of Point := {}
CONST INTERSECT : Point := (+oo,+oo)
xstructure : Heap of (point:Point,to:Point)
with priority point.x
// Insert all end and start points of all line segments
// into the xstructure
FOREACH pq in L DO
xstructure.insert(( p ,q))
xstructure.insert(( q ,p))
// All line segments intersecting the sweep line
ystructure : Heap of (point:Point,to:Point)
with priority point.y
WHILE !xstructure.empty() DO
(p,q) := xstructure.deleteMin()
pred := ystructure.pred(p,q)
succ := ystructure.succ(p,q)
// We have an intersection point
IF q==INTERSECT THEN
result := result \/ {p}
ystructure.delete((p,q))
ppred := ystructure.pred(pred)
ssucc := ystructure.succ(succ)
IF ppred intersects with succ to the right of
the sweep line (p.x,oo)>(p.x,+oo) THEN
xstructure.insert((intersectPoint(ppred,succ),INTERSECT))
IF ssucc intersects with pred to the right of
the sweep line (p.x,oo)>(p.x,+oo) THEN
xstructure.insert((intersectPoint(ssucc,pred),INTERSECT))
ystructure.swap(pred,succ)
CONTINUE WHILE
// We have a line start point
IF p.xq.x THEN
IF pred intersects with succ to the right of
the sweep line (p.x,oo)>(p.x,+oo) THEN
xstructure.insert((intersectPoint(pred,succ),INTERSECT))
ystructure.delete(p,q)
RETURN result
Assume no line segment to be vertical and no two line endpoints or
intersection points have the same xcoordinate.
This algorithm moves a sweep line across the plane. All line segment
start and end points (left and right points) are event points. These
event points are stored left to right in a heap called xstructure.
Furthermore, all intersection points are also event points. Since
these are unknown at the beginning, they are calculated on the fly
as the sweep line moves and added to xstructure. All time, a heap
called ystructure stores all line segments intersecting the sweep
line, sorted from above to below. Now the sweep line moves to the
next event point and different things happen:
* If we meet the start of a line, then this line is added to
ystructure and eventual intersections of this line with its
neighbors in ystructure are recorded as event points in
xstructure
* If we meet the end of a line, then we check whether its
neighbors in the ystructure intersect to the right of the sweep
line. If this is the case, we add the intersection point as a
new event point to the xstructure. Then we delete the line
from the ystructure.

i\ When line j ends, i and k are checked
 / for interection and q is found and
\ / added to the xstructure
j_____ \/.....q
 /\
/ \
 \
k/
 sweep line
* If we met an intersection point, we check whether the upper
line will intersect with the successor of the lower line in the
ystructure. We also check whether the lower line will intersect
with the predecessor of the upper line in the ystructure. We add
the intersection points as event points to the xstructure. Then
we interchange the upper and the lower line in the ystructure.

i\  / When the intersection of i and k is met,
j_____\____/______ i and k swap their order in the
\  / ystructure. It is checked whether the
\/ lower line k meets the predecessor of the
X upper line i, which is j here
/\
/  \
k/ 
sweep line
Be n the number of lines and s the number of intersection points.
As each of the heap operations takes time O(log(n)), the total time
amounts to O((n+s)*log(n)).
Proof of correctness: Be x an intersection point of two line
segments pq, rs. Assume that p is left of q and r is left of s.
All points p,q,r,s will be event points. Assume p is left of r, so p
will be in the ystructure when r is inserted. Assume all intersection
points left of the sweep line have already been found (trivially true
in the beginning).
There are three cases:
* There is no line segment reaching into the polygon p,x,r.
As there is no line in between, p and r will be neighbors
in the ystructure. When r is inserted, the intersection with its
neighbors is checked and so x will be found.
* There are line segments reaching into p,x,r. Be k the end point of
the rightmost of these lines. As soon as k is deleted from the
xstructure, it is also deleted from the ystructure. p and r become
neighbors and their intersection is checked, x is found.
* There is a line segment intersecting with pq left of x and to
the right of the rightmost endpoint of any line reaching into p,x,r.
Then this intersection point will have been found and will be an
event point. Before exchanging the two lines in the ystructure,
eventual intersections with the neighbors are checked and x is
found.
If x is found, it is added as an event point to the xstructure and
reported as an intersection.
Polygon intersection sweeping: A modified version of the line
intersection algorithm, which calculates the intersection of
two convex polygons. The edges of the polygons are the line segments
the algorithm works on. The ystructure can only contain 4 line
segements at a maximum. The sweep line starts at the leftmost point
of both polygons. The xstructure does not contain all line end and
startpoints, but just the endpoints of the line segments stored
in the ystructure. Then there is also a fixed number of possible
intersection points. As a result, the xstructure also contains a
constant number of lines. Every endpoint is at the same time a
startpoint. New endpoints are added on the fly. Hence all heap
accesses take constant time. Be n the number of vertices and s the
number of intersections of the polygon edges, then the total time
is O((n+s)*1) = O(n). Knowing the intersection points in the order
from left to right, the polygon intersection can be calculated.
HalfPlaneIntersectionAlgorithm: The following algorithm, which,
given a set of halfplanes, determines their intersection:
FUNCTION hpintersect(H : Array of HalfPlane)
: Set of LineSegment
IF H==1 THEN RETURN { LineSegment(Line(H[0])) }
n2 := H.length/2
RETURN polygonIntersection_(hpintersect(H[0..n2]),
hpintersect(H[n2..H.length1])
This is a recursive divideandconquerapproach, which relies on
the fact that the polygonintersection sweeping algorithm also works
for "unbounded polynoms", if adjusted appropriately. If n is the number
of halfplanes, then we have log(n) recursive calls, each of which
takes time O(n) for the intersection, resulting in a total time of
O(n*log(n)).
Binary tree: A tree, all nodes of which have 0, 1 or 2 children.
Checktree: The following type
CLASS CheckTree
n : Integer
depth : Integer
tree : array[depth][n] of (mean:Real,checkedLeaves:Integer)
PROCEDURE check(r: Real, b:Boolean)
// Checks a leaf
pos:=0
FOR level:=0 TO depth1 DO
tree[level][pos].checkedLeaves +:= b?1:1
IF r>tree[level][pos].mean THEN pos:=pos*2+1
ELSE pos:=pos*2
FUNCTION inbetween(r,s:Real):Integer
// Returns the number of checked leaves betwen r and s, inclusive
result:=0
pos:=0
FOR level:=0 TO depth2 DO
IF r>tree[level][pos].mean THEN
// Go right
pos:=pos*2+1
ELSE
// Go left, add right checked leaves
pos:=pos*2
result:=result+tree[level+1][pos*2+1].checkedLeaves
IF tree[depth1][pos].mean==r THEN
result := result + tree[depth1][pos].checkedLeaves
FOR level:=0 TO depth2 DO
IF s>tree[level][pos].mean THEN
pos:=pos*2+1
ELSE
pos:=pos*2
result:=resulttree[level+1][pos*2+1].checkedLeaves
RETURN result
FUNCTION makeCheckTree(S : Set of Real)
// Returns a CheckTree for S
t : FakeTree
t.depth:= ceil(log2(S))
t.n:=2^t.depth
// Sort the reals
sorted : array[t.n] of Real := sort(S,INCREASING)
// The sorted numbers are the leaves
FOR i:=0 TO t.n1 DO
result.tree[t.depth1][i]:=(sorted[i],0)
// Create a tree over them
FOR level := t.depth2 TO 0 DO
FOR node := 0 TO 2^level DO
t.tree[level][node].mean := (t.tree[level+1][node*2]+
t.tree[level+1][node*2+1])/2
t.tree[level][node].checkedLeaves := 0
RETURN t
This CheckTreeclass serves only for the following OrthogonalLine
Sweepingalgorithm. It can handle a set of reals S, some of which
are "checked". It sorts all numbers of S and builds a binary tree over
them. Each node of the binary tree stores the mean value of its
children, so that a leaf can be found by passing from the root to the
leaves. Each node also stores the number of checked leaves below it.
If a number is to be checked or unchecked, the corresponding leaf is
found and the number of checked leaves is decreased on the way from
the root to the leaf. If the number of checked reals between two
reals r and s is to be determined, we run down the tree from the root
to r and sum up all checked leaves of all right children, whenever we
chose a left child. Thereby, we count all checked leaves right to r.
Now we do the same again for s and sum up all checked leaves right to
s. The difference is the number of intermediate checked leaves. All
operations take time O(log(S)).
OrthogonalLineSweeping: The following algorithm, which
calculates the number of intersections between a set of horizontal line
segments and a set of vertical line segments:
FUNCTION OlineIntersect(H,V : Set of LineSegment) : Integer
result : Integer := 0
xstructure : PriorityQueue of (from, to:Point) with
priority from.x
// Insert all line points into the xstructure
FOREACH pq in H \/ V DO
xstructure.insert((p,q))
xstructure.insert((q,p))
// All horizontal line segments will be checked
ystructure : makeCheckTree({y  y is a ycoordinate in H})
WHILE !xstructure.empty() DO
(p,q) := xstructure.deleteMin()
// Hline start
IF p.xq.x THEN
ystructure.check(p.y,false)
// Vline
IF p.x==q.x THEN
result := result + ystructure.inbetween(p.y,q.y)
RETURN result
This algorithm moves a sweep line across the plane like the
line intersection sweeping algorithm. All line points are
event points. When a horizontal line is met, its ycoordinate
is stored in the ystructure. This entry is deleted, when the
horizontal line is over. When a vertical line is met, the number of
horizontal lines between start and endpoint of the vertical
line is added up. The time is O((V+H)*log(V+H)).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Linear Programs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
LP, Linear Program: A pair of
* A maximize expression:
c*x
* A set of constraints
d[1]*x =< d'[1]
d[2]*x =< d'[2]
...
where x is a vector variable, the d[i] are constant vectors and the
d'[i] are real constants. The dimension of x is called the dimension
of the linear program. Here, we will deal with twodimensional linear
programs. The following variations are allowed:
* One may also use the relation '>=', since it can be expressed
by multiplying the line by 1 and using '=<'.
* One may also use the relation '=', as it can be expressed by a line
with '=<' and a line with '>='.
* One may express the lines as sums of the form
SUM {d[i][j]*x[j]  i=1..} =< d'[i]
Feasible region of a linear program: The set of all points (x[0],x[1]),
which satisfy the constraints of the linear program.
Pointer of a linear program: The vector (c[1],c[2]). It holds: For any
point p, p+v has a greater value for the maximize expression than p.
The value of a point (x,y) for the maximize expression is
c[1]*x + c[2]*y = (x,y)*(c[1],c[2])
Solution of a linear program: The point, which lies in the feasible
region and the coordinates of which maximize the maximize expression.
If the feasible region is empty, a solution does not exist.
If the feasible region is unbounded in the direction of the pointer,
the solution is infitity. A solution of a linear program can be found
in polynomial time. A typical execution time is O((m+n)*N), where n
is the dimension of the linear program, m is the number of constraints
and N is the number of nonzero components of the vectors d[i].
// Algorithm treated later
Lineaddingalgorithm: The following incremental algorithm, which
solves a linear program:
FUNCTION linProg(H : Array[n] of HalfPlane, pointer : Vector)
: Point
// Find two halfplanes, which bound the feasible region in the
// direction of p, put them to the beginning of H
j := 0
FOR i:=0 TO n1 DO
// If there is no upper bound, we're done
IF i==H.length THEN RETURN +oo*pointer
Be ph a point on the boundary of H[i]
IF not ph+pointer in H[i] THEN
// h bounds feasible region
swap(H[j],H[i])
j := j + 1
IF j==2 THEN BREAK FOREACH
// tentative result: Intersection of these two upper bounds
result : Point := intersection(Line(H[0]),Line(H[1]))
// Store for every halfplane its boundary with the feasible
// region by two points and . We order these points, such
// that a*pointer > b*pointer.
Hboundary : array[n] of (a,b:Point)
Hboundary[0]:= (result, )
Hboundary[1]:= (result, )
// Add other halfplanes, check whether result is not OK
FOR i:=2 TO n1 DO
IF not result in H[i] THEN
// Check intersection with H[i] with all previous ones
bestsofar := oo
Hboundary[i]:=(,
)
FOR j:=0 TO i1 DO
intsct : Point = intersection(Line(H[i]),Line(H[j]))
IF intsct in Hboundary[j].aHboundary[j].b &&
intsct*pointer < result*pointer THEN
bestsofar := intsct*pointer
result := intsct
Hboundary[j]:= (intsct,Hboundary[j].b)
IF intsct*pointer > Hboundary[i].a*pointer THEN
Hboundary[i]:=(intsct,Hboundary[i].b)
ELSE
Hboundary[i]:=(Hboundary[i].a,intsct)
// If we didn't find anything, there is no solution
IF not result in H[i] THEN RETURN <>
RETURN result
This algorithm first finds two halfplanes, which bound the feasible
region in direction of the pointer. If there is none, the problem is
solved. We take the intersection point of the upper bounds as a
tentative result. Then we add the other halfplanes one after the
other. If the tentative result point lies within the newly added
halfplane, we can keep the point. Else, we check all intersections of
the new line with the previous lines, as these may become new results.
We pick the best among them. We only take those intersection points
into account, which lie on a boundary with the feasible region.
In the best case, this algorithm has a runtime of O(n): The lines are
added in such an order, that the seek for intersection points is never
necessary. In the worstcase, it has runtime of O(n^2), because each
of the newly added lines may require a check with all previous lines.
i\ /j For {i,j}, p was the optimum. Now k is
\/......p added. Since p does not lie in the half
^ /\ space blow k, the intersections of k with
 c ___/__\_____k all previous lines have to be checked.
/ \
Shuffling: The following algorithm:
FUNCTION shuffle(A: array) : array
FOR i := A.length1 DOWNTO 1 DO
swap(A[i],A[random(i)]) // random(i) in {0...i1}
RETURN A
This algorithm creates a randomized permutation of the array. The
probability of any given element being at any given position is
1/A.length.
Randomized algorithm: An algorithm, which, for the same inputs, may
do different things.
Expected runtime of a randomized algorithm: Its average runtime for
an infinite number of runs on the same inputs.
Backwardanalysis: The following technique for determining the
runtime of an incremental randomized algorithm: Instead of adding
a constraint and asking what happens with what probability, one
assumes a situation with the first constraints already given and
asks what happens with what probability if the last constraint is
removed.
Randomized lineaddingalgorithm: The linealgorithm, which shuffles
the array of halfplanes first. This shuffling excludes the first
two halfplains, which have been found to be upper bounds. The
shuffling causes the algorithm to run in expected time O(n): If the
halfplains are added in random order, chances are good that
* either the tentative result point already satisfies the newly
added constraint or
* if it does not, there are few previous lines with which to check
intersections.
Proof: We do a backward analysis. Assume the halfplanes H[0]...H[i]
to be given and p to be the tentative result point. Now we remove
H[i] and imagine the situation before H[i] was added. Was the previous
p different from the new p? p always lies on an intersection of two
lines. These lines are said to determine p. If by chosing H[i], we
incidentally chose a line which determined p, then this means p was
different before H[i] was added. What is the probability that by
removing H[i], we remove a line determining p? Only two lines out of
the i lines determine p, hence the probability is 2/(i2). (We
subtract 2 to take into account the first two predetermined lines).
With this probability of 2/(i2), we will have to check all (i1)
previous lines for intersection points with H[i]. The expected
running time is thus
O(n)+ for finding the two upper bound lines
O(SUM {2/(i2)*(i1)  i=3..n}) for the intersection checks
= O(n).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Smallest enclosing Disc & Duality
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Smallest enclosing disc of a set of sites: A disc, which contains
all sites and has the smallest radius of all these discs. This
disc can be computed by the following incremental algorithm:
FUNCTION smallestDisc3(q1, q2, q3 : Point) : Disc
// Computes the smallest disc containing q1, q2, q3
p : Point := (q1+q2+q3)/3
RETURN (radius: q1p, center: p)
FUNCTION smallestDisc2(P : array of Point, q1,q2 : Point) : Disc
// Computes the smallest disc containing P and with
// q1, q2 on its boundary
D : Disc := (radius:q1q2, center:(q1+q2)/2)
FOR i := 0 TO P.length1
IF P[i] is not contained in D THEN
D := smallestDisc3(q1,q2,P[i])
RETURN D
FUNCTION smallestDisc1(P : array of Point, q1 : Point) : Disc
// Computes the smallest disc containing P and with
// q1 on its boundary
D : Disc := (radius:q1P[0], center:(q1+P[0])/2)
FOR i := 1 TO P.length1
IF P[i] is not contained in D THEN
D := smallestDisc2(P[0..i1],P[i])
RETURN D
FUNCTION smallestDisc(P : array of Point) : Disc
P := shuffle(P)
// Preliminary result
D : Disc := (radius:P[0]P[1]/2, center:(P[0]+P[1])/2)
FOR i:=2 TO P.length1 DO
IF P[i] is not contained in D THEN
// D is the disc containing P[0..i1] with P[i]
// on its boundary
D := smallestDisc1(P[0..i1], P[i])
RETURN D
This algorithm starts off with a function to compute the smallest
disc of three sites. The next function computes the smallest
disc for a set of sites, where it is already known that two sites
q1 and q2 lie on its boundary. It just computes the smallest disc
of q1,q2 and the first site of the set. If any other site p is not
in this disc, then the disc around q1,q2,p is taken instead.
Does this work(?) Suppose
q1
p1 p2
q2
First, the circle goes through q1 and q2. It contains p1, so the
tentative circle does not change for p1. Now, p2 is added and the
circle is computed as smallestDisc3(q1,q2,p2), which may grow so
large, that the circle line between q1 and q2 might become
straight, p1 is excluded from the disc! This doubts lemma 1(iii).
The next function computes a smallest enclosing disc for a set of
sites, where it is already known that one site q1 lies on its
boundary. The sites of the set are added incrementally and whichever
site is not in the disc must be on its boundary. In this case, the
previous function is called. The main function does a shuffling and
proceeds according to the same pattern. The runtimes are as follows:
* smallestDisc3: O(1)
* smallestDisc2 with i points: O(i)
* smallestDisc1 with i points:
It has to be calculated how often smallestDisc1 calls
smallestDisc2. smallestDisc2 is called for j1 points, if P[j]
does not lie on the tentative disc. We do a backwards analysis:
We suppose a given state j with a disc containing all sites.
Now we remove the last site P[j]. What is the probability that
the disc was different before? The smallest disc was only
different if P[j] is one of the three sites on the boarder
of the current disc. This is the case with a probability of
3/j. Since in this case, smallestDisc2 will have been called
before, we can compute the running time of smallestDisc1 as
SUM {O(j)*3/j + 1*(13/j)  j=2..i } = O(i)
* smallestDisc with n points: O(n) by a similar analysis
Dual point of a line y=a*xb: The point(a,b).
Notation for the dual point of a line L: L*
Dual line of a point (a,b): The line y=a*xb.
Notation for the dual line of a point P: P*
Selfinverselemma: The fact that (p*)=p.
Proof: (p*)* = ({(x,y)y=p.x*xp.y})* = (px.py) = p
Orderreversinglemma: The fact that iff point p lies above line L, then
point L* lies also above line p*. This also applies analogously to the
inverted situation and if p lies on L.
Proof: Suppose p=(p.x,p.y) lies above L:y=la*xlb
p.y>la*p.xlb
<=> lb > p.x*lap.y
...which means that L* lies above p*
Incidencepreservinglemma: The fact that the lines L1 and L2
intersect at point p, iff the line p* passes through the points L1*
and L2*.
Proof: Suppose lines L1 and L2 intersect in point p.
p is on L1 && p is on L2
<=> L1* is on p* && L2* is on p* (by orderreversinglemma)
Thus p* passes through L1* and L2*.
Collinearity/Coincidencelemma: The fact that three points are
collinear, iff their dual lines intersect in a common point.
If the three points are exactly one above the other, the lemma only
holds if one accepts parallel lines to intersect in (+oo,+oo).
Proof: Suppose lines L1, L2, L3 intersect in point p.
p is on L1 && p is on L2 && p is on L3
<=> L1* is on p* && L2* i on p* && L3* is on p*
(by orderreversinglemma)
Thus, L1*, L2* and L3* are collinear.
Upper envelope of lines: The boundary of intersection of the
halfplanes above these lines.
LowerHull/UpperEnvelopeTheorem: The fact that the lefttoright
order of the lower hull of a set of sites is the same as the
lefttorightorder of the sequence of their dual lines on their
upper envelope.
Proof: Suppose p[i]p[j] is in the lower hull of a set P of sites.
all p lie above p[i]p[j]
<=> all p* lie below (p[i]p[j])* (by orderreversinglemma)
<=> all p* lie below p[i]* /\ p[j]* (by incidencepreservinglemma)
<=> all dual lines below intersection
of p[i]* and p[j]*
This means that p[i]* and p[j]* form two adjacent edges of the
upper envelope. Suppose p[i] is left of p[j]. Then p[i]* will have
a smaller slope than p[j]*. Hence the contribution of p[i]* to the
upper envelope must lie left of the contribution of p[j]*.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
VoronoiDiagrams
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Binary tree structure of type T with priority f:T>Real: A heap of
type T with priority f:T>Real with all its operations, basing
on the following type structure:
CLASS BinaryTree of type T with priority f:T>Real
left, right : BinaryTree of type T with priority f:T>Real
element : T
Balanced binary tree: A binary tree, where no leaf is on a level
smaller than the maximum level  2 (?).
Balanced binary tree structure: A binary tree structure, which with
each delete or insertoperation reorganizes the tree such that
the corresponding binary tree is balanced. As a result, all
operations always take time O(log(n)), if n is the number of
nodes in the tree.
Box with the edge points p and q: The set
p(p.x,q.y) \/ p(p.y,q.x) \/ q(p.x,q.y) \/ q(p.y,q.x).
Noncrossing line segments: Two line segments, which either have no
point in common or just one endpoint.
Curve: A sequence of (possibly infinitely small) noncrossing
line segments, each of which has the startpoint identical to the
endpoint of its predecessor  except for the first one.
Planar graph: An undirected graph, the nodes of which are points
and the links of which are curves between the points. The
curves may only have the nodes in common.
Straight Planar graph: A planar graph, all edges of which are line
segments. There are additional linesegments (edges) allowed, which
start at one node and are unbound in one direction. These lines are
usually dealt with by adding a big box around the graph, such that
these lines can be given by their start point and their
intersection point with the box.
Face of a straight planar graph: A (possibly unbounded) polygon of its
nodes, which does not include other nodes.
Euler's formula: The fact that in a straight planar graph the number of
nodes plus the number of faces minus the number of edges equals
two: nodes+facesedges=2
This implies that the expected degree of any vertex is 6. (?)
Cocircular points: Points, for which there is a circle, such that
they all lie on it.
General position assumption for a set of points P: The assumption
that no 4 points of P are cocircular.
Voronoicell of a site q in a set of sites P: The set of all points
in the plane, which are closer to q than to another site.
V(q) := {x  qx =< p[i]x for all p[i] in P}
A voronoicell is a possibly unbounded convex polygon and given
by the sequence of its vertices.
Voronoidiagram of a set P of sites: The straight planar graph made up
by the set of all Voronoicells of P. The general position assumption
will be made. The following properties hold:
* Each point on an edge between two Voronoicells V(p[i]) and
V(P[j]) has the same distance to the sites p[i] and p[j]
* A point, which has the same distance to three sites is a
vertex of the Voronoidiagram. It is the center of the circle
passing through the three sites with no other site inside the
circle.
* The vertices of the Voronoidiagram all have degree 3.
* The Voronoidiagram has P=:n faces, less than 2*n5
vertices and less than 3*n6 edges.
Notation: Vor(P)
\  . /
\ . __/ .
\___/ . \_______
. / . \___/ .
\___/ . \
/ .  \
/ 
Naive Voronoi algorithm: An algorithm, which computes the
Voronoidiagram of n sites in time O(n^2*log(n)) by computing the
intersection of n halfplanes per site. (?)
Distance of a site to a line: The minimum among the distances of
all points on the line to the site.
Parabola of a site s and a sweepline: The set of all points of the
plane, which have the same distance to s as to the sweepline.
If we have a horizontal sweepline at ycoordinate ly, then the
parabola is given by
{(x,y)  (yly)^2 = (xs.x)^2+(ys.y)^2}
It is assumed that the sweep line moves downwards.
Beach line of sites S and a horizontal sweepline: The set of
all those points on the parabolae, which do not have another
parabola point below them.
Breakpoint of a beach line: A point of a beach line, which lies on two
parabolae. This point has the same distance to two sites. Hence, it
will lie on a Voronoiedge. If it lies on three parabolae at the same
time, then it has the same distance to three points. Hence, it must
be a Voronoivertex. As the sweepline moves down, the breakpoints
will trace the edges of the Voronoidiagram. Whenever three parabolae
meet, the breakpoint becomes a Voronoivertex. This happens exactly,
when the sweep line touches the bottom point of the circle through
these three sites, because in this moment, any point which is to
come further down will have a greater distance to the three sites than
the three sites have among themselves.
\ / x has the same distance to p and q. y has
\ x q y / the same distance to q and r. Both trace a
\ p /\___/\ r / Voronoiedge. If x and y meet and q's arc
\___/ \___/ disappears, a Voronoivertex is found.
SL_______________________
Adjacent breakpoints: Two breakpoints b1 and b2, such there is no
breakpoint with an xcoordinate between b1.x and b2.x.
Arc of a beach line: The set of points on the beach line,
which have an xcoordinate between the xcoordinates of two given
adjacent breakpoints.
Voronoisweeping: The following algorithm, which uses a beach line
to calculate the Voronoidiagram of a set of points:
// This is only a very restricted algorithm (returns only vertices),
// which has to involve plenty of pseudocode due to its complexity.
// Sorry for that...
FUNCTION voronoi(P : Set of Point) : List of Vertex
result: List of Vertex := <>
// This structure stores the arcs of the beachline in the
// leaves of a balanced binary tree. Such an arc is stored
// in 'arc' by the corresponding site with x giving an
// estimate of its xposition
beachLine: Heap of (x:Real,arc:Point)
with priority x := <>
events: PriorityQueue of (p:Point,type:{SITE,VERTEX},
center:Point,arcIds:Set of Real) with the
priority p.y := <>
// Insert all sites into the eventsstructure
FOREACH p in P DO
events.insert((p,SITE,nil,{}))
// Run over first site
(p,_,_) := event.deleteMin()
beachline.insert((p.x,p))
WHILE !events.empty() DO
(p,type,center) := events.deleteMin()
IF type==SITE THEN
beachLine.insert((p.x,p))
// Get my two neighbors
pred:=beachLine.pred((p.x,p))
succ:=beachLine.succ((p.x,p))
// Find the closest arc among them
IF ppred.arc<psucc.arc THEN
// All vertexevents involving pred.arc must die...
// Patch a piece of pred right to p
beachLine.insert(( (p.x+succ.x)/2, pred)
ELSE
// All vertexevents involving succ.arc must die...
// Patch a piece of succ left to p
beachLine.insert(( (p.x+pred.x)/2, succ)
// Find circles with predecessors and successors
predpred:=beachLine.pred(pred)
succsucc:=beachLine.succ(succ)
disc:=smallestDisc3(predpred.arc,pred.arc,p)
// Is this a good one?
IF disc.center.ydisc.radius
beachLine.delete(above)
result := result \/ center
RETURN result
This algorithm works as follows (if it works at all in this
version...): All sites are events. The beach line is kept in a heap
(in the normal algorithm, it is kept in a balanced binary tree). When
the beach line meets a site p2, it patches the beach line as follows:
p1 p1
\____/ ~~~~~~> \_____/
\/
sweep line p2 p2
beachLine: ....p1.... ...p1,p2,p1....
Moreover, all vertex events which have been introduced on basis of the
arc, which is just being split up, have to be removed: The arc can no
longer contribute to any voronoi cell, because p2 is i its way.
Furthermore, it is checked whether the sites of the neighboring arcs
of the new arc make up a circle, which crosses the sweep line. If this
is the case, then the bottom of this circle is added as a vertex event.
If the sweep line reaches such a vertex event, then no further site
below the sweep line can influence the area of the circle. Hence, its
center is declared a vertex. Moreover, the arc directly above this
vertex is deleted, as it contributes no longer to the beach line.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DelaunayTriangulations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Triangle: A polygon with 3 vertices.
Quadrilateral: A polygon with 4 vertices. As for polygons,
quadrilateral vertices are given in clockwise order.
Skinny triangle: A triangle with one of its angles being small.
Triangulation of a set of sites P: A maximal straight planar graph with
* the sites as its vertices
* line segments between the sites as edges
* each bounded face being a triangle
Every triangulation of P has the same number of triangles.
Good triangulation: A triangulation with very few skinny triangles.
Good triangulations improve interpolated function values, if only the
function values of the vertices are known.
Example: One knows the height of certain landscape locations and one
would like to estimate the height of other locations. One can find a
triangulation among the known locations and calculate the heights on
the edges as weighted meanvalues of their end points. This
aproximation will be more realistic with good triangulations.
Angle vector of a triangulation T: The tuple of all angles of all
triangles of T, sorted in ascending order. One defines
(a1,a2,...am) < (b1,b2,...bm) :<=> a1 < b1 
a1==b1 && (a2,...am)<(b2,...bm)
Good triangulations have a large angle vector.
DelaunayTriangulation of a set of sites P: A triangulation of P with
a maximum angle vector. The general position assumption
will be made. The following properties hold:
* The triangulation is a kind of inverse of the VoronoiDiagram: The
edges of the DelaunayTriangulation connect exactly those points of P
separated by one Voronoiedge.
* Three sites p1, p2, p3 are the vertices of a Delaunaytriangle
iff the circle through them contains no site of P in its interior.
* Two sites form a Delaunayedge, iff there is a circle with the two
sites on its boundary and no site in its interior.
Notation: Del(P)
Additional information: The convex hull can also be calculated for
a set of 3dimensional sites. Then the convex hull is a set of
twodimensional polygons, which fit together like the faces of a
football. By the help of a 3dimensional convexhullalgorithm, one
can calculate the DelaunayTriangulation of a 2dimensional set of
points: Just imagine a paraboloid structure flying above the plane.
This is as if someone placed a saladbowl on the plane. Now have all
2dimensional sites fly vertically upwards until they hit the
saladbowl. Calculate the convex hull of these 3dimensional points.
This gives a set of triangles between the sites, sticking on
the saladbowl. Now have these triangles fall back on the plane: This
is the DelaunayTriangulation. Proof: Consider the three sites of a
triangle on the saladbowl and draw a circle through them. There
cannot be another site inside it, because all points interior to the
circle have a greater distance to the saladbowl center than the
points on its boundary. Thus, it does not lie to the interior side of
the triangle, which contradicts the principle of the convex hull.
Illegal edge of a quadrilateral