Background: 20x Going Viral

I think Plants Vs. Zombies Heroes is an underrated game.

It's not lacking in the lore (the Earth is canonically hollow), but it's even richer in meta-lore - you can trace the downfall of PopCap through its update history and concept art, and also see it clearly though the extremely buggy state in which they released the latest update (after five years of silence, overturning years of optimization at the highest level of play).

Anyway, there's this card called Going Viral, which buffs your zombies (you can play as zombies) and shuffles two copies of itself into your deck when played. Combined with cost-reduction effects from Zombology Teacher, you see a lot of hackers running 20x Going Viral + 20x Zombology Teacher decks to secure cheap wins.

Going Viral's play animation

Footage captured by ResearcherFerd for the PvZ wiki

But of course, if your deck gets flooded with Going Virals, you won't be able to draw any Teachers, so you won't have any zombies to buff. And then you will lose.

So, basically, if you're a hacker, you want to draw exactly four Zombology Teachers - the maximum playable on field at one time - and then nothing but Going Viral, to close out the game.


The Puzzle: Jeju's Cursed Cauldron

A witch doctor's cauldron contains 5 blessing orbs and 5 curse orbs. Jeju wants to retrieve the 5 blessing orbs. When reaching into the cauldron he has a uniform probability of retrieving any of the orbs in the cauldron.

When retrieving a blessing orb, the orb will bless him and vanish.
When retrieving a curse orb, the orb will curse him, put two copies of itself in the cauldron, and vanish.

How many times should Jeju expect to be cursed until he retrieves the kk-th blessing for 1k51\leq k\leq5?


The Dipping In of the Fingertips

Letting E[Xk]\mathbb{E}[X_k] be the expected number of curses until the kk-th blessing, we have a straightforward expression for E[X1]\mathbb{E}[X_1]:

E[X1]=x=0xP(X1=x)=x=0x(510+xk=0x15+k10+k)\mathbb{E}[X_1]=\sum_{x=0}^{\infty}x\mathbb{P}(X_1=x)=\sum_{x=0}^{\infty}x(\frac{5}{10+x}\prod_{k=0}^{x-1}\frac{5+k}{10+k})

using the fact that after kk failures the cauldron contains 55 blessings and 5+k5+k curses for a total of 10+k10+k orbs.

The product in the right side (k=0x15+k10+k\prod_{k=0}^{x-1}\frac{5+k}{10+k}) telescopes down to

(56789)(6+s)(7+s)(8+s)(9+s)(10+s)s=(x1)\frac{(5\cdot6\cdot7\cdot8\cdot9)}{(6+s)(7+s)(8+s)(9+s)(10+s)}\Big\rvert _{s=(x-1)}

meaning that the entire right side simplifies to

x=05(56789)x(5+x)(6+x)(7+x)(8+x)(9+x)(10+x)\htmlStyle{color: var(--katex-purple)}{\sum_{x=0}^{\infty}\frac{\htmlStyle{color: var(--katex-base)}{5(5\cdot6\cdot7\cdot8\cdot9)}x}{(5+x)(6+x)(7+x)(8+x)(9+x)(10+x)}}

It's at this point that I briefly got stuck, because this doesn't immediately telescope well. But the official Mathematics Discord weighed in with a beautiful trick - we can express the numerator - that is, xx - in terms of things that telescope:

x=0x(5+x)...(10+x)\htmlStyle{color: var(--katex-purple)}{\sum_{x=0}^{\infty}\frac{x}{(5+x)...(10+x)}}

=x=0(5+x(5+x)...(10+x)5(5+x)...(10+x))=\sum_{x=0}^{\infty}(\frac{5+x}{(5+x)...(10+x)}-\frac{5}{(5+x)...(10+x)})

=x=01(6+x)...(10+x)5x=01(5+x)...(10+x)=\htmlStyle{color: var(--katex-red)}{\sum_{x=0}^{\infty}\frac{1}{(6+x)...(10+x)}}-\htmlStyle{color: var(--katex-blue)}{5\sum_{x=0}^{\infty}\frac{1}{(5+x)...(10+x)}}

where the left term telescopes as

