Mathematical Programming MATH20014: Group projects
March 29, 2022
Outline of the group projects
The group projects for MATH20014 will be organised as follows. From the set of potential
projects detailed below, each of you will choose a favourite project and backup projects in
decreasing order of preference, according to their interest to you. The projects are all of a
similar level of difficulty, but focus on different areas of mathematics, including mechanics,
statistical mechanics, combinatorics, cryptography, number theory and knot theory.
Each group will then consist of four students who have shown interest in the same project;
the groups will be assembled by the team of lecturers.
Format of the group projects
We expect the group project to take the form of a jupyter notebook, similar to the lectures and
tutorials. Your jupyter notebook will contain the different pieces of code, and the same jupyter
notebook will also need to contain explanations, formulas, graphs, context and references for
your project.
As you have seen in the lectures and tutorials, the notebooks understand latex, so typesetting
equations is straightforward. On the other hand, it is difficult to include external images into
them (as opposed to plots generated from your code, which are fine). Note that html image
links break as soon as you move from your own machine, even when they shouldn’t. Therefore,
in addition to the core jupyter notebook, we may also create and submit other forms of content,
including the following.
• Python code (.py) files with, for example, modules or classes that you have created, or
scripts (for example scripts that carry out the time-consuming calculations or animations
required by some projects).
• High-resolution images and movies obtained as output from your program.
• Additional notes with graphs and images, as well as explanations or discussions, as a
compiled latex file (i.e. a pdf), or even as a Word file.
Assessment criteria
The projects here contain a core part (marked [core] or (core) below) which you need to
satisfactorily complete to obtain a passing mark. The other parts are pick-and-choose, and open
ended. A brief overview of the grading criteria is given below. (A more complete summary will
be provided before the Easter break.)
• Correctness of the code, and its ability to perform the computations required in the core
parts of the project.
• Quality of your code, including structure, number and usefulness of comments, and use of
functions, arrays, libraries and other elements introduced in the course.
1
• Scope of your project, i.e. how far the topic has been developed beyond the core parts.
• Quality of the surrounding explanations and analysis, and quality of the typesetting and
English language.
There will also be a peer assessment of your group where you rate the commitment and
performance or your group members. This will take a similar form to the group projects in
Mathematical Investigations in year 1. The projects and peer assessments are due on
Monday, May 9 (i.e. the first day of revision week).
How to collaborate on a computing project
The first step is to make sure to have a shared drive that is accessible to all members of your
team, for example by creating a shared folder, for example on Onedrive or Dropbox. Then create
your Jupyter notebook for the project, and edit it. A key issue is version control, i.e. that
there is only one current version of the project (and not, say, four different branches that four
different people are working on). At the same time, you will want to keep the previous versions
of the notebook, so that potentially important pieces are not overwritten or lost and you can
always undo changes that were a mistake.
The simplest (though not the most efficient) form of version control is to save a copy of the
Jupyter notebook with the current date and your name each time you make an edit, and to
always, always start from the most recently modified file.
If you are keen to explore more formal options for working on a coding project as a team,
please consider using the git software, and optionally the github online code sharing platform
(https://github.com/).
A second very important point: the projects generally start from a basic idea and then
develop this through several stages. It is crucial that you do the beginning pieces first. It is
also crucial that you work on one version of it together. For example, to figure out energy
conservation in planetary motion, you need to have a working integrator of the four ODEs first,
etc. Please do not try to develop pieces of the project individually and then merge, as this is
almost impossible to do due people’s different choices for the structure of the code. You should
begin by having an online four-person meeting where you work through the project and set up
the code structure, followed by regular meetings.
2
1 Planetary motion
We want to investigate how planets move around the sun, and how moons move around the
planets. Everything starts with the Newtonian equations of motion for the positions of each
object (sun, planet, moon) i:
mir¨i = G
∑
j
mimj
rj − ri
|rj − ri|3 , (1)
where themj are the masses of the other objects, rj are their positions andG = 6.6734810−11 m
3
kgs2
is the universal gravitational constant.
Gravitation is a conservative force; that means that the total energy of the system, E = T+V
does not vary in time
E =
1
2
∑
i
mi|r˙i|2 −G
∑
i
∑
j
mimj
|rj − ri| . (2)
It is theoretically not trivial to find a numerical integration scheme that is consistent with
energy conservation, or what is called Hamiltonian dynamics in classical mechanics. If we find
it, such an algorithm is called a symplectic integrator. The standard algorithms that we have
discussed, like the Euler method or the Runge-Kutta algorithm, are not symplectic integrators.
Instead, we are going to use the velocity Verlet algorithm, which is accurate to second order.
Suppose that our equations of motion can be written as standard Newton equations with a force
F(r) which depends on positions only:
r˙ = v (3)
v˙ = F(r)/m = a (4)
Then in the velocity Verlet algorithm, we include the acceleration in our calculations. As in
the Runge-Kutta, we need to compute positions, velocities and accelerations in a precise order:
rk+1 = rk + vkdt+
1
2
ak dt
2 (5)
ak+1 =
F(rk+1)
m
(6)
vk+1 = vk +
1
2
(ak + ak+1) dt (7)
Simulate the motion of at least one planet with the velocity Verlet algorithm. In the points
below, the non-core numbers 4-7 are pick and choose, and you can develop the point as far as
you would like to.
1. [core] While the equations above are written in three dimensions, conservation of angular
momentum makes the orbits of planets planar. All in all, for a single planet, you will be
left with only two degrees of freedom, i.e. two coupled second order ODEs, corresponding
to four coupled first-order ODEs. First simulate only the motion of an earth-like planet
around the sun. The mass of the sun is 1.989 1030kg, and the radius of the earth is on
average 149.59787 million kilometers (i.e. R = 149.59787 109m). Assume that the sun is
fixed and that only the earth moves (a pretty good approximation in practice).
2. [core] Verify that you find a closed orbit for a planet with the earth’s initial position
r0 = (R, 0) and tangential initial velocity v0 = (0, 29.8km/s). Check that the total en-
ergy (kinetic plus potential) remains constant. Analytically derive that this value of v0
is the only one compatible with a (near) circular orbit. You will find the Kepler prob-
lem homework from Mechanics 2 very useful here. Hint: Consider the energy in polar
coordinates.
3
3. [core] Now vary the total energy of the orbit by changing the value, but not the direction,
of v0, between 0 and a maximum value to be determined below. Plot the now elliptical
and hyperbolic trajectories, and the positions of the aphelion/perihelion. For the elliptical
trajectories, calculate the orbital eccentricity (defined as the ratio of the minor to the major
axis). Construct a criterion to distinguish elliptical and hyperbolic trajectories. Determine
the escape speed of an earth-like planet from the sun’s gravitational well. Determine the
corresponding energy, and compare it to the theoretical result, which you can derive from
(T + V )|initial = (T + V )|r=∞.
4. For the elliptical trajectories, test Kepler’s second law, which states that “A line joining
a planet and the Sun sweeps out equal areas during equal intervals of time”. In other
words, if A is the area enclosed between the sun at the centre and the planetary position
at (r, θ) in polar coordinates, the rate of change of area dAdt =
r2
2
dθ
dt is constant. This
result is a consequence of the conservation of the out-of-plane angular momentum in the
equations of motion: rederive Kepler’s second law analytically as well.
5. Investigate the three-body problem where three objects interact with each other gravi-
tationally. For example, consider three sun-like stars with different but not too dissimilar
masses. The three-body problem does not have an analytical solution and almost always
leads to chaotic trajectories (that is, when it doesn’t eject one of the three suns).
6. Simulate the solar system including multiple or all planets (and yes, planets interact with
each other). This is harder than it looks, since the periods of the outer planets are much
longer than the inner planets. We therefore obtain numerically stiff differential equations,
where we need a small time step because of the inner planets, but also a long simulation
time for the outer planets. Several compromises are possible, e.g. only simulating inner or
outer planets, and not simulating long enough to find closed trajectories for outer planets.
7. Simulate the trajectory of a large (at least Jupiter sized) exoplanet in a star system, with
the planet close to the star. Drop the assumption that the star itself is not moving, and
compute the “wobble” of the star’s trajectory itself around the common centre of gravity.
Make an estimate at what distance this wobble would be observable from earth: this is
one of the methods for detecting exoplanets.
4
2 Hybrid Cryptography
The Vigenère cipher is a method of encrypting alphabetic text by the use of a sequence of
shift ciphers based on the letters of a keyword. In this project you will firstly develop a hybrid
encryption system using the Vigenère cipher combined with the the RSA protocol for encryption
of the keyword. You will then design a further hybrid system by adding an extra layer of random
encoding. You will analyse and illustrate the security issues inherent to these systems. (You
may choose to do one of the two extensions suggested in parts 7 and 8.)
Note. The text, Jupyter Notebook, and pdf files mentioned below are available in the GitHub
repository for this project at https://github.com/cmh42/hc.
Remark. The Vigenère keyword is not assumed to be an actual word. In fact it can be of any
length up to the length of the message to be encrypted. Note that we will refer to it as the
Vigenère key or simply as the key if the context is unambiguous.
1. (core) Implement functions to encrypt and decrypt messages using the Caesar cipher and
the Vigenère cipher. Note that your functions only need to handle messages containing
alphabetic characters. Your Caesar cipher functions should be able to perform 26 possible
shifts (including the trivial shift) and your Vigenère cipher should have 26 possible choices
for any character in the key.
Note. You should test your functions using randomly generated Caesar shifts and randomly
generated Vigenère keys. You should organise this so that the reader can perform these
tests. For your tests you might, for example, extract the alphabetic content of two or more
of the message_*.txt files provided for this project. You should implement encryption
and decryption to and from files. Similar comments apply to the tests that you should
apply throughout this project.
2. (core) Implement a function to systematically break the Caesar cipher using letter fre-
quency analysis. (Regarding the latter see the example in get_online_texts.ipynb.)
3. (core) In our presentation of the RSA protocol in Week 9 we use 512-bit (i.e. 154 digits in
decimal) primes p and q. The security of the protocol relies on the fact that that it is VERY
HARD to recover p and q—i.e. to factorise N—from N = p · q if you only know N . To see
that this is indeed the case, and also to see what happens when we allow p and q to be
smaller, your task here is to test the performance of the smallest factor function1 from
lectures. To do this generate primes p, q and input N = p · q to the smallest factor
function. Starting with l = 16 bit primes write an algorithm that shows the average
computation time on input N = p · q for k-bit primes p, q for k = l, l + 1, l + 2, . . . .
Continue this analysis for as long as the outcome is a matter of minutes—e.g. up to 15
minutes. Plot your results and the expected outcomes (extrapolated from your results)
on longer bit lengths. Hence conjecture at what bit length the use of smallest factor
becomes unfeasible.
4. (core) A more efficent way of factoring large integers is via the Pollard rho method.
Write a function that implements the Pollard rho method using the outline given in the
file pollard rho.pdf. Carry out the analysis that you carried out on smallest factor
on your function. Plot the results and, for comparison, include on your graph the results
from above relating to smallest factor. Again conjecture from your analysis at what
bit length the use of your function becomes unfeasible.
1Note the decompose function (which uses smallest factor as a subfunction) is not required here since you
are working under the assumption that N = p · q with p and q prime. Thus once you have found one factor,
which must be either either p or q, and assuming for the sake of argument that it is p, then you can extract the
other factor q directly by dividing N by p.
5
5. (core) Write functions that implement the Hybrid System described below with the en-
cryption and decryption of the message carried out using your Vigenère functions from
above and that of the key being carried out by the RSA functions from lectures. Note
that, since the Vigenère key may be long—for example 200 characters—your system will
need to slice it in to one or more parts (so no slicing if only one part) in preparation for
integer conversion/encoding and RSA encryption with the resulting ciphertext integers
being transmitted as a tuple2.
Hybrid System. Alice generates her private and public key. Bob generates a Vigenère
key and Vigenère encrypts/enciphers his message with this key. Then, after slicing it into
parts (if necessary) he encodes and RSA encrypts his Vigenère key using Alice’s public key
and finally sends both the resulting tuple of ciphertext integers and his Vigenère encrypted
message to Alice. Alice uses her private key to RSA decrypt the tuple of ciphertext integers.
She then converts/decodes the resulting integers to strings and so reconstructs the Vigenère
key. She uses this to Vigenère decrypt/decipher Bob’s message.
6. (core) There are 26 · 25 = 650 many 2-grams3 made up of distinct letters. Redesign your
system by performing a random encoding of each letter of the alphabet in to one or more of
such 2-grams. For example you might disguise the frequency of the letter “e” by randomly
choosing a set of 8 (maybe more) different 2-grams to represent different instances of this
letter4. Your message is now randomly encoded before Vigenère encryption and decoded
after Vigenère decryption. Your encoding information should be recorded in a key which
is then appended to the Vigenère key. As before, the resulting string should be sliced
into parts for integer conversion and RSA encryption before transmission as a tuple of
ciphertext integers. After transmission and RSA decryption of both keys the receiver will
be able to obtain the original message. Discuss and compare the security of this hybrid
system and that of part 5, taking the implications of part 8 below into account.
7. (extension) Applying the idea from part 6, modify your system so that it handles mes-
sages that contain, not only letters, but also numbers, punctuation and white space. (You
may need to use the full set of 650 distinct letter 2-grams.) Use this system to encrypt
and decrypt the entire content of two of the message_*.txt files5.
8. (extension) Implement a function to systematically break the Vigenère cipher. To do this
you will need to first perform a Kasiski style analysis of the positions of repeated n-grams
in the encrypted message to work out the length of the key. You will then reapply the
letter frequency analysis that you developed in part 2 to establish the letters of the key.
9. (after thought) With the knowledge that you gained in parts 3, 4 and 8 you can write
functions to check how small the primes p and q used in the RSA encryption of the key
need to be for the RSA part of the hybrid system developed in part 5 to be as vulnerable
to attack as the Vigenère encryption of the message. For example, with p and q of bit
length 64 will Pollard rho break the key encryption before your Vigenère function breaks
the message encryption directly?
2For example, the RSA protocol with 512 bit primes can safely handle strings of length b(512 − 1)/8c = 63
(where we note that each character is encoded with 8 bits and a 1 is added to the front of the ciphertext integer).
Thus if the Vigenère key is of length 200 we will want to slice it into 4 parts, convert/encode these 4 strings into
integers for RSA encryption, and transmit the resulting 4 ciphertext integers as a 4-tuple.
3An n-gram is a string of n letters. The term n-gram is often used to denote a contiguous sequence (i.e. a
substring) of n letters within a longer string. For example ‘ing’ is a 3-gram in ‘Flying high’.
4Once the set of 2-gram encodings of “e” has been chosen, the choice of which 2 gram to use for which instance
of “e” is made randomly during the encoding process and does not need to be recorded.
5Your system should be able to handle the white space character and all the 32 characters in the constant
string string.punctuation. Thus 26+1+32 = 59 characters in all. Any other characters such as tabs or newlines
can be left unmodified by your system during the whole encryption, transmission and decryption process. (And
you should make sure that your system preserves the case of letters.)
6
3 Aliquot sequences
For a positive integer n the sum of proper divisors function is
s(n) =
∑
{d|n,0d.
It gives the sum of all positive divisors of n, excluding n. Interest in this function goes
back to the Pythagoreans (6th century B.C.E.). Aliquot sequences are the seqeunces formed by
repeatedly applying this function.
Definition. For each n ∈ N the sequence
An = {n, s(n), s2(n), s3(n), . . .}
generated by applying the function s repeatedly is the aliquot sequence starting at n. If sj(n) = 0
for some j, the sequence terminates after j.
There are some interesting seqeunces. For example, s(220) = 284 and s(284) = 220. This
means {220, 284, 220, 284, 220, . . .} is an infinitely looping aliquot sequence. These are called
‘amicable’ numbers. Other examples arise from ‘perfect’ numbers such as 6 which has the
property that s(6) = 6. There are also sequences which terminate at 0, such as {7, 1, 0}. The
same should happen starting at any prime.
Surprisingly little is known about aliquot sequences. It is easy to see that there are three
possible types of aliquot sequences:
1. Those which terminate at zero.
2. Those which enter a loop.
3. Those which continue infinitely but do not contain repeats.
It is not currently known whether there are any of type 3 but equally it is possible that most
aliquot sequences are of this type. The goal of this project is explore these sequences computa-
tionally.
One of the difficulties is that calculating s(n) becomes computationally difficult once n is
large because it involves factoring n into primes. Your project should consist of two main parts:
• Write code to compute some aliquot sequences.
• Use this code to explore questions about these sequences.
The efficiency of your code will determine how effectively you can investigate the sequences.
To do:
1. (core) Write a function to calculate s(n). (Note: In week 9 we’ll see a method that
exploits the fact that the sum of all divisors is a multiplicative function.)
2. (core) For a given k, compute the aliquot sequence starting at k (up to a sensible point).
Practical suggestion: Write your code so that it computes at most the first n terms of
the sequence. Also, write it so that for some i, your code stops computing new terms once
sj(k) > i. At first, you can use lower values but you should aim to get code that runs in
a reasonable time for n = 30 and i = 109.
3. (core) Find a way to detect loops. (There’s a little bit to this – look at the sequence
starting at 562 to see why.)
7
4. (core) For each k < 20000 try to classify it according to the end state of the aliquot
sequence starting at k. It should either terminate at zero, enter a loop or be unknown
(you might want to distinguish between the cases where the calculation was cut short
because you reached term n and those where the sequence exceeded i).
5. (core) If your code takes too long to do (4) for this range for k, revisit steps (1) and (2).
Here are some questions about aliquot sequences to explore:
1. (core) How long can the loop be if an aliquot sequence ends in a loop?
2. (core) If a sequence terminates at zero, how long can it be (relative to the starting n)?
3. Find a nice way of plotting to visualise aliquot sequences which enter loops.
4. Related to perfect numbers are ‘abundant’ numbers with s(n) > n and ‘deficient’ numbers,
where s(n) < n. A looping aliquot sequence should contain some of each. Compare the
number of each up to a fixed value n.
5. For a fixed n, how large can the preimage s−1({n}) be? This should tell you something
about when aliquot sequences starting at different values merge together.
6. Can you make your code faster by precomputing factorisations for small integers (e.g. for
all integers up to 106)? Investigate using the sieve of Eratosthenes to do this efficiently.
8
4 Knot Theory and Recursion
Continued fractions
Let pq be a rational number, where p and q are coprime. A continued fraction for
p
q is a sequence
[a1, a2, . . . an] where each ai ∈ Z so that
p
q
= a1 +
1
a2 +
1
a3 +
1
. . . + 1an
.
You might also see continued fractions for any real number. Rational numbers have finite
continued fraction expansions. Continued fractions can be calculated recursively using the quo-
tients you get in Euclid’s algorithm. They are not unique, for example
[1, 3, 2] = 1 +
1
3 +
1
2
=
9
7
= 2 +
1
−2 + 1
2 +
1
−3
= [2,−2, 2,−3].
To do:
1. (core) Write at least two functions which calculate a continued fraction for a given pq (to
generate different expansions).
2. (core) Write a function that recalculates the rational number from a list [a1, . . . an].
2-bridge links
A knot or link is a embedded loop (for a knot) or several loops (for a link) in 3-dimensional
space. Two knots are equivalent if one can be deformed into the other without passing through
itself. It is usually best to visualise them using knot diagrams. The diagrams in Figure 1 are
called 2-bridge because the pictures each have two local maxima.
Figure 1: The 2-bridge diagrams corresponding to [1, 3, 2] and [2,−2, 2,−3]. They are different
pictures of the same knot.
9
For every rational number pq (with p and q coprime) there is an associated 2-bridge link
S(p, q). Here’s the algorithm for drawing it: Start with four parallel strands. Then take a
continued fraction expansion for pq . For each ai add twists between two strands with ai crossings
– use the rightmost two strands if i is odd and the middle two if i is even. The crossings should
be right-handed (i.e. right strand over left as you go down) if i is odd and ai is positive and
should change orientation if you switch either of these. Then connect all these twist regions
together and connect the top and bottom of the diagram. The pattern for this depends on if
you had an odd or even number of twist regions (but the key is that it shouldn’t immendiately
trivialise the last one). (See Figure 1, [1, Figure 3.9 on p59] or [2, Section 1]. The sign choice
varies depending on where you look but the important thing is to pick one and stick with it.)
To do:
1. (core) Write code that prints a 2-bridge link diagram for a given pq . It is fine for the
diagram to be quite ugly, e.g. built from pieces similar to
if crossing_type == 1:
print("| | \ / ")
print("| | / ")
print("| | / \ ").
2. (core) Use this function help you check these facts about 2-bridge links for three or four
cases:
(a) S(p, q) is a knot if p is odd and has 2 components if p is even;
(b) the link S(p, q) does not depend on the choice of continued fraction, even though the
diagram does;
(c) if n ∈ Z then S(p, q) = S(p, q + np) (i.e. we only need q mod p);
(d) changing the sign of all the crossings in S(p, q) gives you S(p,−q) = S(p, p− q).
3. Investigate the relation between the links S(p, q) and S(p, q′) when qq′ ≡ 1 mod p.
d-invariant
The d invariant of S(p, q) is a collection of p rational numbers6. These can be denoted by
d(S(p, q), i) and when p > 0 they are defined recursively by the rules[3, Proposition 4.8]7:
• For any q and i, d(S(1, q), i) = 0;
• If a%p denotes the integer between 0 and p which is congruent to a mod p then d(S(p, q), i) =
d(S(p, q%p), i%p);
• When p > q > 0 and 0 ≤ i < p+ q,
d(S(p, q), i) =
(2i+ 1− p− q)2 − pq
4pq
− d(S(q, p), i).
Since these are rational, you’ll need to need to be careful about rounding errors. For example
d(S(9, 7)) = [0,
−2
9
,
4
9
, 0,
4
9
,
−2
9
, 0,
−8
9
,
−8
9
].
To do:
6Technically there is a missing step. First we should use the link to determine a 3-dimensional manifold L(p, q)
and then d is an invariant of that manifold.
7The orientation convention in that paper is different, hence the different sign.
10
1. (core) Write a function that calulates d(S(p, q), i).
2. (core) Write a function that computes the list of d invariants for S(p, q).
3. For examples where you know S(p, q1) = S(p, q2) determine how this is reflected in the
d-invariant.
Examine these two questions about the d invariant:
4. (core) Does it completely distinguish between 2-bridge links. In other words, is it true
that whenever S(p, q1) and S(p, q2) are inequivalent links they have differing d invariants?
Check this for a collection of example, such as when p < 100. If possible, try to refine your
code to work for all p < 200 or p < 1000.
5. In the case where p = m2 for an odd number m, it is interesting to look at which knots
S(m2, q) have the property that m of the d-invariants are zero. (This is related whether
the knot is what is called a slice knot – this is something you can detect from a diagram
but it’s hard to computationally.) For m < 105 determine which of the knots S(m2, q)
have this property and compare the numbers m
2
q to the set R in [2, Definition 1.1].
References
[1] P. Cromwell, Knots and Links, Cambridge University Press (2004)
[2] P. Lisca, Lens spaces, rational balls and the ribbon conjecture, Geom. Topol. 11 (2007) 429–
472. (Also at https://arxiv.org/abs/math/0701610.)
[3] P. Oszváth and Z. Szabó, Absolutely graded Floer homologies and intersection
forms for four-manifolds with boundary, Adv. Math. 173 (2003), 179–261. (Also at
http://arxiv.org/abs/math/0110170v2.)
11
5 Percolation
Do an image web search for “critical percolation” for some inspiration.
Consider the square grid {0, . . . , n−1}2 = {(i, j) : i, j ∈ {0, . . . , n−1}}. We call its elements
sites. We declare each site to be yellow with probability p, or else blue with probability 1− p,
where p is a parameter between 0 and 1, and different sites receive independent colours. Two
sites are considered adjacent if they are Euclidean distance exactly 1, so each site is adjacent
to at most 4 others (but fewer if it is at a side or corner of the grid). A path means a finite
sequence of sites in which each consecutive pair is adjacent. We are interested in the probability
Fn(p) = P
(
there is a yellow path connecting the left and right sides of the grid
)
.
1. (core) Produce a visualization of the randomly coloured grid for small values of n and
several choices of p, and sample it several times. This could be as simple as a text repre-
sentation of a matrix, but we recommend using a ‘pcolor’ plot as for the Mandelbrot set
tutorial, or something similar. Get a rough idea of how the probability Fn(p) behaves by
inspecting the visualization manually.
2. (core) Write a function to determine whether there is a yellow path connecting the left
and right sides.
Here is a suggested approach. Iteratively determine the set of sites that are reachable by
yellow paths starting from the left side, as follows. Initially all yellow sites on the left side
are known to be reachable. For every known reachable site x and every adjacent yellow
site y not currently known to be reachable, we can declare y reachable. Keep doing this
until no new reachable sites can be added. Then check whether any site on the right edge
is reachable. Test your code by comparing it with the visualization. It might be helpful to
also label the reachable sites in the visualization. Do not worry about efficiency for now.
Hence estimate Fn(p) for some small n, perhaps 5 or 10, and some choices of p.
3. (core) We now want to investigate larger n. Make your function in (b) more efficient, so
that it runs as quickly as possible. Here are some suggestions. Ensure that each site is not
examined more times than necessary. One approach is to maintain a list of reachable sites
and iterate over it in order, examining the neighbours, and adding any new reachable sites
discovered to the end of the list. To check whether a site is already known to be reachable,
using “if site in list” is not efficient because it iterates over the whole list. Instead,
maintain a separate n-by-n array (or a python set). The search can be stopped as soon
as we reach any site on the right side of the grid.
4. (core) The function has a limit: F (p) = limn→∞ Fn(p) which satisfies
F (p) =
{
0, p < pc
1, p > pc,
where pc is called the critical point. Plot graphs of Fn(p), and use them to estimate pc
as precisely as you can. Even with the optimizations discussed above, expect to run your
code for a few hours at least.
Further options (not necessarily to be followed sequentially):
5. Do the same for rectangular grids of different shapes, e.g. 2n-by-n or 3n-by-2n. The critical
point should be the same.
12
6. Investigate similarly the probability
Gn(p) = P
(
there is a yellow path connecting the centre of the grid to the boundary
)
(where, for simplicity, n can be assumed to be odd) and its limit G(p) = limn→∞Gn(p).
Try to understand what the graph of G looks like.
7. Do the same for the triangular lattice, which consists of sites at the corners of equilateral
triangles, with 6 triangles meeting at each site. This can be conveniently implemented as
a square grid, but with diagonal pairs of sites of the form x and x+ (1, 1) now in addition
considered adjacent. The critical point of site percolation on the triangular lattice has a
very simple exact form. Use your estimates to try and guess it. (You can look it up to
check).
8. For the triangular lattice, G has an asymptotic power law behaviour near pc:
G(pc + ) ≈ const β as ↓ 0.
Estimate the power β. You will need to use the exact value of pc, and even then it is tricky
to get a good estimate. You may want to look up the expected answer.
9. Implement a different method for detecting crossings, in which we explore the boundary
between reachable and non-reachable sites by "wall following”. How many sites typically
need to be explored? How does the answer depend on n and p?
10. Implement the following more sophisticated method, which allows us to effectively sample
all values of p simultaneously. Each site is assigned an independent Uniform[0, 1] random
variable. The model with parameter p is then defined by declaring all those sites with label
less than p yellow. We can compute the minimum p for which a yellow path connecting
the regions of interest exists, by a modification of the iterative scheme used before.
11. (advanced) For the triangular lattice, take the parameter to be exactly p = pc, take a
region in the shape of an equilateral triangle of side length n (now with respect to the
true geometry of the lattice, not the square grid implementation above). Investigate the
probability Tn(r) that the base of the triangle is connected by a yellow path to the right
side within distance rn of the top vertex. In the limit n → ∞ there is a very simple
formula for this. Moreover, it should be the same if the equilateral triangle is rotated by
an arbitrary angle.
12. (advanced) Investigate the following percolation model. We have infinitely many equally
spaced parallel wires. Each wire has breaks according to a Poisson process. Each adjacent
pair of wires has bridges connecting them according to a Poisson process. The critical
point is when the rates of the Poisson processes are equal. At the critical point it should
have the same asymptotic behaviour for crossing probabilities of equilateral triangles as
above, although this has not been proved rigorously. It is also necessary to choose the
correct spacing of wires as a function of the rate of the Poisson process for this to work
out correctly.
13
6 Code Breaking with Statistical Physics
In the appendix there are several secret messages for you to decrypt, in roughly increasing order
of difficulty. Each message is written using the 26 letters of the alphabet together with spaces.
(There are no other symbols such as punctuation). Each message is encrypted using a simple
substitution cipher; that is, a certain (secret) permutation of the 27 characters has been applied
to each character of the message in turn. Your job is to crack the code and discover the message.
Ideally your final answers you should include your best guess at the correctly punctuated and
cased message, as well as the raw decrypted version. Of course, the last part will require human
involvement.
The provided messages are there to motivate and focus you, but are not necessarily the final
destination. The broader goal is to investigate decryption methods.
For your convenience, a notebook containing the messages as strings will be provided on the
Blackboard page.
1. (core) In the first message, the cipher is simply a cyclic shift of the 27 characters. E.g. a
shift of 2 would map A 7→ C, B 7→ D, . . . , Y 7→ space, Z 7→ A, space 7→ B. Find the shift
and decode the message. One simple approach is to try all possibilities and check by hand
which one produces readable text.
2. (core) The other messages each use an arbitrary permutation of the characters, so trying
all the possibilities is no longer practical. Instead we will make use of statistical properties
of text. In preparation for this, find a large body of English text (the larger the better) to
download and analyse. Replace all lower case letters with upper case, and find a way to
get rid of extra characters. One simple approach is to replace every character that is not a
capital letter or a space with a space, and then replace multiple consecutive spaces with a
single space. More nuanced schemes are possible - for instance, apostrophes and hyphens
could be treated differently.
For example, the following code will download the text of the novel Moby Dick and store
it in a single very long string.
import requests
url=’http://www.gutenberg.org/files/2701/2701-0.txt’
moby=requests.get(url).text
To avoid repeatedly spamming someone’s server, it is best to save the text to your own
machine and access it from there:
with open(’C:/data/moby.txt’,’w’,errors=’ignore’) as f:
f.write(moby) # write to file
with open(’C:/data/moby.txt’,’r’,errors=’ignore’) as f:
moby=f.read() # read from file
3. (core) One simple approach is to use character frequencies. In the sample English text,
count the frequency with which each of the 27 characters appears, and sort them by
frequency. Do the same for the coded message, and try the cipher in which the kth most
common character in the message maps to the kth most common letter in English, for
each k.
14
4. (core) The above may not be sufficient to get the full correct answer, but for the longer
messages it should be close. Write a program to make it easy for you to try swapping the
roles of any two characters (of your choice) and immediately see the effect of the decoded
message. Then, starting from the result above, adjust the cipher by hand until is readable.
5. (core) The above method is less likely to work for the later messages, because they are not
long enough to get accurate frequency counts (and some might involve unusual language).
Instead use the following approach based on ideas from statistical physics. In the sample
text, record the frequency (i.e. number of occurrences) of each of the 272 bigrams, i.e.
pairs of consecutive characters such as SH, as well as the 27 characters. For each pair of
characters i, j compute
p(i, j) :=
frequency(ij) + 1
frequency(i)
.
This represents the probability that i is followed by j; the +1 is a ‘fudge factor’ to avoid
zeros. To any potential decoded message m = i1i2 . . . in we assign a score
S(m) =
n−1∑
j=1
log p(ij , ij+1),
which measures how plausible m is as a piece of English text (it represents the logarithm
of the probability of seeing m under a certain “Markov model”). Write a function that
computes S(m) for any message m.
6. (core) Implement this method, called the Metropolis algorithm (an example of a Monte
Carlo method). Start with any guess m at the decrypted message, which might simply
be the encrypted message itself or the result of applying a randomly chosen cipher. Now
repeatedly do the following.
Choose two characters at random and swap their roles in the cipher, to get a new candidate
message m′. If S(m′) > S(m) then replace m with m′. If S(m′) ≤ S(m) then replace m
with m′ with probability
exp
S(m′)− S(m)
T
,
otherwise just keep m.
Repeat this for many steps (e.g. 100000 or 1000000).
Here T > 0 is a parameter, which represents temperature in a physics analogy, and indi-
cates how likely the algorithm is to accept steps that make the score worse (in hopes of
finding a route making it better later). You may need to experiment to find the value of
T that works best. Start by trying T = 1. One approach (called simulated annealing) is
to start with T large (say 10), and gradually reduce it as the algorithm proceeds, ending
at something small (say 0.1), perhaps by reducing it by a constant ratio at each step.
It is useful to print the current guessm every so often (e.g. every 10000 steps), and monitor
its progress.
It may be useful to try running the algorithm multiple times, from different starting points,
and to adjust the parameters, and to make manual adjustments to the cipher at the end
(as before).
Further options (not necessarily to be followed sequentially):
For the extensions below, besides the messages provided, you will likely find it helpful to
produce your own coded messages.
15
7. How can the method be improved? The final few messages will require thinking outside
the box. Does it help to consider trigrams, or just common ones, or whole words? What
if the message is in a unknown language? Can it help to have sample texts in several
languages? Might the person encoding the message have been careless somehow?
8. For the simple frequency analysis approach in steps 3-4, how much manual adjustment
is typically needed, and how can that be quantified? Can you represent the permutation
graphically? How does the answer depend on the length of the message? A proper inves-
tigation will require generating your own messages. Do the answers make sense in terms
of a probabilistic model of what is going on?
9. Can you analyse the progress of the Metropolis algorithm over time, perhaps by plotting
graphs of some quantities? Can this help with adjusting T, and choosing an appropriate
number of steps to run for?
10. Can the basic method be adapted to other ciphers such as transposition ciphers, the
Vigenere cipher, or a one-time pad that itself consists of English text? You would need to
look up what these are (e.g. on Wikipedia), and write your own programs for encrypting
as well as decrypting, so that you have ciphertexts to work with. (Perhaps one member of
your group might choose the text and do the encryption in secret from the others.)
16
Coded messages
1. GPKZFER JRSER EKWIGIWKWVRZ YZRCWMWCRSEVRYWEWISCRGLIGFJWRGIFYISDD
EYRCSEYLSYWRGPKZFEJRVWJ YERGZ CFJFGZPRWDGZSJ QWJRUFVWRIWSVST
C KPRN KZR KJREFKSTCWRLJWRFXRJ YE X USEKR EVWEKSK FER
KJRCSEYLSYWRUFEJKILUKJRSEVRFTAWUKRFI WEKWVRSGGIFSUZRS
DRKFRZWCGRGIFYISDDWIJRNI KWRUCWSIRCFY
USCRUFVWRXFIRJDSCCRSEVRCSIYWRJUSCWRGIFAWUKJ
2. RECFV KUWE VJCRWQFCRCFICYWEZRWKLVJWQF SCRRZ ALVWJLAACFWNATZVW
NFWIZRZT FWELYWVCSTWNRWLVTE NBEWZTWHLRWCLRMWS FWJCWHE WUACHWEZJWR
WHCVVWT WRCCWTELTWECWHLRWQF S NAYVMWCPKZTCYWTECWJ JCATWTELTWEZVT
AWKNGZTTWRWGF LYWGLKUWELYWYZRLQQCLFCYWTEF NBEWTECWY FWJMWK
JFLYCWFNRECYWT WTECWTLGVCWVLZYW NTWLVVWTECWRVZQRW SWQLQCFWK
ATLZAZABWYLAKZABWJCAWZAWSF ATW SWEZJWLAYWTEFCHWEZJRCVSWZAT
WLAWZATFZKLTCWLAYWCVLG FLTCWKLVKNVLTZ AWS FWTH WE NFRWZWHLTKECYWEZJWLRWECWK
ICFCYWRECCTWLSTCFWRECCTW SWQLQCFWHZTEWSZBNFCRWLAYWVCTTCFRWR WK JQVCTCVMWLGR
FGCYWZAWEZRWTLRUWTELTWECWELYWCIZYCATVMWS FB TTCAWJMWQFCRCAKCWR
JCTZJCRWECWHLRWJLUZABWQF BFCRRWLAYWHEZRTVCYWLAYWRLABWLTWEZRWH FUWR
JCTZJCRWECWHLRWQNDDVCYWLAYWH NVYWRZTWS FWV ABWRQCVVRWHZTEWLWSNFF
HCYWGF HWLAYWLWILKLATWCMCWSZALVVMWECWRQFLABWSF JWEZRWKELZFWHZTEWLWKFMW
SWRLTZRSLKTZ AWLAYWHLVUCYWNQWLAYWY HAWTECWF JWFNGGZABWEZRWELAYRWT
BCTECFWTECAWECWHF TCWLWV ABWTCVCBFLJWNQ AWLWKLGVCWS FJWZSWJMWLARHCFWT
WTEZRWZRWLRWZWE QCWM NWHZVVWELICWLWICFMWQFCTTMWKLRCWT WLYYWT WM
NFWK VVCKTZ AWHLTR AWRLZYWECWZWCPQCKTWTELTWHCWRELVVWGCWLGVCWT
WB WY HAWT WA FS VUWT J FF HWLAYWT WTLUCW NFWSFZCAYWR
JCWICFMWYCSZAZTCWACHRWLRWT WTECWRCKFCTW SWEZRWLAA MLAKCWZWK
ASCRRWTELTWZWHLRWSZVVCYWHZTEWKNFZ RZTMWGNTWZWHLRWLHLFCWTELTWE
VJCRWVZUCYWT WJLUCWEZRWYZRKV RNFCRWLTWEZRW HAWTZJCWLAYWZAWEZRW
HAWHLMWR WZWHLZTCYWNATZVWZTWRE NVYWRNZTWEZJWT WTLUCWJCWZAT WEZRWK
ASZYCAKCWGNTWTECFCWHLRWLWYCVLMWZAWTELTWLARHCFZABWTCVCBFLJWLAYWTH WYLMRW
SWZJQLTZCAKCWS VV HCYWYNFZABWHEZKEWE VJCRWQFZKUCYWNQWEZRWCLFRWLTWCICFMWFZABW
SWTECWGCVVW AWTECWCICAZABW SWTECWRCK AYWTECFCWKLJCWLWVCTTCFWSF
JWEZVT AWKNGZTTWLVVWHLRWXNZCTWHZTEWEZJWRLICWTELTWLWV ABWZARKFZQTZ
AWELYWLQQCLFCYWTELTWJ FAZABWNQ AWTECWQCYCRTLVW SWTECWRNAYZLVWECWZAKV
RCYWLWK QMW SWZTWHEZKEWZRWECFCWFCQF YNKCYWCVRZCWQFCQLFCWT
WJCCTWTEMWB YWE VJCRWGCATW ICFWTEZRWBF TCRXNCWSFZCDCWS FWR
JCWJZANTCRWLAYWTECAWRNYYCAVMWRQFLABWT WEZRWSCCTWHZTEWLAWCPKVLJLTZ AW
SWRNFQFZRCWLAYWYZRJLMWEZRWSLKCWHLRWELBBLFYWHZTEWLAPZCTM
17
3. TULOMREMAEKBLRSWCTMWBIHB BTN KKUT NBTULOMRB DWBKSYMBYSWMRDBTULOMRKBJUNNBRMQM
NBKE EUKEUT NBUDPSRY EUSDB ISCEBEOMBLN UDEMAEB DWBEO
EBUDPSRY EUSDBT DBSPEMDBIMBCKMWBESBIRM XBEOMBTULOMRB
PEMRBEOMBWUKTSQMRHBSPBPRMVCMDTHB D NHKUKBIHBEOMB R IBY EOMY EUTU DB
DWBLSNHY EOB NBXUDWUB NKSBXDSJDB KB NXUDWCKBUDBEOMBDUDEOBTMDECRHBDM
RNHB NNBKCTOBTULOMRKBTSCNWBIMBIRSXMDBIHB DBUDPSRYMWB EE TXMRBKCTOBTN
KKUT NBTULOMRKBKEUNNBMDGSHBLSLCN RUEHBESW HBEOSCFOBYSKENHB KBLCZZNMKB
NBXUDWUBJRSEMB BISSXBSDBTRHLESFR LOHBMDEUENMWBRUK N OBPUBUKEUXOR
GB NBYC YY BY DCKTRULEBPSRBEOMBWMTULOMRUDFBTRHLESFR LOUTBYMKK
FMKBJOUTOBWMKTRUIMWBEOMBPURKEBXDSJDBCKMBSPBPRMVCMDTHB D NHKUKB DWBTRHLE
D NHKUKBEMTODUVCMKB DBUYLSRE DEBTSDERUICEUSDBSPBUIDB WN DBJ KBSDBK
YLNMBKUZMBPSRBCKMBSPBPRMVCMDTHB D NHKUK
4. RJLHXTNERXWZEWZSXLESNIFWX TEKJAAS NYSVEWZSECYWKZEKJTVYAW
TWVERSLSELYTTXTNETYPDSLVEDJVH AXVVELJASER VED VXK AAIEWJECJEWZSEK
AKYA WXJTVEWZSEVYPVEDSLCJRVHXEV XCEVJERZSTEWZSEVWSLTEK PSEULSSEJTEPJTC
IEPJLTXTNERSEK AKYA WSCEWZ WERSEVZJYACEASWEWJTTSVEJUER WSLED AA VWEXTE
WEWZSELS LEJUEWZSEOSVVSAEWJEFYVZEWZSEVWSLTECJRTE TCEAXUWEWZSEDJR
5. IZWLXSZF QEFLZTQUVLSZMIHZCIVLUZEDDZIUPZXLYAHLPZCQZLICZECZMIHZC QAO CZ
QMLWLXZC ICZECHZIBBLCECLZTEO CZKLZHCETADICLPZKSZIZBEULIBBDL
6. RVRRKZHKTVYPJCGVCHGVGXECUKRKZGRSKZKWGVGRSPKAGVCHGVGHZVLJCAYOGRZET
KRKZGSKWGTOGSKZJGAVPZOGHVCHOGRPQUYPCLGRSKGAVCQOGJAGSPWGYVHOGAZPKCHGRSKGCOT
SGPCGOKYYJMGQVCGMKGWKKGRSKGTVWRKZGWRZJUKGMSVRGVGNEVKZKGAKYYJM
7. LOJ WUZCWKCTFHHFHTWUZFHQWMZEUW EAWKCWUZCWCHV
8. EYENXRPSURPGUUIRKSUGURSHYY NXRCUUIREPRCHBIRENBRMAOKUGIRKUGURS
NX NXROREHABRFEAUBONRBGUKROHPRSURBGONURENBRPORSUGRJ JURKEIRI
NX NXRORPKEIRJ CGOFSRIENXRIPGEPSIJUWIRENBRGUUAIRISURB
GABRPSUYREMMRMHRFAUEGAWRORKSUNRPSUGURFEYRERWUAARORMOGU
XNRIVHUUAIRPSEPRBENXRSUGRPEJIEAPUUG URO
9. WDYRDYLDQCSLR KTYDPYZ LXSKTYWDYQ U KTYZOITYGIDYCDYJSILPTYLDUDKOTYCIQYKQTP
QPYODYPDRZTYCDQPDQCLDYODTYJL GIDRDQPTYSLH
QKGIDTYCDTYMSKTDLKDTYCSIULKLYODTYADINYZSILYXKNDLYODYB
ODKCSTJSZDYCDYOSMTJILKPDYCDYHSYPDLYHLYJDYIQDYOIDILYRSRDQP
QYDYCDYJSQTJKDQJDYODYTSRRDKOYSIYDP KDQPYZOSQHYTYODTYRDIMODTYO
YJF RMLDYODYPSIPYCSQPYWDYQDP KTYGIIQDYZDPKPDYZ LPKDYDPY
YOKQTDQTKMKOKPDYCIGIDOYWDYLDPSILQ KTYUKPDYRIQKL
10. GYYUXPUXGTBQTDFU
11. ZVVFEKTFZKKO
12. RIZOGS BXTPK YWH JQNE LUC MADFV
18
7 Recursive Animations and The Towers of Hanoi
Animations can be a useful tool to visualise and understand solutions to mathematical problems
that can be represented graphically. In this project you will begin by creating animations of
various recursive problems using the python graphical library pygame, in part by integrating its
use with other tools that you have learnt during the course. You will then tackle the Towers of
Hanoi problem and create pygame animations of solutions to the original problem and two of its
variants. As an extra challenge you can either develop a game type animation of the problem
or otherwise tackle the issues involved in dealing with larger numbers of discs.
Note. The Python, Jupyter Notebook, and pdf files mentioned below are available in the
GitHub repository for this project at https://github.com/cmh42/toh.
1. (core) A basic example animation is to simulate a ball bouncing within a square (box)
in two dimensions. You are given a simple pygame animation function bouncing_ball
in the notebook animation_examples.ipynb. You should further develop this function to
achieve the following. (1) The user should be able to change the original position of the
ball and should be able to change the speed of a ball. (2) The ball should slow down under
the effect of gravity. (3) Two or more balls bounce within the same square (and off one
another).
2. (core) Animating fractal constructions gives us direct insight into how the function gener-
ating the fractal operates. You are given a pygame animation function draw_sierpinski
that draws the Sierpinski triangle. Using this function the user is able to choose the depth
of the triangle and to stop and start the animation. You should develop this function
so that the user is also able to change the speed of the animation. You should also add
colours to the triangle drawing.
3. (core) Fractals can simulate shapes found in the natural word. One simple example of
this is the construction of a recursively defined tree. The tree of depth 1 is just a trunk
with three straight branches. Then given the tree T of depth n the tree of depth n+ 1 is
the tree T where every branch has been replaced by a tree of depth 1. In Figure 1 (from
the left) is such a tree of depth 1, 2, 3 and 4 and then a tree of depth 9 partially and
completely drawn.
Figure 2: Recursive Tree
You should develop a pygame animation function that constructs a similar recursive tree.
The user should be able to choose the depth of the tree and should be able to start and
stop the animation and control its speed. The tree should be coloured with the trunk and
branches being brown and the leaves (i.e. the last level of branches) being green.
4. (core) In Homework 3 we saw that there are interesting Julia sets with parameter jp =
0.7885eai where a is a small non-negative real number. In fact if we allow a to vary in
repeated cycles over [0, 2pi] we are able to create a film like animation of Julia sets with
this form of parameter. The images in Figure 3 are snapshots of such Julia sets.
Your task is to develop a function (using the matplotlib.pyplot tools you developed for
Homework 3) to generate a number (for example 200) of image files of the Julia sets with
parameter jp = 0.7885eai where the numbers a are chosen at equally spaced intervals in
[0, 2pi]. You should then develop a pygame function to display these files as a film like
19
Figure 3: Julia sets
animation. The user should have speed and stop/start control over the animation. (Can
you find a number other that 0.7885 that can be used to create such striking visual effects?)
5. (core) In the tutorial of Week 4 we developed functions to display the Mandelbrot set.
Your task here is to create an animation in which the user is able to repeatedly focus in
and out on any part of the Mandelbrot set. To do this you will want to integrate the use
of a function that generates image files (using matplotlib.pyplot) with a pygame function
to control the animation.
6. (core) In the Towers of Hanoi problem we have three poles and n discs that fit onto the
poles. The discs differ in size and are initially arranged in a stack on one of the poles, in
order from the largest disc (disc n) at the bottom to the smallest disc (disc 1) at the top.
The problem is to move the stack of discs from one pole to another pole while obeying the
following rules.
• Move only one disc at a time.
• Never place a disc on one smaller than it.
This problem can be solved by issuing a sequence/list of instructions for moving the discs
in the appropriate way. We assume that the poles are arranged in a row and that each
instruction to move a disc specifies its number and whether to move it left or right. We
allow wrap around : if a disc is on the left pole an instruction to move left means to wrap
around to the right pole, whereas if a disc is on the right pole an instruction to move right
means to wrap around to the left pole. A solution to the Towers of Hanoi problem for 3
discs is displayed in Figure 4.
Figure 4: Towers of Hanoi with 3 discs
We can solve the problem by recursion. First we move the stack comprising the top n− 1
discs to an empty pole, then we move the largest disc (i.e. the bottom disc of the starting
stack) to the other empty pole. Finally we move the stack of n − 1 discs onto (the pole
with) the largest disc. Figure 5 illustrates this recursive approach.
Your task here is to write a function that implements this idea to create, given input n, the
list of moves required to solve the Towers of Hanoi problem for n discs8. You should then
8To get a better idea of how to do this see the description of the Towers of Hanoi problem in the file
towers_of_hanoi.pdf. Note that this is an excerpt of Section 2.3 of [1]. (Do not be confused however by
the fact that [1] uses its own specialist libraries and functions, such as stdio and stdio.writeln.)
20
Figure 5: Solving the Towers of Hanoi recursively
develop a function that displays a printout of the sequence of configurations of a solution
to the problem for n discs. (To be used for a small number of discs in practice - this is a
warm up for the graphical animation.)
7. (core) Using the function you developed in part 6 to compute the list of moves for the
Towers of Hanoi problem you should now develop a pygame function that gives an an-
imation of the problem. The user should be able to control the initial number of discs,
stopping/starting, and the speed of the animation. The initial number of discs may be any
integer in the interval [1, 16]. (Discs should be moved directly from one pole to another,
i.e. without showing any intermediate motion.)
8. (core) Work out the definition for functions required to solve the two following variants of
the Towers of Hanoi problem. (1) The problem where wrap around is not permitted: you
can only move discs between adjacent poles. (2) The problem where initially there are 2n
discs numbered according to size (as before). The discs with odd number are in a stack on
the left pole, the discs with even number are in a stack on the right pole. The problem is
to move the discs with odd number to the right pole and the discs with even number to the
left pole (still satisfying the two rules given above, with wrap around allowed). Develop
your pygame animation function to show solutions to these two problems.
9. (extensions) You can choose to do (a) and either (b) or (c) of the following.
(a) Develop a version of your pygame animation function from part 7 that moves the
discs smoothly between the poles (i.e. now showing the intermediate motion).
(b) Develop a version where a player can play the problem as a game, by moving the
discs from pole to pole.
(c) For larger numbers of discs video simulation of the Towers of Hanoi becomes unfea-
sible9 and the design of the function that generates the sequence of moves becomes
critical to the time/space needed. Develop a simulation of the problem that can deal
with larger numbers of discs and which samples a limited number of the sequence
of stack configurations. Experiment with the use of different designs of the function
that generates the sequence of moves and provide a time complexity analysis of the
performance of each.
Technical Note. To run the video animations in this project you will need to work locally, i.e.
using your own distribution of python or that of the University. (You cannot do this via the
Noteable server.)
References
[1] R. Sedgewick, K. Wayne, and R. Dondero. Introduction to Programming in Python: An
Interdisciplinary Approach. Addison-Wesley Professional, 2015.
9For example a video simulation of the problem with 27 discs, moving 4 discs per second would take more
than one year.
21