How Many Soldiers Do You Need To Beat The Night King?
Posted on Sun 26 May 2019 in misc
Here's an interesting puzzle from Five Thirty Eight's Riddler section that I worked out in Julia.
At a pivotal moment in an epic battle between the living and the dead, the Night King, head of the army of the dead, raises all the fallen (formerly) living soldiers to join his ranks. This ability obviously presents a huge military advantage, but how big an advantage exactly?
Forget the Battle of Winterfell and model our battle as follows. Each army lines up single file, facing the other army. One soldier steps forward from each line and the pair duels — half the time the living soldier wins, half the time the dead soldier wins. If the living soldier wins, he goes to the back of his army’s line, and the dead soldier is out (the living army uses dragonglass weapons, so the dead soldier is dead forever this time). If the dead soldier wins, he goes to the back of their army’s line, but this time the (formerly) living soldier joins him there. (Reanimation is instantaneous for this Night King.) The battle continues until one army is entirely eliminated.
What starting sizes of the armies, living and dead, give each army a 50-50 chance of winning? "
Initial Thoughts¶
This looked like a textbook random walk problem with the task of finding the right initial parameters that would make both the walks hit the X axis before the other equally likely. The dead army is at obvious advantage and on an average would remain constant whereas the living army would be going down at a steady rate of half a soldier per turn. I was tempted to attack this mathematically right away but decided to simulate first to get an estimate of the answer first. My hunch was n times log2(n)
# Some simulations with 100 living army soldiers and 10 dead army soldiers vstack(hstack([plotrandomwalk(10) for run in 1:3]), hstack([plotrandomwalk(10) for run in 1:3]))
Julia Walks¶
# Played around with Julia for kicks and # I'm a convert now! More than the blinding # speed part, I was impressed with just how # bloody pleasant and expressive it is. And I # haven't even explored the homoiconicity and # multiple dispatch parts. # Evolution function starts with m living soldiers and # n dead soldiers and returns 1 if the living win and # 0 if the dead win function evolution(m,n) while min(m,n) > 0 duel = rand([0,1]) duel == 0 ? (m,n) = (m-1, n+1) : (m,n) = (m, n-1) end m > 0 ? (return 1) : (return 0) end # Simula function runs the Evolution function t number of # times and gives an estimate of the probability of living # army winning the duel. # Takes the number of dead soldiers and a function that computes # the number of living soldiers as a function of the number of # dead soldiers as inputs function simula(n,t,f) return sum([evolution(f(n),n) for trial in range(1,stop=t)])/t end
Let's look at how the probabilities emerge for twice, thrice etc. times the soldiers for different numbers of soldiers.
[[simula(n, 10000, x -> m*x) for m in 1:n] for n in 4:10] [0.2242, 0.3556, 0.4437, 0.4938] [0.1813, 0.3123, 0.3864, 0.4507, 0.5031] [0.1453, 0.2676, 0.3447, 0.4141, 0.4529, 0.4922] [0.115, 0.2208, 0.3144, 0.3755, 0.4221, 0.4524, 0.4893] [0.096, 0.1964, 0.2794, 0.3448, 0.3998, 0.4287, 0.4609, 0.4971] [0.0758, 0.1733, 0.2443, 0.322, 0.3578, 0.4041, 0.4359, 0.4593, 0.4952] [0.0605, 0.1531, 0.2304, 0.2881, 0.3421, 0.3893, 0.4088, 0.4451, 0.4809, 0.4867]
$n\;\implies \;n^2$ looks like a promising bet. Let's look at how this stays for different values of n
[simula(n, 10000, x -> x^2) for n in 1:10:100] 0.496 0.4902 0.4831 0.4817 0.4864 0.4893 0.4781 0.48 0.4807 0.4774
I began to notice a clear trend of the probability falling for larger values of n. Tried out a couple of other functions $n^2\; +\;n,\;n^2\;+\;log(n)$ etc and figured I should try solving it analytically. So I start exploring some simple configurations and notice a neat recurrence relation.
$$P(m,n)\;=\;P(m-1,n+1)/2\:+\:P(m,n-1)/2$$I then strayed into combinatorial methods and thought I got it right but soon noticed an error which made me understand this is equivalent to a partition problem with some cunning constraints.
If the living and dead started with armies of sizes $m$ and $n$ respectively, the living could win in $n,\;n+2,\;n+4,....\;\;n+2(m-1)$ steps, corresponding to the living losing $1,\;2,\;...\;m-1$ times
So, the total probability of the living winning this duel $P(L_{w})\;=\;\sum\limits_{i=0}^{m-1} \;P(L_{w_{n+2i}})$
This now seems pretty straightforward since $P(L_{w_{n+2i}})$ is equal to the number of ways in which $i$ losses can be distributed among $n+2i$ duels divided by $2^{n+2i}$, which is ${{n+2i-2} \choose i}/2^{n+2i}$. ($n+2i-2$ since the living couldn't have lost in the last two duels)
Easy Peasy right? Not so fast!
We are overcounting here since a lot of configurations that wouldn't have been possible are being included.
Consider the following configuration for example that would be counted as part of $m\;=\;6,\;n\;=\;10,\;i\;=\;5$
$WWWWWWWWWWWLLLLLWW$
This would have been counted as one among ${{10+2*5-2} \choose 5}$, but it is clearly not valid since the duel would've ended after the first 10 wins itself.
What we have here is the classical problem of finding the number of integral solutions to $a\:+\:b\;+\;c\;+\;....+\;m\;=\;n$ but with a twisted set of constraints $a\;<\;p,\;b\;<\;f(p),\;c\;<\;g(p+f(p))$etc. There were a couple of threads that I felt were worth exploring using recursion and complimentary counting techniques but eventually decided to programmatically solve using the recurrence relation.
The naive programmatic solution is quite trivial and horribly inefficient as one might guess. The optimized function uses the dynamic programming technique. Was pretty pleased with the solution actually where I iterated along the family of lines $x\:+\:y\:=\:k$ and exploited the boundary conditions $f(m,0)\;=\;1,\;f(0,n)\;=\;0$ You can refer to this excellent write up by Dan Swenson, a math professor at Black Hills State University where it is explained in detail. He even proved that P(undead army wins) = P(between x and (x + y − 1) Heads in (2x + y − 1) flips of a fair coin) and also gave a polynomial approximation for the same - 1.099054669y^2 - 0.5y + 0.5 Very cool!
function wincalc_naive(m,n) (n == 0) && (return 1) (m == 0) ? (return 0) : (return wincalc_naive(m-1,n+1)/2 + wincalc_naive(m,n-1)/2) end liner(i) = [(x,y) for x in 2:i-1 for y in 2:i-1 if x+y == i+1] # In hindsight, I should've worked # with a model of origin at (1,1) # and the boundary conditions # f(1,m) = f(1, m-1)/2 etc. function wincalc_opt(n) q = zeros(Float64, (n+1, n+1)) for index in 1:n+1 q[index,1] = 1.0 end for i in 3:n for (x,y) in liner(i) (q[x,y] == 0) && (q[x,y] = q[x-1,y+1]/2 + q[x,y-1]/2) end end return q[2:n,2:n] end
$n^2$ approximation vs the actual function¶
# The minimum number of living for a given number of dead at which living have a better chance of winning (>.5) Gadfly.set_default_plot_size(18cm, 24cm) z = wincalc_opt(1600) y = [findfirst(x -> x > 0.5, z[:,dead]) for dead in 1:1600 if findfirst(x -> x > 0.5, z[:,dead]) != nothing] x = 1:length(y) plot(layer(x=x, y=y, Geom.line,Geom.point, Theme(default_color="orange")), layer(x=1:40, y=map(x -> x^2,x), Geom.line, Theme(default_color="red")))
Some Helper Functions etc.¶
using Gadfly using LaTeXStrings Gadfly.set_default_plot_size(30cm, 20cm) function randomwalk(o, n) steps = rand([-1, 1], n) path = [o] for epoch in 1:n append!(path, path[length(path)] + steps[epoch]) end return path end function firstcrosssimulation(m,bw) epoch = 0 while m > 0 epoch += 1 m += rand([-1,bw]) end return epoch end function firstcross(m,n,bw) return sum([firstcrosssimulation(m,bw) for trial in 1:n])/n end function evolutionplot(m,n) mpath, npath = [m], [n] while min(m,n) > 0 duel = rand([0,1]) duel == 0 ? (m,n) = (m-1, n+1) : (m,n) = (m, n-1) append!(mpath, m) append!(npath, n) end m > 0 ? (result = 1) : (result = 0) return result, mpath, npath end #### Overcounted Function winlivin(living,dead) = sum([binomial(BigInt(dead + 2*deaths - 2),BigInt(deaths))/(2^(dead + 2*deaths)) for deaths in range(0, stop=living-1)]) function plotrandomwalk(n) result, mpath, npath = evolutionplot(n^2,n) p = plot(layer(x=range(1,stop=length(mpath)), y=mpath, Geom.line, Theme(default_color="orange")), layer(x=range(1,stop=length(npath)), y=npath, Geom.line, Theme(default_color="red"))) return p end
plot(x = range(1,stop=21), y = randomwalk(30, 20), Geom.point, Geom.line)