x=01(6+x)...(10+x)=x=014(10+x)(6+x)(6+x)...(10+x)\htmlStyle{color: var(--katex-red)}{\sum_{x=0}^{\infty}\frac{1}{(6+x)...(10+x)}}=\sum_{x=0}^{\infty}\frac{1}{4}\frac{(10+x)-(6+x)}{(6+x)...(10+x)}

=14x=0(1(6+x)...(9+x)1(7+x)...(10+x))=\frac{1}{4}\sum_{x=0}^{\infty}(\frac{1}{(6+x)...(9+x)}-\frac{1}{(7+x)...(10+x)})

=14(x=01(6+x)...(9+x)x=11(6+x)...(9+x))=\frac{1}{4}(\sum_{x=0}^{\infty}\frac{1}{(6+x)...(9+x)}-\sum_{x=1}^{\infty}\frac{1}{(6+x)...(9+x)})

=146789=\htmlStyle{color: var(--katex-red)}{\frac{1}{4\cdot6\cdot7\cdot8\cdot9}}

and the right term telescopes as

5x=01(5+x)...(10+x)=x=055(10+x)(5+x)(5+x)...(10+x)\htmlStyle{color: var(--katex-blue)}{5\sum_{x=0}^{\infty}\frac{1}{(5+x)...(10+x)}}=\sum_{x=0}^{\infty}\frac{5}{5}\frac{(10+x)-(5+x)}{(5+x)...(10+x)}

=x=0(1(5+x)...(9+x)1(6+x)...(10+x))=\sum_{x=0}^{\infty}(\frac{1}{(5+x)...(9+x)}-\frac{1}{(6+x)...(10+x)})

=(x=01(5+x)...(9+x)x=11(5+x)...(9+x))=(\sum_{x=0}^{\infty}\frac{1}{(5+x)...(9+x)}-\sum_{x=1}^{\infty}\frac{1}{(5+x)...(9+x)})

=156789=\htmlStyle{color: var(--katex-blue)}{\frac{1}{5\cdot6\cdot7\cdot8\cdot9}}

condensing the original expression down to

5(56789)(146789156789)5(5\cdot6\cdot7\cdot8\cdot9)(\htmlStyle{color: var(--katex-red)}{\frac{1}{4\cdot6\cdot7\cdot8\cdot9}}-\htmlStyle{color: var(--katex-blue)}{\frac{1}{5\cdot6\cdot7\cdot8\cdot9}})

=5(541)=54=5(\frac{5}{4}-1)=\frac{5}{4}

As a sanity check, if each curse orb added only one copy of itself to the cauldron when drawn, Jeju should expect to be cursed E[Geom(510)]=1\mathbb{E}[\text{Geom}(\frac{5}{10})]=1 time(s) before being blessed. Each curse orb adding two copies of itself to the cauldron should only be a little bit worse than that. So it checks out!

That's a pretty good start, isn't it?


Crunching Nested Telescopes

Seeing the success of that first step, the next thing I did was carry that approach forward. Letting Δxk=xkxk1\Delta x_k = x_k - x_{k-1}:

E[X2]=x1=0E[X2X1=x1]P(X1=x1)\mathbb{E}[X_2]=\sum_{x_1=0}^{\infty}\mathbb{E}[X_2|X_1=x_1]\mathbb{P}(X_1=x_1)

=x1=0(Δx2=0Δx2P(X2=x1+Δx2X1=x1))P(X1=x1)=\sum_{x_1=0}^{\infty}(\sum_{\Delta x_2=0}^{\infty}\Delta x_2\mathbb{P}(X_2=x_1+\Delta x_2|X_1=x_1))\mathbb{P}(X_1=x_1)

=x1=0(Δx2=0Δx2410+x1+kk=0Δx214+x1+k10+x1+k)5(56789)(5+x1)...(10+x1)=\sum_{x_1=0}^{\infty}(\sum_{\Delta x_2=0}^{\infty}\Delta x_2\frac{4}{10+x_1+k}\prod_{k=0}^{\Delta x_2-1}\frac{4+x_1+k}{10+x_1+k})\frac{5(5\cdot6\cdot7\cdot8\cdot9)}{(5+x_1)...(10+x_1)}

