Competitive Programming in Python
Competitive Programming in Python
Competitive Programming in Python
in Python
128 Algorithms to Develop Your Coding Skills
CHRISTOPH DÜRR
CNRS, Sorbonne University
JILL-JÊNN VIE
Inria
www.cambridge.org
Information on this title: www.cambridge.org/9781108716826
DOI: 10.1017/9781108591928
© Cambridge University Press 2021
Translation from the French language edition:
Programmation efficace - 128 algorithmes qu’il faut avoir compris et codés en Python au cour de sa vie
By Christoph Dürr & Jill-Jênn Vie
Copyright © 2016 Edition Marketing S.A.
www.editions-ellipses.fr
All Rights Reserved
This publication is in copyright. Subject to statutory exception
and to the provisions of relevant collective licensing agreements,
no reproduction of any part may take place without the written
permission of Cambridge University Press.
First published 2021
Printed in the United Kingdom by TJ Books Limited, Padstow Cornwall
A catalogue record for this publication is available from the British Library.
Library of Congress Cataloging-in-Publication Data
Names: Dürr, Christoph, 1969– author. | Vie, Jill-Jênn, 1990– author. |
Gibbons, Greg, translator. | Gibbons, Danièle, translator.
Title: Competitive programming in Python : 128 algorithms to develop your
coding skills / Christoph Dürr, Jill-Jênn Vie ; translated by Greg
Gibbons, Danièle Gibbons.
Other titles: Programmation efficace. English
Description: First edition. | New York : Cambridge University Press, 2020.
| Includes bibliographical references and index.
Identifiers: LCCN 2020022774 (print) | LCCN 2020022775 (ebook) |
ISBN 9781108716826 (paperback) | ISBN 9781108591928 (epub)
Subjects: LCSH: Python (Computer program language) | Algorithms.
Classification: LCC QA76.73.P98 D8713 2020 (print) | LCC QA76.73.P98
(ebook) | DDC 005.13/3–dc23
LC record available at https://lccn.loc.gov/2020022774
LC ebook record available at https://lccn.loc.gov/2020022775
ISBN 978-1-108-71682-6 Paperback
Cambridge University Press has no responsibility for the persistence or accuracy of
URLs for external or third-party internet websites referred to in this publication
and does not guarantee that any content on such websites is, or will remain,
accurate or appropriate.
Contents
Preface page ix
1 Introduction 1
1.1 Programming Competitions 1
1.2 Python in a Few Words 5
1.3 Input-Output 13
1.4 Complexity 17
1.5 Abstract Types and Essential Data Structures 20
1.6 Techniques 28
1.7 Advice 37
1.8 A Problem: ‘Frosting on the Cake’ 39
2 Character Strings 42
2.1 Anagrams 42
2.2 T9—Text on 9 Keys 43
2.3 Spell Checking with a Lexicographic Tree 46
2.4 Searching for Patterns 48
2.5 Maximal Boundaries—Knuth–Morris–Pratt 49
2.6 Pattern Matching—Rabin–Karp 56
2.7 Longest Palindrome of a String—Manacher 59
3 Sequences 62
3.1 Shortest Path in a Grid 62
3.2 The Levenshtein Edit Distance 63
3.3 Longest Common Subsequence 65
3.4 Longest Increasing Subsequence 68
3.5 Winning Strategy in a Two-Player Game 70
4 Arrays 72
4.1 Merge of Sorted Lists 73
4.2 Sum Over a Range 74
4.3 Duplicates in a Range 74
4.4 Maximum Subarray Sum 75
vi Contents
5 Intervals 82
5.1 Interval Trees 82
5.2 Union of Intervals 85
5.3 The Interval Point Cover Problem 85
6 Graphs 88
6.1 Encoding in Python 88
6.2 Implicit Graphs 90
6.3 Depth-First Search—DFS 91
6.4 Breadth-First Search—BFS 93
6.5 Connected Components 94
6.6 Biconnected Components 97
6.7 Topological Sort 102
6.8 Strongly Connected Components 105
6.9 2-Satisfiability 110
10 Trees 171
10.1 Huffman Coding 172
10.2 Lowest Common Ancestor 174
10.3 Longest Path in a Tree 178
10.4 Minimum Weight Spanning Tree—Kruskal 179
11 Sets 182
11.1 The Knapsack Problem 182
11.2 Making Change 184
11.3 Subset Sum 185
11.4 The k-sum Problem 187
13 Rectangles 200
13.1 Forming Rectangles 200
13.2 Largest Square in a Grid 201
13.3 Largest Rectangle in a Histogram 202
13.4 Largest Rectangle in a Grid 204
13.5 Union of Rectangles 205
13.6 Union of Disjoint Rectangles 212
16 Conclusion 245
16.1 Combine Algorithms to Solve a Problem 245
16.2 For Further Reading 245
16.3 Rendez-vous on tryalgo.org 246
material of this text. Thanks to all those who proofread the manuscript, especially
René Adad, Evripidis Bampis, Binh-Minh Bui-Xuan, Stéphane Henriot, Lê Thành
Dũng Nguyễn, Alexandre Nolin and Antoine Pietri. Thanks to all those who improved
the programs on GitHub: Louis Abraham, Lilian Besson, Ryan Lahfa, Olivier Marty,
Samuel Tardieu and Xavier Carcelle. One of the authors would especially like to thank
his past teacher at the Lycée Thiers, Monsieur Yves Lemaire, for having introduced
him to the admirable gem of Section 2.5 on page 52.
We hope that the reader will pass many long hours tackling algorithmic problems
that at first glance appear insurmountable, and in the end feel the profound joy when
a solution, especially an elegant solution, suddenly becomes apparent.
Finally, we would like to thank Danièle and Greg Gibbons for their translation of
this work, even of this very phrase.
Attention, it’s all systems go!
1 Introduction
You, my young friend, are going to learn to program the algorithms of this book,
and then go on to win programming contests, sparkle during your job interviews,
and finally roll up your sleeves, get to work, and greatly improve the gross national
product!
Mistakenly, computer scientists are still often considered the magicians of modern
times. Computers have slowly crept into our businesses, our homes and our machines,
and have become important enablers in the functioning of our world. However, there
are many that use these devices without really mastering them, and hence, they do not
fully enjoy their benefits. Knowing how to program provides the ability to fully exploit
their potential to solve problems in an efficient manner. Algorithms and programming
techniques have become a necessary background for many professions. Their mastery
allows the development of creative and efficient computer-based solutions to problems
encountered every day.
This text presents a variety of algorithmic techniques to solve a number of classic
problems. It describes practical situations where these problems arise, and presents
simple implementations written in the programming language Python. Correctly
implementing an algorithm is not always easy: there are numerous traps to avoid and
techniques to apply to guarantee the announced running times. The examples in the
text are embellished with explanations of important implementation details which
must be respected.
For the last several decades, programming competitions have sprung up at every
level all over the world, in order to promote a broad culture of algorithms. The prob-
lems proposed in these contests are often variants of classic algorithmic problems,
presented as frustrating enigmas that will never let you give up until you solve them!
a web interface, where it is compiled and tested against instances hidden from the
public. For some problems the code is called for each instance, whereas for others the
input begins with an integer indicating the number of instances occurring in the input.
In the latter case, the program must then loop over each instance, solve it and display
the results. A submission is accepted if it gives correct results in a limited time, on the
order of a few seconds.
Figure 1.1 The logo of the ICPC nicely shows the steps in the resolution of a problem.
A helium balloon is presented to the team for each problem solved.
To give here a list of all the programming competitions and training sites is quite
impossible, and such a list would quickly become obsolete. Nevertheless, we will
review some of the most important ones.
ICPC The oldest of these competitions was founded by the Association for
Computing Machinery in 1977 and supported by them up until 2017. This
contest, known as the ICPC, for International Collegiate Programming Contest,
is organised in the form of a tournament. The starting point is one of the regional
competitions, such as the South-West European Regional Contest (SWERC),
where the two best teams qualify for the worldwide final. The particularity of
this contest is that each three-person team has only a single computer at their
disposal. They have only 5 hours to solve a maximum number of problems
among the 10 proposed. The first ranking criterion is the number of submitted
solutions accepted (i.e. tested successfully against a set of unknown instances).
The next criterion is the sum over the submitted problems of the time between
the start of the contest and the moment of the accepted submission. For each
erroneous submission, a penalty of 20 minutes is added.
There are several competing theories on what the ideal composition of a
team is. In general, a good programmer and someone familiar with algorithms
is required, along with a specialist in different domains such as graph theory,
dynamic programming, etc. And, of course, the team members need to get along
together, even in stressful situations!
For the contest, each team can bring 25 pages of reference code printed in an
8-point font. They can also access the online documentation of the Java API and
the C++ standard library.
Google Code Jam In contrast with the ICPC contest, which is limited to students
up to a Master’s level, the Google Code Jam is open to everyone. This more
recent annual competition is for individual contestants. Each problem comes in
1.1 Programming Competitions 3
general with a deck of small instances whose resolution wins a few points, and
a set of enormous instances for which it is truly important to find a solution
with the appropriate algorithmic complexity. The contestants are informed of
the acceptance of their solution for the large instances only at the end of the
contest. However, its true strong point is the possibility to access the solutions
submitted by all of the participants, which is extremely instructive.
The competition Facebook Hacker Cup is of a similar nature.
Prologin The French association Prologin organises each year a competition
targeted at students up to twenty years old. Their capability to solve algorithmic
problems is put to test in three stages: an online selection, then regional
competitions and concluding with a national final. The final is atypically an
endurance test of 36 hours, during which the participants are confronted with a
problem in Artificial Intelligence. Each candidate must program a “champion”
to play a game whose rules are defined by the organisers. At the end of the
day, the champions are thrown in the ring against each other in a tournament to
determine the final winner.
The website https://prologin.org includes complete archives of past
problems, with the ability to submit algorithms online to test the solutions.
France-IOI Each year, the organisation France-IOI prepares junior and senior
high school students for the International Olympiad in Informatics. Since
2011, they have organised the ‘Castor Informatique’ competition, addressed at
students from Grade 4 to Grade 12 (675,000 participants in 2018). Their website
http://france-ioi.org hosts a large number of algorithmic problems (more
than 1,000).
Chinese online judges Several training sites now exist in China. They tend to
have a purer and more refined interface than the traditional judges. Nevertheless,
sporadic failures have been observed.
poj http://poj.org
tju http://acm.tju.edu.cn (Shut down since 2017)
zju http://acm.zju.edu.cn
Modern online judges Sphere Online Judge http://spoj.com and Kattis
http://open.kattis.com have the advantage of accepting the submission
of solutions in a variety of languages, including Python.
spoj http://spoj.com
kattis http://open.kattis.com
zju http://acm.zju.edu.cn
Other sites
codechef http://codechef.com
codility http://codility.com
gcj http://code.google.com/codejam
prologin http://prologin.org
slpc http://cs.stanford.edu/group/acm
Throughout this text, problems are proposed at the end of each section in rela-
tion to the topic presented. They are accompanied with their identifiers to a judge
site; for example [spoj:CMPLS] refers to the problem ‘Complete the Sequence!’ at
the URL www.spoj.com/problems/CMPLS/. The site http://tryalgo.org contains
links to all of these problems. The reader thus has the possibility to put into practice
the algorithms described in this book, testing an implementation against an online
judge.
The languages used for programming competitions are principally C++ and Java.
The SPOJ judge also accepts Python, while the Google Code Jam contest accepts
many of the most common languages. To compensate for the differences in execution
speed due to the choice of language, the online judges generally adapt the time limit
to the language used. However, this adaptation is not always done carefully, and it is
sometimes difficult to have a solution in Python accepted, even if it is correctly written.
We hope that this situation will be improved in the years to come. Also, certain judges
work with an ancient version of Java, in which such useful classes as Scanner are
not available.
Accepted Your program provides the correct output in the allotted time. Congrat-
ulations!
Presentation Error Your program is almost accepted, but the output contains
extraneous or missing blanks or end-of-lines. This message occurs rarely.
Compilation Error The compilation of your program generates errors. Often,
clicking on this message will provide the nature of the error. Be sure to compare
the version of the compiler used by the judge with your own.
Wrong Answer Re-read the problem statement, a detail must have been over-
looked. Are you sure to have tested all the limit cases? Might you have left
debugging statements in your code?
Time Limit Exceeded You have probably not implemented the most efficient
algorithm for this problem, or perhaps have an infinite loop somewhere. Test
your loop invariants to ensure loop termination. Generate a large input instance
and test locally the performance of your code.
Runtime Error In general, this could be a division by zero, an access beyond the
limits of an array, or a pop() on an empty stack. However, other situations can
also generate this message, such as the use of assert in Java, which is often not
accepted.
The taciturn behaviour of the judges nevertheless allows certain information to be
gleaned from the instances. Here is a trick that was used during an ICPC / SWERC
contest. In a problem concerning graphs, the statement indicated that the input con-
sisted of connected graphs. One of the teams doubted this, and wrote a connectivity
test. In the positive case, the program entered into an infinite loop, while in the negative
case, it caused a division by zero. The error code generated by the judge (Time Limit
Exceeded ou Runtime Error) allowed the team to detect that certain graphs in the input
were not connected.
The programming language Python was chosen for this book, for its readability and
ease of use. In September 2017, Python was identified by the website https://
stackoverflow.com as the programming language with the greatest growth in high-
income countries, in terms of the number of questions seen on the website, notably
thanks to the popularity of machine learning.1 Python is also the language retained
for such important projects as the formal calculation system SageMath, whose critical
portions are nonetheless implemented in more efficient languages such as C++ or C.
Here are a few details on this language. This chapter is a short introduction to
Python and does not claim to be exhaustive or very formal. For the neophyte reader
we recommend the site python.org, which contains a high-quality introduction as
well as exceptional documentation. A reader already familiar with Python can profit
1 https://stackoverflow.blog/2017/09/06/incredible-growth-python/
6 Introduction
enormously by studying the programs of David Eppstein, which are very elegant and
highly readable. Search for the keywords Eppstein PADS.
Python is an interpreted language. Variable types do not have to be declared, they
are simply inferred at the time of assignment. There are neither keywords begin/end
nor brackets to group instructions, for example in the blocks of a function or a loop.
The organisation in blocks is simply based on indentation! A typical error, difficult
to identify, is an erroneous indentation due to spaces used in some lines and tabs
in others.
Data Structures
The principal complex data structures are dictionaries, sets, lists and n-tuples. These
structures are called containers, as they contain several objects in a structured manner.
Once again, there are functions dict, set, list and tuple that allow the conversion of
an object into one of these structures. For example, for a string s, the function list(s)
returns a list L composed of the characters of the string. We could then, for example,
replace certain elements of the list L and then recreate a string by concatenating the ele-
ments of L with the expression ‘’.join(L). Here, the empty string could be replaced
by a separator: for example, ‘-’.join([’A’,’B’,’C’]) returns the string “A-B-C”.
Lists The indices of a list start with 0. The last element can also be accessed with
the index −1, the second last with −2 and so on. Here are some examples of
operations to extract elements or sublists of a list. This mechanism is known as
slicing, and is also available for strings.
1.2 Python in a Few Words 7
A loop in Python is written either with the keyword for or with while. The nota-
tion for the loop for is for x in S:, where the variable x successively takes on the
values of the container S, or of the keys of S in the case of a dictionary. In contrast,
while L: will loop as long as the list L is non-empty. Here, an implicit conversion of
a list to a Boolean is made, with the convention that only the empty list converts to
False.
At times, it is necessary to handle at the same time the values of a list along with
their positions (indices) within the list. This can be implemented as follows:
For a dictionary, the following loop iterates simultaneously over the keys and
values:
>>> n = 5
>>> squared_numbers = [x ** 2 for x in range(n + 1)]
>>> squared_numbers
[0, 1, 4, 9, 16, 25]
nb_occurrences = {}
for letter in my_string:
nb_occurrences[letter] = 0
range(k, n) from k to n − 1
range(k, n, 2) from k to n − 1 two by two
range(n - 1, -1, -1) from n − 1 to 0 (−1 excluded) in decreasing order.
def all_pairs(L):
n = len(L)
for i in range(n):
for j in range(i + 1, n):
yield (L[i], L[j])
math This module contains mathematical functions and constants such as log,
sqrt, pi, etc. Python operates on integers with arbitrary precision, thus there
is no limit on their size. As a consequence, there is no integer equivalent to
represent −∞ or +∞. For floating point numbers on the other hand, float(’-
inf’) and float(’inf’) can be used. Beginning with Python 3.5, math.inf (or
from math import inf) is equivalent to float(’inf’).
fractions This module exports the class Fraction, which allows computations
with fractions without the loss of precision of floating point calculations. For
example, if f is an instance of the class Fraction, then str(f) returns a string
similar to the form “3/2”, expressing f as an irreducible fraction.
bisect Provides binary (dichotomous) search functions on a sorted list.
heapq Provides functions to manipulate a list as a heap, thus allowing an element
to be added or the smallest element removed in time logarithmic in the size of
the heap; see Section 1.5.4 on page 22.
string This module provides, for example, the function ascii_lowercase, which
returns its argument converted to lowercase characters. Note that the strings
themselves already provide numerous useful methods, such as strip, which
removes whitespace from the beginning and end of a string and returns the
result, lower, which converts all the characters to lowercase, and especially
split, which detects the substrings separated by spaces (or by another separa-
tor passed as an argument). For example, “12/OCT/2018”.split(“/”) returns
[“12”, “OCT”, “2018”].
Packages
One of the strengths of Python is the existence of a large variety of code packages.
Some are part of the standard installation of the language, while others must be
imported with the shell command pip. They are indexed and documented at the web
site pypi.org. Here is a non-exhaustive list of very useful packages.
tryalgo All the code of the present book is available in a package called tryalgo
and can be imported in the following manner: pip install tryalgo.
collections To simplify life, the class from collections import Counter can
be used. For an object c of this class, the expression c[x] will return 0 if x is not
1.2 Python in a Few Words 11
a key in c. Only modification of the value associated with x will create an entry
in c, such as, for example, when executing the instruction c[x] += 1. This is
thus slightly more practical than a dictionary, as is illustrated below.
>>> c = {} # dictionary
>>> c[’a’] += 1 # the key does not exist
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
KeyError: ’a’
>>> c[’a’] = 1
>>> c[’a’] += 1 # now it does
>>> c
{’a’: 2}
>>> from collections import Counter
>>> c = Counter()
>>> c[’a’] += 1 # the key does not exist, so it is created
Counter({’a’: 1})
>>> c = Counter(’cowboy bebop’)
Counter({’o’: 3, ’b’: 3, ’c’: 1, ’w’: 1, ’y’: 1, ’ ’: 1, ’e’: 1, ’p’: 1})
The collections package also provides the class deque, for double-ended
queue, see Section 1.5.3 on page 21. With this structure, elements can be added
or removed either from the left (head) or from the right (tail) of the queue. This
helps implement Dijkstra’s algorithm in the case where the edges have only
weights 0 or 1, see Section 8.2 on page 126.
This package also provides the class defaultdict, which is a dictionary that
assigns default values to keys that are yet in the dictionary, hence a generalisa-
tion of the class Counter.
See also Section 1.3 on page 13 for an example of reading a graph given as
input.
numpy This package provides general tools for numerical calculations involving
manipulations of large matrices. For example, numpy.linalg.solve solves
a linear system, while numpy.fft.fft calculates a (fast) discrete Fourier
transform.
While writing the code of this book, we have followed the norm PEP8, which
provides precise recommendations on the usage of blanks, the choice of names for
variables, etc. We advise the readers to also follow these indications, using, for exam-
ple, the tool pycodestyle to validate the structure of their code.
12 Introduction
A = [1, 2, 3]
B = A # Beware! Both variables refer to the same object
A = [1, 2, 3]
B = A[:] # B becomes a distinct copy of A
The notation [:] can be used to make a copy of a list. It is also possible to make a
copy of all but the first element, A[1 :], or all but the last element, A[: −1], or even
in reverse order A[:: −1]. For example, the following code creates a matrix M, all of
whose rows are the same, and the modification of M[0][0] modifies the whole of the
first column of M.
A square matrix can be correctly initialised using one of the following expressions:
The module numpy permits easy manipulations of matrices; however, we have cho-
sen not to profit from it in this text, in order to have generic code that is easy to translate
to Java or C++.
Ranges
Another typical error concerns the use of the function range. For example, the follow-
ing code processes the elements of a list A between the indices 0 and 9 inclusive, in
order.
To process the elements in descending order, it is not sufficient to just swap the
arguments. In fact, range(10, 0, -1)—the third argument indicates the step—is the
list of elements with indices 10 (included) to 0 (excluded). Thus the loop must be
written as:
1.3 Input-Output
If you wish to save the output of your program to a file called test.out, type:
A little hint: if you want to display the output at the same time as it is being written
to a file test.out, use the following (the command tee is not present by default in
Windows):
The inputs can be read line by line via the command input(), which returns the
next input line in the form of a string, excluding the end-of-line characters.3 The
module sys contains a similar function stdin.readline(), which does not suppress
the end-of-line characters, but according to our experience has the advantage of being
four times as fast!
If the input line is meant to contain an integer, we can convert the string with the
function int (if it is a floating point number, then we must use float instead). In
the case of a line containing several integers separated by spaces, the string can first be
cut into different substrings using split(); these can then be converted into integers
with the method map. For example, in the case of two integers height and width to be
read on the same line, separated by a space, the following command suffices:
import sys
If your program exhibits performance problems while reading the inputs, our expe-
rience shows that a factor of two can be gained by reading the whole of the inputs with
a single system call. The following code fragment assumes that the inputs are made
up of only integers, eventually on multiple lines. The parameter 0 in the function
os.read means that the read is from standard input, and the constant M must be an
upper bound on the size of the file. For example, if the file contains 107 integers,
each between 0 and 109 , then as each integer is written with at most 10 characters
and there are at most 2 characters separating the integers (\n and \r), we can choose
M = 12 · 107 .
3 According to the operating system, the end-of-line is indicated by the characters \r, \n, or both, but this
is not important when reading with input(). Note that in Python 2 the behaviour of input() is
different, so it is necessary to use the equivalent function raw_input().
1.3 Input-Output 15
import os
3
paris tokyo 9471
paris new-york 5545
new-york singapore 15344
nb_edges = int(input())
g = defaultdict(dict)
for _ in range(nb_edges):
u, v, weight = input().split()
g[u][v] = int(weight)
# g[v][u] = int(weight) # For an undirected graph
def readint():
return int(stdin.readline())
def readarray(typ):
return list(map(typ, stdin.readline().split()))
def readmatrix(n):
M = []
for _ in range(n):
row = readarray(int)
assert len(row) == n
M.append(row)
return M
if __name__ == "__main__":
n = readint()
A = readmatrix(n)
B = readmatrix(n)
C = readmatrix(n)
print(freivalds(A, B, C))
Note the test on the variable __name__. This test is evaluated as True if the file
containing this code is called directly, and as False if the file is included with the
import keyword.
Problem
Enormous Input Test [spoj:INTEST]
To display numbers with fixed precision and length, there are at least two pos-
sibilities in Python. First of all, there is the operator % that works like the func-
tion printf in the language C. The syntax is s % a, where s is a format string, a
character string including typed display indicators beginning with %, and where a
consists of one or more arguments that will replace the display indicators in the format
string.
>>> i_test = 1
>>> answer = 1.2142
>>> print("Case #%d: %.2f gigawatts!!!" % (i_test, answer))
Case #1: 1.21 gigawatts!!!
The letter d after the % indicates that the first argument should be interpreted as an
integer and inserted in place of the %d in the format string. Similarly, the letter f is
used for floats and s for strings. A percentage can be displayed by indicating %% in
the format string. Between the character % and the letter indicating the type, further
numerical indications can be given. For example, %.2f indicates that up to two digits
should be displayed after the decimal point.
Another possibility is to use the method format of a string, which follows the
syntax of the language C#. This method provides more formatting possibilities and
is in general easier to manipulate.
Finally, beginning with Python 3.6, f-strings, or formatted string literals, exist.
In this case, the floating point precision itself can be a variable, and the formatting
is embedded with each argument.
>>> precision = 2
>>> print(f"Case #{testCase}: {answer:.{precision}f} gigawatts!!!")
Case #1: 1.21 gigawatts!!!
1.4 Complexity
Landau Symbols
The complexity of an algorithm is, for example, said to be O(n2 ) if the execution time
can be bounded above by a quadratic function in n, where n represents the size or some
parameter of the input. More precisely, for two functions f ,g we denote f ∈ O(g)
if positive constants n0,c exist, such that for every n ≥ n0 , f (n) ≤ c · g(n). By an
abuse of notation, we also write f = O(g). This notation allows us to ignore the
multiplicative and additive constants in a function f and brings out the magnitude and
form of the dependence on a parameter.
Similarly, if for constants n0,c > 0 we have f (n) ≥ c · g(n) for every n ≥ n0 , then
we write f ∈ (g). If f ∈ O(g) and f ∈ (g), then we write f ∈ (g), indicating
that f and g have the same order of magnitude of complexity. Finally, if f ∈ O(g)
but not g ∈ O(f ), then we write f ∈ o(g)
Complexity Classes
If the complexity of an algorithm is O(nc ) for some constant c, it is said to be polyno-
mial in n. A problem for which a polynomial algorithm exists is said to be polynomial,
and the class of such problems bears the name P. Unhappily, not all problems are
polynomial, and numerous problems exist for which no polynomial algorithm has
been found to this day.
One such problem is k-SAT: Given n Boolean variables and m clauses each
containing k literals (a variable or its negation), is it possible to assign to each variable
a Boolean value in such a manner that each clause contains at least one literal with the
value True (SAT is the version of this problem without a restriction on the number of
variables in a clause)? The particularity of each of these problems is that a potential
solution (assignment to each of the variables) satisfying all the constraints can be
verified in polynomial time by evaluating all the clauses: they are in the class NP
(for Non-deterministic Polynomial). We can easily solve 1-SAT in polynomial time,
hence 1-SAT is in P. 2-SAT is also in P; this is the subject of Section 6.9 on page 110.
However, from 3-SAT onwards, the answer is not known. We only know that solving
3-SAT is at least as difficult as solving SAT.
It turns out that P ⊆ NP—intuitively, if we can construct a solution in polynomial
time, then we can also verify a solution in polynomial time. It is believed that
P = NP, but this conjecture remains unproven to this day. In the meantime,
researchers have linked NP problems among themselves using reductions, which
transform in polynomial time an algorithm for a problem A into an algorithm for a
problem B. Hence, if A is in P, then B is also in P: A is ‘at least as difficult’ as B.
The problems that are at least as difficult as SAT constitute the class of problems
NP-hard, and among these we distinguish the NP-complete problems, which are
defined as those being at the same time NP-hard and in NP. Solve any one of these
in polynomial time and you will have solved them all, and will be gratified by eternal
recognition, accompanied by a million dollars.4 At present, to solve these problems
in an acceptable time, it is necessary to restrict them to instances with properties
4 www.claymath.org/millennium-problems/p-vs-np-problem
1.4 Complexity 19
that can aid the resolution (planarity of a graph, for example), permit the program
to answer correctly with only a constant probability or produce a solution only close
to an optimal solution.
Happily, most of the problems encountered in programming contests are
polynomial.
If these questions of complexity classes interest you, we recommend Christos H.
Papadimitriou’s book on computational complexity (2003).
For programming contests in particular, the programs must give a response within
a few seconds, which gives current processors the time to execute on the order of tens
or hundreds of millions of operations. The following table gives a rough idea of the
acceptable complexities for a response time of a second. These numbers are to be
taken with caution and, of course, depend on the language used,5 the machine that
will execute the code and the type of operation (integer or floating point calculations,
calls to mathematical functions, etc.).
1000000 O(n)
100000 O(n log n)
1000 O(n2 )
The reader is invited to conduct experiments with simple programs to test the time
necessary to execute n integer multiplications for different values of n. We insist here
on the fact that the constants lurking behind the Landau symbols can be very large,
and at times, in practice, an algorithm with greater asymptotic complexity may be
preferred. An example is the multiplication of two n × n matrices. The naive method
costs O(n3 ) operations, whereas a recursive procedure discovered by Strassen (1969)
costs only O(n2.81 ). However, the hidden multiplicative constant is so large that for the
size of matrices that you might need to manipulate, the naive method without doubt
will be more efficient.
In Python, adding an element to a list takes constant time, as does access to an
element with a given index. The creation of a sublist with L[i : j ] requires time
O(max{1,j − i}). Dictionaries are represented by hash tables, and the insertion of
a key-value pair can be done in constant time, see Section 1.5.2 on page 21. However,
this time constant is non-negligible, hence if the keys of the dictionary are the integers
from 0 to n − 1, it is preferable to use a list.
Amortised Complexity
For certain data structures, the notion of amortised time complexity is used. For exam-
ple, a list in Python is represented internally by an array, equipped with a variable
length storing the size. When a new element is added to the array with the method
append, it is stored in the array at the index indicated by the variable length, and then
5 Roughly consider Java twice as slow and Python four times as slow as C++.
20 Introduction
the variable length is incremented. If the capacity of the array is no longer sufficient,
a new array twice as large is allocated, and the contents of the original array are copied
over. Hence, for n successive calls to append on a list initially empty, the time of each
call is usually constant, and it is occasionally linear in the current size of the list.
However, the sum of the time of these calls is O(n), which, once spread out over each
of the calls, gives an amortised time O(1).
We will now tackle a theme that is at the heart of the notion of efficient programming:
the data structures underlying our programs to solve problems.
An abstract type is a description of the possible values that a set of objects can take
on, the operations that can be performed on these objects and the behaviour of these
operations. An abstract type can thus be seen as a specification.
A data structure is a concrete way to organise the data in order to treat them
efficiently, respecting the clauses in the specification. Thus, we can implement an
abstract type by one or more different data structures and determine the complexity in
time and memory of each operation. Thereby, based on the frequency of the operations
that need to be executed, we will prefer one or another of the implementations of an
abstract type to solve a given problem.
To program well, a mastery of the data structures offered by a language and its
associated standard library is essential. In this section, we review the most useful data
structures for programming competitions.
1.5.1 Stacks
A stack (see Figure 1.2) is an object that keeps track of a set of elements and pro-
vides the following operations: test if the stack is empty, add an element to the top
of the stack (push), access the element at the top of the stack and remove an ele-
ment (pop). Python lists can serve as stacks. An element is pushed with the method
append(element) and popped with the method pop().
stack
queue
deque
Figure 1.2 The three principal sequential access data structures provided by Python.
the case for all objects implementing the method __len__; for these, if x: behaves
exactly like if len(x) > 0:. All of these operations execute in constant time.
1.5.2 Dictionaries
A dictionary allows the association of keys to values, in the same manner that
an array associates indices to values. The internal workings are based on a hash
table, which uses a hash function to associate the elements with indices in an array,
and implements a mechanism of collision management in the case where different
elements are sent to the same index. In the best case, the read and write operations
on a dictionary execute in constant time, but in the worst case they execute in linear
time if it is necessary to run through a linear number of keys to handle the collisions.
In practice, this degenerate case arrives only rarely, and in this book, we generally
assume that accesses to dictionary entries take constant time. If the keys are of the
form 0,1, . . . ,n − 1, it is, of course, always preferable for performance reasons to use
a simple array.
Problems
Encryption [spoj:CENCRY]
Phone List [spoj:PHONELST]
A concrete simulation [spoj:ACS]
1.5.3 Queues
A queue is similar to a stack, with the difference that elements are added to the end
of the queue (enqueued) and are removed from the front of the queue (dequeued).
A queue is also known as a FIFO queue (first in, first out, like a waiting line), whereas
a stack is called LIFO (last in, first out, like a pile of plates).
In the Python standard library, there are two classes implementing a queue. The
first, Queue, is a synchronised implementation, meaning that several processes can
access it simultaneously. As the programs of this book do not exploit parallelism,
the use of this class is not recommended, as it entails a slowdown because of the
use of semaphores for the synchronisation. The second class is deque (for double-
ended queue). In addition to the methods append(element) and popleft(), which,
respectively, add to the end of the queue and remove from the head of the queue,
deque offers the methods appendleft(element) and pop(), which add to the head of
the queue and remove from the end of the queue. We can thus speak of a queue with
two ends. This more sophisticated structure will be found useful in Section 8.2 on
page 126, where it is used to find the shortest path in a graph the edges of which have
weights 0 or 1.
We recommend the use of deque—and in general, the use of the data structures
provided by the standard library of the language—but as an example we illustrate here
how to implement a queue using two stacks. One stack corresponds to the head of the
queue for extraction and the other corresponds to the tail for insertion. Once the head
22 Introduction
stack is empty, it is swapped with the tail stack (reversed). The operator __len__ uses
len(q) to recover the number of elements in the queue q, and then if q can be used
to test if q is non-empty, happily in constant time.
class OurQueue:
def __init__(self):
self.in_stack = [] # tail
self.out_stack = [] # head
def __len__(self):
return len(self.in_stack) + len(self.out_stack)
2 3
4 5 6 7
8 9 10 11 12
0 1 2 3 4 5 6 7 8 9 10 11 12
Implementation details
The structure contains an array heap, storing the heap itself, as well as a dictionary
rank, allowing the determination of the index of an element stored in the heap. The
principal operations are push and pop. A new element is inserted with push: it is added
as the last leaf in the heap, and then the heap is reorganised to respect the heap order.
The minimal element is extracted with pop: the root is replaced by the last leaf in the
heap, and then the heap is reorganised to respect the heap order, see Figure 1.4.
The operator __len__ returns the number of elements in the heap. The existence of
this operator permits the inner workings of Python to implicitly convert a heap to a
Boolean and to perform such conditional tests as, for example, while h, which loops
while the heap h is non-empty.
The average complexity of the operations on our heap is O(log n); however, the
worst-case complexity is O(n), due to the use of the dictionary (rank).
5 7
13 9 8 30
20 17 11 12 15
9 7
13 11 8 30
20 17 15 12
Figure 1.4 The operation pop removes and returns the value 2 of the heap and replaces it by the
last leaf 15. Then the operation down performs a series of exchanges to place 15 in a position
respecting the heap order.
1.5 Abstract Types and Essential Data Structures 25
class OurHeap:
def __init__(self, items):
self.heap = [None] # index 0 will be ignored
self.rank = {}
for x in items:
self.push(x)
def __len__(self):
return len(self.heap) - 1
def pop(self):
root = self.heap[1]
del self.rank[root]
x = self.heap.pop() # remove last leaf
if self: # if heap is not empty
self.heap[1] = x # move the last leaf
self.rank[x] = 1 # to the root
self.down(1) # maintain heap order
return root
The reorganisation is done with the operations up(i) and down(i), which are called
whenever an element with index i is too small with respect to its parent (for up)
or too large for its children (for down). Hence, up effects a series of exchanges of
a node with its parents, climbing up the tree until the heap order is respected. The
action of down is similar, for an exchange between a node and its child with the
smallest value.
Finally, the method update permits the value of a heap element to be changed. It
then calls up or down to preserve the heap order. It is this method that requires the
introduction of the dictionary rank.
26 Introduction
1.5.5 Union-Find
Definition
Union-find is a data structure used to store a partition of a universe V and that allows
the following operations, also known as requests in the context of dynamic data
structures:
1.5 Abstract Types and Essential Data Structures 27
• find(v) returns a canonical element of the set containing v. To test if u and v are
in the same set, it suffices to compare find(u) with find(v).
• union(u,v) combines the set containing u with that containing v.
Application
Our principal application of this structure is to determine the connected components
of a graph, see Section 6.5 on page 94. Every addition of an edge will correspond to
a call to union, while find is used to test if two vertices are in the same component.
Union-find is also used in Kruskal’s algorithm to determine a minimal spanning tree
of a weighted graph, see Section 10.4 on page 179.
1. When traversing an element on the way to the root, we take advantage of the
opportunity to compress the path, i.e. we link directly to the root that the elements
encountered.
2. During a union, we hang the tree of smaller rank beneath the root of the tree of
larger rank. The rank of a tree corresponds to the depth it would have if none of
the paths were compressed.
8 12 10
7 2 4 11 9 8 12 2 4 11
3 6 1 7 3 6 1 9 10
5 5
Figure 1.5 On the left: The structure union-find for a graph with two components {7,8,12}
and {2,3,4,5,6,9,10,11}. On the right: after a call to find(10), all the vertices on the path
to~the root now point directly to the root 5. This accelerates future calls to find for these
vertices.
28 Introduction
class UnionFind:
def __init__(self, n):
self.up_bound = list(range(n))
self.rank = [0] * n
Problem
Havannah [gcj:2013round3B]
1.6 Techniques
1.6.1 Comparison
In Python, a comparison between n-tuples is based on lexicographical order. This
allows us, for example, to find, at the same time, the largest value in an array along
with the corresponding index, taking the largest index in the case of equality.
use of a dictionary, where n is the number of words given and k is the maximal
length of a word.
def majority(L):
count = {}
for word in L:
if word in count:
count[word] += 1
else:
count[word] = 1
# Using min() like this gives the first word with
# maximal count "for free"
val_1st_max, arg_1st_max = min((-count[word], word) for word in count)
return arg_1st_max
1.6.2 Sorting
An array of n elements in Python can be sorted in time O(n log n). We distinguish two
kinds of sorts:
• a sort with the method sort() modifies the list in question (it is said to sort ‘in
place’) ;
• a sort with the function sorted() returns a sorted copy of the list in question.
Given a list L of n distinct integers, we would like to determine two integers
in L with minimal difference. This problem can be solved by sorting L, and then
running through L to select the pair with the smallest difference. Note the use of
minimum over a collection of pairs of integers, compared in lexicographical order.
Hence, valmin will contain the smallest distance between two successive elements
of L and argmin will contain the corresponding index.
def closest_values(L):
assert len(L) >= 2
L.sort()
valmin, argmin = min((L[i] - L[i - 1], i) for i in range(1, len(L)))
return L[argmin - 1], L[argmin]
Sorting n elements requires (n log n) comparisons between the elements in the
worst case. To be convinced, consider an input of n distinct integers. An algorithm
must select one among n! possible orders. Each comparison returns one of two values
(greater or smaller) and thus divides the search space in two. Finally, it requires at
most log2 (n! )
comparisons to identify a particular order, giving the lower bound
(log(n! )) = (n log n).
Variants
In certain cases, an array of n integers can be sorted in time O(n), for example when
they are all found between 0 and cn for some constant c. In this case, we can simply
30 Introduction
scan the input and store the number of occurrences of each element in an array count
of size cn. We then scan count by increasing index, and write the values from 0 to cn
to the output array, repeating them as often as necessary. This technique is known as
counting sort; a similar variant is pigeonhole sort.
Problems
Spelling Lists [spoj:MIB]
Yodaness Level [spoj:YODANESS]
Example—Intersection of Intervals
Given n intervals [i ,ri ) for i = 0, . . . ,n − 1, we wish to find a value x included
in a maximum number of intervals. Here is a solution in time O(n log n). We sort
the endpoints, and then sweep them from left to right with an imaginary pointer x.
A counter c keeps track of the number of intervals whose beginning has been seen,
but not yet the end, hence it counts the number of intervals containing x.
Note that the order of processing of the elements of B guarantees that the right
endpoints of the intervals are handled before the left endpoints of the intervals, which
is necessary when dealing with intervals that are half-open to the right.
def max_interval_intersec(S):
B = ([(left, +1) for left, right in S] +
[(right, -1) for left, right in S])
B.sort()
c = 0
best = (c, None)
for x, d in B:
c += d
if best[0] < c:
best = (c, x)
return best
Problem
Back to the future [prologin:demi2012]
local criterion. A formal notion exists that is related to combinatorial structures known
as matroids, and it can be used to prove the optimality or the non-optimality of greedy
algorithms. We do not tackle this here, but refer to Kozen (1992) for a very good
introduction.
For some problems, an optimal solution can be produced step by step, making a
decision that is in a certain sense locally optimal. Such a proof is, in general, based on
an exchange argument, showing that any optimal solution can be transformed into the
solution produced by the algorithm without modifying the cost. The reader is invited
to systematically at least sketch such a proof, as intuition can often be misleading and
the problems that can be solved by a greedy algorithm are in fact quite rare.
An example is making change with coins, see Section 11.2 on page 184. Suppose
you work in a store, and dispose of an infinity of coins with a finite number of distinct
values. You must make change for a client for a specific sum and you wish to do it
with as few coins as possible. In most countries, the values of the coins are of the
form x · 10i with x ∈ {1,2,5} and i ∈ N, and in this case we can in a greedy manner
make change by repeatedly selecting the most valuable coin that is not larger than
the remaining sum to return. However, in general, this approach does not work, for
example with coins with values 1, 4 and 5, and with a total of 8 to return. The greedy
approach generates a solution with 4 coins, i.e. 5 + 1 + 1 + 1, whereas the optimal is
composed of only 2 coins: 4 + 4.
Application
Suppose there are n tasks to be assigned to n workers in a bijective manner, i.e. where
each task must be assigned to a different worker. Each task has a duration in hours and
each worker has a price per hour. The goal is to find the assignment that minimises the
total amount to be paid.
x0 ≤ xk
x0 (yj − yi ) ≤ xk (yj − yi )
32 Introduction
x0 yj − x0 yi ≤ xk yj − xk yi
x0 yj + xk yi ≤ x0 yi + xk yj .
By repeating the argument on the vector x with x0 removed and on y with yj removed,
we find that the product is minimal when i → yπ(i) is decreasing.
Problems
Minimum Scalar Product [gcj:2008round1A]
Minimal Coverage [timus:1303]
F (0) = 0
F (1) = 1
F (i) = F (i − 1) + F (i − 2).
This number gives, for example, the number of possibilities to climb an n-step
stairway by taking one or two steps at a time. A naive implementation of F by a
recursive function is extremely inefficient, since for the same parameter i, F (i) is
computed several times, see Figure 1.6.
def fibo_naive(n):
if n <= 1:
return n
return fibo_naive(n - 1) + fibo_naive(n - 2)
1.6 Techniques 33
F4
F2 F3
F0 F1 F1 F2
F0 F1 F2 F3 F4
F0 F1
Figure 1.6 On the left, the exhaustive exploration tree of the recursive computation of the
Fibonacci number F4 . On the right, the acyclic oriented graph of the value dependencies of
dynamic programming, with many fewer nodes.
Try this function on a not-very-large value, n = 32, and you will see that the code
can already take more than a second, even on a powerful machine; the computation
time increases exponentially with n.
A solution by dynamic programming simply consists of storing in an array of size
n+1 the values of F (0) to F (n), and filling it in increasing order of the indices. Hence,
at the moment we compute F (i), the values F (i − 1) and F (i − 2) will already have
been computed; they are the last two values entered into the array.
def fibo_dp(n):
mem = [0, 1]
for i in range(2, n + 1):
mem.append(mem[-2] + mem[-1])
return mem[-1]
However, the above code uses n + 1 memory locations, whereas 2 are sufficient,
those corresponding to the last two Fibonacci numbers computed.
def fibo_dp_mem(n):
mem = [0, 1]
for i in range(2, n + 1):
mem[i % 2] = mem[0] + mem[1]
return mem[n % 2]
dictionary, to avoid the combinatorial explosion of the calls. In other words, adding
@lru_cache(maxsize=None) transforms any old naive recursive implementation into
a clever example of dynamic programming. Isn’t life beautiful?
@lru_cache(maxsize=None)
def fibo_naive(n):
if n <= 1:
return n
return fibo_naive(n - 1) + fibo_naive(n - 2)
Problem
Fibonacci Sum [spoj:FIBOSUM]
{} 0 empty set
{i} 1 << i this notation represents 2i
{0,1, . . . ,n − 1} (1 << n) − 1 2n − 1 = 20 + 21 + · · · + 2n−1
A∪B A|B the operator | is the binary or
A∩B A&B the operator & is the binary and
(A\B) ∪ (B\A) AˆB the operator ˆ is the exclusive or
A⊆B A & B == A test of inclusion
i∈A (1 << i)& A test of membership
{min A} −A & A this expression equals 0 if A is empty
The justification of this last expression is shown in Figure 1.7. It is useful, for
example, to count in a loop the cardinality of a set. There is no equivalent expression
to obtain the maximum of a set.
To illustrate this encoding technique, we apply it to a classic problem.
7 The limit comes from the fact that integers in Python are, in general, coded in a machine word, which
today is usually 63 bits plus a sign bit, for a total of 64 bits.
1.6 Techniques 35
We wish to obtain the minimum of the set {3,5,6} encoded by an integer. This integer is
23 + 25 + 26 = 104.
64+32+8=104 01101000
complement of 104 10010111
-104 10011000
-104&104 = 8 00001000
The result is 23 encoding the singleton {3}.
def three_partition(x):
f = [0] * (1 << len(x))
for i, _ in enumerate(x):
for S in range(1 << i):
f[S | (1 << i)] = f[S] + x[i]
for A in range(1 << len(x)):
for B in range(1 << len(x)):
if A & B == 0 and f[A] == f[B] and 3 * f[A] == f[-1]:
return (A, B, ((1 << len(x)) - 1) ^ A ^ B)
return None
For another utilisation of this technique, see the resolution for numbers in the TV
contest ‘Des Chiffres et des Lettres’ (1965), presented as ‘Countdown’ (1982) in the
UK, Section 15.6 on page 243.
space is reduced either to [,m] or to [m + 1,h]. Note that because of the rounding
below in the computation of m, the second interval is never empty, nor the first for
that matter. The search terminates after log2 (n)
iterations, when the search interval
is reduced to a singleton.
Libraries
Binary search is provided in the standard module bisect, so that in many cases
you will not have to write it yourself. Consider the case of a sorted array tab of n
elements, where we want to find the insertion point of a new element x. The function
bisect_left(tab, x, 0, n) returns the first index i such that tab[i] ≥ x.
Continuous Domain
This technique can equally be applied when the domain of f is continuous and we
seek the smallest value x0 with f (x) = 1 for every x ≥ x0 . The complexity will then
depend on the precision required for x0 .
maximal value. The number of iterations necessary is again logarithmic, on the order
of log3/2 n.
Invert a Function
Let f be a continuous and strictly monotonic function. Thus there exists an inverse
function f −1 , again monotonic. Imagine that f −1 is much easier to compute than f ,
in this case we can use it to compute f (x) for a given value of x. Indeed, it suffices to
find the smallest value y such that f −1 (y) ≥ x.
Example—Filling Tasks
Suppose there are n tanks in the form of rectangular blocks at different heights, inter-
connected by a network of pipes. We pour into this system a volume V of liquid, and
wish to know the resulting height of the liquid in the tanks.
Problems
Fill the cisterns [spoj:CISTFILL]
Egg Drop [gcj:eggdrop]
1.7 Advice
Here is a bit of common sense advice to help you rapidly solve algorithmic problems
and produce a working program. Be organised and systematic. For this, it is important
to not be carried away with the desire to start hacking before all the details are clear.
It is easy to leap into the implementation of a method that can never work for a reason
38 Introduction
that would not have escaped you if you had taken a bit more time to ponder before
starting to pound on the keyboard.
When possible, in your program, separate the input of the instance from the com-
putation of the solution. Be nice to yourself and systematically add to the comments
the name of the problem, if possible with its URL, and specify the complexity of
your algorithm. You will appreciate this effort when you revisit your code at a later
date. Above all, stay coherent in your code and reuse the terms given in the problem
statement to highlight the correspondence. There is little worse than having to debug
a program whose variable names are not meaningful.
iteration is a classic error. Take, for example, a program solving a graph problem
whose input consists of the number of instances, followed by each instance,
beginning with two integers: the number of vertices n and the number of edges
m. These are followed by two arrays of integers A,B of size m, encoding the
endpoints of the edges for each instance. Imagine that the program encodes
the graph in the form of an adjacency list G and for every i = 0, . . . ,m − 1,
adds B[i] to G[A[i]] and A[i] to G[B[i]]. If the lists are not emptied before
the beginning of each input of an instance, the edges accumulate to form a
superposition of all the graphs.
Debug
Make all your mistakes now to have the proper reflexes later.
Generate test sets for the limit cases (Wrong Answer) and for large instances
(Time Limit Exceeded or Runtime Error).
Explain the algorithm to a team-mate and add appropriate comments to the pro-
gram. You must be capable of explaining each and every line.
Simplify your implementation. Regroup all similar bits of code. Never go wild
with copy-and-paste (note from the translators—one of the most observed
sources of buggy programs from novice programmers!).
Take a step back by passing to another problem, and come back later for a fresh
look.
Compare the environment of your local machine with that of the server where
your code will be tested.
See [icpcarchive:8154].
Iskander the Baker is decorating a huge cake by covering the rectangular sur-
face of the cake with frosting. For this purpose, he mixes frosting sugar with lemon
juice and beetle juice, in order to produce three kinds of frosting: yellow, pink and
white. These colours are identified by the numbers 0 for yellow, 1 for pink and 2
for white.
To obtain a nice pattern, he partitions the cake surface into vertical stripes of width
A1,A2, . . . ,An centimeters, and horizontal stripes of height B1,B2, . . . ,Bn centime-
ters, for some positive integer n. These stripes split the cake surface into n × n
rectangles. The intersection of vertical stripe i and horizontal stripe j has colour
number (i + j ) mod 3 for all 1 ≤ i,j ≤ n, see Figure 1.8. To prepare the frosting,
Iskander wants to know the total surface in square centimeters to be coloured for each
of the three colours, and asks for your help.
Input
The input consists of the following integers:
40 Introduction
Limits
The input satisfies 3 ≤ n ≤ 100 000 and 1 ≤ A1, . . . ,An,B1, . . . ,Bn ≤ 10 000.
Output
The output should consist of three integers separated by single spaces, representing
the total area for each colour 0, 1 and 2.
Solution
When you see the upper bound on n, you want to find a solution that runs in time O(n)
or possibly O(n log n). This rules out the naive solution which loops over all n2 grid
cells and accumulates their areas in variables corresponding to each colour. But with
a little thought about the problem, we can make the following observation. Permuting
columns or rows preserves the total area of each colour. Hence, we can reduce to the
n = 3 case, by simply summing the values of each colour class A3k , A3k+1 and A3k+2 .
Then the answer per colour class is just the sum of the areas of three rectangles.
Hence, a solution could be as short as the following one. The tricky part relies in
not mixing up the colours.
1.8 A Problem: ‘Frosting on the Cake’ 41
input() # n
A = cat(read_ints())
B = cat(read_ints())
print("{} {} {}".format(B[2] * A[0] + B[0] * A[2] + B[1] * A[1],
B[2] * A[1] + B[0] * A[0] + B[1] * A[2],
B[2] * A[2] + B[0] * A[1] + B[1] * A[0]))
2 Character Strings
The resolution of problems involving character strings (or simply strings) rapidly
became an important domain in the study of algorithms. These problems occur, of
course, in text processing, such as spell checking and the search for factors (sub-
strings) or for more general patterns. With the development of bioinformatics, new
problems arose in the alignment of DNA chains. This chapter presents a selection of
string-processing algorithms that we consider important.
A character string could be represented internally by a list of symbols, but, in
general, we use the Python native type str, which behaves essentially like a list. The
characters can be encoded with two bytes for a Unicode encoding, but are usually
encoded with only one byte, using the ASCII encoding: each integer between 0 and
127 corresponds to a distinct character, and the codes are organised in such a way that
the successive symbols ‘0’ to ‘9’, ‘a’ to ‘z’ and ‘A’ to ‘Z’ are consecutive. Thereby, if
a string s contains only capitals, we can recover the rank of the ith letter by computing
ord(s[i])-ord(’A’). Conversely, the j th capital letter of the alphabet numbered
from~0—is obtained with chr(j+ord(’A’)).
When we speak of a factor (or substring) of a string, we require the characters to
be consecutive, in contrast with the more general notion of a subsequence.
2.1 Anagrams
Definition
A word w is an anagram of a word v if a permutation of the letters transforming w
into v exists. Given a set of n words of length at most k, we would like to detect all
possible anagrams.
input: below the car is a rat drinking cider and bending its elbow while this thing
is an arc that can act like a cat which cried during the night caused by pain in its
bowel1
output: {bowel below elbow}, {arc car}, {night thing}, {cried cider}, {act cat}
1 Believe it or not: the authors did not consume cider in order to produce this sample input.
2.2 T9—Text on 9 Keys 43
Complexity
The proposed algorithm solves this problem in time O(nk log k) on average, and in
O(n2 k log k) in the worst case, due to the use of a dictionary.
Algorithm
The idea is to compute a signature for each word, so that two words are anagrams of
each other if and only if they have the same signature. This signature is simply a new
string made up of the same letters sorted in lexicographical order.
The data structure used is a dictionary associating with each signature the list of
words with this signature.
Problem
Anagram verifier [spoj:ANGRAM]
input: 2665687
output: bonjour
Application
Mobile telephones with keys offer an interesting input mode, sometimes called T9.
The 26 letters of the alphabet are distributed over the keys 2 to 9, as in Figure 2.1.
To input a word, it suffices to input the corresponding sequence of digits. However, as
several words could begin with the same digits, a dictionary must be used to propose
the most probable word. At any moment, the telephone displays the prefix of the most
probable word corresponding to the sequence of digits entered.
44 Character Strings
Definition
The first part of the problem instance is a dictionary, composed of pairs (m,w) where
m is a word over the alphabet of 26 lower-case letters, and w is a weight of the
importance of the word. The second part is a series of sequences of digits from 2
to 9. For each sequence s, the word in the dictionary of maximal weight is to be
displayed. A word m corresponds to s if s is a prefix of the sequence t obtained from
m by replacing each letter by the corresponding digit, according to the correspondence
table given in Figure 2.1. For example, ‘bonjour’ corresponds to 26, but also to 266
or 2665687.
ABC DEF
1 2 3
GHI JKL MNO
4 5 6
PQRS TUV WXYZ
7 8 9
* 0 #
Algorithm
The complexity is O(nk) for the initialisation of the dictionary, and O(k) for each
request, where n is the number of words in the dictionary and k an upper bound on the
length of the words.
In a first phase, we compute for each prefix p of a word in the dictionary the
total weight of all the words with prefix p. This weight is stored in a dictionary
total_weight. In a second phase, we store in a dictionary prop[seq] the prefix to
be proposed for a given sequence seq. A scan over the keys in total_weight allows
the determination of the prefix with greatest weight.
A principal ingredient is the function code_word, which for a given word returns
the corresponding sequence of digits.
To improve readability, the implementation below is in O(nk 2 ).
2.2 T9—Text on 9 Keys 45
t9 = "22233344455566677778889999"
# abcdefghijklmnopqrstuvwxyz mapping on the phone
def letter_to_digit(x):
assert ’a’ <= x <= ’z’
return t9[ord(x) - ord(’a’)]
def code_word(word):
return ’’.join(map(letter_to_digit, word))
def predictive_text(dic):
# total_weight[p] = total weight of words having prefix p
total_weight = {}
for word, weight in dic:
prefix = ""
for x in word:
prefix += x
if prefix in total_weight:
total_weight[prefix] += weight
else:
total_weight[prefix] = weight
# prop[s] = prefix to display for s
prop = {}
for prefix in total_weight:
code = code_word(prefix)
if (code not in prop
or total_weight[prop[code]] < total_weight[prefix]):
prop[code] = prefix
return prop
Problem
T9 [poj:1451]
46 Character Strings
Application
How should the words of a dictionary be stored in order to implement a spell checker?
For a given word, we would like to quickly find a dictionary word close to it in the
sense of the Levenshtein edit distance, see Section 3.2 on page 63. If we store the
dictionary words in a hash table, then we lose all proximity information between
words. It is better to store them in a lexicographic tree, also known as a prefix tree
or a trie (pronounced like ‘try’).
Definition
A trie is a tree that stores a set of words. The edges of a node towards its children are
labelled by distinct letters. Each word in the dictionary then corresponds to a path from
the root to a node in the tree. The nodes are marked to distinguish those corresponding
to words in the dictionary from those that are only strict prefixes of such words, see
Figure 2.2.
Spell Checking
With such a structure, it is easy to find a dictionary word that is at a distance dist from
a given word, for the Levenshtein edit distance defined in Section 3.2 on page 63. It
suffices to simulate the edit operations at each node, and invoke recursive calls to the
search with the parameter dist - 1.
a p
s o r a p
r e s or re
t e s t t e s t
Figure 2.2 A lexicographic tree for the words as, port, pore, pré, près, prêt (but without the
accents). In this figure, the label of each edge is indicated in the node of the child. The circles
with solid boundaries mark the nodes that correspond to words in the dictionary. On the right,
we show a Patricia trie for the same dictionary.
2.3 Spell Checking with a Lexicographic Tree 47
class TrieNode:
def __init__(self): # each node will have
self.is_word = False # 52 children -
self.s = {c: None for c in ascii_letters} # most will remain empty
Variant
A more complex structure exists that merges nodes as long as they have a single child.
A node is thus labelled with a word, rather than by a simple letter, see Figure 2.2. This
structure, optimal in memory and in traversal time, is known as a Patricia trie (see
Morrison, 1968).
Problem
Spell checker [icpcarchive:3872]
Definition
Given a string s of length n and a pattern t of length m, we want to find the first index
i such that t is a factor of s at the position i. The response should be −1 if t is not a
factor of s.
Complexity
O(n + m) (see Knuth et al., 1977).
Naive Algorithm
This consists of testing all possible alignments of t under s and for each alignment i
verifying character by character if t corresponds to s[i..i + m − 1]. The complexity
is O(nm) in the worst case. The following example illustrates the comparisons
performed by the algorithm for an example. Each line corresponds to a choice of
i, and indicates the characters of the ‘pattern motif’ that match, or an × for a
difference.
l a l o p a l a l a l i
0 l a l ×
1 ×
2 l ×
3 ×
4 ×
5 ×
6 l a l a
Observe that after handling a few i, we already know about a good portion of the
string s. We could use this information to skip the comparison of t[0] with s[i] in
the above example. This observation is extended to a study of the boundaries of the
2.5 Maximal Boundaries—Knuth–Morris–Pratt 49
entrée: abracadabra
sortie: abra abra
Definition
The boundary of a word w is a word that is at the same time a strict prefix and a strict
suffix of w, where by strict we mean that the length of the boundary must be strictly
less than the length of w. The maximal boundary of w is its longest boundary, and
is denoted β(w). For example, abaababaa has for boundaries abaa, a and the empty
word ε, hence β(abaababaa) = abaa. See Figure 2.3 for an illustration.
w: a b a a b a b a a
w: a b a a b a b a a
maximal boundary
Figure 2.3 Intuitively, the maximal boundary represents the longest overlap of a word with
itself: take two copies of the word w one over the other, and slide the second to the right until
the letters of the two words coincide. The portion of w corresponding to this overlap is then its
maximal boundary.
indices: 0 1 2 3 4 5 6 7 8
w: a b a a b a b a a
f0=0
f1=0
f2=1 a
f3=1 a
f4=2 a b
f5=3 a b a
f6=2 a b
f7=3 a b a
f8=4 a b a a
Figure 2.4 The array f of boundary lengths illustrated with an example. For a given string w, fi
is the largest value k such that the suffix wi−k+1 . . . wi is equal to the prefix w0 . . . wk−1 .
Key Observation
The relation of being a boundary is transitive: if v is a boundary of a boundary b of
u, then v is also a boundary of u. Moreover, if v is a boundary of u shorter than β(u),
then v is a boundary of β(u), see Figure 2.5. Hence, the iterated application of β to
a word w generates all the boundaries of w. For example, with w = abaababa, we
have β(w) = aba, then β(β(w)) = a and finally β(β(β(w))) = ε. This will prove very
useful for our algorithm.
ux u x
β(u) ? β(u) x
β(ux) v x v x
Figure 2.5 Visual explanation of the proof: every boundary of uwi can be written vwi where v is
a boundary of u, hence knowing the lengths of the boundaries of u allows the identification of
the longest boundary of uwi .
of u. Thus, to find the longest boundary of uwi it suffices to verify whether for each
boundary v of u in order of decreasing length, the word vwi is a boundary of uwi , see
Figure 2.5. We already know that vwi is a suffix of uwi and that v is a prefix of u,
hence it suffices to compare the two letters wk =? wi where k = |v|: indeed, wk is the
letter following the prefix v in uwi . How can we iterate over the boundaries v of u in
order of decreasing size? We simply successively apply the function β to the word u,
until we hit the empty word. Note that to perform this test, we only need to know the
lengths of the boundaries v of u, which have been stored in the array f .
Illustration
We use a window that exposes a letter wi of the first copy of w and a letter wk of the
second copy. Three cases must be considered:
1. If the two characters are equal, then we know that w0 · · · wk is the longest bound-
ary of w0 · · · wi , and set fi = k + 1, since |w0 · · · wk | = k + 1. The win-
dow is then shifted one step to the right in order to process the next prefix, see
Figure 2.6; this corresponds to an increment of both i and k.
2. If not, in the case wk = wi and k > 0, we move on to the next smaller boundary
β(w0 . . . wk−1 ), of size fk−1 . This boils down to shifting to the right the second
copy of w, until its contents coincide with the first along the whole of the left part
of the window, see Figure 2.6.
3. The final case corresponds to k = 0 and w0 = wi . At this point, we know that the
maximal boundary of w0 . . . wi is the empty word. We thus set fi = 0.
a b a a b a b a a
a b a a b a b a a
k
a b a a b a b a a
a b a a b a b a a
a b a a b a b a a
a b a a b a b a a
Figure 2.6 An example of the execution of the Knuth–Morris–Pratt algorithm. When the
window reveals two identical characters, we set fi = k + 1 and shift to the right, which comes
down to incrementing i and k. However, when the window exposes two different characters,
the bottom word must be shifted to the right, so that its prefix determined by the window is of
length fk−1 .
52 Character Strings
Code
All these cases can be captured in just a few lines of code. An outer while loop looks
for the boundary of a suitable prefix (case 2 above). We exit this loop for one of
two reasons. First, if k = 0 (no non-empty boundary was found), then only a single
character remains to test: wk =? wi . If yes, then fi = 1 (case 1); if not, fk = 0
(case 3). Second, we exit because a non-empty boundary is found where wk = wi and
fi = k + 1: we increment k and then assign k to fi . We thus ensure that at each start
of the loop for, the value of k corresponds to the length of the longest boundary of the
preceding prefix.
def maximum_border_length(w):
n = len(w)
f = [0] * n # init f[0] = 0
k = 0 # current longest border length
for i in range(1, n): # compute f[i]
while w[k] != w[i] and k > 0:
k = f[k - 1] # mismatch: try the next border
if w[k] == w[i]: # last characters match
k += 1 # we can increment the border length
f[i] = k # we found the maximal border of w[:i + 1]
return f
Complexity
This algorithm consists of a while loop nested in a for loop, which suggests a
quadratic complexity. However, the behaviour of the algorithm should be considered
using the example illustrated above. For each comparison wk =? wi in the algorithm,
either the word on the bottom is shifted to the right or the window is shifted to the
right. Each of these movements can only be executed at most |w| times, which shows
the linear complexity in |w|: we speak of amortised complexity, as the long iterations
are on average compensated by other shorter iterations.
The most important application of the maximal boundary algorithm is the search for
the first occurrence of a pattern t within a string s whose length is greater than |t|. The
2.5 Maximal Boundaries—Knuth–Morris–Pratt 53
naive approach consists of testing all the possible positions of t with respect to s. More
precisely, for each position i = 0 . . . n − m in s we test if the substring s[i : i + m]
is equal to t. This equality test involves comparing for every j = 0 . . . m − 1 the
character s[i + j ] with t[j ]. These two nested loops have complexity O(nm).
The Knuth–Morris–Pratt algorithm selects a letter # occurring neither in t nor in s
and we consider the lengths fi of the maximal boundaries of the prefixes of w = t#s.
Note that these boundaries can never be longer than t, due to the character #, hence
fi ≤ |t|. However, if ever fi = |t|, then we have found an occurrence of t in s. In
the positive case, the answer to the problem is the index i − 2|t|: the length of t is
subtracted twice from i, once to arrive at the start of the boundary and another time
to obtain the index in s rather than in w, see Figure 2.7. Note that this algorithm
is a bit different from the classic presentation of the Knuth–Morris–Pratt algorithm,
composed of a pre-computation of the array f of boundaries of t followed by a very
similar portion seeking the maximal alignments of t with the prefixes of s.
t#s: a b a a # a b a b a a b a a
f: 0 0 1 1 0 1 2 3 2 3 4 2 3 4
a b a a
a b a a
Figure 2.7 An execution of the maximal boundary search algorithm. Each occurrence of t in s
corresponds to a boundary of length |t| in the array f of lengths of maximal boundaries.
Note
All standard libraries of common programming languages propose a function that
looks for a pattern needle in a string haystack. However, the function provided by
Java 10 has worst-case complexity (nm), which is astonishingly expensive. We
invite the reader to measure, for their favourite language, the execution time as a
function of n of the search for the pattern a n b in the string a 2n , and protest if it is not
linear.
Given a word x, we would like to identify the largest integer k such that x can be
written as zk for a word z. If k > 1, then we say that x is not primitive.
We can use our Swiss army knife of the computation of boundaries to solve this
problem: if x can be written as y , then all the y p for p = 0, . . . , − 1 are boundaries
of x. It remains to prove that y −1 is the maximal boundary of x. So, if n is the length
of u and if n − fn divides n, the word is not primitive and the largest value of k we
seek is n/(n − fn ), see Figure 2.8.
u a b a a b a a b a
n− n
β(u) a b a a b a
a b a a b a
Figure 2.8 Knowing the maximal boundary of a periodic word allows the determination of its
smallest period, here aba.
Proof
Suppose that w = zk for k maximal and β(w) can be written as zk−1 q where qb = z
for a certain non-empty b. First, note that z = qb = bq (see Figure 2.9), hence
|b| ≤ |q|, otherwise zk−1 b would be a larger boundary of w. Hence, b is a boundary of
z and bz = zb = bqb. Thus, bw = wb, meaning that at the same time b is a boundary
of w and w is a boundary of bw. We now prove the key observation: as long as b is
smaller that w, then b is a boundary of w. We already know this is true for = 1,
and if it is true for b , then bb is a boundary of bw = wb; hence b+1 and w are
boundaries of bw. This means that either b+1 is a boundary of w, or w is a boundary
of b+1 . Let L be the largest integer for which bL is a boundary of w: then w = bL r =
rbL for a word r strictly smaller than b. If r is non-empty, then we have found a
boundary larger than the maximal boundary zk−1 q, a contradiction. Hence, r is empty
and since |b| < |z|, the factorisation bL is suitable for L > k, again a contradiction.
See Figure 2.9
2.5 Maximal Boundaries—Knuth–Morris–Pratt 55
def powerstring_by_border(u):
f = maximum_border_length(u)
n = len(u)
if n % (n - f[-1]) == 0: # does the alignment shift divide n ?
return n // (n - f[-1]) # we found a power decomposition
return 1
w
z z z q b
q z z z
z z z q
Figure 2.9 If w = z4 = z3 qb and we suppose that the maximal boundary of w is z3 q, then the
boundary property allows us to see that z = bq = qb.
Note that a shorter implementation exists using the substring search function inte-
grated with Python, see Figure 2.10.
def powerstring_by_find(u):
return len(x) // (x + x).find(x, 1)
a b a a b a a b a a b a a b a a b a
a b a a b a a b a
Figure 2.10 Detection of the first non-trivial position of x in xx allows the identification of its
smallest period, here 3.
Another classic problem consists in detecting whether two words x and y are
conjugate, i.e. if they can be written as x = uv and y = vu for words u and v. In
the affirmative case, we would like to find the decomposition of x and y minimising
the length of u. This boils down to simply looking for the first occurrence of y in the
word xx.
Problems
Find the maximal product of string prefixes [codility:carbo2013]
A Needle in the Haystack [spoj:NHAY]
Power strings [kattis:powerstrings]
Period [spoj:PERIOD]
56 Character Strings
Complexity
The expected time is O(n + m), but the worst case complexity is O(nm).
Algorithm
The Rabin–Karp algorithm 1987 is based on a completely different idea from
Knuth–Morris–Pratt. To find a pattern t in a large string s, we slide a window of
length len(t) over s and verify if the content of this window is equal to t. Since a
character-by-character test is very costly, we instead maintain a hash value for the
contents of the current window and compare it to the hash value of the search string t.
If ever the hash values coincide, we proceed to the expensive character-by-character
test, see Figure 2.11. To obtain an interesting complexity, it is necessary to efficiently
update the hash value as the window shifts. We thus use what is known as a rolling
hash function.
If the hash function, with values in {0,1, . . . ,p − 1}, is well-chosen, we would
expect a collision—i.e. when two distinct strings u,v with the same size, selected
uniformly, give the same hash value—to occur with probability on the order of 1/p.
In this case, the mean complexity of the algorithm is O(n + m + m/p). Our imple-
mentation uses p on the order of 256 , hence in practice the complexity is O(n + m),
but in the worst case, it is O(nm).
Figure 2.11 The idea of the Rabin–Karp algorithm is to first compare the hash values between t
and a window on s before performing a costly character-by-character comparison.
character amounts to subtracting the first term, shifting the characters to the left
corresponds to a multiplication by 128 and modifying the last character is done by
adding a term. Consequently, shifting the window on s, updating the hash value and
comparing it with that of t takes constant time.
Note in the code below the addition of the value DOMAIN * PRIME (which is of
course 0 mod p) to ensure that the calculations remain in the non-negative integers.
This is not strictly necessary in Python but is required in languages such as C++, where
the calculation of modulo can return a negative number.
Next, the implementation of the Rabin–Karp algorithm itself begins with the com-
putation of the hash values of t and of the first window on s, followed by a loop over
all the factors of s, until a match is found.
hash_s = 0
# Now check for matching factors in s
for i in range(k): # preprocessing
hash_s = (DOMAIN * hash_s + ord(s[i])) % PRIME
for i in range(len(s) - k + 1):
if hash_s in pos: # is this signature in s?
for j in pos[hash_s]:
if matches(s, t, i, j, k):
return (i, j)
if i < len(s) - k:
hash_s = roll_hash(hash_s, ord(s[i]), ord(s[i + k]), last_pos)
return None
Problem
Longest Common Substring [spoj:LCS]
input: babcbabcbaccba
output: abcbabcba
Definition
A word s is a palindrome if the first character of s is equal to the last, the second is
equal to the next-to-last and so on.
The problem of the longest palindrome consists in determining the longest factor
that is a palindrome.
60 Character Strings
Complexity
This problem can be solved in quadratic time with the naive algorithm, in time
O(n log n) with suffix arrays, but in time O(n) with Manacher’s algorithm (1975),
described here.
Algorithm
First, we transform the input s by inserting a separator # around each character and
by adding sentinels ^ and $ around the string. For example, abc is transformed into
^#a#b#c#$. Let t be the resulting string. This allows us to process in an equivalent
manner palindromes of both even and odd length. Note that with this transformation,
every palindrome begins and ends with the separator #. Thus, the two ends of a
palindrome have indices with the same parity, which simplifies the transformation
of a solution on the string t into one on the string s. The sentinels avoid special care
for the border cases.
The word nonne contains a palindrome of length 2 (nn) and one of length 3 (non). Their
equivalents in
|-----|
^#n#o#n#n#e#$
|---|
The output of the algorithm is an array p indicating for each position i, the largest
radius r such that the factor from i − r to i + r is a palindrome. The naive algorithm
is the following: for each i, we initialise p[i] = 0 and increment p[i] until we find the
longest palindrome t[i − p[i], . . . ,i + p[i]] centred on i.
c d
j i
Figure 2.12 Manacher’s algorithm. Having already computed p for the indices < i, we wish to
compute p[i]. Suppose there is a palindrome centred on c of radius d − c with d maximal, and
let j be the mirror image of i with respect to c. By symmetry, the palindrome centred on j of
radius p[j ] must be equal to the word centred on i, at least up to the radius d − i.
Consequently, p[j ] is a lower bound for the value p[i].
Manacher’s improvement concerns the initialisation of p[i]. Suppose we already
know a palindrome centred on c with radius r, hence terminating on the right at
2.7 Longest Palindrome of a String—Manacher 61
d = c + r. Let j be the mirror image of i with respect to c, see Figure 2.12. There is
a strong relation between p[i] and p[j ]. In the case where i + p[j ] is not greater than
d, we can initialise p[i] by p[j ]. This is a valid operation, as the palindrome centred
on j of radius p[j ] is included in the first half of the palindrome centred on c and of
radius d − c; hence it is also found in the second half.
After having computed p[i], we must update c and d, to preserve the invariant that
they code a palindrome with d − c maximal. The complexity is linear, since each
comparison of a character is responsible for an incrementation of d.
def manacher(s):
assert set.isdisjoint({’$’, ’^’, ’#’}, s) # Forbidden letters
if s == "":
return (0, 1)
t = "^#" + "#".join(s) + "#$"
c = 1
d = 1
p = [0] * len(t)
for i in range(2, len(t) - 1):
# -- reflect index i with respect to c
mirror = 2 * c - i # = c - (i-c)
p[i] = max(0, min(d - i, p[mirror]))
# -- grow palindrome centered in i
while t[i + 1 + p[i]] == t[i - 1 - p[i]]:
p[i] += 1
# -- adjust center if necessary
if i + p[i] > d:
c = i
d = i + p[i]
(k, i) = max((p[i], i) for i in range(1, len(t) - 1))
return ((i - k) // 2, (i + k) // 2) # extract solution
Application
A man is walking around town, and his smartphone registers all of his movements.
We recover this trace and seek to identify a certain type of trajectory during the day,
notably round trips between two locations that use the same route. For this, we extract
from the trace a list of street intersections and check it for palindromes.
Problems
Longest Palindromic Substring [spoj:LPS]
Casting Spells [kattis:castingspells]
3 Sequences
Definition
Consider a grid of (n + 1) × (m + 1) cells, labelled (i,j ) with 0 ≤ i ≤ n and
0 ≤ j ≤ m. The cells are linked by weighted arcs: the predecessors of a cell (i,j )
are (i − 1,j ),(i − 1,j − 1) and (i,j − 1), with the exception of the cells in the first
column and first row (see Figure 3.1). The goal is to find the shortest path from (0,0)
to (n,m).
1 2 3 1 1 1 1
0 1 2 1 0 2 1
Figure 3.1 Shortest path in a directed grid, illustrated by the heavy arcs.
Algorithm in O(nm)
As this graph is directed and acyclic (abbreviated DAG, for directed acyclic graph),
we can solve this problem by dynamic programming, by computing the distances from
3.2 The Levenshtein Edit Distance 63
(0,0) to all of the cells (i,j ) in a certain order. First of all, the distances to the cells of
the first row and first column are easy to compute, since a unique path towards these
cells exists. Next, we compute the distances to the cells (i,j ) with 1 ≤ i ≤ n,1 ≤
j ≤ m in lexicographical order. Each distance from (0,0) to an (i,j ) is the minimum
of three alternatives, which can be determined in constant time, as the distances to the
predecessors of (i,j ) have already been computed.
Variants
Several classic problems can be reduced to this simple problem, as described in the
following sections.
Problem
Philosophers Stone [spoj:BYTESM2]
Definition
Given two strings x,y, we wish to know how many operations (insertion, deletion, or
substitution) are required to transform x into y. This distance is used, for example, in
the Unix command diff, which displays line-by-line a minimal number of operations
to transform one file into another.
Algorithm in O(nm)
For n = |x|,m = |y|, we present an algorithm in O(nm) using dynamic programming,
see Figure 3.2. We construct an array A[i,j ] giving the distance between the prefix of
x of length i and the prefix of y of length j . To begin, we initialise A[0,j ] = j and
A[i,0] = i (the distances from the empty string to the various prefixes of x and y).
Then, in general, when i,j ≥ 1, there are three possibilities for the last letters of the
prefixes. Either xi is deleted, or yj is inserted (at the end), or xi is replaced by yj (if
they are not already equal). This gives the following recurrence formula, where match
is a Boolean function which returns 1 if its arguments are identical:
⎧
⎨ A[i − 1,j − 1] + match(xi ,yj )
⎪
A[i,j ] = min A[i,j − 1] + 1
⎪
⎩
A[i − 1,j ] + 1.
This function encodes the cost of a letter substitution, which can be adjusted if desired.
For example, the cost could depend on the distance between letters on the keyboard.
64 Sequences
L A D A
0 1 1 1 2 1 3 1 4
A 1 1 1 0 1 1 1 0 1
1 1 1 1 1 1 2 1 3
U 1 1 1 1 1 1 1 1 1
2 1 2 1 2 1 2 1 3
D 1 1 1 1 1 0 1 1 1
3 1 3 1 3 1 2 1 3
I 1 1 1 1 1 1 1 1 1
4 1 4 1 4 1 3 1 3
Figure 3.2 Shortest path in a directed graph, determining the edit distance between two words.
Sequence of Operations
As well as computing the edit distance between the strings, it would be nice to know
the operations necessary to transform x into y. This can be done as in the search
for the shortest path in a graph. By examining the distances of its predecessors, we
can determine a choice that realises a minimal distance to a node. Thus we can trace
the path from the node (n,m) back to (0,0), and along this path produce the list of
operations corresponding to the optimal edits. At the end, it suffices to reverse this list.
Implementation Details
Note that in the description of the dynamic program, the dimension of the array A is
(n + 1) × (m + 1), whereas the strings are of length n and m. This is because A takes
into account all prefixes of x and y, including the empty string . So x has n + 1
prefixes, and y has m + 1.
3.3 Longest Common Subsequence 65
To Go Further
Algorithms with better performance in practice have been proposed. For example, if
an upper bound s is known for the edit distance, then the above dynamic program can
be restricted to entries in A at a distance of at most s from the diagonal, to obtain a
complexity of O(s min{n,m}) (Ukkonen, 1985).
Problems
Edit distance [spoj:EDIST]
Advanced Edit Distance [spoj:ADVEDIST]
Definition
Let be a set of symbols. For two sequences s,x ∈ , s is said to be a subsequence
of x if indices i1 < . . . < i|s| exist such that xik = sk for every k = 1, . . . ,|s|. Given
two sequences x,y ∈ , we wish to find a maximal length sequence s ∈ , which
is a subsequence of both x and y.
66 Sequences
Another way to see the problem involves matchings, see Section 9.1 on page 139.
We seek a maximum matching between equal letters of x and y such that the links
between the matches do not cross each other, see Figure 3.3.
Application
The program ‘diff’ reports differences between two files, expressed as a minimal list
of line changes to bring either file into agreement with the other (Hunt and McIlroy,
1976). To do this, it is necessary to identify the largest common subsequence of lines
between both files
This problem also arises in bioinformatics, for example when aligning two chains
of nitrogenous bases.
A B C D A E
A E D B E A
Figure 3.3 The problem of the longest common subsequence can be seen as a maximal
alignment without crossings between equal elements of the two given sequences. For example,
an additional alignment of the letters B would cross the alignment of the letters D.
Algorithm in O(nm)
We present an algorithm in O(nm), where n = |x|,m = |y|. The idea is to compute
for every 0 ≤ i ≤ n and 0 ≤ j ≤ m the longest common subsequence of the prefixes
x1 . . . xi and y1 . . . yj . This gives n·m subproblems. The optimal solution for (i,j ) can
be computed from the solutions for (i − 1,j ),(i,j − 1), and (i − 1,j − 1), in constant
time. We can thus solve the problem for (n,m) in time O(nm). The algorithm is based
on the following observation.
Key Observation
Let A[i,j ] be a longest common subsequence of x1 . . . xi and y1 . . . yj . If i = 0 or
j = 0, then A[i,j ] is empty. If xi = yj , then necessarily one of the letters xi ,yj
is not matched in an optimal solution and A[i,j ] is the longest sequence between
A[i − 1,j ] and A[i,j − 1]. If xi = yj , then an optimal solution exists that makes
a match between these letters and A[i,j ] is A[i − 1,j − 1] · xi . Here, the symbol
‘·’ represents concatenation, and by maximum we mean the sequence of maximal
length.
Implementation
In the simplified implementation given below, the variable A[i][j] does not represent
the sequence A[i,j ], but instead represents its length.
3.3 Longest Common Subsequence 67
In Practice
The algorithm BLAST (for Basic Local Alignment Search Tool) is widely used, but it
does not always guarantee to produce an optimal solution.
Problems
Longest Common Substring [spoj:LCS]
Longest Common Subsequence [spoj:LCS0]
68 Sequences
Definition
Given a sequence of n integers x, find a strictly increasing subsequence of s of
maximal length.
Application
Consider a straight road leading to the sea, with houses built along the road. Each
house has several floors and enjoys a sea view if all the houses between it and the sea
have fewer floors. We wish to give each house a sea view. What is the minimal number
of houses that must be demolished to reach this goal (see Figure 3.4)?
Figure 3.4 Demolish the least number of houses so that the rest have a sea view.
such that b[k − 1] < xi < b[k], then we can improve the sequence of length k − 1 by
appending xi , and obtain a better subsequence of length k, since it ends with a smaller
element. If xi is greater than all the elements of b, we can extend b with the element xi .
This is the only possible improvement, and the search for the index k can be done by
dichotomy in time O(log |b|), where |b| is the size of b.
Implementation Details
The subsequences are represented by linked lists with the aid of arrays h and p. The
heads of the lists are held in h, such that b[k] = x[h[k]], while the predecessor of an
element j is p[j ]. The lists are terminated by the constant None, see Figure 3.5.
index 0 1 2 3 4
b − 1 2 4 5
h
x 3 1 4 1 5 9 2 6 5 4 5 3 9 7 9
p
Figure 3.5 Computation of the longest increasing subsequence. In grey the input x being
processed. For the prefix considered, the sequences so far are (),(1),(1,2),(1,2,4),(1,2,4,5).
After processing the element 3, the sequence (1,2,4) is dethroned by (1,2,3).
def longest_increasing_subsequence(x):
n = len(x)
p = [None] * n
h = [None]
b = [float(’-inf’)] # - infinity
for i in range(n):
if x[i] > b[-1]:
p[i] = h[-1]
h.append(i)
b.append(x[i])
else:
# -- binary search: b[k - 1] < x[i] <= b[k]
k = bisect_left(b, x[i])
h[k] = i
b[k] = x[i]
p[i] = h[k - 1]
# extract solution in reverse order
q = h[-1]
s = []
while q is not None:
s.append(x[q])
q = p[q]
return s[::-1] # reverse the list to obtain the solution
70 Sequences
Problem
Easy Longest Increasing Subsequence [spoj:ELIS]
Definition
Consider a game over a stack of positive integers played by two players (in this
instance, both players are female, in order to avoid convoluted grammar), see
Figure 3.6. Player 0 starts. If the stack is empty, she loses. If not, then if the top
of the stack contains an integer x, she has the choice between popping one element
or popping x times. The latter option is authorised only if the stack contains at least x
elements. Then it is the turn of Player 1 to play, and so on. Suppose you are given a
stack P with n integers; the question is whether Player 0 has a winning strategy, i.e. if
she can play in a way that wins no matter what Player 1 chooses to do.
4 3 2 1 3 2 4 1
Complexity
The complexity is linear using dynamic programming.
that Player 0 has a winning strategy if she manages in one shot to place Player 1 in a
position where she can never win.
The base case is G[0] = False, and the inductive case for i > 0 is
G[i − 1] ∨ G[i − P [i]] if P [i] ≤ i
G[i] =
G[i − 1] otherwise.
By a simple scan in a linear number of steps, the array G can be populated, and the
problem is solved by returning G[n − 1].
4 Arrays
Arrays figure as one of the most important of the elementary data types. For many
simple problems, no other data structure is required. In Python, arrays are stored in
the data type called list. The choice of this name might be confusing, because Python
lists have nothing to do with linked lists, which are standard objects in C++ (std::list )
and Java (LinkedList ). Typically, a linked list allows the deletion or insertion of a
given element at a given position in constant time. However, inserting into a Python
list requires building in linear time a new list consisting of a prefix, the new element
and a suffix. The elements of a Python list are indexed from 0. This must be taken into
consideration when reading or displaying data, if ever the numbering of the elements
begins at 1 in the format of the inputs or outputs.
In what follows, we use ‘array’ as the generic language-independent term, and ‘list’
when referring to the specific Python data structure.
An array can be used to represent a binary tree in a very simple manner. The
parent of a node at index i is found at i/2, its left child at index 2i and its right
child at index 2i + 1. The root is at index 1 and the element at index 0 is simply
ignored.
This section deals with classic problems on arrays, and presents data structures to
perform operations on intervals of indices known as ranges, for example, the compu-
tation of a minimal value within a range. Two of the sections describe dynamic data
structures used to provide efficient operations of element modification and queries
over such ranges.
Recall: In Python, the last element of a list t can be accessed with the index -1.
A reversed copy of t is obtained with t[:: −1]. A new element is added with the
method append. In this sense, Python lists are similar to the class ArrayList
of Java or the class vector of C++.
Definition
Given two sorted lists x and y, produce a sorted list z, containing the elements of
x and y.
Application
This operation is useful for the merge sort algorithm. To sort a list, we divide it into
two portions of equal size (up to a difference of one element if the length is odd). The
two portions are then sorted recursively, and the sorted results are merged using the
procedure described below. This algorithm has an optimal complexity of O(n log n)
comparisons.
17
10 12
7
1 3 4 6
8
11 15
18
Figure 4.1 Merging two lists works a bit like a zipper with minimum selection instead of
alternation.
Problem
Mergesort [spoj:MERGSORT]
Definition
Each request is of the form of an interval of indices [i,j ) and must return the sum of
the elements of t between the indices i (included) and j (excluded).
Definition
Each request is in the form of an interval of indices [i,j ) and we must either find an
element x of t appearing at least twice within this interval or announce that all the
elements are distinct.
index 0 1 2 3 4 5 6 7 8 9
t a b a a c d b a b c
p −1 −1 0 2 −1 −1 1 3 6 4
q −1 −1 0 2 2 2 2 3 6 6
Figure 4.2 The algorithm for finding duplicates on an example. In the interval [3,6) all elements
of t are distinct because q[5] < 3. However, the interval [4,10) contains duplicates, for
example, b at the positions q[9] and p[q[9]].
4.5 Query for the Minimum of a Range—Segment Tree 75
Problem
Unique Encryption Keys [icpcarchive:5881]
Definition
This statistical problem seeks to find, for an array of values t, the maximum of t[i] +
t[i + 1] + · · · + t[j ] over every pair of indices i,j with i ≤ j .
Algorithm in O(n)
This solution by dynamic programming was found by Jay Kadane in 1984. For each
index j we look for the maximal sum t[i] + · · · + t[j ] over all indices 0 ≤ i ≤ j .
Denote A[j ] as this value. This value must be either t[j ] itself, or made up of t[j ] plus
a sum of the form t[i] + · · · + t[j − 1], which is itself maximal. Hence, A[0] = t[0],
as this is the only option, and for j ≥ 1, we have the recurrence relation A[j ] =
t[j ] + max{A[j − 1],0}.
Variants
This problem can be generalised to matrices. Given a matrix M of dimension n × m,
a range of row indices [a,b] and a range of column indices [i,j ] define a rectangle
[a,b] ⊗ [i,j ] within the matrix. We wish to find the rectangle with the largest sum.
For this, it suffices to loop over all pairs of row indices [a,b],a ≤ b, and for each pair
generate an array t of length m, such that t[i] = M[a,i] + · · · + M[b,i]. This array
can be obtained in time O(m) from the array of the previous iteration corresponding
to the pair [a,b − 1], by simply adding the bth row of M. With the aid of the above
algorithm, we can find in time O(m) a range of columns [i,j ] which describes a
rectangle [a,b] ⊗ [i,j ] of maximal sum bounded by the rows a and b. The largest
sum seen when looping over a and b gives the maximal sum for M. This approach
produces a solution with complexity O(n2 m).
Problems
Maximum Sum Sequences [spoj:MAXSUMSQ]
Largest Rectangle [codechef:LARGEST]
Definition
We wish to maintain a data structure storing an array t of n values and allowing the
following operations:
Variant
With only minor changes we can also determine the index of the minimum element
along with its value.
Figure 4.3 The tree used in the structure to determine the minimum in a range.
The analysis terminates with the observation that the children of an overlapping node
are always either full or empty, and that only the root can be a strict node.
class RangeMinQuery:
def __init__(self, t, INF=float(’inf’)):
self.INF = INF
self.N = 1
while self.N < len(t): # find size N
self.N *= 2
self.s = [self.INF] * (2 * self.N)
for i in range(len(t)): # store values of t
self.s[self.N + i] = t[i] # in the leaf nodes
for p in range(self.N - 1, 0, -1): # fill inner nodes
self.s[p] = min(self.s[2 * p], self.s[2 * p + 1])
A weighted variant of this structure is put to use in order to efficiently calculate the
area of a union of rectangles, see Section 13.5 on page 205.
Problem
Negative Score [spoj:RPLN]
Definition
We wish to maintain a data structure t storing n values and allowing the following
operations for a given index a ∈ {0, . . . ,n − 1}:
78 Arrays
• update t[a].
• calculate t[0] + · · · + t[a].
For technical reasons, internally we work with an array T whose indices vary from
1 to n, with t[a] = T [a + 1]. Hence, in Python, T [0] is ignored, and the length of T
will actually be n + 1!. The first statement in the access methods to our structure is a
shift of the index, in order to hide this technical detail from the users.
This problem can be solved using a segment tree, as described in the previous sec-
tion. By replacing the constant ∞ by 0 and the operation min by an addition, we obtain
a structure that solves the problem with a complexity O(log n) per request. The struc-
ture described in this section has a similar performance, but is quicker to implement.
Variant
With a slight modification, this structure can also be used to perform the following
operations in time O(log n).
• add a value across the whole of an interval t[a],t[a + 1], . . . ,t[b] for given indices
a,b.
• request the value of t[a] for a given index a.
For this, it suffices to use a Fenwick tree storing an array t , with t[a] = t [0] +
· · · + t [a]. Hence, reading t[a] reduces to calculating the sum of a prefix in t , and
adding a value to a suffix in t involves modifying a value in t . Hence, adding a value
to an interval of indices in t comes down to modifying two values in t at the endpoints
of the interval.
1 We usually use subscript 2 to denote the binary representation, but it is omitted if no confusion is
possible.
4.6 Query the Sum over a Range—Fenwick Tree 79
• The left neighbour of the index i = x10 is j = x00 . The interval I (j ) is the one
k k
16
T 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Figure 4.4 Example of a Fenwick tree. The relation ‘is parent of’ can be read vertically from
top to bottom. For example, 8 is the parent of 4, 6 and 7. An index i is left neighbour of j if the
interval I (i) touches I (j ) on the left. For example, 4 is the left neighbour of 5 and 6.
80 Arrays
class Fenwick:
def __init__(self, t):
self.s = [0] * (len(t) + 1) # create internal storage
for a, v in enumerate(t):
self.add(a, v) # initialize
# variante:
def intervalAdd(self, a, b, val):
self.add(a, +val)
self.add(b + 1, -val)
Definition
Given a sequence x of n elements and an integer k, we wish to determine all the
intervals [i,j ), maximal by inclusion, such that xi , . . . ,xj −1 include exactly k distinct
elements, see Figure 4.5.
1 2 1 3 3 2 1 2 3
i j
Application
A cache is a rapid-access memory which is placed in front of a slower memory. The
address space of the slow memory is partitioned into equal-sized blocks, known as
4.7 Windows with k Distinct Elements 81
pages. The cache has limited capacity and can only hold k pages. The memory access
of a processor over time forms a sequence x whose elements are the pages requested.
When a requested page is in the cache, access to it is very fast, otherwise we speak
of a cache miss and the page must be loaded into the cache from slow memory in
place of another page (one that is used infrequently or has not been used in a long
time). Looking for intervals in x containing exactly k distinct elements reduces to
finding the intervals of time during which no cache misses have occurred, under the
hypothesis of a favourable initial cache configuration.
Algorithm in O(n)
The idea is to run over the sequence x with two cursors i and j that define a window.
We maintain the number of distinct elements of the set xi , . . . ,xj −1 in a variable dist,
with the aid of a counter of occurrences occ for each element. When dist exceeds the
parameter k, we advance i, otherwise we advance j . The following implementation
returns an iterator, and can be used by another function to individually process each
interval. As each of the two cursors i and j only advance, we perform at most 2n
operations, hence the linear time complexity.
Variant
Given a ring of symbols (a list of characters looping back on itself), what is the length
of the smallest interval containing at least one instance of each symbol? This is indeed
a variant, as it suffices to run through the list concatenated with itself and determine
the smallest interval containing k distinct elements, where k is the total number of
symbols. A solution with two cursors again has linear complexity.
Problem
Épiphanie [prologin:2011:epiphanie]
5 Intervals
Definition
The problem consists of storing n given intervals in a structure in order to rapidly
answer queries of the following form: for a given value p, what is the list of all the
intervals containing p? We suppose that all the intervals are of the half-open form
[l,h), but the structure can be adapted to other forms.
center
1 5
3 7
4 5
0 3 6 9
7 10
9 10
Choice of center
For the binary tree to be balanced, we can choose center as the median value of the
left endpoints of the intervals to be stored. Then half the intervals will be stored in the
right subtree, guaranteeing a logarithmic depth. In a manner similar to the behaviour
of the sort algorithm quicksort, if the centre is chosen randomly among the left
endpoints, the performance will be similar in expectation.
Complexity
The construction of the tree takes time O(n log n), and the processing of a query
costs O(log n + m), where the logarithmic term is due to the binary search in the
sorted lists.
Implementation Details
In our implementation, the intervals are represented by n-tuples whose first two ele-
ments contain the endpoints of an interval. Other elements can be added to transport
supplementary information.
The binary search is done via the function bisect_right(t,x), which returns i
such that t[j ] > x if and only if j > i. Note that you should not loop in the array
by_high with the aid of the sublist selection by_high[i:], as the creation of the sublist
is linear in len(by_high), thus increasing the complexity to O(log n + n) instead of
O(log n + m), where m is the size of the list returned.
The arrays contain pairs (value, interval), thus we search with bisect_right the
insertion point of an element x of the form (p,(∞,∞)).
84 Intervals
class _Node:
def __init__(self, center, by_low, by_high, left, right):
self.center = center
self.by_low = by_low
self.by_high = by_high
self.left = left
self.right = right
def interval_tree(intervals):
if intervals == []:
return None
center = intervals[len(intervals) // 2][0]
L = []
R = []
C = []
for I in intervals:
if I[1] <= center:
L.append(I)
elif center < I[0]:
R.append(I)
else:
C.append(I)
by_low = sorted((I[0], I) for I in C)
by_high = sorted((I[1], I) for I in C)
IL = interval_tree(L)
IR = interval_tree(R)
return _Node(center, by_low, by_high, IL, IR)
Variant
If the goal is only to determine the number of intervals containing a given value,
and not the list, then a much simpler solution exists with a sweep line algorithm. We
sweep the intervals contained in a node from left to right, advancing from one end to
the other. At any given moment, we maintain the number of intervals open (i.e. the
right endpoint has not yet been seen), and thus produce an array containing pairs of
5.3 The Interval Point Cover Problem 85
the form (x,k), such that for two successive pairs (x,k),(x ,k ), the number of intervals
containing p is k for x ≤ p < x .
Definition
Given a set S of n intervals, we wish to determine their union, in the form of an ordered
list L of disjoint intervals, with the property I ∈S = I ∈L .
Algorithm
The complexity is O(n log n) using a sweep line algorithm. We sweep the endpoints
of the intervals from left to right. At any given moment, we maintain in nb_open the
number of intervals open (i.e. the right endpoint has not yet been seen). When this
number becomes zero, we add to the solution a new interval [last,x), where x is the
current position of the sweep line and last is the last sweep position where nb_open
became positive.
Implementation Details
Note the order in which the endpoints of the intervals are processed. This order is
correct when the intervals are closed or half-open. For open intervals, the end of the
interval (x,y) must be handled before the beginning of the interval (y,z).
def intervals_union(S):
E = [(low, -1) for (low, high) in S]
E += [(high, +1) for (low, high) in S]
nb_open = 0
last = None
retval = []
for x, _dir in sorted(E):
if _dir == -1:
if nb_open == 0:
last = x
nb_open += 1
else:
nb_open -= 1
if nb_open == 0:
retval.append((last, x))
return retval
Application
We need to install a minimal number of antennas along a straight beach in order to
provide coverage to a number of small offshore islands (as small as individual points!).
86 Intervals
An antenna can cover all of the islands within a given radius r (too bad for islands
further out than r), which is the same for all the antennas.
Observation
The intersection of the beach with a circle of radius r drawn around an island defines
an interval on the beach that must contain an antenna. The problem thus reduces to the
following, see Figure 5.2.
Definition
For n given intervals, find a minimal set of points S such that each interval intersects S.
Algorithm
The complexity is O(n log n) using the sweep line technique. We process the intervals
in increasing order of their right endpoint. At every step, we maintain a solution S
for the intervals already seen, which minimises |S| and in case of equality maximises
max S.
The algorithm is simple. If for an interval [l,r] we have l ≤ max S, then we do
nothing, otherwise we add r to S. The idea is that in any case we must cover [l,r],
and by choosing the largest value that covers it, we increase the chance to additionally
cover subsequent intervals.
To be convinced of the optimality, let Sn be the solution constructed for the set of
intervals I1, . . . ,In already processed. Obviously, S1 is optimal for I1 . Suppose Sn is
optimal, and consider In+1 = [l,r]. Either l ≤ max Sn , in which case Sn is already an
optimal solution for In+1 , or l > max Sn , and |Sn+1 | = |Sn | + 1. In the latter case,
sea
beach
Figure 5.2 How many antennas are needed to cover all the islands? This problem reduces to the
interval point cover problem.
5.3 The Interval Point Cover Problem 87
suppose Sn+1
is an optimal solution for I1, . . . ,In+1 . Then Sn+1
\ {max Sn+1 } must
cover I1, . . . ,In , and hence, by the optimality of Sn , |Sn+1 | = |Sn | + 1 = |Sn+1 |.
def interval_cover(I):
S = []
# sort by right endpoints
for start, end in sorted(I, key=lambda v: v[1]):
if not S or S[-1] < start:
S.append(end)
return S
Problem
Radar Installation [onlinejudge:1193]
6 Graphs
A simple manner to encode a graph is to identify the n vertices by the integers from
0 to n − 1. However, it is common to number the vertices starting from 1 when
displaying a graph, or when reading from an input file. Thus, remember to add or
subtract 1 from the indices when displaying or reading a graph.
# adjacency list
G = [[1,2],[0,2,3],[0,1],[1]]
0
# adjacency matrix
1 3 G = [[0,1,1,0],
[1,0,1,1],
2 [1,1,0,0],
[0,1,0,0]]
The algorithms presented in this book generally work with adjacency lists.
In a directed graph, at times we require two structures G_out, G_in, containing the
arcs originating from and terminating on each vertex. However, rather than storing the
arc, we store the endpoints of the arcs. Hence, for each vertex u, G_out[u] contains
the list of vertices v for each outgoing arc (u,v) and G_in[u] the list of vertices v for
each incoming arc (v,u).
A simple way to store labels on the vertices and the edges is to use complementary
arrays indexed by the vertices or matrices indexed by pairs of vertices. In this way
the structure of the encoding of the graph G itself is not affected, and G can be
used without modification in the portions of the code that are not concerned with the
labels.
Certain implementations of graph algorithms require the vertices to be simply
the integers between 0 and n − 1, whereas at times the vertices need to be iden-
tified by a name or by a more complex but immutable object, such as a string
or an n-tuple of integers. To provide an interface between these code fragments,
a small class can be written, translating between the index of a vertex and the
complex object that it represents. Such a class would contain a dictionary name2node
relating the name of a vertex to its index, and an array node2name providing
the inverse function. If G is an instance of the class Graph, then the expression
len(G) should return the number of its vertices and G[u] the adjacency list of
the vertex with index u. Then the class behaves exactly like an adjacency list,
and in addition allows vertices and edges to be added using the names of the
vertices.
90 Graphs
class Graph:
def __init__(self):
self.neighbors = []
self.name2node = {}
self.node2name = []
self.weight = []
def __len__(self):
return len(self.node2name)
At times, the graph is given implicitly, for example, in the form of a grid, where the
vertices are the cells of the grid and the edges correspond to the adjacent cells, as in a
labyrinth. Another example of a graph given implicitly is with combinatorial objects,
where an arc corresponds to a local modification.
Exit
Figure 6.2 A configuration of the game Rush Hour.
vectors reachable in one move. The encoding of the configuration given in Figure 6.2
is the following.
For example, if orient[i] = 0, then the vehicle i occupies all the cells (x,y) for
coovar[i] ≤ x < coovar[i] + length[i] and y = coofix[i]. If orient[i] = 1, then
the vehicle i occupies all the cells (x,y) for x = coofix[i] and coovar[i] ≤ y <
coovar[i] + length[i].
Problem
Ricochet Robots [kattis:ricochetrobots]
Definition
Depth-first search is an exploration of a graph that begins at a given vertex and
recursively explores its neighbours.
Complexity
The complexity in time is O(|V | + |E|).
Application
The principal application is to discover all the nodes accessible starting from a given
vertex of a graph. This technique is the basis for several algorithms that we will tackle
further on, such as the detection of biconnected components or topological sort, see
Sections 6.6 on page 97 and 6.7 on page 102.
Implementation Details
In order to avoid repeatedly exploring the neighbours of a vertex, we mark the vertices
already visited, using a Boolean array.
92 Graphs
Improved Implementation
The above recursive implementation cannot handle very large graphs, as the call
stack is limited. In Python, the function setrecursionlimit allows this limit to be
somewhat exceeded, but, in general, to not more than a depth of some thousands of
recursive calls. To overcome this problem and for better efficiency, we propose an iter-
ative implementation. The stack to_visit contains the vertices discovered but not yet
processed.
Problems
ABC Path [spoj:ABCPATH]
A Bug’s Life [spoj:BUGLIFE]
6.4 Breadth-First Search—BFS 93
Definition
Rather than exploring as far as possible starting from the current node (depth-first
search), we can enumerate the nodes of a graph by their increasing distance from an
initial node (breadth-first search, abbreviated BFS).
Key Observation
We must process the nodes by increasing distance from an initial node, hence we need
a data structure to efficiently maintain this order. A queue is a good choice: if for every
vertex extracted from the head of the queue we add its neighbours to the end of the
queue, it is easy to see that at a given moment, it only contains nodes at a distance d
at the head and d + 1 at the tail, since as long as vertices at distance d remain at the
head, there can only be vertices at distance d + 1 added to the queue.
Implementation Details
The principal interest of breadth-first traversal is to determine the distances away from
a given starting node in a non-weighted graph. Our implementation computes these
distances as well as the predecessors in the shortest path tree. The array of distances
is also used to mark the vertices visited.
Problem
Hike on a Graph [spoj:HIKE] Prime Path [spoj:PPATH]
94 Graphs
Definition
A subset A ⊂ V of a graph is called a connected component if for every pair of
vertices u,v in A, there exists a path from u to v. We could, for example, wish to know
the number of connected components of a graph. A graph is said to be connected if it
has exactly one connected component.
Application
Consider a photo taken from above of a die on a table. We would like to easily
determine the number of pips on the exposed face. For this, we posterise1 the image
with a greyscale threshold to obtain an image in black and white, so that each pip
corresponds to a connected component in the image, see Figure 6.3.
Figure 6.4 shows the logo of the British electronic music group Clean Bandit in
ASCII art. It can be seen as a graph on a grid whose vertices are the sharps and two
vertices are linked by an edge if and only if they touch, horizontally or vertically. This
graph is made up of four connected components.
........................................ ........................................
...................#........#####....... ...................1........22222.......
..................###......#######...... ..................111......2222222......
.................#####....#########..... .................11111....222222222.....
................#######....#######...... ................1111111....2222222......
...............#########....#####....... ...............111111111....22222.......
..............###########............... ..............11111111111...............
.............#############.............. .............1111111111111..............
............###############............. ............111111111111111.............
.............#############.............. .............1111111111111..............
..............###########............... ..............11111111111...............
.........#.....#########................ .........3.....111111111................
........###.....#######...#########..... ........333.....1111111...444444444.....
.......#####.....#####....#########..... .......33333.....11111....444444444.....
......#######.....###.....#########..... ......3333333.....111.....444444444.....
.....#########.....#......#########..... .....333333333.....1......444444444.....
....###########...........#########..... ....33333333333...........444444444.....
........................................ ........................................
Figure 6.4 The state of the grid before and after the execution of the algorithm on the graph of
Clean Bandit.
Each cell containing a sharp is visited only once, so the algorithm is O(|V |), i.e. it
is linear in the number of vertices.
96 Graphs
Complexity
The complexity of this method is slightly worse than that of our first approach;
nevertheless, it becomes necessary when we have a graph whose number of edges
evolves in time and we wish to know at each moment the number of connected
components.
def nb_connected_components(graph):
n = len(graph)
uf = UnionFind(n)
nb_components = n
for node in range(n):
for neighbor in graph[node]:
if uf.union(node, neighbor):
nb_components -= 1
return nb_components
Problems
The Ant [spoj:ANTTT]
Lego [spoj:LEGO]
6.6 Biconnected Components 97
Input
Output
Application
Given a graph with costs of sabotage for each vertex and each edge, we wish to
determine a unique vertex or a unique edge of minimal cost which disconnects
the graph. Note the difference with the problem of finding a minimal set of edges
which disconnect the graph, the minimum cut problem, described in Section 9.8 on
page 162.
Definition
Let G be an undirected connected graph.
Figure 6.5 An edge between two cut vertices (vertices shown in bold) is not always a cut edge
and the endpoints of a cut edge (edges shown in bold) are not always cut vertices.
A biconnected component S also has the property that for every pair of vertices
s,t ∈ S two paths from s to t exist, which pass by distinct vertices. Note that the
biconnected components are defined by a partition of the edges and not by a partition
of the vertices. Indeed, a vertex can belong to several biconnected components, see
vertex 5 in Figure 6.8.
Given an undirected graph, the goal is to decompose it into biconnected
components.
Complexity
Linear by a depth-first traversal (see Hopcroft and Tarjan, 1973).
• a tree arc if v is seen for the first time during the processing of u. These arcs form
a spanning tree constructed by the traversal, also known as the DFS tree;
• an inverse tree arc if (v,u) is a tree arc;
• a back arc if v has already been seen and is an ancestor of u in the DFS tree;
• a forward arc if v has already been seen and is a descendant of u in the DFS tree.
For a depth-first traversal of a directed graph, an additional type of arc exists, a
cross arc, which goes to a vertex already seen, but is neither ancestor nor descen-
dant. As in this section we only consider non-directed graphs, we can ignore this
type of arc.
6.6 Biconnected Components 99
num
/min
u C tree arc
4/4
back arc
B
forward arc
F
B v C 9/5 cross arc
5/4 C
F A
B
6/4
8/5
B 7/6
1 2 5 7 10 11
4 3 6 8 9
Figure 6.7 The vertices are labelled by the number in their order of processing, stored in
dfs_num. The line style is solid for tree arcs, in dashes for back arcs and in dots for forward
arcs. For every tree arc (u,v) an inverse tree arc exists (v,u), not shown for readability.
100 Graphs
Key Observation
The values of dfs_low allow us to determine the cut vertices and edges.
1. A vertex u at the root of the DFS tree is a cut vertex if and only if it has at least
two children in the tree.
Note that each child v satisfies dfs_low[v] ≥ dfs_num[u].
2. A vertex u which is not the root is a cut vertex if and only if it has at least one
child v with dfs_low[v] ≥ dfs_num[u].
3. An edge (u,v)—up to a swap of u and v—is a cut edge if and only if (u,v) is a
tree arc and dfs_low[u] ≥ dfs_num[v].
Figure 6.8 The vertices are labelled with dfs_num followed by dfs_low. The cut vertices and
edge are highlighted in bold.
To determine the biconnected components, it thus suffices to apply the above defini-
tions in order to detect when a new biconnected component starts.
Implementation Details
The array parent contains the predecessor in the DFS tree of each vertex and also
allows the identification of the roots of the DFS trees. For each vertex u, we count in
critical_children[u] the number of children v in the tree such that dfs_low[v] ≥
dfs_num[u]. This information allows us at the end of the traversal to determine the
cut vertices and edges. The value dfs_low[v] is updated for every back arc detected
emanating from v. At the end of the process, this value is propagated up to the parent
vertex in the DFS tree.
6.6 Biconnected Components 101
Problem
Police Query [spoj:POLQUERY]
102 Graphs
5 3 2
4 3 1 5 4 2
Input Output
Definition
Given a directed graph G(V ,A), we want to order the vertices according to a rank
function r in such a way that for every arc (u,v), we have r(u) < r(v).
Application
The graph could represent a set of tasks where the arcs encode the dependencies
u → v ⇔ ‘u must be executed before v’. We are interested in an execution order
compatible with the dependencies.
First of all, a few remarks:
• Several topological sorts for a given graph can exist. For example, if the sequence
s is a topological sort of G1 and the sequence t a topological sort of G2 , then, for
example, the sequences st and ts are topological sorts of the graph formed by the
union of G1 and G2 .
• A graph containing a cycle never admits a topological sort: each vertex of a cycle
would require the prior processing of all the others before it (including itself!).
• A graph without cycles admits a topological sort. See below.
Complexity
The execution time is linear in the size of the input.
• either v was seen before u, in which case v was already handled (if not, this would
mean that u is accessible from v, hence the graph contains a cycle) and the inverse
topological order is respected;
• or v is reached from u, in which case u will be handled after v, hence the inverse
topological order is again respected.
The DFS traversal described previously must be adapted, as it does not allow
determination of when the processing of a vertex is complete.
6.7 Topological Sort 103
The following implementation uses an array seen, which is set to −1 for all vertices
that have not yet been seen, and then to the number of direct descendants that have
already been seen. When this counter is equal to the number of children of the node,
then the processing of the node is finished, and the node is added to the sequence
order. This sequence contains an inverse topological order, which must be reversed at
the end of the algorithm.
def topological_order_dfs(graph):
n = len(graph)
order = []
times_seen = [-1] * n
for start in range(n):
if times_seen[start] == -1:
times_seen[start] = 0
to_visit = [start]
while to_visit:
node = to_visit[-1]
children = graph[node]
if times_seen[node] == len(children):
to_visit.pop()
order.append(node)
else:
child = children[times_seen[node]]
times_seen[node] += 1
if times_seen[child] == -1:
times_seen[child] = 0
to_visit.append(child)
return order[::-1]
Greedy Algorithm
An alternative solution is based on the input degree of the vertices. Consider an acyclic
graph. Intuitively, we can begin by adding to the solution sequence the nodes without
predecessors in an arbitrary order, trimming them from the graph, then adding to
the sequence the new nodes without predecessors in an arbitrary order, and so on.
This process terminates because an acyclic graph always contains a vertex without
predecessors, and the deletion of nodes preserves the absence of cycles.
def topological_order(graph):
V = range(len(graph))
indeg = [0 for _ in V]
for node in V: # compute indegree
for neighbor in graph[node]:
indeg[neighbor] += 1
Q = [node for node in V if indeg[node] == 0]
order = []
while Q:
node = Q.pop() # node without incoming arrows
order.append(node)
for neighbor in graph[node]:
indeg[neighbor] -= 1
if indeg[neighbor] == 0:
Q.append(neighbor)
return order
104 Graphs
Applications
Given an acyclic graph and two vertices s,t, we might like to count the number of
paths from s to t, or discover the longest path when there are weights on the arcs.
A solution in linear time consists of performing a topological sort and executing a
dynamic program on the nodes in that order.
For example, the dynamic program P [s] = 0,P [v] = 1 + maxu P [u] computes
the longest path from s to t where the maximisation is made over all the arcs (u,v)
terminating on v.
Problems
Topological Sorting [spoj:TOPOSORT]
Project File Dependencies [spoj:PFDEP]
counters become 0) are added to the set of their colour. The algorithm is in linear time,
as the work carried out on each arc is in constant time.
Definition
A subset A ⊂ V of a directed graph is called a strongly connected component if for
every pair (u,v) of vertices of A a path in A exists that links u to v. Note that in this
case, a path in A exists that links v to u, see Figure 6.9.
Figure 6.9 An example of the partition of a graph into strongly connected components.
Observation
The component (or reduced) graph is defined as the result of the contraction of the
strongly connected components into ‘super vertices’. This graph is acyclic, as each
cycle of a directed graph is included in exactly one strongly connected component.
Complexity
The following algorithms execute in linear time in the size of the graph (|V | + |E|).
Tarjan’s Algorithm
Tarjan’s algorithm (1972) consists of a single depth-first search, numbering the ver-
tices in the order of the beginning of their processing. It also places the vertices
encountered into a stack waiting of nodes waiting to be grouped in a strongly con-
nected component. Once a component is detected, it is trimmed from the graph and
removed from the stack waiting. This means that all the arcs entering into an already
discovered component are ignored.
A depth-first search enters into a component C by a first vertex v0 , and then explores
all the vertices and eventually also the components reachable from C. Vertex v0 is
processed after all its descendants in the DFS tree. By construction, at the moment of
the final processing of v0 , clearly the stack waiting contains above v0 all the vertices
106 Graphs
Note that this value differs from the value dfs_low[v] described in Section 6.6 on
page 97, where there was no notion of waiting vertices and where dfs_low[v] could
take on the value ∞.
Consider the case of a component C without any outgoing arcs, and a vertex v ∈ C.
Let A be the DFS subtree with root v, and suppose that a back arc exists that leaves A
towards a vertex u, predecessor of v. As C does not have any outgoing arcs, u must
belong to C and v is not the representative of C. In this case, we have dfs_min[v]
< dfs_num[v]. If a back arc leaving A does not exist, then it contains a strongly
connected component, and since it is included in C, this tree covers C, and v is its
representative. This situation is detected by dfs_min[v] == dfs_num[v].
Implementation Details
In parallel with the stack waiting, the algorithm maintains a Boolean array waits,
which allows the presence in the stack to be tested in constant time. This way, it is
easy to trim a component by setting the corresponding elements to False.
The algorithm returns a list of lists, each containing the vertices of one of
the components. Note that the components are determined in reverse topological
order. This will be useful when we come to the resolution of a 2-SAT formula
later on.
6.8 Strongly Connected Components 107
def tarjan_recursif(graph):
global sccp, waiting, dfs_time, dfs_num
sccp = []
waiting = []
waits = [False] * len(graph)
dfs_time = 0
dfs_num = [None] * len(graph)
def dfs(node):
global sccp, waiting, dfs_time, dfs_num
waiting.append(node) # new node is waiting
waits[node] = True
dfs_num[node] = dfs_time # mark visit
dfs_time += 1
dfs_min = dfs_num[node] # compute dfs_min
for neighbor in graph[node]:
if dfs_num[neighbor] is None:
dfs_min = min(dfs_min, dfs(neighbor))
elif waits[neighbor] and dfs_min > dfs_num[neighbor]:
dfs_min = dfs_num[neighbor]
if dfs_min == dfs_num[node]: # representative of a component
sccp.append([]) # make a component
while True: # add waiting nodes
u = waiting.pop()
waits[u] = False
sccp[-1].append(u)
if u == node: # until representative
break
return dfs_min
Iterative Version
The iterative version of the algorithm is required for very large graphs, larger than,
say, 100,000 vertices. Here, a counter times_seen allows us at the same time to
mark the vertices encountered and to memorise the number of neighbours already
considered.
108 Graphs
def tarjan(graph):
n = len(graph)
dfs_num = [None] * n
dfs_min = [n] * n
waiting = []
waits = [False] * n # invariant: waits[v] iff v in waiting
sccp = [] # list of detected components
dfs_time = 0
times_seen = [-1] * n
for start in range(n):
if times_seen[start] == -1: # initiate path
times_seen[start] = 0
to_visit = [start]
while to_visit:
node = to_visit[-1] # top of stack
if times_seen[node] == 0: # start process
dfs_num[node] = dfs_time
dfs_min[node] = dfs_time
dfs_time += 1
waiting.append(node)
waits[node] = True
children = graph[node]
if times_seen[node] == len(children): # end of process
to_visit.pop() # remove from stack
dfs_min[node] = dfs_num[node] # compute dfs_min
for child in children:
if waits[child] and dfs_min[child] < dfs_min[node]:
dfs_min[node] = dfs_min[child]
if dfs_min[node] == dfs_num[node]: # representative
component = [] # make component
while True: # add nodes
u = waiting.pop()
waits[u] = False
component.append(u)
if u == node: # until repr.
break
sccp.append(component)
else:
child = children[times_seen[node]]
times_seen[node] += 1
if times_seen[child] == -1: # not visited yet
times_seen[child] = 0
to_visit.append(child)
return sccp
Kosaraju’s Algorithm
A different algorithm was proposed by Kosaraju (1979). Its complexity is also linear,
and in practice quite comparable to Tarjan’s algorithm; however, it is perhaps easier
to understand.
The principle consists in performing a first depth-first search, and then a sec-
ond on the graph with the orientation of the arcs reversed. Formally, denote AT :=
{(v,u)|(u,v) ∈ A} the result of the reversal of the arcs. The algorithm is then composed
of the following steps:
6.8 Strongly Connected Components 109
1. Perform a depth-first search of G(V ,A) and store as f [v] the time of the
completion of processing of the vertex v.
2. Perform a depth-first search of G(V ,AT ), selecting roots v in order of decreasing
f [v].
Each tree encountered in the second traversal is a strongly connected component.
The validity of the algorithm is based on the idea that if we associate with each
strongly connected component C the integer F (C) := maxu∈C fu , then F gives
a topological order on the graph induced by the strongly connected components in
G(V ,AT ). Hence, each tree generated by the second traversal remains in the interior
of a component, since the only outgoing arcs lead to components already traversed.
Implementation Details
The list sccp (for strongly connected component) contains the list of strongly con-
nected components.
def reverse(graph):
rev_graph = [[] for node in graph]
for node, _ in enumerate(graph):
for neighbor in graph[node]:
rev_graph[neighbor].append(node)
return rev_graph
def kosaraju(graph):
n = len(graph)
order = []
sccp = []
kosaraju_dfs(graph, range(n), order, [])
kosaraju_dfs(reverse(graph), order[::-1], [], sccp)
return sccp[::-1] # follow inverse topological order
110 Graphs
Problem
Capital City [spoj:CAPCITY]
6.9 2-Satisfiability
Definition
We are given n Boolean variables. A literal is a variable or the negation of a variable.
A clause is a disjunction (OR) of several literals, so that a clause is satisfied (i.e. , made
true) if at least one of its literals is true. A formula is a conjunction (AND) of several
clauses, so that a formula is satisfied if all its clauses are satisfied. The goal is to see
whether there exists an assignment to the variables that satisfies the formula.
A formula is said to be of class 2-SAT if each clause contains at most two literals.
One of the fundamental results in complexity theory is that we can verify in linear
time if a 2-SAT formula is satisfiable, whereas in general (for a 3-SAT formula, for
example), in the worst case, a polynomial time algorithm is not known (it is NP-hard).
a b c ⎧
⎪
⎪a∨b
⎨
b∨c
⎪
⎪ a∨c
⎩
a ∨ c.
a b c
Figure 6.10 Graph obtained from an instance of 2-SAT. There are two strongly connected
components. The bottom one points to the top one, hence assigning ‘False’ to all the literals on
the bottom and ‘True’ to all on top will satisfy this formula.
Complexity
The complexity is linear, i.e. time O(n + m) for n and m the number of variables and
clauses (see Aspvall et al., 1979).
Key Observation
We can easily show that a 2-SAT formula is not satisfiable if for some x there is a path
from x to x and a path from x to x in the implication graph. What is surprising is that
the converse is also true!
As the implications are transitive, clearly if x ⇒ x ⇒ x is implied by the formula,
then it cannot be satisfied, since each of the assignments x = False or x = True
leads to a contradiction False ⇒ True. Conversely, suppose that each variable is in
a different strongly connected component from its negation.
As the directed graph is symmetric, we can consider pairs of strongly connected
components, where each contains all the negations of literals of the other. More-
over, all the literals of a strongly connected component must have the same Boolean
value. From this graph, it is thus possible to determine an assignment satisfying the
formula: it suffices to assign True to all the literals of a component without any
outgoing arcs. There is forcibly one, as the component graph is acyclic. Then we
assign False to the opposing component, trim the two components from the graph and
recurse.
To perform this efficiently, note that the algorithm presented previously to deter-
mine the strongly connected components finds them in reverse topological order. It
thus suffices to traverse the components in this order, assign True to each component
that does not yet have a value and assign False to the opposite component.
Implementation Details
We encode the literals by integers, such that +1, . . . , + n encode the n variables
and −1, . . . , − n encode their negations. A clause is encoded by a pair of inte-
gers and a formula as a list of clauses. Each literal corresponds to one of the 2n
nodes of the associated directed graph. These nodes are numbered from 0 to 2n − 1,
where 2i encodes the variable xi+1 and 2i + 1 encodes xi+1 . Here is the correspond-
ing code.
112 Graphs
def two_sat(formula):
# num_variables is the number of variables
num_variables = max(abs(clause[p])
for p in (0, 1) for clause in formula)
graph = [[] for node in range(2 * num_variables)]
for x_idx, y_idx in formula: # x_idx or y_idx
graph[_vertex(-x_idx)].append(_vertex(y_idx)) # -x_idx => y_idx
graph[_vertex(-y_idx)].append(_vertex(x_idx)) # -y_idx => x_idx
sccp = tarjan(graph)
comp_id = [None] * (2 * num_variables) # each node’s component ID
assignment = [None] * (2 * num_variables)
for component in sccp:
rep = min(component) # representative of the component
for vtx in component:
comp_id[vtx] = rep
if assignment[vtx] is None:
assignment[vtx] = True
assignment[vtx ^ 1] = False # complementary literal
for i in range(num_variables):
if comp_id[2 * i] == comp_id[2 * i + 1]:
return None # insatisfiable formula
return assignment[::2]
Problems
Manhattan [kattis:manhattan]
Soldiers on Parade [spoj:SOPARADE]
7 Cycles in Graphs
Several classic problems concern cycles in graphs, whether they refer to geographical
displacements or to anomalies in a dependency graph. The simplest problems consist
of detecting the existence of cycles, the existence of cycles with negative weight or
the identification of a minimal total weight or minimal mean weight cycle.
Other problems are concerned with exploring the whole of a graph to find paths
traversing each edge exactly once (Eulerian path) or, when this is not possible, at least
once (Chinese postman problem). These problems are polynomial, whereas determin-
ing a cycle that visits each vertex exactly once (Hamiltonian cycle) is NP-hard.
8
4 7
1 5
3 6
2
Input Output
Application
You are in Kaliningrad—formerly Königsberg—and want to know if it is possible to
take a walk that crosses each bridge exactly once and finishes back at the starting
point. It is this situation that led Leonhard Euler to study the problem in 1759.
114 Cycles in Graphs
Definition
Given a connected graph G(V ,E) such that each vertex has even degree, we wish to
find a cycle in this graph that traverses each edge exactly once. The same problem can
be posed for directed and strongly connected graphs such that the out-degree of each
vertex is equal to its in-degree.
0 1 2 3 4
Figure 7.1 Hierholzer’s algorithm running on an example. The first cycle detected is
0 → 1 → 2 → 0. Then along this cycle the vertices are explored to discover further attached
cycles. The vertices 0 and 1 have been processed and are pushed onto the stack P = [0,1].
It remains to process the remaining vertices of this cycle, namely Q = [0,2]. The algorithm
then explores a cycle starting from vertex 2 (start), which was popped from the stack Q, now
reduced to [0]. The exploration discovers the path 2 → 3 → 4 whose vertices are pushed onto
the stack R = [3,4]. The extreme point of the path is node = 4. The final path is
0 → 1 → 2 → 3 → 4 → 2 → 0.
For the algorithm to be in linear time, the search for the vertex v must be efficient.
For this, we decompose the cycle C into two portions P and Q. The vertices in P are
those that have no more allowable edges; they will be part of the final tour. As long as
Q is non-empty, we remove the vertex v from the head of Q and add it to the end of P ,
hence advancing the cycle. Next, we attempt to insert here a cycle passing through v.
For this, if there is an allowable edge incident to v, we find by simple exploration a
cycle R starting from v and returning to v. The cycle R is then added to the head of Q.
The algorithm is in linear time, as each vertex is considered only once.
Implementation Details
We first describe the implementation for directed graphs. To ease the manipulation of
the data, P ,Q and R are represented by stacks, encoded by lists. The current cycle
7.1 Eulerian Tour 115
consists of the stack P , the vertex start, followed by the mirror list Q. The list R
describes a path, which once completed as a cycle will be attached to the current cycle
at the vertex start. To rapidly find an allowable outgoing arc of a vertex, we store
in next[node] the number of outgoing arcs of a vertex node that have already been
explored. It suffices to increment this counter when we traverse the arc towards the ith
neighbour of node, where i=next[node].
def eulerian_tour_directed(graph):
P = [] # resulting tour
Q = [0] # vertices to be explored, start at 0
R = [] # path from start node
next_ = [0] * len(graph) # initialize next_ to 0 for each node
while Q:
start = Q.pop() # explore a cycle from start node
node = start # current node on cycle
while next_[node] < len(graph[node]): # visit all allowable arcs
neighbor = graph[node][next_[node]] # traverse an arc
next_[node] += 1 # mark arc traversed
R.append(neighbor) # append to path from start
node = neighbor # move on
while R:
Q.append(R.pop()) # add to Q the discovered cycle R
P.append(start) # resulting path P is extended
return P
The variant of this algorithm for non-directed graphs is a bit more subtle. Once an
arc (u,v) is traversed, we must mark the arc (v,u) as forbidden. We have chosen to
store in seen[v] the set of neighbours u of v such that (v,u) is no longer allowed. For
more efficiency, v does not need to be added to seen[u] at the time of traversal of the
arc (u,v), since at this same moment the counter next[u] is incremented and this arc
will no longer be considered.
def eulerian_tour_undirected(graph):
P = [] # resulting tour
Q = [0] # vertices to be explored, start at 0
R = [] # path from start node
next_ = [0] * len(graph) # initialize next_ to 0 for each node
seen = [set() for _ in graph] # mark backward arcs
while Q:
start = Q.pop() # explore a cycle from start node
node = start # current node on cycle
while next_[node] < len(graph[node]): # visit all allowable arcs
neighbor = graph[node][next_[node]] # traverse an arc
next_[node] += 1 # mark arc traversed
if neighbor not in seen[node]: # not yet traversed
seen[neighbor].add(node) # mark backward arc
R.append(neighbor) # append to path from start
node = neighbor # move on
while R:
Q.append(R.pop()) # add to Q the discovered cycle R
P.append(start) # resulting path P is extended
return P
116 Cycles in Graphs
Problem
Goldbach graphs [spoj:GOLDG]
Application
In 1962, the mathematician Meigu Guan (or Kwan Mei-Ko) worked as a postman
during the Chinese Cultural Revolution. He then posed the problem of finding the
shortest cycle in a graph that visits each edge at least once. This is exactly the problem
that a postman needs to solve to deliver the mail in each neighbourhood of a city.
Definition
Suppose we are given a connected undirected graph G(V ,E). The goal is to find a
cycle in the graph that visits each edge at least once. When all the vertices are of even
degree, this boils down to determining an Eulerian cycle, the problem treated in the
preceding section. Note that in the original formulation of this problem, the graph was
weighted.
Complexity
The algorithm executes in time O(n3 ), where n is the number of nodes of the graph
(see Edmonds and Johnson, 1973).
Algorithm
The principle of the algorithm consists of working with a multigraph, i.e. a graph that
may have several edges between the same two vertices. The idea is to add edges in
order to render the graph Eulerian, and then to produce a Eulerian cycle.
The edges added must help link vertices of odd degree in order to make the graph
Eulerian. As we wish to add as few edges as possible, these extra edges will form a
collection of paths that will match edges of odd degree.
The central problem thus consists of computing a perfect matching in the complete
graph of vertices in V of odd degree, where the weight of an edge (u,v) is equal to the
distance between u and v in G.
The computation of all the distances costs O(n3 ) by the Floyd–Warshall algo-
rithm. Then a perfect minimum weight matching can be computed in time O(n3 ) with
Gabow’s algorithm (1976). However, this sophisticated algorithm is beyond the scope
of the present text.
7.3 Cycles with Minimal Ratio of Weight to Length—Karp 117
Problem
Free Tour [spoj:FTOUR]
Definition
Given a weighted directed graph, the goal is to find a cycle C that minimises the mean
weight over the edges, i.e. that minimises the following
e∈C w(e)
.
|C|
Application
Imagine a system modelled by a graph whose vertices represent states and whose arcs
are the possible transitions between states. Each transition is labelled with its resource
consumption. At each step in time, the system is in a particular state and must evolve
to another state by following an outgoing arc. The goal is to minimise the long-term
resource consumption, and an optimal solution is a minimal mean weight cycle.
where the minimum is taken over all incoming arcs (u,v) of weight wuv . Let d[v] =
mink∈N d[k][v] be the distance from the source s to v.
Key Observation
Let C be a cycle of minimal mean weight. Its weight λ satisfies the following
d[n][v] − d[k][v]
λ = min max . (7.1)
v k=0,...,n−1 n−k
118 Cycles in Graphs
To see this, it suffices to make a few observations, beginning with the case λ = 0.
We must show that the right-hand side of (7.1) is zero, or equivalently,
min max d[n][v] − d[k][v] = 0.
v k=0,...,n−1
Since λ = 0, the graph does not have a negative weight cycle. For each vertex v, a
shortest path exists from the source s to v which is acyclic, and hence
d[v] = min d[k][v].
k=0,...,n−1
Consequently,
max d[n][v] − d[k][v] = d[n][v] − d[v].
k=0,...,n−1
It remains to exhibit a vertex v for which the equality d[n][v] = d[v] holds. Let u be
a vertex of C. As there is no negative weight cycle, a simple path P exists that is of
minimal weight from the source s to u. We complete P by as many copies of C as
required to obtain a path P from s to u with length at least n. As C has weight zero,
P is again a minimal weight path to u. Let P be a prefix of P of length exactly n
and let v be the last vertex of P . P is also a shortest path from the source to v, see
Figure 7.2. Hence d[n][v] = d[v] for this vertex, which proves the equality (7.1) for
the case λ = 0.
s P u C
Figure 7.2 Illustration of the proof for the case λ = 0. P is a minimal weight path from the
source to u, and adding copies of the cycle C with zero weight results in a minimal weight
path; hence every vertex encountered on the cycle is itself reached by a minimal weight path.
Implementation Details
The matrix dist contains the distances as described above. There is also a variable
prec which stores the predecessors on the shortest paths. Once these matrices are
populated, the second step consists of finding the pair v,k optimising the expression
(7.1), and a third step extracts a cycle. In the case where there are no cycles reachable
from the source, the function returns None.
Problem
Word Rings [spoj:WORDRING]
120 Cycles in Graphs
Definition
Let G be a directed graph, with two weights on its arcs—costs c and times t. The times
are non-negative, while the costs are arbitrary. The goal is to find a cycle minimising
the ratio between the total cost and the total time. This is known as the tramp steamer
problem.
Application
The captain of a cargo vessel wants to determine a route that maximises his profit.
For this, he has available a complete graph covering different ports, where each (u,v)
is labelled by the travel time between u and v as well as the profit he can make by
purchasing merchandise at the port u and selling it at the port v. The goal is to find a
cycle that maximises the ratio of total profit to total time.
To detect a cycle C with a ratio strictly better than δ, it suffices to detect a negative
cycle in the graph where the arcs a are weighted by c(a) − δt(a). Hence, this test can
be used in a binary search on δ to obtain a solution with the required precision. In the
case of integer costs and times, a precision of 1/ a t(a) suffices to solve the problem
exactly. The complexity is then O(log( t(a))·|V |·|A|), where V is the set of vertices
and A the set of arcs.
As initially we do not have an a priori upper bound on the optimum δ∗ , we begin
the test with δ = 1. As long as the test fails, we multiply δ by 2, to finally obtain a
δ where the test is positive, thus providing upper and lower bounds on the optimum.
In the case where the initial test is positive, we proceed in the same way, but instead
divide by 2 until the test fails.
Definition
Given a graph with weights on its arcs, we would like to compute a shortest path
beginning at a given source and visiting each vertex exactly once: such a path is said
to be Hamiltonian.
7.6 Full Example: Menu Tour 121
Complexity
The solution presented, using dynamic programming, executes in time O(|V |2 2|V | ).
Algorithm
As the associated decision problem is NP-complete, we propose an acceptable solution
for when the number of vertices is small, on the order of twenty or so. Suppose that
the vertices are numbered from 0 to n − 1, where for ease of notation we take n − 1
as the source.
For each set S ⊆ {0,1, . . . ,n − 2}, assign O[S][v] as the minimal weight of a path
from the source to v passing exactly once by all the vertices of S and terminating on v,
with v ∈ S.
For the base case, O[∅][v] is simply the length of the arc from n−1 to v. Otherwise,
for non-empty S, O[S][v] is the minimum of
O[S\{u}][u] + wuv
over all the vertices u ∈ S, where wuv is the weight of the arc (u,v).
Problems
Small TSP [spoj:STSP]
Collecting Beepers [kattis:beepers]
See [spoj:MENUTOUR].
Rachid is in Barcelona for the first time, and wants to enjoy a truly fine dinner. The
restaurant guides suggest that the ultimate dinner consists of C specialities, numbered
from 1 to C, tasted in the specified order.
The streets of Barcelona form a regular grid: east–west oriented streets cross-
ing north–south oriented streets. There are R restaurants in the city and they are
located at the crossings. Walking from crossing (i1,j1 ) to crossing (i2,j2 ) takes exactly
|i1 − i2 | + |j1 − j2 | minutes, where |x| denotes the absolute value of x. Here, (i,j )
means the ith street from west to east and the j th street from south to north.
Unfortunately, none of the restaurants offer all C dishes. If a restaurant k offers
dish c, then it costs P [k,c] euros. Otherwise, the value P [k,c] = 0 indicates that
the item is not offered. Rachid has B euros that he can spend on his dinner. He
would like to choose a sequence of restaurants so that he can have his ultimate dinner
without exceeding the available budget, while minimising the total travel time between
restaurants. The tour can start and end at an arbitrary crossing, and can visit the
same restaurant several times.
122 Cycles in Graphs
Example
4 5
-2- 8-9
3
--3
2
1
1-- -97
In this example, there are three menu items, whose prices are displayed in order, for
every restaurant (a ‘-’ corresponds to a dish that is not offered, i.e. when P [k,c] = 0).
If Rachid has a budget of 9 euros, the optimal tour consists of the restaurants 1–4–3
for a cost of 6 euros and a total travel time of 12 minutes. There is a shorter tour of 2
minutes consisting of the restaurant sequence 1–2, but its cost is 17 euros, exceeding
the available budget of 9 euros.
Input
The input begins with a line consisting of the integers C,R,B, separated by a single
space. Then R lines follow. The kth line describes the kth restaurant and consists of
2 + C integers separated by a single space, namely i[k], j [k], P [k,1], . . . ,P [k,C],
where (i[k],j [k]) defines the location of the restaurant.
Limits
• 1 ≤ C ≤ 20;
• 1 ≤ R ≤ 100;
• 0 ≤ B ≤ 100;
• for every 1 ≤ k ≤ R:
– 1 ≤ i[k] ≤ 1 000;
– 1 ≤ j [k] ≤ 1 000;
– for every 1 ≤ c ≤ C, 0 ≤ P [k,c] ≤ 40.
7.6 Full Example: Menu Tour 123
Output
The output should consist of a single integer y: the minimum total travel time of an
optimal restaurant tour for Rachid. If there is no feasible tour, the value −1 should be
returned.
359 12
11100
31097
62003
35020
65809
Solution
This problem sounds very much like a travelling salesman tour, but the order in which
the dishes need to be consumed simplifies it a lot. It can be computed by dynamic
programming. For any dish number c, restaurant r and budget b, let A[c,r,b] be the
minimum time for a tour that covers dishes from 0 to c, ends in restaurant r and costs
exactly b euros. By convention denote A[c,r,b] = +∞ if there is no such tour.
The base cases are for c = 0. For any restaurant r, if P [r,0] > 0 the tour could
start in r and spend P [r,0] on the first dish (numbered 0), without any travel time.
Hence, A[0,r,b] = 0 for b = P [r,0] in that case, and A[0,r,b] = +∞ otherwise.
The induction step for c > 0 and any r,b follows the following recursion
+∞ if b < P [r,c]
A[c,r,b] =
mini A[c − 1,i,b − P [r,c]] + dist[i][r] otherwise,
where the dist[i][r] denotes the distance between the restaurants i and r. These
distances should be precomputed in order for the algorithm to be efficient. There are
O(NCB) variables, each is the minimisation over O(N ) alternatives, hence the total
complexity is O(N 2 CB), which is acceptable given the upper bounds.
8 Shortest Paths
A classic problem on graphs consists of finding a shortest path between two vertices,
a source s and a destination v. For the same cost we can obtain the shortest paths
between the source s and all possible destinations v ; this is why the algorithms pre-
sented in this chapter solve this more general problem of the shortest paths from a
unique source in a directed graph.
The length of a path is defined as the sum of its arc weights, and the distance from
s to v is defined as the length of a shortest path between s and v. To simplify the
presentation, in general we show only how to compute the distances. To obtain a path
realising this distance, it suffices to maintain, in addition to the array of distances, an
array of predecessors. Hence, if dist[v] is the distance from s to v with dist[v] =
dist[u] + w[u][v] for a vertex u, we store prec[v] = u in the array of predecessors.
By following the predecessors back to the source s, we can determine a shortest path
in reverse order, leading from the source to a given destination.
The shortest paths exhibit a composition property, which is the key to the different
shortest path algorithms. Richard Bellman called this the Principal of Optimality,
which lies at the heart of problems of dynamic programming. Consider a path P from
s to v (also known as an s − v path), passing by a vertex u, see Figure 8.1. It is thus the
concatenation of a path P1 from s to u with a path P2 from u to v. The length of P is
the sum of the lengths of P1 and P2 . Thus, if P is a shortest path from s to v, then P1
must be a shortest path from s to u and P2 a shortest path from u to v. This is obvious,
as the replacement of P1 by a shorter path would result in a shorter path from s to v.
P1 P2
s u v
Figure 8.1 A shortest path from s to v passing through u is the concatenation of a shortest path
from s to u with a shortest path from u to v.
8.1 Composition Property 125
1
5 ∞
3 2 3
1 8 7 2
1 2 6 ∞ ∞
1 3 2 3
4 2 2 1
s 0 4 6 ∞ ∞
Figure 8.2 Example of the colouring of vertices by Dijkstra’s algorithm. Each vertex v is
labelled with dist[v], its current distance from the source, and the arcs specified by prec are
highlighted in bold.
Key Observation
Which grey vertex should we select to colour black? The choice is made with a vertex
v with the smallest dist[v]. Why is this the correct choice? Consider an arbitrary s −v
path P . Since s is black and v is grey, a first grey vertex u exists on this path. Hence, P
can be decomposed into an s − u path P1 and a u − v path P2 , where P1 contains only
intermediate black vertices, excluding u. By the choice of v, dist[u] ≥ dist[v],
and by the hypothesis that the weights of the arcs are non-negative, the length of P2
is also non-negative. Hence, the length of P is at least dist[v]. As the choice of P
was arbitrary, this shows that the shortest s − v path is of length dist[v]. Colouring
v black is thus valid, and to maintain the invariant, for each vertex v reachable from
v by an arc, we must colour v in grey if this is not already the case, and dist[v] +
w[v][v’] becomes a new candidate for the distance dist[v’].
of dist. This is the implementation chosen for Dijkstra’s algorithm. The removal
of the vertex with the smallest distance can thus be done in time logarithmic in the
number of vertices.
If the graph has a special structure, then a simpler structure could possibly be used.
For example, in the case where the weights of the arcs are in {0,1}, there will only be
two types of grey vertices, those with a distance d and those with a distance d +1. This
priority queue can then be implemented in a simpler manner using a double-ended
queue, abbreviated deque. This queue contains the list of grey vertices in priority order,
and it suffices to remove from the left in constant time a vertex v to be coloured black.
The neighbours of v are added to the queue on the left or on the right, according
to whether the weight on the corresponding arc is 0 or 1, see 8.2. Finally, all the
operations of manipulation of the queue can be done in constant time, and we shave
off a logarithmic factor compared to Dijkstra’s algorithm.
If the graph is even simpler, with unit weights on all of its arcs—in a manner of
speaking a non-weighted graph—then this double-ended queue can be replaced by
a simple queue, and managed using the algorithm of breadth-first search described
previously in Section 6.4 on page 93.
Definition
Given a graph whose arcs are weighted by 0 or 1, and a source vertex s, we wish to
compute the distances from s towards the other vertices.
Application
You are provided with a map of an N ×M rectangular labyrinth containing obstructing
walls and would like to exit while demolishing a minimum number of walls. The
labyrinth can be seen as a directed graph whose arcs (from one cell to an adjacent
cell) are weighted either by 0 (free access to the cell) or by 1 (access to a cell blocked
by a wall). We seek to demolish as few walls as possible, which comes down to finding
the shortest path in this graph from the entrance to the exit.
Algorithm
We return to the structure of a generic shortest path algorithm. At any moment, the
vertices of the graph are divided into three groups: black, grey and white.
8.3 Graphs with Non-negative Weights—Dijkstra 127
The algorithm maintains a double-ended queue containing all the grey vertices, as
well as eventually some vertices that were grey at the time of their insertion, but that
since have become black. The queue has the property that for a certain value x, all
the black vertices vb satisfy dist[vb ] = x, and up to a certain position all the grey
vertices vg0 satisfy dist[vg0 ] = x, whereas those after satisfy dist[vg1 ] = x + 1.
As long as this queue is non-empty, the algorithm extracts a vertex v from the
head of the queue, which hence has a minimal value dist[v]. If this vertex is
already black, there is nothing to do. Otherwise, it is coloured black. In order to
maintain the invariant, we must then add certain neighbours v of v to the queue. Let
= dist[v] + w[v][v ]. If v is already black, or if dist[v’] ≤ , there is no reason
to add v to the queue. Otherwise, v can be coloured grey, dist[v’] is lowered to
and v is added to the head or tail of the queue depending on whether w[v][v ] = 0 or
w[v][v ] = 1.
Problem
KATHTHI [spoj:KATHTHI]
Definition
Given a directed graph with non-negative weights on its arcs, we seek the shortest path
between a source and a destination.
128 Shortest Paths
Complexity
A brutal implementation is in O(|V |2 ), while the use of a priority queue lowers the
complexity to O(|E| log |V |). By using a Fibonacci priority queue, we can obtain a
complexity of O(|E| + |V | log |V |) (see Fredman and Tarjan, 1987); however, the
improvement is perhaps not worth the implementation effort.
Algorithm
We continue with the formalism introduced in Section 8.1 on page 125. Dijkstra’s
algorithm maintains a set of vertices S for which a shortest path from the source has
already been computed. Hence, S is precisely the set of black vertices. Initially, S
contains only the source itself. In addition, the algorithm maintains a tree of shortest
paths rooted at the source and covering S. We assign prec[v] as the predecessor of v
in this tree and dist[v] as the distance computed, see Figure 8.3.
2 1
4
Figure 8.3 Dijkstra’s algorithm captures the vertices via balls of increasing radius centred on
the departure vertex. There are no vertices at distance 1, thus the algorithm leaps directly to the
closest vertex, at distance 2, via the priority queue.
We then examine the edges at the boundary of S. More precisely, we consider the
arcs (u,v) with u in S and v not in S. These arcs define paths composed of a shortest
path from the source to u followed by the arc (u,v). This path has a weight that we
assign as priority to the arc (u,v).
As explained in Section 8.1 on page 125, the algorithm extracts from the priority
queue an arc (u,v) of minimal priority, thus defining a shortest s − v path. We can then
add the arc (u,v) to the tree of shortest paths, add v to S and iterate.
Priority Queues
The principal ingredient of this algorithm is a priority queue, which is a data structure
containing a set of elements, allowing the addition of elements and the extraction of
the one with the smallest priority.
Such a structure is, in general, implemented with a heap (here a min-heap, see
Section 1.5.4 on page 22), where these operations have logarithmic cost in the size
of the set.
8.3 Graphs with Non-negative Weights—Dijkstra 129
A Slight Improvement
For more efficiency, we can ignore arcs already known to not lead to a shortest path,
and not store them in the priority queue. For this, we update dist[v] whenever an arc
to v is added to the queue. Hence, at the moment of considering an arc, we can decide
if it improves the shortest path to v known up to this point.
Implementation Details
The code presented allows the computation of the shortest paths to each and every
destination, if a particular destination is not specified when the function is invoked.
For each outgoing arc (u,v), we store in the queue the pair (d,v), where d is the
length of the path associated with this arc.
It is possible for a vertex to be found several times in the queue with different
weights (priorities). However, once the first occurrence of this vertex is extracted, it
will be coloured black, and the further occurrences will be ignored at the time of their
extraction.
Variant
The implementation of Dijkstra’s algorithm can be slightly simplified with the use of
a priority queue allowing the priority of an element to be changed, such as the queue
presented in Section 1.5.4 on page 22. Instead of storing the arcs in the queue, we store
the vertices, in fact all the vertices of the graph, and the priority of a vertex v is set
to dist[v]. In practice, the queue contains pairs of the form (dist[v], v), and we
compare them lexicographically. When a shorter path to v is discovered, we update the
pair corresponding to v with the new shorter distance. The advantage of this variant is
that we can skip marking the vertices black. For each vertex v, the queue contains only
130 Shortest Paths
a single pair with v, and when (dist[v], v) is extracted from the priority queue, we
know that we have discovered a shortest path towards v.
Problems
Mice and Maze [spoj:MICEMAZE]
Almost Shortest Path [spoj:SAMER08A]
Dragon Maze [gcj:China2014roundB]
Potholers [spoj:POTHOLE]
Definition
For this problem, negative arc weights are allowed. If a cycle of negative total weight
is reachable from the source, the distance from the source to the vertices of the cycle
is −∞. Indeed, by traversing this cycle an arbitrary number of times, we obtain a path
from the source to a vertex in the cycle of negative weight with negative length and
arbitrarily large magnitude. This bothersome situation is detected by the algorithm
presented.
Complexity
Dynamic programming provides a solution in O(|V | · |E|); see Bellman (1958); Ford
(1956) as well as Moore (1959), who independently discovered the same algorithm.
Algorithm
The central operation is the relaxation of the distances (see Figure 8.4), consisting
for each arc (u,v) in testing if the use of this arc could diminish the distance from
8.4 Graphs with Arbitrary Weights—Bellman–Ford 131
the source to v. The value du + wuv is thus a candidate for the distance dv . This is
performed with two nested loops; the inner loop relaxes the distances via each arc,
and the outer loop repeats this procedure a certain number of times. We show that
after k iterations of the outer loop, we have computed for each vertex v the length of a
shortest path from the source to v, using at most k arcs. This is true for k = 0, and if
we assign dk [v] as this length for every k = 1, . . . ,|V | − 1, then
dk+1 [v] = min dk [u] + wuv .
u:(u,v)∈E
9 6
9 9
2 2
0 4 0 4
4 4
Figure 8.4 An example of relaxation via the greyed arc. Following this edge shortens the path to
the destination vertex of the arc. The vertices are labelled with the distance from the source
vertex, which has label 0. The shortest paths are shown in bold.
Problem
Negative Cycles [spoj:NEGCYC]
132 Shortest Paths
Definition
Given a graph with weights on its arcs, we want to compute the shortest path between
each pair of nodes. Once again, the problem is well-posed only in the absence of cycles
with negative weight, and the algorithm detects this exception.
Complexity
The Floyd–Warshall algorithm executes in O(n3 ), where n is the number of vertices;
see Floyd (1962); Warshall (1962) as well as Roy (1959), who independently discov-
ered the same algorithm.
Algorithm
The distances between the vertices are computed by dynamic programming, see
Figure 8.6. For each k = 0,1, . . . ,n a square matrix Wk is computed, containing
in Wk [u][v] the length of the shortest path from u to v using only the intermediate
vertices of index strictly inferior to k, where the vertices are numbered from 0 to n−1.
Hence, for k = 0, the matrix W0 contains only the weights of the arcs, and +∞ for the
elements (u,v) not having an arc (u,v). The updating is simple using the composition
principle: a shortest path from u to v passing through a vertex k is composed of
a shortest path from u to k and a shortest path from k to v. We thus compute, for
k,u,v ∈ {0, . . . ,n − 1}:
Wk+1 [u][v] = min{Wk [u][v],Wk [u][k] + Wk [k][v]}.
It is indeed the same k as index of the vertices and as position in the array: for the
computation of Wk+1 [u][v], we are allowed to consider the paths passing through the
vertex k.
n−1
(+1) n−1 () ()
auv = buk × ckv Wuv = min Wuk + Wkv
k=0
k=0
Figure 8.5 The Floyd–Warshall algorithm can be seen as a multiplication of matrices on the
tropical algebra (R, min ,+). However, we cannot employ the usual algorithms of rapid
multiplication since we are only operating in a semiring.
2 2
4 k 4 k
v v
u u
9 9
Figure 8.6 An example of relaxation via the vertex k. Passing through k shortens the path
from u to v.
8.6 Grid 133
Implementation Details
The implementation maintains a unique matrix W , which plays the role of the succes-
sive Wk ; it modifies the given weight matrix. It detects the presence of negative cycles
and returns False in this case.
def floyd_warshall(weight):
V = range(len(weight))
for k in V:
for u in V:
for v in V:
weight[u][v] = min(weight[u][v],
weight[u][k] + weight[k][v])
for v in V:
if weight[v][v] < 0: # negative cycle found
return True
return False
8.6 Grid
Problem
Given a rectangular grid, where certain cells are inaccessible, find a shortest path
between the entrance and the exit.
134 Shortest Paths
The Algorithm
This problem is easily solved with a breadth-first search on the graph determined by
the grid. However, rather than explicitly constructing the graph, it is simpler to perform
the search directly on the grid. The grid is given in the form of a 2-dimensional array
containing the character # for the inaccessible cells and a space for the accessible cells.
Our implementation uses this same array to mark the vertices already visited, and
hence avoids an additional data structure. The cells visited will then contain symbols
‘→’, ‘←’, ‘↓’, ‘↑’, indicating the predecessor on the shortest path from the source.
Variants
The proposed implementation is easy to modify to permit diagonal steps on the grid,
between cells sharing only a corner. A hexagonal grid can be treated in a similar
manner by transforming it into a grid with square cells, but with a particular neigh-
bourhood. Figure 8.7 illustrates this transformation.
−→ −→
Problem
A Famous Grid [spoj:SPIRALGR]
8.7 Variants
Given the importance of the shortest path problem, we present a number of classic
variants.
Application
This problem is given as an exercise in Cormen et al. (2010, Ex. 24-6, p. 618). We
wish to determine a route from home to work that starts by ascending and finishes by
descending, in order to first get some exercise and then cool down. For this, we model
the city as a graph, where the vertices represent the intersections with their elevation,
and the edges represent the roads with their lengths. Using the previous algorithm, we
can compute for each vertex v the distance from the source using only ascending arcs,
and at the same time compute the distance from v to the target using only descending
arcs. The answer is the sum of these distances minimised over all vertices v.
where S \ u is the set S with the vertex u removed, and w[u][v] is the weight of the arc
(u,v). The base case is then
w[u][v] if the arc (u,v) exists
D[∅][v] =
−∞ otherwise.
Problems
Almost Shortest Path [spoj:SAMER08A]
Flowery Trails [kattis:flowerytrails]
9 Matchings and Flows
Traditionally, the heart of combinatorial optimisation lies in the matching and flow
problems. These problems are related and exist with several variants. The principle
underlying most of these algorithms is the repeated improvement of a solution, ini-
tially empty, leading to an optimal solution.
Suppose that in the bipartite graph of Figure 9.1 we wish to determine a perfect
matching, i.e. a pairing of each vertex on the left with a distinct vertex on the right.
The edges of the graph indicate the authorised pairings. If we begin by matching u0 to
v0 , then we immediately find ourselves blocked. To obtain a perfect matching, we are
forced to backtrack and undo this pairing. This principle is illustrated in the following
sections.
u0 v0 u0 v0 u0 v0
u1 v1 u1 v1 u1 v1
The flow algorithms presented here require that for every arc (u,v), the reverse arc
(v,u) also exists. Thus, beforehand they invoke the following function, which as nec-
essary completes the graph with reverse arcs of capacity 0. For every arc (u,v), it tests
if u is present in the adjacency list of v. This executes with complexity O(|E| · |V |),
negligible in our cases.
Matchings
bipartite without weights O(|V | · |E|) augmenting paths
bipartite with weights O(|V |3 ) Kuhn-Munkres
bipartite with preferences O(|V |2 ) Gale-Shapley
Flows
bounded capacities O(|V | · |E| · C) Ford-Fulkerson
bounded capacities O(|V | · |E| · log C) binary blocking flows
arbitrary capacities O(|V | · |E|2 ) Edmonds-Karp
arbitrary capacities O(|V |2 · |E|) Dinic
Cuts
arbitrary graph O(|V |2 · |E|) Dinic
planar graph O(|E| log |V |) Dijkstra
u0 v0 u0 v0
v1 u1 v1
u1
u2 v2 u2 v2
u3 v3 u3 v3
Input Output
Application
Assigning n courses to n classrooms, assigning n tasks to n workers—the applications
of the problem of maximum bipartite matching are widespread. A weighted version
of this problem exists that is described in Section 9.2 on page 145. As the applications
are most of the time modelled with bipartite graphs, we limit ourselves here to this
case, but the problem can also be solved for general graphs.
Definition
Let G(U,V ,E) be a bipartite graph with E ⊆ U ×V . A matching is a set M ⊆ E such
that no two edges in M have the same endpoint. Clearly, any subset of a matching is
also a matching. A matching M is said to be maximal if there is no matching M that
140 Matchings and Flows
Key Observation
A natural approach to tackle an optimisation problem is to begin with a naive solu-
tion and attempt to improve it iteratively until eventually it becomes optimal. To
understand how we might improve a matching, take a look at the symmetric differ-
ence between matchings. Given a blue matching M (shown as solid lines) and a red
matching M (shown as dashed lines), the symmetric difference M ⊕ M is defined
as M\M ∪ M \M. It is composed of cycles and paths of alternating colours, red and
blue, along the edges. By a counting argument, we see that if |M | > |M|, then an
alternating path P always exists that begins and ends with a red edge, see Figure 9.2.
M¢
M M¢
Figure 9.2 The symmetric difference M ⊕ M is made up of an alternating cycle, and a path P
augmenting for M and two paths augmenting for M.
With respect to M, this path P starts and ends on a free vertex and alternates
between edges that are not in M and edges in M. Such a path is said to be augmenting,
as the difference M ⊕ P is a matching with one edge more than M.
9.1 Maximum Bipartite Matching 141
Technical Details
In the following implementation,1 the matching is encoded by an array match associ-
ating with each vertex v ∈ V the matching vertex in U , where match[v] = None if v
is free.
The graph is represented by an array E such that E[u] is the list of neighbours
of u. The sets of vertices U,V must be given in the form [0,1,...,|U|-1] and
[0,1,...,|V|-1].
Note that this implementation uses Python’s lazy evaluation of boolean expressions,
a technique common to almost all programming languages. This means that in the OR
expression, the function augment is called if and only if the first term is false.
def max_bipartite_matching(bigraph):
n = len(bigraph) # same domain for U and V
match = [None] * n
for u in range(n):
augment(u, bigraph, [False] * n, match)
return match
1 Do you understand why only the vertices of V are marked when visited?
142 Matchings and Flows
Other Algorithms
An algorithm from Hopcroft and Karp (1972) allows the resolution of the maximum
√
matching problem in a bipartite graph in time O( |V | · |E|). Its principle is to dis-
cover several augmenting paths with the same traversal, and to privilege the short
√
augmenting paths. Alt et al. (1991) found an algorithm in O(|V |1.5 |E|/ log |V |)
which is interesting for dense graphs (with many edges). However, these algorithms
are complicated to implement and not that much faster in practice than the algorithm
presented in this section.
Application
Consider a city the map of which consists of n north–south streets and n east–west
streets, crossing in n2 intersections. We would like to monitor a number of the most
important intersections. A camera can monitor all of the intersections of any given
street. What is the minimum number of cameras that need to be installed in order to
monitor the given set of critical intersections?
Kőnig’s Theorem
Since each edge in a matching must have at least one endpoint in S, the size of the
maximal matching is a lower bound on the minimal size of the smallest covering.
Kőnig’s theorem states that there is, in fact, equality between these optimal values.
The proof is constructive and provides an algorithm for finding a minimal covering
from a maximal matching M of the graph G(U,V ,E). Let Z be the set of vertices
unmatched by M, and that are in U . We add to Z all the vertices reachable from Z by
alternating paths, see Figure 9.3. We define the set
S = (U \Z) ∪ (V ∩ Z).
The construction of Z implies that for each edge (u,v) ∈ M, if v ∈ Z, then u ∈ Z
also. Conversely, since u was not initially in Z with the free vertices of U , u was added
to Z thanks to a matching edge; hence v ∈ Z.
This implies that for each matching edge (u,v) ∈ M, either both its endpoints are
in Z or neither of them are. Hence, every matching edge has exactly one endpoint in
S, and |S| ≥ |M|.
We now show that S covers all of the edges of the graph. Let (u,v) be an edge of
the graph. If u ∈ Z, then the edge is covered, so suppose that u ∈ Z. If (u,v) is not
a matching edge, then the maximal character of Z means that v must also be in Z,
and hence in S. If (u,v) is in M, then by the preceding observation v ∈ Z, and hence
9.1 Maximum Bipartite Matching 143
u1 u2 u3 u4 u5 u6 u7 u8
v1 v2 v3 v4 v5 v6 v7 v8
v ∈ S. This shows that S is a covering of the vertices, implying |S| ≤ |M|, hence we
can conclude that |S| = |M|.
Problems
Crime Wave [icpcarchive:5584]
Dungeon of Death [spoj:QUEST4]
Muddy Fields [spoj:MUDDY]
u1 u2 u3 u4 u5 u6 u7 u8
v1 v2 v3 v4 v5 v6 v7 v8
A Solution
First, observe that if we can diminish the height of a pile without getting caught, then
we can shrink it down to a single crate. So—which piles are ripe for the taking? To
fix the notation, the height of a row or column is defined as the maximal height of
its piles.
Consider a stack of height h > 1 at the intersection of a row and column the heights
of which are strictly greater than h. Then, we can safely make off with all but one of
the crates in this stack, as there are other piles in the same row and same column that
will be higher.
Now, fix a height h > 1 and consider the bipartite graph the vertices of which
are the rows on one hand, and the columns on the other, and for which there is an
edge for each stack of height exactly h. We ignore the isolated vertices in this graph
and mark the vertices corresponding to the rows and columns with maximal height
exactly h. The non-isolated unmarked vertices have a value strictly greater than h. We
must find a minimal set S of edges in this graph covering all the marked vertices.
These edges correspond to piles that we must preserve in order to go undetected
by the cameras.
If an edge is between two unmarked vertices, it will not be part of an optimal
solution S, since it cannot be used to cover a marked vertex. Now consider two
edges (u,v) and (u,v ) with u,v marked and v unmarked. Then, an optimal solution
exists without the edge (u,v ), as we can always exchange it with the edge (u,v).
Hence, for each marked vertex u in this graph, if it has at least one marked
neighbour, then we can suppress all the edges incident to u leading to unmarked
neighbours.
Similarly, if a marked vertex exists without unmarked neighbours, we can choose
an arbitrary incident edge for the optimal solution S. The result of these rules is the
classic problem of the covering by edges in a bipartite graph, since all the remaining
edges are between marked vertices. Note that from a practical point of view, it is
not necessary to construct such a graph for each height h; instead, we can construct
the union of these graphs and invoke one single time an algorithm to find a minimal
covering by edges.
9.2 Maximal-Weight Perfect Matching—Kuhn–Munkres 145
v0 v1 v2 v3 v0 v1 v2 v3
u0 3 3 1 2 u0 3 3 1 2
u1 0 4 -8 3 u1 0 4 -8 3
u2 0 2 3 0 u2 0 2 3 0
u3 -1 1 2 -5 u3 -1 1 2 -5
input output
Application
In a teaching department there are n courses to distribute amongst n teachers. Each
teacher has provided preferences in the form of weights for each course. The problem
does not need normalisation, as some teachers might have more influence than others.
No one ever said that life had to be fair! The goal is to find a bijection between
courses and teachers maximising the total weights of the preferences over all of the
assignments. This is the problem of perfect maximum-weighted matching in a bipar-
tite graph.
Definition
Let G(U,V ,E) be a bipartite graph with weights w : E → R on the edges. Without
loss of generality, we can suppose |U | = |V | and that the graph is complete, i.e. E =
U × V . The goal is to find a perfect matching M ⊆ E maximising the weight
e∈M w(e) (also called profit).
Variants
A variant of this problem is to compute a perfect matching of minimal cost. In this
case, it suffices to change the signs of the weights and then solve the problem of
maximal profit. If |U | > |V |, it suffices to add new vertices to V , linked to all the
vertices of U with weight 0. There is then a correspondence between the perfect
matchings of this new graph with the maximum matchings of the original graph. If
the graph is not complete, it suffices to complete it with edges of weight −∞; these
will never be selected for an optimal solution.
Complexity
The Kuhn–Munkres algorithm, also known as the Hungarian method, has complexity
in time of O(|V |3 ). This name is a reference to the Hungarian mathematicians Denes
Kőnig and Jenő Egerváry, who inspired the American Harold Kuhn to discover this
method in 1955. James Munkres then proved, in 1957, its polynomial complexity (see
Kuhn, 1955; Munkres, 1957). In 2006, it was found that the German mathematician
Carl Gustav Jacob Jacobi had, in fact, discovered this algorithm in 1840, published
posthumously only in 1890, see Ollivier (2009a,b)!
146 Matchings and Flows
Algorithm
This is an algorithm belonging to the class of what is known as primal-dual algo-
rithms, so named because they are modelled using linear programming. The descrip-
tion given here is quite simple, and will not use the vocabulary of linear programming.
The principal idea consists of considering a related problem, that of minimal valid
labels. The labels are weights on the vertices, and they are valid if for every edge
(u,v)
(u) + (v) ≥ w(u,v). (9.1)
The dual problem consists of finding valid labels with minimal total weight. By sum-
ming the inequality over all of the edges (u,v) of a matching we easily see that the
total weight of any valid labelling is equal to or greater than the weight of any perfect
matching.
The labels define a set E made up of all the edges (u,v) such that
(u) + (v) = w(u,v).
The graph restricted to these edges is called the equality subgraph.
If ever we could find valid labels and a perfect matching M ⊆ E then the value
of E is equal to the value of |M|. Since the value of E bounds above that of
any perfect matching, this shows that M is of maximal weight.
The goal is thus to produce such a couple ,M. At any moment, the algorithm has
valid labels as well as a matching M ⊆ E . Then repeatedly, the algorithm will
attempt to augment the matching M, or if this is not possible, to improve the labels
(diminish their total value).
3 u0 v0 0 3 u0 v0 0 3 u0 v0 0 3 u0 v0 0
4 u1 v1 0 3 u1 v1 1 u1 v1 3 u1 v1 1
AU AU
3 1
AV AV
2 u2 v2 1 1 u2 v2 2 1 u2 v2 2 1 u2 v2 2
1 u3 v3 0 0 u3 v3 0 0 u3 v3 0 0 u3 v3 0
Figure 9.5 The steps of the Kuhn–Munkres algorithm. (1) The bipartite graph described by the
weight matrix given at the start of the section, with the labels of the vertices, the edges of the
equality subgraph and in thick lines the edges of a matching. The edges of a maximal
alternating tree rooted at u3 are shown in dotted lines. This tree covers the sets of vertices
AU ,AV and can no longer be extended. (2) The improvement of the labels reveals a new
(u1,v3 ). (3) The alternating tree includes the free vertex v3 . (4) The matching can be extended
along the path from u3 to v3 .
Progress
This procedure has the advantage that an additional edge between AU and V \AV will
be introduced into the equality graph. This is notably the case of the edge determining
the minimal slack , since its slack then becomes 0. Note that other edges of this
graph might disappear, but this is of no importance for this algorithm.
Initialisation
To initiate the algorithm, we begin with a valid labelling and an empty match-
ing M. To keep things simple, we choose (v) = 0 for every v ∈ V , and (u) =
maxv∈V w(u,v) for every u ∈ U .
148 Matchings and Flows
Implementation Details
The vertices of U are encoded by the integers 0 to n − 1, as are those of V . This
leads us to encode in two arrays lu and lv, one for each set U,V . Similarly, mu
and mv hold the matching, where mu[u] = v, mv[u] = v if u ∈ U is matched with
v ∈ V.
The free vertices are indicated by mu[u] = None or mv[v] = None. Finally, the
Boolean arrays au and av indicate whether a vertex is covered by an alternating tree.
Implementation Details
This time, we would like to know for each vertex in V if it belongs to the tree, and
if so, what is its predecessor. As the tree is not constructed by a depth-first search,
this information must be stored in an array. It will allow us to climb back to the
root when we need to update the matching. We assign Av[v] as the predecessor of
v ∈ AV , and set Av[v] = None for every vertex v ∈ V \AV outside of the tree. The
name of the array begins with a capital letter to distinguish it from the Boolean array
encoding AU .
The arrays slackVal and slackArg will be represented by a single array slack,
containing the pairs (val, arg). This implementation accepts graphs G(U,V ,E)
with |U | = |V |.
150 Matchings and Flows
Problem
Selfish Cities [spoj:SCITIES]
input output
Definition
We are given n white points and n black points in the plane. We suppose that these
are non-degenerate, in the sense that no three of the points are collinear. The goal is to
152 Matchings and Flows
find a perfect matching between the white points and the black points, such that if we
link all of the matched pairs with line segments, none of these intersect.
Key Observation
Imagine for an instant that the points are matched but some segments cross each
other. Consider two matched pairs u − v and u − v the segments of which cross,
see Figure 9.6. If we change their assignments to u − v and u − v, their segments no
longer intersect. Does this operation improve matters somewhat?
Note that the total number of crossings can increase with this ‘un-crossing’, as in
Figure 9.6. Hence, this is not the correct measure of progress. However, the total length
of the segments has diminished! This observation leads us to a first algorithm.
An algorithm from Matousek, Lo and Steiger (1994) computes this line in time
O(n). However, it is complex to implement.
How does this help us? With the guarantee that there are no three collinear points,
we obtain the following property. Either n is even, so that the separating line contains
none of the points, or n is odd and the separating line contains exactly one white
point and one black point. In the latter case, it is correct to match them. And in any
case, we can iterate on the two independent sub-instances separated by the line. This
recursive decomposition has a depth O(log2 n), hence the final algorithm is indeed in
O(n log n).
Definition
Suppose there are n women and n men. Each man ranks the women in order of
preference, and the women rank the men in the same way. A marriage is a perfect
matching in a complete bipartite graph between the men and the women. It is called
stable if a man i and a woman j do not exist such that the man prefers j to his
current wife, and the woman prefers i to her current husband. The goal is to compute
a stable marriage from the 2n stated preferences. The solution is not necessarily
unique.
Complexity
The Gale–Shapley algorithm (1962) solves this problem in time O(n2 ).
Algorithm
The algorithm starts with no married couples. Then, as long as there are single men,
the algorithm selects a bachelor i and the woman j that i likes best among the women
who have not yet been considered. The algorithm then seeks to marry i with j . The
wedding will take place if j is still single, or if j is married to a man k but prefers i
to k. In this case, k is summarily divorced and sent off to continue partying with the
group of single men.
Analysis
For the complexity, every couple i,j is considered at most once by the algorithm,
and the work to do for each couple takes constant time. For the validity, it suffices to
observe that as the algorithm proceeds, (1) a given woman is married to men that she
prefers more and more, while (2) a given man is married to women that he prefers less
and less. Hence, for a proof by contradiction of the algorithm’s validity, suppose that
at the end a man i married to a woman j and a woman j married to a man i exist such
that i prefers j to j and j prefers i to i . By observation (2), at some point in time
the algorithm examined the couple i,j . However, by (1), the algorithm never married
i with j . This means that at the moment the couple (i,j ) was considered, j must have
154 Matchings and Flows
been married with a man k that she preferred to i. This contradicts the fact that she
ends up being married to a man she likes less than i.
Implementation Details
In our implementation, the men are numbered from 0 to n − 1, as are the women. The
input consists of two arrays. The array men contains for each man his list of n preferred
women, in decreasing order. The array women contains the women’s corresponding
lists of preferences. The array women is first transformed into a array rank containing
for each woman j and man i the rank of preference that j has for i. For example, if
rank[j][i] = 0, then i is the preferred man of j , and if rank[j][i’] = 1, then i is
the second on the list for j , etc.
Finally, the array spouse associates with each woman the man with whom she is
currently2 married. The queue singles contains the single men, and for each man i,
current_suitor[i] indicates the index of the next woman to be wooed in his list of
preferences.
Problem
Stable Marriage Problem [spoj:STABLEMP]
2 Note from the translator: the divorce lawyers must be rubbing their hands with glee!
9.5 Maximum Flow by Ford–Fulkerson 155
1 4 3 10 t
10 2 8 5 10
s 10 2 9 4
input
1 4 3 9 t
10 6 5 10
s 9 2 9 4
output
Application
This algorithm can be used to find the bottlenecks in the evacuation of passengers
from an aircraft, or the resistance of an electrical network. In fact, numerous problems
can be modelled as a maximum flow in a graph.
Much of the pioneering work in this domain was done in the 1950s, in the context
of the Cold War. The original paper by Ford and Fulkerson was in fact a memorandum
of the RAND corporation, and was motivated by . . . studies of the Soviet Railway
system3 ! A top-secret report, ‘Fundamentals of a Method for Evaluating Rail Net
Capacities’, unclassified only in 1999, describes how this algorithm, and the related
minimum-cut algorithms (see Section 9.8), could be used to evaluate the capacity of
the Soviet railway system, and the optimal way to disable it! A fascinating article on
this topic is Schrijver (2002).
Definition
Formally, we are given a directed graph G(V ,A) with capacities on the arcs c : A →
R+ and two distinct vertices: a source s ∈ V and a sink t ∈ V . Without loss of
generality, we can suppose that for each arc (u,v) the arc (v,u) is also present in the
graph, since we can always add arcs of capacity zero without modifying the optimal
solution.
A flow is a function f : A → R, satisfying that:
∀(u,v) ∈ A : f (u,v) = −f (v,u), (9.2)
while respecting the capacities
∀e ∈ A : f (e) ≤ c(e), (9.3)
and respecting the conservation of flows at each vertex, other than at the source and
the sink:
The interpretation of f (u,v) > 0 is that f (u,v) units are flowing through the arc
(u,v). However, if f (u,v) < 0, then no flow is traversing the arc (u,v). The skew
symmetric condition (9.2) has been imposed to simplify the notation in (9.4).
The value of the flow is the quantity emanating from the source, v f (s,v). The
goal is to find the flow of maximal value.
Non-Directed Graphs
In a non-directed graph, we replace each edge (u,v) by the two arcs (u,v) and (v,u)
with the same capacity.
Implementation Details
In our implementation, the graph is described with an adjacency list; however, to ease
the manipulation of the flows, the flow is represented by a two-dimensional matrix F .
The procedure augment attempts to augment the flow along a path from u to target
with a value at most val. The procedure returns the amount of the augmentation of the
flow in the case of success and zero otherwise. It looks for this path via a depth-first
search and for this uses the marker visit.
9.5 Maximum Flow by Ford–Fulkerson 157
Problems
Potholers [spoj:POTHOLE]
Maximum Flow [spoj:FLOW]
158 Matchings and Flows
Observation
The Ford–Fulkerson and Goldberg–Rao algorithms have the annoying property that
their complexities depend on the given capacities. Happily, a small transformation
exists that renders the complexity independent of C (see Edmonds and Karp, 1972).
Surprisingly, applying the same algorithm, but first to the shortest augmenting
paths, results in a complexity that is independent of the maximal capacity. The idea
is that the length of the shortest augmenting path increases strictly for at most |E|
iterations. This is due to the following observations:
• Let Lf be the level graph of G, where s is at level 0, all the vertices reachable
from s by a non-saturated arc4 form the level 1, and so on. It is hence an acyclic
subgraph of the residual graph.
• An arc (u,v) is admissible if its flow is strictly below its capacity and if v is on the
next level after u.
• A shortest path from s to t in the residual graph is necessarily a path in the level
graph. If we augment the flow along this path, one of its arcs must become saturated.
• An augmentation along a path can render certain arcs non-saturated, but only arcs
towards lower levels. Hence, the arcs that become admissible cannot decrease the
distance from s to v counted in number of arcs. By symmetry, the distance from v
to t cannot decrease either.
• After at most |E| + 1 iterations, an arc (u,v) and its inverse (v,u) must have been
saturated. This proves that at some point the vertex u has changed level, and by the
preceding observation the distance from s to t has strictly increased.
• As there are only n levels, there are at most |V | · |E| iterations in all.
• The search for a shortest path is done with breadth-first search in time O(|V |). The
total complexity is thus O(|V | · |E|2 ).
Implementation Details
The shortest augmenting paths are found with a breadth-first search (BFS), with the
help of a queue Q. The array augm_path has two functions. First, it serves to mark
the vertices already visited by the BFS traversal.5 Then, rather than simply marking
the vertices with a Boolean, we store the predecessor on the path from the source
to the vertex in question. Thus, we can trace back to the source by following the
values of augm_path to augment the flow along this path. The array A contains for
each visited vertex v the minimal residual capacity along the path from the source
to v. With this array, we can determine by how much the flow can be augmented along
the augmenting path found.
4 An arc is saturated if the flow equals the capacity of the arc.
5 Do you understand why it is necessary to mark augm_path[source]?
9.7 Maximum Flow by Dinic 159
The idea of Dinic’s algorithm, discovered at the same time as the preceding algorithm,
is the following. Rather than seeking a set of augmenting paths one by one, up until
the distance between s and t increases, we find such a flow with a single traversal.
This reduces the complexity to O(|V |2 · |E|) (see Dinitz, 2006).
Imagine a function dinic(u,val) that attempts to pass a flow from u to t in the
level graph. The restriction is that the flow is not greater than val. The function then
returns the value of this flow, which can be improved along the path from s to u, by
the successive invocations that led to u. More precisely, to push a flow with value
val from u to t, we visit all the neighbours v of u in the level graph, and recursively
attempt to pass a maximum flow from v to t. The sum is then the maximum flow that
we can push from u to t.
160 Matchings and Flows
The function dinic(u,val) then detects if no additional flow can pass from u to t,
and t thus becomes inaccessible from u. In this case, it removes u from the level graph,
simply by indicating a dummy level −1. Hence, the subsequent calls will no longer
attempt to pass flow through u. Clearly, after O(n) iterations, s and t are disconnected,
and we compute a new level graph, see Figure 9.7.
0 1 2 3 4 [level]
1 4/4 3 4/9 t
s 0/9 2 0/9 4
1 4/4 3 4/9 t
s 4/9 2 4/9 4
1 3 4/9 t
0/2 0/5
s 4/9 2 4/9 4
Figure 9.7 The evolution of the level graph during the execution of Dinic’s algorithm. Only the
arcs from one level to the next are shown, labelled by f/c where f is the flow and c the
capacity. The dashed lines show the saturated, hence non-passable, arcs, while the dotted lines
show the arcs that will later become part of the level graph. After having augmented the flow
of value 4 along the path s − 2 − 4 − t, the vertex t becomes disconnected. It is thus necessary
to reconstruct the level graph, shown on the bottom.
Implementation Details
The array level contains the level of a vertex in the level graph. This graph is
recomputed as soon as the sink is no longer accessible from the source. The function
9.7 Maximum Flow by Dinic 161
dinic_step pushes as much flow as possible from u to the sink, without exceeding
the given limit. For this, it pushes a maximum of flow to its neighbours v in the level
graph, and accumulates in val the quantity of flow already pushed. When no further
flow can be sent, the vertex v can be removed from the level graph, by setting its level
to None.
Problem
Fast Maximum Flow [spoj:FASTFLOW]
Application
In a country, there are several roads linking various cities. A railway company wants
to attract customers to its new train line connecting city s with city t, and wishes to
place billboard advertisements along certain of these roads, so that every trip from s
to t must traverse one of these roads. The goal is to minimise the number of billboards
for this task.
Definition
An instance of this problem consists of a directed graph G(V ,A) with two distinct
vertices s,t and a cost for each arc c : A → R+ . An s − t cut is a set S ∈ V
containing s but not t. Equivalently, the cut S is sometimes identified by the arcs
leaving S, i.e. the arcs (u,v) with u ∈ S and v ∈ S, see Figure 9.8. The value of the
cut is the total cost of these arcs. The goal is to find the cut with minimal cost.
b 4 e
7
s 2 2 1
6 1
S a 3 c 8 d 7 t
Figure 9.8 A minimum s − t cut S of cost 3 + 2 + 2 + 1 = 8. The arcs leaving S are drawn with
solid lines.
The Max-Flow-Min-Cut theorem states that the value of the maximum flow is equal
to the value of the minimum cut.
This theorem was shown independently by Elias, Feinstein and Shannon (1956) on
one hand and by Ford and Fulkerson (1956) on the other. The proof follows from a
sequence of several simple observations.
1. For a flow f the quantity of the flow leaving a cut S, i.e. f (S) = u∈S,v∈S
f (u,v), is the same for every S. This is simply because for any w ∈ S,w = t, we
have (by (9.4) and then by (9.2))
9.9 s − t Minimum Cut for Planar Graphs 163
= f (u,v) = f (S ∪ w).
u∈S∪w,v∈S∪w
2. By (9.3), we have f (S) ≤ c(S), i.e. the flow leaving S is never greater than that
of S. This already proves half of the theorem—the value of the maximum flow is
at most the value of the minimum cut.
3. Now, if for a given flow f , we have f (S) < c(S) for every cut S, then an
augmenting path exists. For this, we start with S = {s} and A = ∅. Since
f (S) < c(S), an arc (u,v) exists with u ∈ S,v ∈ S,f (u,v) < c(u,v). We add (u,v)
to A and v to S, and repeat while t ∈ S. During this process, A is a spanning tree
of S rooted at the source, using only non-saturated arcs. Eventually, A will reach
the sink t and contain an augmenting path.
Algorithm
This observation leads us to an algorithm to solve the minimum cut problem. First,
we compute a maximum flow. Then, in the residual graph, we determine the set S of
vertices v reachable from s by arcs with positive residual capacity. Then, S does not
contain t, as the existence of an augmenting path contradicts the optimality of the flow.
The maximal character of S means that every arc leaving S is saturated by the flow, and
its residual capacity is zero. Hence, S is a cut of minimal value.
Planar Graphs
The s − t minimum cut problem can be solved more efficiently in the case of a planar
graph6 where the embedding in the plane is given. To simplify the presentation, we
suppose that the graph is a planar grid. This grid is made up of vertices linked by
arcs, which decompose the plane into areas known as faces. Imagine that the grid is
rectangular, that the source is the vertex on the bottom left and the sink is the vertex
on the top right.
Dual Graphs
In the dual graph, each face becomes a vertex, and there are two additional vertices s
and t . The vertex s represents the left and top boundary of the grid and the vertex t
the right and bottom boundary. In this new graph, two vertices are connected if the
6 A graph is planar if it is possible to draw it in the plane without edges crossing each other. Formally, a
mapping of the vertices to some points must exist (embedding) such that the arcs (segments) intersect
only at their endpoints.
164 Matchings and Flows
corresponding faces are adjacent. The edge between the dual vertices has the same
weight as the primal edge at the contact line between the faces.
Key Observation
Every s − t path in the dual graph of length w corresponds to an s − t cut in the
primal graph with the same value w, and vice-versa.
3 1 2 t s’
3 1 2
6 2 8 10
6 2 8 10
1 7 2
1 7 2
4 3 1 4 4 3 1 4
6 4 1
6 4 1
5 2 1 3
5 2 1 3
7 6 2
s 7 6 2 t’
Figure 9.9 On the left, the primal graph. On the right, the dual graph: a minimal-length
s − t -path corresponds to an s − t minimum cut.
Variants
An important problem is to find a minimum cut which disconnects the graph, i.e. to
find the smallest s − t cut over all pairs s,t of vertices. This problem can be solved
with the Stoer–Wagner algorithm in time O(|V | · |E| + |V |2 log |V |), which is more
interesting than solving the (|V |2 ) s − t cuts for each pair s,t of vertices, see (Stoer
and Wagner, 1997).
Problem
Landscaping [kattis:landscaping]
two vertices s,t to the graph, and then link s to all the black vertices and t to all the
white vertices. Then, an s − t cut in this graph is a set S ∪ {s} and the value of the cut
is exactly the cost of S as described above, see Figure 9.10.
s t
Figure 9.10 From left to right: the given graph; the instance of the min-cut problem with its
optimal solution shaded; and the corresponding solution of cost $3, where two vertices are
coloured white and one edge has endpoints with distinct colours.
A generalisation of the flow problem is the transport problem, where several sources
and sinks exist. Here, each vertex v has a value dv , representing the supply when it is
positive and the demand when it is negative. The instance must satisfy v∈V dv = 0 in
order to have a solution. The goal is to find a flow transporting the quantities proposed
to the customers, hence a flow respecting the capacities and such that for each vertex
v the flow entering minus the flow leaving is equal to dv . This problem can easily be
reduced to a flow problem by adding a source and a sink vertex and linking the source
to all vertices with a supply, linking all vertices with a demand to the sink and setting
the appropriate capacities to these new arcs.
There is a generalisation of the flow problem, where arcs have costs per unit flow in
addition to capacities. For example, if an arc e has a flow fe = 3 and a cost we = 4, the
cost generated by the flow on this arc is 12. The goal is then to find a flow minimising
the total cost. The flow can either be a maximum flow or a flow respecting all the
demands in the case of a transport problem. The latter is known as a minimal cost
transport problem.
For this problem, algorithms exist that are similar to the Kuhn–Munkres Hun-
garian algorithm (see Ahuja et al., 1993). A first approach could be to start with
a maximum flow, and then try to repeatedly augment the flow along negative-cost
cycles.
Two interesting links exist between the maximum bipartite matching problem and the
maximum flow problem.
166 Matchings and Flows
7 In principle, fractional flows are possible, but all reasonable flow algorithms produce an integer flow if
the capacities are integers.
9.12 Width of a Partial Order—Dilworth 167
c− c+
+
a
d− g−
g+
+
d
e− e+ h−
c +
b
f−
a d g
b e h
s t
f i i−
f+
Figure 9.11 Reduction of the unit-capacity s − t flow problem to a maximal bipartite matching
problem. The maximum flow and corresponding maximum matching are depicted with solid
lines.
8 8
5 5
6 6
4 7 4 7
2 2
3 3
1 1
Input8 Output
Definition
Suppose we are given a partial order, i.e. an acyclic and transitive directed graph
G(V ,A). A chain in G is a directed path and by extension the set of vertices it covers.
A chain decomposition is a set of chains such that every vertex is covered by exactly
one chain. Its size is defined as the number of its chains. An antichain in G is a set S
168 Matchings and Flows
of vertices such that u,v ∈ S with (u,v) ∈ A does not exist. The width of the partial
order is defined as the size of the largest antichain. The problem is to compute this
width.
An important relation between chains and antichains is given by Dilworth’s theo-
rem (1950):
The size of the largest antichain equals the size of the smallest chain decomposition.
Showing that the size of the largest antichain is at most the size of the small-
est chain decomposition is straightforward. Assume that there is an antichain longer
than a smallest decomposition in chains. Then, according to the pigeonhole princi-
ple, at least two elements of the antichain should be in the same chain, which is a
contradiction.
It is possible to prove the other direction either by induction on the size of the graph,
or as a reduction to Kőnig’s theorem (see Section 9.1 on page 142): we construct a
bipartite graph H (V ,V ,E) such that V is a copy of V and for each arc (u,v) ∈ A in
the original graph, there is (u,v ) ∈ E, where v is a copy of v.
Let M be a maximum matching in the new graph H . By Kőnig’s theorem, a set S
of vertices exists that touches every edge of H , and |M| = |S|.
Every matching in H corresponds to a decomposition in chains: following the edges
of the matching, we can recover chains in G. Every unmatched node in V corresponds
to an endpoint of a chain. Hence, the size of the decomposition is the number of
endpoints |V | − |M|, which is minimum because M is maximum.
We now want to build a longest antichain from S. If we select the nodes v in V
such that neither v nor v is in S, then we have at least |V | − |S| elements. If there was
an edge between any pair (v,v ) of these elements, then either v or v would be in S,
which is a contradiction. Thus, the selected elements form an antichain, and we have
proved Dilworth’s theorem.
8 8 8’ 8 8’ 8
7 7 7’ 7 7’ 7
6 6 6’ 6 6’ 6
5 5 5’ 5 5’ 5
4 4 4’ 4 4’ 4
3 3 3’ 3 3’ 3
2 2 2’ 2 2’ 2
1 1 1’ 1 1’ 1
Figure 9.12 From left to right: a partial order G, the associated bipartite graph H with a
covering of the vertices S marked in grey, a maximum matching in H , the decomposition into
minimum chains in G associated with an antichain marked in grey.
a precise location in the city at a certain time, in order to reach a given destination.
We suppose that all the source–destination travel times are known. We can thus
decide if the same taxi can honour the reservation j after having completed the
reservation i. In this case, we denote i j , and the relation is a partial order
on the day’s reservations. The solution to the problem is simply the width of this
partial order.
Generalisation
If the graph is weighted, we can look for a minimal decomposition into chains which
also minimises the total weight of the arcs. This problem can be solved in time
O(|V |3 ) by reduction to the problem of minimal-cost perfect matching in a bipartite
graph.
Implementation Details
Our implementation returns an array p encoding an optimal partition, where p[u] is
the index of the chain containing the vertex u, where the chains are numbered starting
with 0. The input is a square matrix M indicating for each pair of vertices u,v the cost
of the arc (u,v) or None if there is no such arc. This implementation supposes that the
vertices are given in topological order, otherwise a topological sort must be done first,
see 6.7. Hence, if the arc (u,v) exists, then u < v.
170 Matchings and Flows
def dilworth(graph):
n = len(graph)
match = max_bipartite_matching(graph) # maximum matching
part = [None] * n # partition into chains
nb_chains = 0
for v in range(n - 1, -1, -1): # in inverse topological order
if part[v] is None: # start of chain
u = v
while u is not None: # follow the chain
part[u] = nb_chains # mark
u = match[u]
nb_chains += 1
return part
Problems
Stock Charts [gcj:2009round2C]
Taxi Cab Scheme [icpcarchive:3126] [poj:2060]
10 Trees
Rooted trees are combinatorial structures which appear naturally when considering
objects with a structure of recursive decomposition. This is the case, for example, for
classifications, hierarchies or genealogies. The recursive nature of trees is an invita-
tion to recursive algorithms, and the key to an algorithm is to find the appropriate
method of exploration. In this section, we review certain classic problems concerning
trees.
Formally, a tree is a connected acyclic graph. One of its vertices can be designated
as the root, and in this case we generally speak of a rooted tree. This root provides
the tree with an orientation from parent to child. Starting from a vertex and climbing
up the links, we arrive at the root. The vertices without children are called leaf nodes.
A tree with n vertices contains exactly n − 1 edges. For proof, it suffices to observe
that if we repeatedly remove a leaf from the tree, along with its incident edge, we end
up with an isolated vertex, i.e. a tree with 1 vertex and 0 edges.
The vertices of a rooted tree are partitioned into levels. The root is the unique vertex
at level 0, its child nodes are at level 1 and so on. The largest non-empty level number
defines the depth of the tree, which is also the maximum distance from the root to a
leaf. The subtree rooted at vertex v consists of all vertices and edges reachable from
the root of the original tree only through vertex v. A disjoint union of trees is called a
forest, as the reader might have guessed.
There are numerous dynamic data structures based on trees, such as binary red-
black search trees or interval trees. These structures invoke rebalancing operations on
the trees in order to guarantee logarithmic-time insert, delete and query operations.
However, in programming competitions, the inputs are given only once, and thus it
is often possible to skip the insert/delete operations by directly constructing balanced
structures.
A tree can be represented by essentially one of two methods. The first is the clas-
sic representation of a graph with an adjacency list, which does not distinguish any
particular vertex as the root of the tree. Another representation commonly used is that
of an array of predecessors. For a rooted tree (often with vertex 0 as the root), each
vertex other than the root has a unique predecessor, encoded in the array. Based on the
context, one or the other of the representations will be best suited, while translating
between them can easily be done in linear time.
172 Trees
Definition
Let be a finite alphabet. We denote by ∗ the set of all words on . Let S ⊆ ∗
be a finite set of words. A binary code for an alphabet of n characters is a function
c : → {0,1}∗ such that no word c(a) of the code is a prefix of another word c(b) for
a,b ∈ S. This code applied to each of the letters transforms a word from S into a word
in {0,1}∗ , and the property for the prefixes allows an unambiguous decoding of the
original. In general, we would like the encoding to be as short as possible, hence the
more frequent letters should be encoded by shorter words. Formally, given a frequency
function f : S → R+ , we seek a code minimising the cost
f (a) · |c(a)|.
a∈S
we begin with a forest, where each word from S is a tree consisting of a single node,
labelled by its frequency. Then, as long as there are multiple trees, we merge the two
trees with minimal frequency by attaching them under a new root, and label the new
root with the sum of the frequencies of the two subtrees, see Figure 10.1. The arcs
connecting the root to the merged subtrees are labelled, respectively and arbitrarily,
with 0 and 1. By an exchange argument, we can show that the encoding thus produced
is optimal (see Huffman, 1952).
To manipulate the trees, we place them in a data structure allowing efficient inser-
tion of elements and retraction of the minimal element: a priority queue. These opera-
tions have a logarithmic cost in the number of objects in the structure. In general, this
structure is implemented with a heap. In Python, such a structure can be found in the
module heapq.
The elements stored in the priority queue are pairs (fA,A) where A is a binary
tree and fA is the total frequency of the letters stored under A. A tree is encoded in
two ways. A letter a represents a tree consisting of a unique leaf (and root) a. A tree
composed of a left subtree and a right subtree r is represented by the pair (,r),
which must be a list instead of a tuple to avoid deep copies of the structure.
def huffman(freq):
forest = [] # build forest with singletons
for item in freq:
heappush(forest, (freq[item], item))
while len(forest) > 1:
(f_l, left) = heappop(forest) # merge two trees
(f_r, right) = heappop(forest)
heappush(forest, (f_l + f_r, [left, right]))
code = {} # build code from unique tree
extract(code, forest[0][1], [])
return code
Problems
Huffman Trees [poj:1261]
Variable Radix Huffman Encoding [spoj:VHUFFM]
174 Trees
A D
B M C J
H I L G K E F O
Input: B, G
Output: A
Definition
Given a tree with n nodes, we would like to answer, in time O(log n), queries for the
lowest common ancestor or LCA of two given nodes u,v, i.e. the unique node u such
that u and v are contained in the subtree rooted at u and with no direct descendant of
u having this property.
Implementation Details
We suppose that the tree is given in the form of an array parent, indicating the
parent node parent[u] for each node u ∈ {0,1, . . . ,n − 1} in the tree. We also
suppose that a parent has a smaller index than its children, and that the root is the
vertex 0.
10.2 Lowest Common Ancestor 175
9 15 31 25 20
A B C D E
24
0 1
9 15 31 25 20
A B C D E Letter Frequency Code
Input Output
44 A 7 00
0 1
B 7 01
C 7 10
20 24 D 7 11
E 0 1 A 40 1
B 5 01
9 15 31 25 C 2 001
A B C D D 1 000
44
0 1
20 24 56
E 0 1 0 1
9 15 25 31
A B D C
100
0 1
44 56
0 1 0 1
20 24 25 31
E 0 1 D C
9 15
A B
Figure 10.1 Construction of a Huffman code. Each node is labelled by the sum of the
frequencies of the letters at the leaves of the subtrees rooted at this node. On the right, two
examples of Huffman codes resulting from two different inputs.
176 Trees
class LowestCommonAncestorShortcuts:
def __init__(self, prec):
n = len(prec)
self.level = [None] * n # build levels
self.level[0] = 0
for u in range(1, n):
self.level[u] = 1 + self.level[prec[u]]
depth = log2ceil(max(self.level[u] for u in range(n))) + 1
self.anc = [[0] * n for _ in range(depth)]
for u in range(n):
self.anc[0][u] = prec[u]
for k in range(1, depth):
for u in range(n):
self.anc[k][u] = self.anc[k - 1][self.anc[k - 1][u]]
Implementation Details
The input to this implementation is the graph in the form of an adjacency list. It makes
no hypothesis as to the numbering of the vertices. Consequently, the trace dfs_trace
contains not only the vertices but the couple (depth, vertex). As the input could be very
large, the DFS traversal is implemented iteratively using a stack to_visit. The array
next indicates for each vertex how many descendants have already been explored.
2 9
3 5
4 6 7 8
1 2 3 4 3 2 5 6 5 7 5 8 5 2 1 9 1
Figure 10.2 Reduction of the problem of lowest common ancestor to the problem of the
minimum over a range of the trace of a depth-first traversal. The query illustrated is LCA(3,7),
which returns the node 2.
class LowestCommonAncestorRMQ:
def __init__(self, graph):
n = len(graph)
dfs_trace = []
self.last = [None] * n
to_visit = [(0, 0, None)] # node 0 is root
succ = [0] * n
while to_visit:
level, node, father = to_visit[-1]
self.last[node] = len(dfs_trace)
dfs_trace.append((level, node))
if succ[node] < len(graph[node]) and \
graph[node][succ[node]] == father:
succ[node] += 1
if succ[node] == len(graph[node]):
to_visit.pop()
else:
neighbor = graph[node][succ[node]]
succ[node] += 1
to_visit.append((level + 1, neighbor, node))
self.rmq = RangeMinQuery(dfs_trace, (float(’inf’), None))
def query(self, u, v):
lu = self.last[u]
lv = self.last[v]
if lu > lv:
lu, lv = lv, lu
return self.rmq.range_min(lu, lv + 1)[1]
178 Trees
Variant
This structure also permits the computation of the distance between two nodes of the
tree, as the shortest path must pass through the lowest common ancestor.
Problem
Lowest Common Ancestor [spoj:LCA]
Definition
Given a tree, we want to compute the longest path in this tree.
Complexity
The following algorithm is linear in the size of the tree.
The program can dispense with tests on the number of children by using -1 as the
default value. Note that it is not necessary to sort the children of a node to obtain the
two children maximising their b-value.
Some Observations
Let r be an arbitrary vertex of the tree and u a vertex at the maximal distance from r.
Then, a longest path exists in the tree with endpoint u. To be convinced, suppose we
have a longest path with endpoints v1,v2 . Let u be a vertex on the path from r to u
and v a vertex on the path from v1 to v2 that minimises the distance from u to v . If
the two paths intersect, we can choose any vertex in the intersection, see Figure 10.3.
Assign d as the distance in the tree. By the choice of u we have
This implies d(u,v ) = 0, i.e. u = v , and d(v,v1 ) = d(v,u) , thus the path from v2
to u is also a longest path in the tree.
r u
u
v
v1 v2
Figure 10.3 An exchange argument shows that if u is a vertex at a maximal distance from r then
a longest path exists starting from u.
Variant
We wish to remove the fewest possible edges in a tree so that in the remaining forest
no path is longer than R. To solve this problem, it suffices to remove the critical edges
determined in the dynamic programming solution above.
Consider the processing of a vertex v with its children u1, . . . ,ud ordered such that
b[u1 ] ≤ . . . ≤ b[ud ]. For every i from d down to 1, we must remove the edge (v,ui )
if b[ui ] + 1 > R or if i > 1 and b[ui−1 ] + 2 + b[ui ] > R.
Problems
Labyrinth [spoj:LABYR1]
Longest path in a tree [spoj:PT07Z]
SpeedCameras [codility:calcium2015]
Definition
Given an undirected connected graph, we would like to find a set of edges T such
that every pair of vertices is linked by a path through these edges. In this case, T is
called a spanning tree. The edges have weights, and we want the total weight of our
180 Trees
set of edges to be minimal. Note that such a set is acyclic—and hence a tree—since
the deletion of an edge of a cycle preserves the connectedness. We thus seek what is
known as a minimum weight spanning tree, see Figure 10.4.
Figure 10.4 Minimum weight spanning tree in a complete graph of 256 vertices where the
weight of each edge is the Euclidean distance separating its endpoints.
Application
The edges of a graph are given with a cost (or weight) w, and we need to purchase at
least cost some edges to render the resulting subgraph connected. Thus, we must find
a set T ⊂ E such that G(V ,T ) is connected with e∈T w(e) minimal.
Implementation Details
To visit the edges in order of increasing weight, we create a list of couples (weight,
edge). This list is sorted lexicographically on the couples, then processed in order. We
use a union-find structure (see Section 1.5.5 on page 26) in order to maintain a forest
and efficiently detect if adding an edge forms a cycle.
10.4 Minimum Weight Spanning Tree—Kruskal 181
Prim’s Algorithm
Another algorithm for this problem exists, namely Prim’s algorithm. It works in a
similar manner to Dijkstra’s algorithm. It maintains for a set of vertices S a priority
queue Q containing all the edges leaving S. Initially, S contains a single arbitrary
vertex u. Then, while S does not contain all the vertices, an edge (u,v) with minimal
weight is extracted from Q with u ∈ S,v ∈ S. The vertex v is then added to S and Q
updated accordingly. The complexity of this algorithm is again O(|E| log |E|).
Problems
Cobbled streets [spoj:CSTREET]
Minimum Spanning Tree [spoj:MST]
11 Sets
This chapter extends the chapter on sequences. It illustrates how a single method,
specifically dynamic programming, allows the resolution of a large variety of prob-
lems. We begin with the presentation of two classic problems—the knapsack and
making change.
Definition
Given n objects with weights p0, . . . ,pn−1 and values v0, . . . ,vn−1 , and given a knap-
sack with a capacity C, where C is an integer, the problem is to find a subset of
the objects with maximal total value, the total weight of which does not exceed the
capacity C. This problem is NP-hard.
Key Observation
For i ∈ {0, . . . ,n − 1} and c ∈ {0, . . . ,C}, assign Opt[i][c] as the largest value obtain-
able among objects with index 0 to i without their weight exceeding the capacity c. For
the base case i = 0, we have Opt[0][c] = 0 if p0 > c and Opt[0][c] = v0 otherwise.
For larger values of i, at most two possible choices appear for the object of index i:
we can take it or we can leave it. In the first case, the available capacity is reduced by
pi . We thus have the relation:
⎧
⎨ Opt[i − 1][c − pi ] + vi case where we take the object,
Opt[i][c] = max provided c ≥ pi
⎩
Opt[i − 1][c] case where we leave the object.
Algorithm in O(nC)
An algorithm with such a complexity is said to be pseudo-polynomial.1 A Boolean
matrix Sel is maintained in parallel with the dynamic programming matrix Opt.
1 A pseudo-polynomial algorithm is one that is polynomial in the value of the inputs but not in their size.
Here, adding an extra bit to the size of C doubles the running time.
11.1 The Knapsack Problem 183
Figure 11.1 A representation of the array Opt. The computation of each element requires at
most two elements in the preceding row, including the one immediately above. This is a
special case of the problem of the longest path in a grid, described in Section 3.1 on page 62.
It allows us to remember the choices made leading to the values stored in Opt. Once
these matrices are populated following the recurrence formula described above, a
traversal of the elements in reverse order allows the extraction from the matrix Sel
the set of elements providing the optimal total value.
Problems
The Knapsack Problem [spoj:KNAPSACK]
Knapsack [spoj:KNPSACK]
184 Sets
Suppose we need to make change for a customer, of value R—the difference between
the purchase price and the amount tendered. We have coins or banknotes at our dis-
posal, of respective value x0, . . . ,xn−1 cents (or whatever!). The problem thus con-
sists of determining the existence of a positive linear combination of x0, . . . ,xn−1
that equals R. You might laugh, as this seems so simple, but in Burma (now known
as Myanmar), there used to be banknotes of 15, 25, 35, 45, 75 and 90 kyats (see
Figure 11.2).
For this problem, banknotes or coins of value xi can contribute multiple times to
the sum (or at least until the merchant runs out . . . ). We compute a Boolean array b,
with b[i,s] indicating whether a combination exists of the coins indexed from 0 to i,
of value exactly s. Clearly, b[−1,0] = T rue, corresponding to the empty coin set,
and b[−1,s] = F alse for all s > 0. As a solution either does contain a coin of value
xi or does not, we have the recursion b[i,s] = b[i,s − xi ] ∨ b[i − 1,s]. The imple-
mentation below avoids the use of index i on b by processing its elements in the right
order.
Variant
In the case where a solution exists, we can try to determine a solution using a minimal
number of banknotes or coins.
Observation
You might think that it suffices to consider the currency by decreasing value and
provide at each step the largest banknote or coin that does not exceed the remaining
sum due. This realises a minimal amount of currency for most monetary systems.
11.3 Subset Sum 185
However, if we consider a system with values 1, 3, 4 and 10 and wish to make change
for a sum of 6, the greedy approach results in a solution with the three values 4+1+1,
whereas the unique optimal solution is 3 + 3.
Algorithm in O(nR)
Let A[i][m] be the minimal number of coins to return for an amount 0 ≤ m ≤ R using
the monetary system x0, . . . ,xi−1 , with A[i][m] = ∞ if there is no such solution. We
can derive a recurrence relation similar to that obtained for the knapsack problem: for
any amount m, A[0][m] is m/x0 if x0 divides m and ∞ otherwise. For i = 1, . . . ,n−1
we have the relation:
⎧
⎨ A[i][m − xi ] + 1 case where we can use the coin,
A[i][m] = min provided m ≥ xi
⎩
A[i − 1][m] case where we cannot use the coin.
input output
Definition
Given n positive integers x0, . . . ,xn−1 , we would like to know if a subset exists the
sum of which is exactly a given value R. This problem is NP-hard. This problem has
the flavour of the change-making problem, with the twist that every coin can be used
at most once.
Algorithm in O(nR)
We maintain a Boolean array indicating for each index i and every 0 ≤ s ≤ R whether
a subset of the integers x0,x1, . . . ,xi exists the sum of which is s.
Initially, the entries in this array are true only for the index 0 and the empty set.
Then, for every i ∈ {0, . . . ,n − 1} and every s ∈ {0, . . . ,R}, we can produce a subset
of sum s with the integers x0, . . . ,xi , if and only if a subset exists of x0, . . . ,xi−1 of
sum either s or s − xi . We proceed as for change-making, with a slightly different
recursion, namely b[i,s] = b[i − 1,s − xi ]|b[i − 1,s].
Again using a carefully chosen loop order for s, our implementation avoids the use
of index i on the array b.
186 Sets
b[s] is true and s is as close as possible to ( xi )/2. Hence, for every a = 0,1,2, . . .
and d = +1, − 1, consider b[( xi )/2 + a · d].
Problem
Boat Burglary [spoj:BURGLARY]
Definition
Given a set of n integers x0, . . . ,xn−1 , we want to know if it is possible to select k of
them the sum of which is zero. Without loss of generality, we assume that no element
has value 0, because a solution using 0 is also a solution using only k − 1 elements; a
simpler problem.
Application
For k = 3, this problem is important for discrete geometry, as numerous classic
problems can be reduced to the 3-sum problem. For example, given n triangles, do
they entirely cover another target triangle? Or, given n points, does a line exist that
intersects at least three of these points? For all these problems, algorithms exist in
O(n2 ). This is essentially optimal, as it is conjectured that an algorithm does not exist
with time complexity O(n2−(1) ) (Gajentaan and Overmars, 1995). For larger values
of k, this problem is important for cryptography.
Algorithm in O(nk−1 )
First, an observation for the case k = 2—the problem reduces to deciding if i,j with
xi = −xj exists. If x is sorted, a parallel run-through—as in the merge of two sorted
lists—allows the resolution of the problem (see Section 4.1 on page 73). Otherwise,
we can store the input in a hash table and then for each xi see if −xi can be found in
the table.
For k = 3, we propose an algorithm in time O(n2 ). We start by sorting x. Then,
for each xj , we look to see if the lists x + xj and −x have an element in common.
Such an element would be of the form xi + xj and −xk with xi + xj + xk = 0, with
i,j,k necessarily being distinct indices because of the absence of zero in the input.
This approach can be generalised for larger values of k, but is less efficient than the
following algorithm.
Implementation Details
To avoid taking several times the same indices, we store in A,B not only the sums but
also couples consisting of each sum as well as the indices realising it. Thus, we can
verify for each couple in A and B with sum 0 if the indices are disjoint.
Problem
4 Values whose Sum is 0 [poj:2785]
12 Points and Polygons
The main object at the heart of computational geometry is the point, representing a
position in some space. A number of classic problems concerning sets of points in the
plane are presented in this chapter.
Points in the plane are naturally represented by their pair of coordinates. An impor-
tant basic operation on points is the orientation test, see Figure 12.1: given three points
a,b,c, are they aligned, i.e. all contained in some line, and if not, will the transforma-
tion a → b → c move the points counter-clockwise (to the left) or clockwise (to the
right)? Formally, we look at the sign of the z-component of the vector cross product
! bc.
ab× ! If this is positive, then the transformation moves the points counter-clockwise,
if negative then clockwise, and if zero, the points are aligned.
Figure 12.1 The test left_turn(a,b,c) returns True for these points.
Input Output
Definition
Given a set of n points, the convex hull is defined as the smallest convex polygon
that encloses all of the points of the set. This polygon is also the polygon of smallest
perimeter containing all of these points.
Imagine a wooden board covered with nails. A rubber band enclosing these nails
represents the convex hull of the nail positions.
Implementation Details
The lower part bot of the convex hull is obtained in the same manner. The result is
the concatenation of the two lists, reversing the list top in order to obtain the points
of the convex hull in the ‘normal’ order. This is defined as the reverse direction of the
hands of an analogue clock (a quaint object used by preceding generations to keep
track of time, extremely accurate twice each day if perchance it was not wound up).
Note that the first and last elements of the lists are identical, and would therefore
be duplicated in the resulting concatenation. Hence, it is important to remove these
endpoints.
In order to simplify the code presented, the point p is added to the lists top and bot
only after the deletion of elements that would make the sequence non-convex.
def andrew(S):
S.sort()
top = []
bot = []
for p in S:
while len(top) >= 2 and not left_turn(p, top[-1], top[-2]):
top.pop()
top.append(p)
while len(bot) >= 2 and not left_turn(bot[-2], bot[-1], p):
bot.pop()
bot.append(p)
return bot[:-1] + top[:0:-1]
Problems
Doors and Penguins [spoj:DOORSPEN]
Build the Fence [spoj:BSHEEP]
Input
The first line consists of the integers N and R, separated with a space, where N is
Jacques-Édouard’s age in hours. Then N lines follow, each of them consisting of the
two integer coordinates xi and yi of the ith candle in nanometers, separated with a
space.
Limits
• 3 ≤ N ≤ 2 · 105 ;
• 10 ≤ R ≤ 2 · 108 ;
• for 1 ≤ i ≤ N , xi2 + yi2 ≤ R 2 ;
• all points have distinct coordinates.
Output
Print the value W as a floating point number. An additive or multiplicative error of
10−5 is tolerated: if y is the answer, any number either within [y − 10−5 ;y + 10−5 ]
or within [(1 − 10−5 )y;(1 + 10−5 )y] is accepted.
3 10 7.0710678118654755
00
10 0
0 10
Solution
Without loss of generality, an optimal solution consists of a strip which touches two
of the given points on one border of the strip and another single point on the opposite
border of the strip, see Figure 12.2. Clearly all three points belong to the convex hull
of the given input points and without loss of generality, the two points on the first
border of the strip are consecutive on the convex hull. With this observation done,
the algorithm is quite easy to find. First compute the convex hull. Then loop over all
successive pairs of points pi ,pi+1 along the hull, and at every step maintain a point
pk which is furthest from the line going through the points pi ,pi+1 . This means that
while the pairs of points loop in counter-clockwise order over the convex hull, the
opposite point also moves in counter-clockwise order over the convex hull. But when
the pair makes one step on the hull, the opposite point might remain the same or take
several steps. Nevertheless, the overall complexity of this part of the algorithm is only
O(n) in total. The minimum distance observed for all of these triplets pi ,pi+1,pk is
the solution to the problem.
The bottleneck in the overall running time is the computation of the convex hull,
which takes time O(n log n), acceptable complexity given the upper bounds of the
problem.
12.2 Measures of a Polygon 193
pi
pi+1
pk
Figure 12.2 The structure of an optimal solution. The points of the convex hull are
shown in white.
Figure 12.3 The number of integer points on the boundary (in black) is 4 the number of integer
points inside (in white) is 4, and the area is 4 + 4/2 − 1 = 5.
Area
We can compute the area A in linear time via the formula
n−1
A= (xi yi+1 − xi+1 yi ),
2
i=0
where the index i + 1 is taken modulo n. The ith term corresponds to the signed area
of the triangle (0,pi ,pi+1 ), whose sign depends on the orientation of the triangle.
Removing this triangle reduces the problem to the computation of the area of a
polygon having one less point on the boundary, hence the area of a polygon can be
expressed as a sequence of additions and subtractions of areas of triangles.
def area(p):
A = 0
for i, _ in enumerate(p):
A += p[i - 1][0] * p[i][1] - p[i][0] * p[i - 1][1]
return A / 2.
Input Output
12.3 Closest Pair of Points 195
Problem
Toil for Oil [spoj:OIL]
Input Output
Application
In a campground, tents are pitched arbitrarily, and each is occupied by someone with
a radio. We would like to impose a common threshold for the noise level not to be
exceeded, giving a minimum distance between tents, so that nobody is bothered by
the music of their neighbours.
Definition
Given n points p1, . . . ,pn , we want to find the pair of points pi ,pj minimising the
Euclidean distance between pi and pj .
•
•
•
p
◦
• •
• •
•
•
Figure 12.4 Each cell of the grid contains at most one point. At the time a new point p is
considered, it suffices to measure its distance with the points contained in the neighbouring
cells (white).
see Figure 12.4. If a pair of points at distance d < d is discovered, the procedure is
restarted from the beginning with a new grid of step d /2.
Complexity
We suppose that the access to G takes constant time, as does the computation of the
cell containing a given point. The key argument is that if the points given as input are
chosen in a uniformly random order, then when the ith point is processed (3 ≤ i ≤ n),
we improve the distance d with probability 1/(i − 1). Hence, the expected complexity
is on the order of ni=3 i/(i − 1), hence linear in n.
Implementation Details
To compute the cell associated with a point (x,y) in the grid with a given step, it
suffices to divide each coordinate by step and then round down. In Python, the integer
division by the operator // returns the required value, but as a floating point number.
We must convert it to integer type in order to loop over the neighbouring cells, as
the function range admits only integer parameters. In other languages, particular
attention must be paid to negative coordinates. For example, in C++, the expression
(int)(-1.1 / 2) returns 0 rather than the -1 required, and hence the function floor4
of the standard math library should be used.
4 The floor of a real number x is the greatest integer equal to or less than x.
12.3 Closest Pair of Points 197
def closest_points(S):
shuffle(S)
assert len(S) >= 2
p = S[0] # start with distance between
q = S[1] # first two points
d = dist(p, q)
while d > 0: # distance 0 cannot be improved
r = improve(S, d)
if r: # distance improved
d, p, q = r
else: # r is None: could not improve
break
return p, q
Problems
Closest Point Pair [spoj:CLOPPAIR]
Closest Triplet [spoj:CLOSEST]
198 Points and Polygons
Input
Output not simple simple
Definition
A polygon is rectilinear if its edges alternate between vertical and horizontal orien-
tations. It is simple if its edges do not cross each other. The goal is to test whether a
given rectilinear polygon is simple.
Observation
If a polygon is rectilinear, then each of its vertices has a vertical neighbour as well as
a horizontal neighbour among those that precede or follow it. Hence, the vertices of
the polygon can be of type ‘left’ or ‘right’, and also of type ‘top’ or ‘bottom’; a simple
test of the neighbours suffices to label them, see Figure 12.5.
Algorithm
The algorithm, based on the sweep line technique, is of complexity O(n log n).
left-top
right-bottom
left-bottom
The algorithm will sweep the vertices of the polygon in lexicographical order of the
coordinates, and maintain in a structure S the set of y-coordinates of the left vertices
the right neighbours of which have not yet been visited. Initially, S is empty, and for
each vertex met, we process it as follows:
y-coordinate y crosses the current vertical segment (x,y ) − (x,y), and hence once
again the polygon is not simple.
Complexity
For good performance, the following operations on S must be performed efficiently:
Implementation Details
Rather than working with the coordinates of the vertices, which are not bounded, we
consider their rank. Let y0 < · · · , < yk be the list of distinct y-coordinates of vertices
(k ≤ n/2 since the polygon is rectilinear). Hence, yi ∈ S if and only if t[i] = −1. To
determine if S contains an element yj in the interval [yi ,yk ], it suffices to determine if
the minimum of t between t[i + 1] and t[k − 1] is −1.
def is_simple(polygon):
n = len(polygon)
order = list(range(n))
order.sort(key=lambda i: polygon[i]) # lexicographic order
rank_to_y = list(set(p[1] for p in polygon))
rank_to_y.sort()
y_to_rank = {y: rank for rank, y in enumerate(rank_to_y)}
S = RangeMinQuery([0] * len(rank_to_y)) # sweep structure
last_y = None
for i in order:
x, y = polygon[i]
rank = y_to_rank[y]
# -- type of point
right_x = max(polygon[i - 1][0], polygon[(i + 1) % n][0])
left = x < right_x
below_y = min(polygon[i - 1][1], polygon[(i + 1) % n][1])
high = y > below_y
if left: # y does not need to be in S yet
if S[rank]:
return False # two horizontal segments intersect
S[rank] = -1 # add y to S
else:
S[rank] = 0 # remove y from S
if high:
lo = y_to_rank[below_y] # check S between [lo + 1, rank - 1]
if (below_y != last_y or last_y == y or
rank - lo >= 2 and S.range_min(lo + 1, rank)):
return False # horiz. & vert. segments intersect
last_y = y # remember for next iteration
return True
13 Rectangles
Numerous problems involving geometrical figures are concerned with rectangles, for
example when they arise from drawing house plans or from computer graphics. Often,
the rectangles are rectilinear, i.e. with their sides parallel to the axes, thus easing their
handling. An important algorithmic technique in geometry is the sweep line, already
described in Section 13.5 on page 205.
Definition
Given a set S of n points in the plane, we would like to determine all of the rectangles
with their four corners in S. These rectangles do not necessarily have to be rectilinear.
Algorithm
The complexity is O(n2 + m) where m is the size of the solution. The key observa-
tion is that the two pairs of diagonally opposed corners of a rectangle have a com-
mon signature. This consists of their midpoint and the distance between them. It
then suffices to store in a dictionary all pairs of points p,q with key (c,d) where
c = (p + q)/2 is the midpoint between p and q and d = |q − p| is the distance.
The couples having the same signature are those forming the rectangles of S, see
Figure 13.1.
Figure 13.1 A common signature between pairs of opposing corners: the diagonals intersect at
their midpoint and are of the same length.
13.2 Largest Square in a Grid 201
Implementation Details
Often in programming contests the input consists of integers. For improved efficiency,
one should in these cases avoid the use of floating point numbers whenever possible.
Hence, in this implementation, when calculating c, the division by 2 is skipped, as is
the square root for the calculation of d.
def rectangles_from_points(S):
answ = 0
pairs = {}
for j, _ in enumerate(S):
for i in range(j): # loop over point pairs (p,q)
px, py = S[i]
qx, qy = S[j]
center = (px + qx, py + qy)
dist = (px - qx) ** 2 + (py - qy) ** 2
signature = (center, dist)
if signature in pairs:
answ += len(pairs[signature])
pairs[signature].append((i, j))
else:
pairs[signature] = [(i, j)]
return answ
Problem
Rectangle [spoj:HMSRECT]
Definition
Given a black and white image in the form of an n × m array of pixels, we would like
to determine the largest completely black square, see Figure 13.2.
(i,j)
Figure 13.2 The black square of maximal size k having for bottom right-hand corner (i,j )
contains three black squares of size k − 1 finishing on the corners (i,j − 1), (i − 1,j − 1) and
(i − 1,j ).
202 Rectangles
A[i,j ]
0 if the cell (i,j ) is white,
=
1 + min{A[i − 1,j ],A[i − 1,j − 1],A[i,j − 1]} otherwise.
Problem
Making Chess Boards [spoj:CT101CC]
Definition
Given a histogram in the form of an array of non-negative integers x0, . . . ,xn−1 , the
goal is to place a rectangle of maximal area within this histogram. This means finding
an interval [,r) that maximises the area (r − ) × h, where the height h is h =
min≤i<r xi .
Application
At the bottom of the Atlantic Ocean lies a telecommunications cable linking Europe
and North America. The technical characteristics of this cable vary with time, as
a function of the movements of the ocean, fluctuations of temperature, etc. Hence,
at each moment there is a maximal bandwidth possible for transmission defining a
function max bandwidth vs. time. Changing the transmission bandwidth for a connec-
tion is possible but requires an exchange between the terminals that prohibits any
communication for a short time. Imagine that you know in advance the maximal
bandwidths possibles for each of the 60 × 24 minutes of a day. You would like to find
an interval of time and a bandwidth allowing the transmission of a maximal quantity
of information without a disconnection. This problem reduces to finding a rectangle
with maximal area within a histogram.
13.3 Largest Rectangle in a Histogram 203
10
12
15
16
Figure 13.3 As long as the histogram is increasing, the corners are pushed onto the stack.
When the histogram decreases, the higher corners must be popped, the corresponding
candidate rectangles tested and then the last popped corner is pushed again but with the
smaller height.
def rectangles_from_histogram(H):
best = (float(’-inf’), 0, 0, 0)
S = []
H2 = H + [float(’-inf’)] # extra element to empty the queue
for right, _ in enumerate(H2):
x = H2[right]
left = right
while len(S) > 0 and S[-1][1] >= x:
left, height = S.pop()
# first element is area of candidate
rect = (height * (right - left), left, height, right)
if rect > best:
best = rect
S.append((left, x))
return best
204 Rectangles
Problems
Largest Rectangle in a Histogram [spoj:HISTOGRA]
Galerie d’art [prologin:2012:galerie]
Application
Given a building site with trees scattered about the lot, we wish to find a maximal
rectangular area on which to build a house without cutting any trees.
Definition
Given a black and white image in the form of an array of n × m pixels, we want to
determine the largest completely black rectangle. Here, a rectangle is the intersection
between a set of contiguous rows and a set of contiguous columns, see Figure 13.4.
We now introduce a simple problem for which a plethora of interesting solutions exist.
Definition
Given n rectilinear rectangles, we would like to compute the area of their union. The
same techniques could be used to compute the perimeter or the number of connected
regions.
Key Observation
A first step is to make the complexity of our algorithms independent of the coordinates
of the rectangles: these coordinates span at most 2n points on each axis. Therefore,
they lie on a grid made up of O(n2 ) cells, which are either entirely covered by a
rectangle or disjoint from every rectangle, see Figure 13.5.
Figure 13.5 Grid with O(n) rows and O(n) columns. Each cell is either included in the union of
rectangles or is disjoint from it.
206 Rectangles
Algorithm in O(n3 )
A first simple solution consists then of determining a Boolean array indicating for
each cell if it is contained in the union of rectangles. The area of the union can then
be computed by a simple sum of the areas of the cells. The work performed for each
rectangle is O(n2 ), which leads to the announced complexity of O(n3 ).
def union_rectangles_naive(R):
X = set() # set of all x coordinates in the input
Y = set() # same for y
for x1, y1, x2, y2 in R:
assert x1 <= x2 and y1 <= y2
X.add(x1)
X.add(x2)
Y.add(y1)
Y.add(y2)
j_to_x = list(sorted(X))
i_to_y = list(sorted(Y))
# X and Y partition space into a grid
area = 0
for j in range(len(j_to_x) - 1): # loop over columns in grid
x1 = j_to_x[j]
x2 = j_to_x[j + 1]
for i in range(len(i_to_y) - 1): # loop over rows
y1 = i_to_y[i] # (x1,...,y2) is the grid cell
y2 = i_to_y[i + 1]
if rectangles_contains_point(R, x1, y1):
area += (y2 - y1) * (x2 - x1) # cell is covered
return area
def union_rectangles(R):
events = []
for x1, y1, x2, y2 in R: # initialize events
assert x1 <= x2 and y1 <= y2
events.append((y1, OPENING, x1, x2))
events.append((y2, CLOSING, x1, x2))
current_intervals = Counter()
area = 0
previous_y = 0 # arbitrary initial value,
# ok, because union_intervals is 0 at first event
for y, offset, x1, x2 in sorted(events): # sweep top down
area += (y - previous_y) * union_intervals(current_intervals)
previous_y = y
current_intervals[x1, x2] += offset
return area
Algorithm in O(n2 )
There is an easy trick to shave off the logarithmic factor. The bottleneck in the func-
tion union_intervals is sorting the endpoints of the given intervals. To avoid this
operation, we consider the set of all x-coordinates in the input (left or right borders of
the rectangles) and use it to divide the horizontal axis into O(n) elementary segments.
We maintain an array of integers indicating for each segment how many rectangles
cover it. Then adding or removing a new interval takes O(n) time, as we must modify
the counters of segments covered by the current interval. Computing the total length
spanned by their union also takes O(n), so the total complexity is O(n2 ).
13.5 Union of Rectangles 209
def union_rectangles_fast(R):
X = set() # set of all x coordinates in the input
events = [] # events for the sweep line
for x1, y1, x2, y2 in R:
assert x1 <= x2 and y1 <= y2
X.add(x1)
X.add(x2)
events.append((y1, OPENING, x1, x2))
events.append((y2, CLOSING, x1, x2))
# array of x coordinates in left to right order
i_to_x = list(sorted(X))
# inverse dictionary maps x coordinate to its rank
x_to_i = {xi: i for i, xi in enumerate(i_to_x)}
# nb_current_rectangles[i] = number of rectangles intersected
# by the sweepline in interval [i_to_x[i], i_to_x[i + 1]]
nb_current_rectangles = [0] * (len(i_to_x) - 1)
area = 0
length_union_intervals = 0
previous_y = 0 # arbitrary initial value,
# because length is 0 at first iteration
for y, offset, x1, x2 in sorted(events):
area += (y - previous_y) * length_union_intervals
i1 = x_to_i[x1]
i2 = x_to_i[x2] # update nb_current_rectangles
for j in range(i1, i2):
length_interval = i_to_x[j + 1] - i_to_x[j]
if nb_current_rectangles[j] == 0:
length_union_intervals += length_interval
nb_current_rectangles[j] += offset
if nb_current_rectangles[j] == 0:
length_union_intervals -= length_interval
previous_y = y
return area
The operation change is called whenever the sweep line encounters the bottom side
of a rectangle (d = 1) or its top side (d = −1). The method cover() returns the total
length spanned by the current intervals.
210 Rectangles
The data structure consists of a binary tree, where each node p is responsible for a
range of indices I in t (and hence also in L). It possesses three attributes:
The root of the tree is node 1 and the result of cover() is stored in s[1]. The array t
is stored implicitly in the attribute c. The value t[j ] is the sum of c[p] over all of the
nodes p along the path from the root to the leaf node corresponding to j . As is the case
for the structure for the minimum over a range (minimum range query), the update can
be done in logarithmic time. See Figure 13.6 on the next page for an example of usage
of such a segment tree.
The complexity of O(n log n) is justified by the required sort of the events.
class CoverQuery:
def __init__(self, L):
assert L != [] # L is assumed sorted
self.N = 1
while self.N < len(L):
self.N *= 2
self.c = [0] * (2 * self.N) # --- covered
self.s = [0] * (2 * self.N) # --- score
self.w = [0] * (2 * self.N) # --- length
for i, _ in enumerate(L):
self.w[self.N + i] = L[i]
for p in range(self.N - 1, 0, -1):
self.w[p] = self.w[2 * p] + self.w[2 * p + 1]
def cover(self):
return self.s[1]
0/13/8
0/6/5 0/7/3
c/w/s
1/1/1
0/1/0
0/1/0
0/1/0
0/2/0 1/2/2 0/0/2 1/3/3
t 1 2 1 0 0 0 1 0
sweep line
def union_rectangles_fastest(R):
if R == []: # segment tree would fail on an empty list
return 0
X = set() # set of all x coordinates in the input
events = [] # events for the sweep line
for Rj in R:
(x1, y1, x2, y2) = Rj
assert x1 <= x2 and y1 <= y2
X.add(x1)
X.add(x2)
events.append((y1, OPENING, x1, x2))
events.append((y2, CLOSING, x1, x2))
i_to_x = list(sorted(X))
# inverse dictionary
x_to_i = {i_to_x[i]: i for i in range(len(i_to_x))}
L = [i_to_x[i + 1] - i_to_x[i] for i in range(len(i_to_x) - 1)]
C = CoverQuery(L)
area = 0
previous_y = 0 # arbitrary initial value,
# because C.cover() is 0 at first iteration
for y, offset, x1, x2 in sorted(events):
area += (y - previous_y) * C.cover()
i1 = x_to_i[x1]
i2 = x_to_i[x2]
C.change(i1, i2, offset)
previous_y = y
return area
212 Rectangles
Definition
Given n disjoint rectilinear rectangles, we wish to determine all pairs of adjacent
rectangles.
Application
This can serve to compute the perimeter of the union, for example, by subtracting
from the total of the perimeters of the rectangles the length of the contacts between
all the adjacent rectangles. This last value must be taken with a factor of 2, as a
unit of contact removes a unit of perimeter from each rectangle. Another application
is to determine the connected components formed by the adjacent rectangles, see
Figure 13.7.
2 3
5
4
6
Figure 13.7 Detect the adjacency between the rectangles (2,4),(2,5),(3,5) and (4,5).
Variant
If in the problem statement, we consider two rectangles to be adjacent if even only a
corner is shared, then it suffices to change the priorities and, for two identical points,
first process the bottom corners and then the top corners.
Problem
City Park [icpcarchive:6889]
14 Numbers and Matrices
14.1 GCD
Definition
Given two integers a,b we seek the largest integer p such that a and b can be expressed
as integer multiples of p; this is their greatest common divisor (GCD).
The calculation of the GCD can be implemented recursively to be very efficient.
A mnemonic trick: from the second iteration on, we arrange to always have the second
argument smaller than the first: a mod b < b.
Definition
For two integers a and b, we would like to find two integers u and v such that au+bv =
d where d is the GCD of a and b.
This calculation is based on an observation similar to the above. If a = qb + r, then
au + bv = d corresponds to (qb + r)u + bv = d, or bu + rv = d for
u = qu + v u = v
⇔ .
v =u v = u − qv
Variant
Certain problems involve the calculation of extremely large numbers, and conse-
quently require a response modulo a large prime number p in order to test if the
solution is correct. Since p is prime, we can easily divide by an integer a non-multiple
of p: a and p are relatively prime, hence their Bézout coefficients satisfy au+pv = 1,
hence au is equal to 1 modulo p and u is the inverse of a. Hence, to divide by a, we
can instead multiply by u (mod p).
An alternative is to calculate
Pascal’s triangle by dynamic programming, which
could be interesting if nk needs to be calculated for numerous pairs n,k.
216 Numbers and Matrices
Definition
Given a,b we wish to calculate a b . Again, as the result risks being enormous, we are
often asked to perform the calculation modulo a given integer p, but this does not
change the nature of the problem.
Algorithm in O(log b)
The naive approach performs b − 1 multiplications by a. However, we can rapidly cal-
k k k+1
culate the powers of a of the form a 1,a 2,a 4,a 8, . . . using the relation a 2 ·a 2 = a 2 .
A trick consists of combining these values according to the binary decomposition of b.
Example:
a 13 = a 8+4+1
= a8 · a4 · a1.
Variant
This technique can also be applied to matrix multiplications. Let A be a matrix and b a
positive integer. Rapid exponentiation allows the calculation of Ab with only O(log b)
matrix multiplications.
Problems
The last digit [spoj:LASTDIG]
Tiling a grid with dominoes [spoj:GNY07H]
14.5 Prime Numbers 217
Implementation Details
The proposed implementation skips marking 0, 1 as well as the multiples of 2. This is
why the iteration over p is done with a step of 2, in order to test only the odd numbers.
def eratosthene(n):
P = [True] * n
answ = [2]
for i in range(3, n, 2):
if P[i]:
answ.append(i)
for j in range(2 * i, n, i):
P[j] = False
return answ
y = factor[y] · x
plays the role of factor[y] in the expression y = p · x. The algorithm is correct, since
at the time of processing x, all prime numbers less than or equal to x have been found,
as well as the smallest factor of x. As every number y between 2 and n−1 is examined
exactly once, the complexity of the algorithm is O(n).
def gries_misra(n):
primes = []
factor = [0] * n
for x in range(2, n):
if not factor[x]: # no factor found
factor[x] = x # meaning x is prime
primes.append(x)
for p in primes: # loop over primes found so far
if p > factor[x] or p * x >= n:
break
factor[p * x] = p # p is the smallest factor of p * x
return primes, factor
Problems
The Spiral of Primes [icpcarchive:2120]
Prime or Not [spoj:PON]
Bazinga! [spoj:DCEPC505]
Definition
Given an expression obeying a certain fixed syntax, we would like to construct the
syntax (parse) tree or evaluate the value of the expression, see Figure 14.1.
2 /
* 2
3 4
Figure 14.1 The syntax tree associated with the expression 2 + (3 × 4)/2.
Approaches
In general, this problem is solved by decomposing the stream containing the expres-
sion into a series of lexical tokens, using a tool called a scanner or lexical analyzer.
14.6 Evaluate an Arithmetical Expression 219
Then, as a function of the lexical tokens identified, we construct the syntax tree accord-
ing to a grammar, using a tool known as a parser.
However, if the syntax is similar to that of arithmetical expressions, an evaluation
method using two stacks can be used, which is much simpler to implement. One stack
contains the values and the other the operators. Each value encountered is pushed
as-is onto the value stack. For each operator p encountered, before pushing it onto the
operator stack, the following must be performed: as long as the top element q of the
operator stack has priority at least as high as p, we pop q, as well as the last two a,b,
and then push onto the value stack the result of the expression a q b (the result of the
operation q on the operands a,b).
In this manner, the evaluation of operators is delayed until the priority of the oper-
ators forces the evaluation, see Figure 14.2.
A special case is needed to distinguish the binary minus operator from the unary
minus sign for integer constants. Usually, it suffices to test if the previous token was
an operator or a closing parenthesis.
Spreadsheet Example
Consider a spreadsheet where each cell can contain either a value or an arithmetical
expression. This expression can be composed of numerical constants and identifiers
of cells, linked by the operators - + * / and parentheses.
We represent an expression by an integer, a character string or a triplet composed
of an operator and two operands. The numerical evaluation of an expression is realised
by the following recursive function, where cell is a dictionary, associating the names
of cells to their contents.
Figure 14.2 Example of the evaluation of the expression 2 + 3 × 4/2 with the marker ‘ ; ’ as
‘end of expression’.
220 Numbers and Matrices
The syntax tree is constructed following the method described above. Note the
particular handling of the parentheses. A left parenthesis is always pushed onto
the operator stack without doing anything else. Meeting a right parenthesis provokes
the popping of the stack and the construction of the expressions, until the top of the
operator stack corresponds to the matching left parenthesis; this is in turn popped. To
completely empty the stacks at the end of the processing, we add to the stream a token
‘ ; ’ as end of expression, and assign it the minimal priority.
def arithm_expr_parse(line_tokens):
vals = []
ops = []
for tok in line_tokens + [’;’]:
if tok in PRIORITY: # tok is an operator
while (tok != ’(’ and ops and
PRIORITY[ops[-1]] >= PRIORITY[tok]):
right = vals.pop()
left = vals.pop()
vals.append((left, ops.pop(), right))
if tok == ’)’:
ops.pop() # this is the corresponding ’(’
else:
ops.append(tok)
elif tok.isdigit(): # tok is an integer
vals.append(int(tok))
else: # tok is an identifier
vals.append(tok)
return vals.pop()
14.7 System of Linear Equations 221
Beware of a Trap!
An error often committed during the implementation of this code is to write:
which has the undesirable effect of inverting the left and right values of the expression,
because of the order of evaluation of the arguments of the method append.
Problems
Boolean Logic [spoj:BOOLE]
Arithmetic Expressions [spoj:AREX]
Cells [spoj:IPCELLS]
Definition
A system of linear equations is composed of n variables and m linear equations.
Formally, a matrix A of dimension n × m is given, along with a column vector b
of size m. The goal is to find a vector x such that Ax = b.
the forces (including gravitational) cancel out for every ball. This is again a problem
of resolution of a system of linear equations.
Algorithm in O(n2 m)
If A were the identity matrix, then the solution of the system would be b. We are going
to diagonalise A to obtain a form as close as possible to this ideal situation.
For this, we apply transformations that preserve the solutions of the system, such as
changing the order of variables, exchanging two equations, multiplying one equation
by a constant and adding one equation to another.
In order to simplify these manipulations and not modify the parameters A and b,
we store the system of equations in a matrix S, which is composed of a copy of A
along with an extra column containing b and an extra row containing the indices of
the variables. Thus, when the columns are exchanged in S, it is possible to remember
to which variable each column corresponds.
The invariant is the following: after k iterations, the first k columns of S are all
zero, except on the diagonal where there are 1’s, see Figure 14.4.
To obtain the invariant, we divide the kth row of S—hence the kth equation—by
S[k,k] if it is non-zero, which brings a 1 to the index (k,k). Then we subtract the kth
equation from every equation i = k, but multiplied by the factor S[i,k], bringing a 0
to all of the indices (i,k) for every i = k, as required by the invariant.
14.7 System of Linear Equations 223
⎛ ⎞ ⎛ ⎞
1 0 0 · · · ·
⎜ 0 1 0 · · · ⎟ ⎜ · ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 0 1 · · · ⎟ ⎜ · ⎟
A=⎜
⎜
⎟ b=⎜
⎟ ⎜
⎟
⎟
⎜ 0 0 0 · · · ⎟ ⎜ · ⎟
⎝ 0 0 0 · · · ⎠ ⎝ · ⎠
0 0 0 · · · ·
Implementation Details
As we work here with floating point numbers, with rounding errors possible for each
calculation, it is necessary to permit a slight tolerance when testing equality to zero.
When subtracting the kth row from the ith row, we must assign the coefficient S[i][k]
to a variable fact, as the operation will change the value of S[i][k]. To calculate
the solution exactly, we could work with fractions instead of floating point numbers.
In this case, it is important to reduce the fractions at each iteration, otherwise the
numerators and denominators of the solution could contain an exponential number of
digits. Without the normalisation, the Gauss–Jordan algorithm is not polynomial, but
happily the Python class Fraction reduces the fractions at each operation.
Variants
If the matrix is sparse, i.e. it contains only a small number of non-zero elements per
row and column, then the execution time can be reduced from O(n2 m) to O(n).
224 Numbers and Matrices
GJ_ZERO_SOLUTIONS = 0
GJ_SINGLE_SOLUTION = 1
GJ_SEVERAL_SOLUTIONS = 2
Problems
Ars Longa [icpcarchive:3563]
Linear Equation Solver [spoj:LINQSOLV]
14.8 Multiplication of a Matrix Sequence 225
Figure 14.5 What parentheses should be added to multiply a sequence of matrices with the least
total number of operations?
Definition
Given n matrices M1, . . . ,Mn , where the ith matrix has ri rows and ci columns, with
ci = ri+1 for every 1 ≤ i < n. We would like to calculate the product M1 M2 · · · Mn
with as few operations as possible. By associativity, it is possible to place the paren-
theses in different ways. We have the following equality by associativity of the matrix
product, but the number of operations could differ, based on the sizes of the various
sub-products:
(((M1 M2 )M3 )M4 = M1 (M2 (M3 (M4 ))) = (M1 M2 )(M3 M4 ).
To multiply two matrices Mi and Mi+1 , we use the standard algorithm which executes
ri ci ci+1 numerical multiplications. The goal is to find how to place the parentheses so
as to multiply with least cost, see Figure 14.5.
Algorithm in O(n2 )
The recurrence formula1 is simple. The last multiplication multiplies the result of
M1 · · · Mk with the result of Mk+1 · · · Mn for a certain 1 ≤ k < n. Let opt(i,j ) be the
minimal cost for calculating Mi · · · Mj . Then opt(i,i) = 0 and for i < j
opt(i,j ) = min (opt(i,k) + opt(k + 1,j ) + ri ck cj . (14.1)
i≤k<j
If we wish to compute both the cost of the optimal order and the order itself, we
must store side by side with the matrix opt a matrix arg, containing the index k which
realises the minimisation (14.1). This is done in the following implementation. The
1 The complexity concerns only the computation of the optimal placement of parentheses, and not the
operation of matrix multiplication itself.
226 Numbers and Matrices
def matrix_mult_opt_order(M):
n = len(M)
r = [len(Mi) for Mi in M]
c = [len(Mi[0]) for Mi in M]
opt = [[0 for j in range(n)] for i in range(n)]
arg = [[None for j in range(n)] for i in range(n)]
for j_i in range(1, n): # loop on i, j of increasing j - i = j_i
for i in range(n - j_i):
j = i + j_i
opt[i][j] = float(’inf’)
for k in range(i, j):
alt = opt[i][k] + opt[k + 1][j] + r[i] * c[k] * c[j]
if opt[i][j] > alt:
opt[i][j] = alt
arg[i][j] = k
return opt, arg
def matrix_chain_mult(M):
opt, arg = matrix_mult_opt_order(M)
return _apply_order(M, arg, 0, len(M)-1)
An even better algorithm exists in O(n log n), see (Hu and Shing, 1984), but we do
not present it here.
Problems
Mixtures [spoj:Mixtures]
The Safe Secret [kattis:safesecret]
Sweet and Sour Rock [spoj:ROCK]
15 Exhaustive Search
For some combinatorial problems, no polynomial time algorithm is known. The usual
last-gasp (desperate!) solution is then to exhaustively explore the space of poten-
tial solutions. The word combinatorial here refers to structures which are combined
from simpler structures. Examples include the construction of trees by a combination
of subtrees, or the assembly of tilings with small tiles. Exhaustive exploration then
consists of exploring the implicit tree of possible constructions in order to detect a
solution. The nodes of the tree represent partial constructions, and when they can no
longer be extended to a solution, for example by the violation of some constraint, the
exploration backs up and examines alternate branches. This technique is known as
backtracking.
We will illustrate these principles with a simple example.
Definition
Suppose there is a rectangular grid, surrounded by a border with two openings in
the top row, one on the left and one on the right. Certain cells of the grid contain
double-faced mirrors, which can be oriented in two possible directions, diagonal or
anti-diagonal, see Figure 15.1. The goal is to orient the mirrors in such a way that
a laser beam entering by the left opening exits by the right opening. The laser beam
traverses the empty cells horizontally or vertically, but when it hits a mirror, it is
deflected by 90 degrees to the right or to the left, depending on the orientation of the
mirror. If the beam hits the boundary of the grid, it is absorbed.
Figure 15.1 An orientation of the mirrors permitting the laser beam to traverse the grid.
vertices (a,d) and (b, − d) are present in the graph and linked by an edge, where −d
denotes the opposite direction to d. However, if the direction d from a mirror a leads
to the wall, then there is no vertex (a,d) present in the graph. The construction is
completed with a vertex per opening, linked to the first accessible mirror, if such exists.
Figure 15.2 Reduction of the laser problem (left) to a perfect matching problem (right). The
matching edges are shown with solid lines.
We can verify that every perfect matching in this graph corresponds to a solution
of the laser problem and vice-versa. We distinguish the rectilinear edges linking the
mirrors, and the diagonal edges linking different vertices associated with the same
mirror. The interpretation is as follows. A rectilinear edge in the matching indicates
that the beam will not follow this path, whereas a diagonal edge, say between vertices
(a,d) and (a,e), indicates that the beam can be reflected from the mirror a by the
directions d et e.
The graph thus produced is not bipartite. To compute a maximal matching we must
thus resort to Jack Edmonds’ ‘flowers and petals’ algorithm (1965), with a complexity
of O(n2 ), where n is the number of mirrors. However, this is quite long and difficult
15.1 All Paths for a Laser 229
Algorithm
We propose a solution by exhaustive exploration. For each mirror, we store a state
that can be one of three types: the two orientations and a type ‘without orientation’.
Initially, the mirrors are all without orientation. Then the route of the laser beam is
simulated, from its entrance by the left opening.
• When an oriented mirror is hit, the beam is reflected as a function of the orientation.
• When the beam hits a mirror without orientation, two recursive calls are made, one
for each possible orientation of the mirror. If one of these calls finds a solution, it
is returned. If neither call finds a solution, the mirror is returned to the unoriented
position and a failure code is returned.
• If the beam hits the boundary, the recursive procedure terminates with a failure
code. The exploration then backtracks and considers other possible orientations.
• Finally, if ever the opening on the top right is reached, the solution is returned.
Implementation Details
The n mirrors whose positions are given as input are numbered from 0 to n − 1. Two
dummy mirrors are placed at the openings and numbered n and n + 1. The program
produces an array L with the coordinates and the indices of the mirrors.
The four possible directions of the beam are represented by the integers from 0 to 3,
and the two orientations by 0,1. An array reflex of dimension 4 × 2 indicates the
change of direction made by the beam if it arrived in a given direction on a mirror
with a given orientation.
In a pre-computation, the successive mirrors of a same row or a same column are
linked together, with the help of the array next_. For a mirror i and a direction d, the
element next_[i][d] indicates the index of the next mirror the beam will hit when it
leaves i in the direction d. This entry is set to None if the beam is sent to the boundary.
To populate the array next_, the array L is scanned first of all row by row and then
column by column. Note the sort function that reverses the indices row and column
for the lexicographic sort.
The variables last_i, last_r, last_c contain the information concerning the
last mirror seen in the scan. If the last mirror is in the same row as the current mirror
(for the scan row by row), then the two must be linked by the indications of their
respective references in the array next_.
230 Exhaustive Search
# directions
UP = 0
LEFT = 1
DOWN = 2
RIGHT = 3
# orientations None:? 0:/ 1:\
Definition
The problem of exact cover concerns a set of points U , called the universe, and a
collection of subsets of U , S ⊆ 2U , see Figure 15.3. A set A ⊆ U is said to cover
x ∈ U if x ∈ A. The goal is to find a selection of sets of S, i.e. a collection S ∗ ⊆ S,
covering each element of the universe exactly once.
In terms of a characteristic matrix, the problem can be described as follows. On
input, we are given a binary matrix M, the columns of which represent the elements
of the universe, and the rows represent the sets of the collection S. The element "x,A#
of the matrix is 1 if and only if x ∈ A. On output, we must generate a set of rows S ∗
such that the matrix restricted to S ∗ contains exactly one 1 in each column.
Applications
The game Sudoku can be seen as an exact cover problem, see Section 15.4 on
page 238. This is also the case for the tiling problem, wherein a set of tiles must be
232 Exhaustive Search
used to cover an m × n grid without overlap. Each tile must be used exactly once, and
each cell of the grid must be covered by exactly one tile. The cells and the tiles thus
form the elements of the universe, whereas the placements of individual tiles constitute
the sets.
1 This choice reduces the number of branches at this point, and hopefully the size of the search tree.
15.2 The Exact Cover Problem 233
S2
0 1 2
S3
3
set S1 S4
S1
6 4 5
S0
S5
S2
S2
0 1 2
0 1 2 S3
S3
3
3
item 3 set S3 S4
S4
S1
6 4 5
S1
6 4 5
S0
S0 S5
S5
S2
0 1 2
S3
3
set S5 S4
S1
6 4 5
S0
S5
Figure 15.3 Three subproblems resulting from the choice of the element e = 3. Note that the
third subproblem does not have a solution, as the element 0 cannot be covered.
For each column, we also require a header cell, which is included in the vertical
linking, in order to gain access to the column. At the initialisation of the structure, the
headers are stored in an array col indexed by the column numbers. Finally, we have a
particular cell h, included in the horizontal linking of the column headers, in order to
access them. This cell does not use the fields U,D.
Each cell possesses two additional fields S, C, the role of which depends on the
type of cell. For the matrix cells, S contains the row number and C contains the column
header. For the column header cells, S contains the number of 1’s in the column, and
the field C is ignored. The cell h ignores both of these fields.
234 Exhaustive Search
h 0 1 2 3 4 5 6
2 2 2 3 2 2 3
3,0 3,3
4,1 4,6
⎛ ⎞
0 0 1 0 1 1 0
⎜ 1 0 0 1 0 0 1 ⎟
⎜ ⎟
⎜ 0 1 1 0 0 1 0 ⎟
⎜ ⎟
⎜ 1 0 0 1 0 0 0 ⎟
⎜ ⎟
⎝ 0 1 0 0 0 0 1 ⎠
0 0 0 1 1 0 1
h 0 1 2 3 4 5 6
2 2 2 1 2 2 2
3,0 3,3
4,1 4,6
Figure 15.4 A binary matrix M, above its encoding and below the result of covering the
column 0. The lists are circular, and the links leaving an edge of the image re-enter by the
opposite edge.
15.2 The Exact Cover Problem 235
class Cell:
def __init__(self, horiz, verti, S, C):
self.S = S
self.C = C
if horiz:
self.L = horiz.L
self.R = horiz
self.L.R = self
self.R.L = self
else:
self.L = self
self.R = self
if verti:
self.U = verti.U
self.D = verti
self.U.D = self
self.D.U = self
else:
self.U = self
self.D = self
def hide_verti(self):
self.U.D = self.D
self.D.U = self.U
def unhide_verti(self):
self.D.U = self
self.U.D = self
def hide_horiz(self):
self.L.R = self.R
self.R.L = self.L
def unhide_horiz(self):
self.R.L = self
self.L.R = self
The Links
The idea of the algorithm of dancing links is due to Hitotsumatsu and Noshita in 1979
and was described by Donald Knuth (2000). To remove element c from a doubly
linked list, it suffices to change the targets of the pointers of its neighbours, see
Figure 15.5. To reinsert it in the list, it suffices to perform the reverse operations in
the opposite order.
We use this observation to easily remove and restore a column of the matrix.
The operation of covering a column c consists of removing it from the horizontal
list of column headers, and for each row r with Mrc = 1, removing from the
vertical lists all the cells corresponding to positions (r,c ) in M with Mrc = 1.
We must make sure to keep up to date the counter S in the header of the cell c ,
decrementing it.
236 Exhaustive Search
Figure 15.5 The operation hide removes element c from a doubly linked list. By preserving the
pointers of c, it is then easy to reinsert it at its initial position.
def uncover(c):
assert c.C is None
i = c.U
while i != c:
j = i.L
while j != i:
j.C.S += 1 # one more entry in this column
j.unhide_verti()
j = j.L
i = i.U
c.unhide_horiz()
The Search
In the search procedure, we simply scan all the columns to find a column minimising
the number of elements, hence with minimal counter S. We could accelerate this
search by using a priority queue for the columns, but this would make the procedure
of covering a column more costly. The following function writes the solution in an
array sol if the function returns True.
15.3 Problems 237
15.3 Problems
Sudoku [spoj:SUDOKU]
Making Jumps [spoj:MKJUMPS]
238 Exhaustive Search
15.4 Sudoku
Several naive methods succeed in efficiently solving classic instances of Sudoku (see
Figure 15.6). However, for Sudoku grids of size 16 × 16, the dancing links algorithm
seems to be the right tool.
7 2
5 9
3 8
4 5
3 8
2 3
2 5
6 3
1 9
Figure 15.6 A Sudoku grid. The goal is to fill in the free squares so that each column, each row
and each 3 × 3 block contains all the integers from 1 to 9.
Modelisation
How can we model the problem of Sudoku as an exact cover problem? There are
four types of constraints. Each cell must contain a value. In each row of the Sudoku
grid, every value must appear exactly once. The same holds for the columns and the
blocks. There are thus four types of elements in the universe: the couples row-column,
the couples row-value, the couples column-value and the couples block-value. These
elements constitute the universe U of the instance "U,S# of the exact cover problem.
Then, the sets of S consist of assignments: triplets row-column-value, with each
assignment covering exactly four elements of the universe.
The pleasant aspect of this representation is that it makes abstraction of the notions
of rows, columns, blocks or values, giving a problem much easier to manipulate.
Sudoku grids typically come filled with a few predefined values. How can those be
encoded? Our method consists of adding a new element e to the universe and a new
set A to S that is the unique set containing e. Every solution must then include A. It
then suffices to place in A all elements of the universe already covered by the initial
assignments.
Encoding
The sets of the instance "U,S# of the exact cover correspond to assignments, and
our implementation of the dancing links algorithm returns the solution in the form
of an array with the indices of the selected sets. Hence, to recover the assignments
corresponding to the indices, we must specify an encoding. We chose to encode the
assignment of a value v to the square in row r and column c by the integer 81r +9c +v
(replace the coefficients by 256 and 16 for the case of a 16 × 16 Sudoku).
Similarly, the elements of the universe are encoded by integers. For example, the
couple row-column r,c is encoded by 9r + c, the couple row-value r,v by 81 + 9r + v
and so on.
15.4 Sudoku 239
N = 3 # global constants
N2 = N * N
N4 = N2 * N2
# sets
def assignment(r, c, v): return r * N4 + c * N2 + v
def sudoku(G):
global N, N2, N4
if len(G) == 16: # for a 16 x 16 sudoku grid
N, N2, N4 = 4, 16, 256
e = N * N4
universe = e + 1
S = [[rc(a), rv(a), cv(a), bv(a)] for a in range(N4 * N2)]
A = [e]
for r in range(N2):
for c in range(N2):
if G[r][c] != 0:
a = assignment(r, c, G[r][c] - 1)
A += S[a]
sol = dancing_links(universe, S + [A])
if sol:
for a in sol:
if a < len(S):
G[row(a)][col(a)] = val(a) + 1
return True
return False
240 Exhaustive Search
Problems
Easy sudoku [spoj:EASUDOKU]
Sudoku [spoj:SUDOKU]
Applications
Certain problems, lacking structure, require a naive solution by testing one by one
each element in the space of potential solutions. It is thus sometimes necessary to
explore all permutations of a given array.
S E N D
+ M O R E
= M O N E Y
where a distinct digit has to be assigned to each letter, in such a way that each word
becomes a number without leading zero and the equation of the addition is correct.
To solve this problem by naive enumeration, it suffices to construct an array tab =
“@@DEMNORSY” composed of the letters of the problem, completed with enough @ as
necessary to be of size 10. Then, there is a correspondence between the permutations
of this array and the assignments of distinct digits to the letters, by assigning to each
letter its rank in the array.
The interesting part is then the enumeration of all permutations of an array, and this
is the subject of this section.
Definition
Given an array t of n elements, we want to determine the next permutation of t, in
lexicographical order, or to signal that t is already maximal.
Key Observation
To permute t into the next lexicographical permutation, we need to preserve the
longest possible prefix and exchange only elements in the suffix.
initial array 0 2 1 6 5 2 1
choice of pivot 0 2 [1] 6 5 2 1
swap 0 2 [2] 6 5 [1] 1
reverse suffix 0 2 2 [1 1 5 6]
final array 0 2 2 1 1 5 6
Sorting the suffix in increasing order essentially comes down to reversing its ele-
ments since, initially, it was in decreasing order.
def next_permutation(tab):
n = len(tab)
pivot = None # find pivot
for i in range(n - 1):
if tab[i] < tab[i + 1]:
pivot = i
if pivot is None: # tab is already the last perm.
return False
for i in range(pivot + 1, n): # find the element to swap
if tab[i] > tab[pivot]:
swap = i
tab[swap], tab[pivot] = tab[pivot], tab[swap]
i = pivot + 1
j = n - 1 # invert suffix
while i < j:
tab[i], tab[j] = tab[j], tab[i]
i += 1
j -= 1
return True
242 Exhaustive Search
The solution to the problem of word additions can then be coded as follows.
Problems
The Next Permutation [spoj:NEXT]
Great Swerc [spoj:SWERC14A]
15.6 Le Compte est Bon 243
This problem originated in a television game show ‘Des chiffres et des lettres’, broad-
cast in France since 1972. The UK version was known as ‘Countdown’ (1982).
Input
Integers x0, . . . ,xn−1 and b for n quite small, say n ≤ 20.
Output
An arithmetical expression as close as possible to b, using each of the integers at
most once, and using the operations +, -, *, / as often as necessary. Subtraction is
allowed only if its result is positive, and division is allowed only if the result is without
a remainder.
Algorithm
The algorithm uses exhaustive search and dynamic programming. In a dictionary E,
we associate2 with S ⊆ {0, . . . ,n − 1} all the values obtained by using each input xi
for i ∈ S at most once. More precisely, E[S] is in its turn a dictionary associating each
obtained value with the corresponding expression.
For example, for x = (3,4,1,8) and S = {0,1}, the dictionary E[S] contains the
pairs composed of a key x and a value e, where e is an expression formed from the
inputs x0 = 3,x1 = 4 and x the value of e. Hence, E[s] contains the mapping pairs
1 → 4 − 3,3 → 3,4 → 4,7 → 3 + 4,12 → 3 ∗ 4.
To compute E[S], we loop over all of the partitions of S into two non-empty sets L
and R. For each value vL in E[L] obtained by an expression eL and each value vR in
E[R] obtained by an expression eR , we can form new values to be stored in E[S]. In
particular, vL + vR is attained by the expression eL + eR .
The complexity in the worst case is worse than exponential and is delicate to
analyse. First of all, we must count the number of binary trees with at most n
leaves, and then multiply by the number of possibilities to assign these operations
to internal nodes of the trees. Nonetheless, with the restrictions on subtraction and
division, in practice for values of n in the tens, the observed complexity remains
acceptable.
Implementation Details
Particular attention must be paid to the order in which the sets S are processed. By
respecting the increasing order of their size, we can ensure that the sets E[L] and
E[R] have already been determined.
In conclusion, here a few final remarks before we let you loose into the wilderness.
At times, the problem may require you to find an optimal value. Can this value be
found by dichotomy? If it lies between 0 and N and if the condition verifies a certain
property, this adds only a factor log N to the complexity. This is why it is always prof-
itable to know how to efficiently code a binary search, see Section 1.6.7 on page 35.
Try your skills on this problem, one of the most difficult:
• Flow algorithms are studied from top to bottom in Network Flows: Theory, Algo-
rithms, and Applications by R. K. Ahuja, T. L. Magnanti and J. B. Orlin, Pearson,
2013.
• A good introduction to algorithms for problems in geometry is Computational
Geometry: Algorithms and Applications by M. de Berg, O. Cheung, M. van Kreveld
and M. Overmars, Springer Verlag 2008.
• If you would like to learn more about tricky manipulations with the binary repre-
sentation of integers, the following is a delightful reference: Hacker’s Delight by
Henry S. Warren, Jr, Addison-Wesley, 2013.
• A very good introduction to programming in Python is Think Python: How to Think
Like a Computer Scientist by Allen Downey, O’Reilly Media, 2015.
• Other highly appreciated books include Python Essential Reference by David
M. Beazley, Pearson Education, 2009 and Python Cookbook by Brian K. Jones,
O’Reilly Media, 2013.
• Two texts more adapted to the preparation for programming contests are Com-
petitive Programming by Steven and Felix Halim, Lulu, 2013 and The Algorithm
Design Manual by Steven S. Skiena, Springer Verlag, 2009.
Throughout this work, the references point to the bibliography at the end of the
book. It includes books but also research articles, which are in general less accessible
and destined for a public of specialists. Many of these articles can be accessed freely
on the Internet, or found in a university Computer Science library.
This book is accompanied by a website, tryalgo.org, where the code of all the Python
programs described here can be found, as well as test data, training exercises and
solutions that, of course, must not be read before having attempted the problems. The
package tryalgo is available on GitHub and PyPI under the MIT licence, and can be
installed under Python 2 or 3 via pip install tryalgo.
Happy hacking!
Debugging tool
If you are stymied when trying to solve a problem, speak to this rubber duck. Explain
to it precisely and in detail your approach, go over and comment line-by-line your
code, and you are guaranteed to find the mistake.
References
Ahuja, Ravindra K., Magnanti, Thomas L., and Orlin, James B. 1993. Network Flows: Theory,
Algorithms, and Applications. Prentice Hall.
Alt, Helmut, Blum, Norbert, Mehlhorn, Kurt, and Paul, Markus. 1991. Computing a Maximum
√
Cardinality Matching in a Bipartite Graph in Time O(n1.5 m/ log n). Information Process-
ing Letters, 37(4), 237–240.
Andrew, Alex M. 1979. Another Efficient Algorithm for Convex Hulls in Two Dimensions.
Information Processing Letters, 9(5), 216–219.
Aspvall, Bengt, Plass, Michael F., and Tarjan, Robert Endre. 1979. A Linear-time Algorithm for
Testing the Truth of Certain Quantified Boolean Formulas. Information Processing Letters,
8(3), 121–123.
Bellman, Richard. 1958. On a Routing Problem. Quarterly of Applied Mathematics, 16, 87–90.
Cormen, Thomas H., Leiserson, Charles E., Rivest, Ronald L., and Stein, Clifford. 2010.
Algorithmique. Dunod.
Dilworth, R. P. 1950. A Decomposition Theorem for Partially Ordered Sets. Annals of
Mathematics, 51(1), 161–166.
Dinitz, Yefim. 2006. Dinitz’ Algorithm: The Original Version and Even’s Version. Berlin,
Heidelberg: Springer Berlin Heidelberg. Pages 218–240.
Edmonds, Jack. 1965. Paths, Trees, and Flowers. Canadian Journal of Mathematics, 17(3),
449–467.
Edmonds, Jack, and Johnson, Ellis L. 1973. Matching, Euler Tours and the Chinese Postman.
Mathematical Programming, 5(1), 88–124.
Edmonds, Jack, and Karp, Richard M. 1972. Theoretical Improvements in Algorithmic
Efficiency for Network Flow Problems. Journal of the ACM, 19(2), 248–264.
Elias, Peter, Feinstein, Amiel, and Shannon, Claude. 1956. A Note on the Maximum Flow
through a Network. IRE Transactions on Information Theory, 2(4), 117–119.
Fenwick, Peter M. 1994. A New Data Structure for Cumulative Frequency Tables. Software:
Practice and Experience, 24(3), 327–336.
Floyd, Robert W. 1962. Algorithm 97: Shortest Path. Communications of the ACM, 5(6), 345–.
Ford, Lester Randolph. 1956. Network Flow Theory. Tech. rept. P-923. RAND Corporation.
Ford, L. R., and Fulkerson, D. R. 1956. Maximal Flow through a Network. Canadian Journal
of Mathematics, 8, 399–404.
Fredman, Michael L., and Tarjan, Robert Endre. 1987. Fibonacci Heaps and Their Uses in
Improved Network Optimization Algorithms. Journal of the ACM, 34(3), 596–615.
Freivalds, Rūsin, š. 1979. Fast Probabilistic Algorithms. Pages 57–69 of: Mathematical Founda-
tions of Computer Science 1979. Springer.
Gabow, Harold N. 1976. An Efficient Implementation of Edmonds’ Algorithm for Maximum
Matching on Graphs. Journal of the ACM, 23(2), 221–234.
References 249
Gajentaan, Anka, and Overmars, Mark H. 1995. On a Class of O(n2 ) Problems in Computa-
tional Geometry. Computational Geometry, 5(3), 165–185.
Gale, David, and Shapley, Lloyd S. 1962. College Admissions and the Stability of Marriage.
American Mathematical Monthly, 9–15.
Goldberg, Andrew V., and Rao, Satish. 1998. Beyond the Flow Decomposition Barrier. Journal
of the ACM, 45(5), 783–797.
Gries, David, and Misra, Jayadev. 1978. A Linear Sieve Algorithm for Finding Prime Numbers.
Communications of the ACM, 21(12), 999–1003.
Hierholzer, Carl, and Wiener, Chr. 1873. Über die Möglichkeit, einen Linienzug ohne
Wiederholung und ohne Unterbrechung zu umfahren. Mathematische Annalen, 6(1),
30–32.
Hitotumatu, Hirosi, and Noshita, Kohei. 1979. A Technique for Implementing Backtrack
Algorithms and Its Application. Information Processing Letters, 8(4), 174–175.
Hopcroft, John, and Karp, Richard. 1972. An n5/2 Algorithm for Maximum Matchings in
Bipartite Graphs. SIAM Journal on Computing, 2(4), 225–231.
Hopcroft, John, and Tarjan, Robert. 1973. Algorithm 447: Efficient Algorithms for Graph
Manipulation. Communications of the ACM, 16(6), 372–378.
Hougardy, Stefan. 2010. The Floyd–Warshall Algorithm on Graphs with Negative Cycles.
Information Processing Letters, 110(8-9), 279–281.
Hu, T. C., and Shing, M. T. 1984. Computation of Matrix Chain Products. Part II. SIAM Journal
on Computing, 13(2), 228–251.
Huffman, David A. 1952. A Method for the Construction of Minimum-Redundancy Codes.
Proceedings of the IRE, 40(9), 1098–1101.
Karp, Richard M. 1978. A Characterization of the Minimum Cycle Mean in a Digraph. Discrete
Mathematics, 23(3), 309–311.
Karp, Richard M., and Rabin, M. O. 1987. Efficient Randomized Pattern-Matching Algorithms.
IBM Journal of Research and Development, 31(2), 249–260.
Knuth, Donald E. 2000. Dancing Links. arXiv preprint cs/0011047.
Knuth, Donald E., Morris, Jr, James H., and Pratt, Vaughan R. 1977. Fast Pattern Matching in
Strings. SIAM Journal on Computing, 6(2), 323–350.
Kosaraju, S. Rao. 1979. Fast Parallel Processing array Algorithms for Some Graph Problems
(Preliminary Version). Pages 231–236 of: Proceedings of the Eleventh Annual ACM
Symposium on Theory of Computing. ACM.
Kozen, Dexter C. 1992. The Design and Analysis of Algorithms. Springer Verlag.
Kruskal, Joseph B. 1956. On the Shortest Spanning Subtree of a Graph and the Traveling
Salesman Problem. Proceedings of the American Mathematical Society, 7(1), 48–50.
Kuhn, Harold W. 1955. The Hungarian Method for the Assignment Problem. Naval Research
Logistics Quarterly, 2, 83–97.
Lo, Chi-Yuan, Matoušek, Jiří, and Steiger, William. 1994. Algorithms for Ham-Sandwich Cuts.
Discrete & Computational Geometry, 11(1), 433–452.
Manacher, Glenn. 1975. A New Linear-Time On-Line Algorithm for Finding the Smallest Initial
Palindrome of a String. Journal of the ACM, 22(3), 346–351.
Moore, Edward F. 1959. The Shortest Path through a Maze. Pages 285–292 of: Proceedings
of the International Symposium on Switching Theory 1957, Part II. Harvard Univ. Press,
Cambridge, Mass.
Morrison, Donald R. 1968. PATRICIA - Practical Algorithm To Retrieve Information Coded in
Alphanumeric. Journal of the ACM, 15(4), 514–534.
250 References
Munkres, James. 1957. Algorithms for the Assignment and Transportation Problems. Journal
of the Society for Industrial and Applied Mathematics, 5(1), 32–38.
Ollivier, François. 2009a. Looking for the Order of a System of Arbitrary Ordinary Differential
Equations. Applicable Algebra in Engineering, Communication and Computing, 20(1),
7–32.
Ollivier, François. 2009b. The Reduction to Normal Form of a Non-normal System of
Differential Equations. Applicable Algebra in Engineering, Communication and Computing,
20(1), 33–64.
Papadimitriou, Christos H. 2003. Computational Complexity. John Wiley and Sons Ltd.
Roy, Bernard. 1959. Transitivité et connexité. Comptes rendus de l’Académie des Sciences, 249,
216–218.
Schrijver, Alexander. 2002. On the History of the Transportation and Maximum Flow Problems.
Mathematical Programming, 91(3), 437–445.
Stoer, Mechthild, and Wagner, Frank. 1997. A Simple Min-Cut Algorithm. Journal of the ACM,
44(4), 585–591.
Strassen, Volker. 1969. Gaussian Elimination Is Not Optimal. Numerische Mathematik, 13(4),
354–356.
Tarjan, Robert. 1972. Depth-First Search and Linear Graph Algorithms. SIAM Journal on
Computing, 1(2), 146–160.
Ukkonen, Esko. 1985. Finding Approximate Patterns in Strings. Journal of Algorithms, 6(1),
132–137.
Warshall, Stephen. 1962. A Theorem on Boolean Matrices. Journal of the ACM, 9(1), 11–12.
Yang, I-Hsuan, Huang, Chien-Pin, and Chao, Kun-Mao. 2005. A Fast Algorithm for Com-
puting a Longest Common Increasing Subsequence. Information Processing Letters, 93(5),
249–253.
Index