This is starting to look hairy. Maybe it's a better idea to pause for now, and look for a better approach.

In other news, here's a picture of the scrap paper on which I crunched the nested telescopes.

Picture of one entire sheet of calculations ending in 25/12

oof ow my wrists

In the top right, there's a Δx2=2512\Delta x_2=\frac{25}{12}, meaning that E[X2]=54+2512=103\mathbb{E}[X_2]=\frac{5}{4}+\frac{25}{12}=\frac{10}{3}.

Still... Let's try not to do any more of that.


And Then There Were Nine (Capital Greek Letters)

Instead of doing any more layers of these nested telescopes, let's just do all the layers in one go by solving the expectation in generality. Let e(b,c)e(b,c) be the expected number of curses drawn before drawing 1 blessing from a cauldron starting with bb blessings and cc curses. Then E[Xk]=c=0e(6k,c+5)P(Xk1=c)\mathbb{E}[X_k]=\sum_{c=0}^{\infty}e(6-k,c+5)\mathbb{P}(X_{k-1}=c), and

e(b,c)=n=0nbb+c+nk=0n1c+kb+c+ke(b,c)=\sum_{n=0}^{\infty}n\frac{b}{b+c+n}\prod_{k=0}^{n-1}\frac{c+k}{b+c+k}

=n=0nbb+c+nk=0b1c+kk=n1(b1)n1b+c+k=\sum_{n=0}^{\infty}n\frac{b}{b+c+n}\frac{\prod_{k=0}^{b-1}c+k}{\prod_{k=n-1-\left(b-1\right)}^{n-1}b+c+k}

=n=0nbk=0b1c+kk=nbnb+c+k=\sum_{n=0}^{\infty}n\frac{b\prod_{k=0}^{b-1}c+k}{\prod_{k=n-b}^{n}b+c+k}

=(bk=0b1(c+k))(n=0nk=nbnb+c+k)=(b\prod_{k=0}^{b-1}(c+k))(\htmlStyle{color: var(--katex-purple)}{\sum_{n=0}^{\infty}\frac{n}{\prod_{k=n-b}^{n}b+c+k}})

So far so good. Now we do the fractional decomposition

n=0nk=nbnb+c+k=n=0(b+c+n)(b+c)k=nbnb+c+k\htmlStyle{color: var(--katex-purple)}{\sum_{n=0}^{\infty}\frac{n}{\prod_{k=n-b}^{n}b+c+k}}=\sum_{n=0}^{\infty}\frac{(b+c+n)-(b+c)}{\prod_{k=n-b}^{n}b+c+k}

=n=0((b+c+n)k=nbnb+c+k(b+c)k=nbnb+c+k)=\sum_{n=0}^{\infty}(\frac{(b+c+n)}{\prod_{k=n-b}^{n}b+c+k}-\frac{(b+c)}{\prod_{k=n-b}^{n}b+c+k})

=n=01k=nbn1b+c+k(b+c)n=01k=nbnb+c+k=\htmlStyle{color: var(--katex-red)}{\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-1}b+c+k}}-\htmlStyle{color: var(--katex-blue)}{(b+c)\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n}b+c+k}}

On the left:

n=01k=nbn1b+c+k=n=01b1((b+c+n1)(c+n)k=nbn1b+c+k)\htmlStyle{color: var(--katex-red)}{\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-1}b+c+k}}=\sum_{n=0}^{\infty}\frac{1}{b-1}(\frac{(b+c+n-1)-(c+n)}{\prod_{k=n-b}^{n-1}b+c+k})

=1b1(n=01k=nbn2b+c+kn=01k=nb+1n1b+c+k)=\frac{1}{b-1}(\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-2}b+c+k}-\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b+1}^{n-1}b+c+k})

=1b1(n=01k=nbn2b+c+kn=11k=nbn2b+c+k)=\frac{1}{b-1}(\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-2}b+c+k}-\sum_{n=1}^{\infty}\frac{1}{\prod_{k=n-b}^{n-2}b+c+k})

=1b1(1k=b2b+c+k)=\frac{1}{b-1}(\frac{1}{\prod_{k=-b}^{-2}b+c+k})

=1b1(1k=0b2c+k)=1b1(c+b1k=0b1c+k)=\frac{1}{b-1}(\frac{1}{\prod_{k=0}^{b-2}c+k})=\htmlStyle{color: var(--katex-red)}{\frac{1}{b-1}(\frac{c+b-1}{\prod_{k=0}^{b-1}c+k})}

On the right:

(b+c)n=01k=nbnb+c+k=n=0(b+c)b((b+c+n)(c+n)k=nbnb+c+k)\htmlStyle{color: var(--katex-blue)}{(b+c)\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n}b+c+k}}=\sum_{n=0}^{\infty}\frac{(b+c)}{b}(\frac{(b+c+n)-(c+n)}{\prod_{k=n-b}^{n}b+c+k})

=(b+c)b(n=01k=nbn1b+c+kn=01k=nb+1nb+c+k)=\frac{(b+c)}{b}(\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-1}b+c+k}-\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b+1}^{n}b+c+k})

=(b+c)b(n=01k=nbn1b+c+kn=11k=nbn1b+c+k)=\frac{(b+c)}{b}(\sum_{n=0}^{\infty}\frac{1}{\prod_{k=n-b}^{n-1}b+c+k}-\sum_{n=1}^{\infty}\frac{1}{\prod_{k=n-b}^{n-1}b+c+k})

=(b+c)b(1k=b1b+c+k)=(b+c)b(1k=0b1c+k)=\frac{(b+c)}{b}(\frac{1}{\prod_{k=-b}^{-1}b+c+k})=\htmlStyle{color: var(--katex-blue)}{\frac{(b+c)}{b}(\frac{1}{\prod_{k=0}^{b-1}c+k})}

The sum of these two terms is

=1b1(c+b1k=0b1c+k)(b+c)b(1k=0b1c+k)=\htmlStyle{color: var(--katex-red)}{\frac{1}{b-1}(\frac{c+b-1}{\prod_{k=0}^{b-1}c+k})}-\htmlStyle{color: var(--katex-blue)}{\frac{(b+c)}{b}(\frac{1}{\prod_{k=0}^{b-1}c+k})}

=(1k=0b1c+k)(c+b1b1b+cb)=(\frac{1}{\prod_{k=0}^{b-1}c+k})(\frac{c+b-1}{b-1}-\frac{b+c}{b})

=(1k=0b1c+k)(b(c+b1)(b+c)(b1)b(b1))=\htmlStyle{color: var(--katex-purple)}{(\frac{1}{\prod_{k=0}^{b-1}c+k})(\frac{b(c+b-1)-(b+c)(b-1)}{b(b-1)})}

We bring back the bk=0b1(c+k)b\prod_{k=0}^{b-1}(c+k) coefficient from earlier to yield

(bk=0b1(c+k))(1k=0b1c+k)(b(c+b1)(b+c)(b1)b(b1))(b\prod_{k=0}^{b-1}(c+k))\htmlStyle{color: var(--katex-purple)}{(\frac{1}{\prod_{k=0}^{b-1}c+k})(\frac{b(c+b-1)-(b+c)(b-1)}{b(b-1)})}

=(bc+b2b)(b2+cbbc)b1=cb1=\frac{(bc+b^2-b)-(b^2+cb-b-c)}{b-1}=\frac{c}{b-1}

which is surprisingly simple - and, even more surprisingly, linear in cc. So we can derive, using linearity of expectation, that

E[Xk]=xk1=0E[XkXk1=xk1]P(Xk1=xk1)\mathbb{E}[X_k]=\sum_{x_{k-1}=0}^{\infty}\mathbb{E}[X_k\rvert X_{k-1}=x_{k-1}]\mathbb{P}(X_{k-1}=x_{k-1})

=xk1=0(xk1+E[ΔxkXk1=xk1])P(Xk1=xk1)=\sum_{x_{k-1}=0}^{\infty}(x_{k-1}+\mathbb{E}[\Delta x_k\rvert X_{k-1}=x_{k-1}])\mathbb{P}(X_{k-1}=x_{k-1})

=xk1=0(xk1+e(6k,5+xk1))P(Xk1=xk1)=\sum_{x_{k-1}=0}^{\infty}(x_{k-1}+e(6-k, 5+x_{k-1}))\mathbb{P}(X_{k-1}=x_{k-1})

=E[Xk1]+E[e(6k,5+Xk1)]=\mathbb{E}[X_{k-1}]+\mathbb{E}[e(6-k, 5+X_{k-1})]

=E[Xk1]+e(6k,5+E[Xk1])=\mathbb{E}[X_{k-1}]+e(6-k, 5+\mathbb{E}[X_{k-1}])

=E[Xk1]+5+E[Xk1]6k1=\mathbb{E}[X_{k-1}]+\frac{5+\mathbb{E}[X_{k-1}]}{6-k-1}

=E[Xk1]+5+E[Xk1]5k=\mathbb{E}[X_{k-1}]+\frac{5+\mathbb{E}[X_{k-1}]}{5-k}

(or if you prefer, 5+(6k)E[Xk1]5k\frac{5+(6-k)\mathbb{E}[X_{k-1}]}{5-k}.)

This matches up with what we crunched earlier: taking E[X0]\mathbb{E}[X_0] to be 0, we see that E[X1]=5+5(0)4=54\mathbb{E}[X_1]=\frac{5+5(0)}{4}=\frac{5}{4} and E[X2]=5+4543=103\mathbb{E}[X_2]=\frac{5+4\frac{5}{4}}{3}=\frac{10}{3}, as we would hope. And as an added bonus, the denominator of E[Xk]\mathbb{E}[X_k] (if present) looks to be always equal to 5k5-k. ^_^

So now it's extremely easy to get values for E[X3]=152\mathbb{E}[X_3]=\frac{15}{2} and E[X4]=20\mathbb{E}[X_4]=20.

And then we hit k=5k=5, and the denominator becomes 0.


In Retrospect, The Sandpiles Were Suspiciously High

Unfortunately for Jeju, this isn't a mistake: E[X5]=\mathbb{E}[X_5]=\infty. To get an intuition as to why, let's consider a cauldron starting with 1 blessing and 1 curse. Then

E[X1]=x1=0x1P(X1=x1)=x1=0x1(12+x1k=0x111+k2+k)\mathbb{E}[X^\prime_1]=\sum_{x^\prime_1=0}^{\infty}{x^\prime_1}\mathbb{P}(X^\prime_1=x^\prime_1)=\sum_{x^\prime_1=0}^{\infty}x^\prime_1(\frac{1}{2+x^\prime_1}\prod_{k=0}^{x^\prime_1-1}\frac{1+k}{2+k})

=x1=0x12+x112+(x11)=x1=0x1(1+x1)(2+x1)=\sum_{x^\prime_1=0}^{\infty}\frac{x^\prime_1}{2+x^\prime_1}\frac{1}{2+(x^\prime_1-1)}=\sum_{x^\prime_1=0}^{\infty}\frac{x^\prime_1}{(1+x^\prime_1)(2+x^\prime_1)}

the partial sums of which x1=0n\sum_{x^\prime_1=0}^{n} are clearly growing according to the harmonic numbers Hn=x=1n1xH_n=\sum_{x=1}^{n}\frac{1}{x}, so this infinite summation diverges. Unlucky! Infinitely many curses for you!

GIF of sand pouring down on a man's head with the caption 'the pharaoh's curse'

Taken from KnowYourMeme


The Wind Blows From East to West

I originally came up with this question on a whim and asked it to DeepSeek, who proceeded to pretend it knew what it was on about for about ten minutes before giving me a partially correct answer.

DeepSeek correctly computes E[X1] but fails on E[X5]

DeepSeek's best effort

In the process, it brought up something called the Pólya urn.

The Pólya urn is a "rich get richer" kind of model. It does the opposite of "drawing balls without replacement": when you draw a ball from the urn, that ball is replaced twice over for the next draw. Its equilibrium state is unstable and is almost surely not converged to (for some physics-adjacent notion of equilibrium, at least).

In fact, letting the two colors of balls in the urn be red and blue, the proportion of red balls in the urn converges to a random variable on [0,1][1].

The base Pólya urn is (from what I've read, I didn't take the martingales course) closely related to the beta distribution and its multivariate extension, the Dirichlet distribution[2]. Using these, we can answer questions about things like the stopping times for ball-drawing processes on Pólya urns containing more than two colors of balls.

But Jeju's cauldron only follows Pólya's model on the curse orbs. It draws the blessing orbs without replacement. Hence, it's only semi-Pólya.

For the similar mixed-replacement urn (where blessings vanish but curses are replaced), Wikipedia makes note of the recursion for P(n,k)=P(k blessings in n draws)P(n,k)=\mathbb{P}(k\ \text{blessings in }n\ \text{draws})[3],

P(n,k)=b(k1)c+b(k1)P(n1,k1)+cc+bkP(n1,k)P(n,k)=\frac{b-(k-1)}{c+b-(k-1)}P(n-1,k-1)+\frac{c}{c+b-k}P(n-1,k)

For our semi-Pólya urn, since the number of curses in the cauldron after nn draws should be equal to c+(nk)c+(n-k), we would have something more like

P(n,k)=(b(k1))P(n1,k1)c+(n1(k1))+b(k1)+(c+(n1k))P(n1,k)(c+(n1k))+bkP(n,k)=\frac{(b-(k-1))P(n-1,k-1)}{c+(n-1-(k-1))+b-(k-1)}+\frac{(c+(n-1-k))P(n-1,k)}{(c+(n-1-k))+b-k}

=(b(k1))P(n1,k1)c+b+n12k+(c+(n1k))P(n1,k)c+b+n12k=\frac{(b-(k-1))P(n-1,k-1)}{c+b+n-1-2k}+\frac{(c+(n-1-k))P(n-1,k)}{c+b+n-1-2k}

I'll think I'll leave it as an exercise to check whether or not P(n,k)=0P(n,k)=0 for all k>bk>b for both recurrences.

In any case, for the generalized semi-Pólya urn where colors b0...bnb_0...b_n drawn without replacement have counts x0,t...xn,tx_{0,t}...x_{n,t} after draw tt and colors c0...cnc_0...c_n drawn with Pólya duplication have counts y0,t...ym,ty_{0,t}...y_{m,t} after draw tt, we would expect that

Xt=ixi,t and Yt=iyi,tX_t=\sum_i x_{i,t}\text{ and }Y_t=\sum_i y_{i,t}

should evolve according to the two-color semi-Pólya urn, and within those color groups we should see evolution according to their respective multivariate generalizations (the multivariate hypergeometric and multinomial Dirichlet distributions respectively; and yes, "multinomial Dirichlet distribution" for Dirichlet distribution being "multivariate beta distribution" is indeed "multinomial multivariate beta distribution". Phew.)

In conclusion, regardless of the exact distribution of blessing types, we should expect to never (in no finite amount of time) fully empty the cauldron of its blessings. If you spin it juuust right, maybe you could convince someone that that's a good thing and fleece them for a respectable sum.

Well, you can't sell a cursed cauldron that you don't have. Don't ask me to put you in contact with any witch doctors that might make one for you. I only know one, and he's not taking commissions at the moment.


Further Questions

Telescope-crunching tricks. Is there any way to improve on mdoern partial fraction algorithms using telescope-crunching? My instinct says no - solving a system of linear equations probably perfectly mirrors the process of isolating telescopic terms. But maybe we can say something about the special case of Pochhammer denominators?

The finite state machurn. A Turing machine may behave differently when it reads the same tape symbol while in different internal states. Hence, the finite state machurn, and its associated mapping function δ::Q×CNC×Q\delta::Q\times C\rightarrow\mathbb{N}^{|C|}\times Q, where pulling a ball of color cCc\in C while in state qQq\in Q adds back in the specified count and colors' worth of balls and transitions to state qQq'\in Q. This (from a cursory search) seems like it would add a notion of convergence to probabilistic automata.


  1. Daniel Ahlberg's notes on Pólya's urn ↩︎

  2. Wikipedia page for Pólya's urn ↩︎

  3. Wikipedia page for urn problems ↩︎