Detailed mathematical approach

Problem Summary

Each round costs \(t\) units. The payout is \(2^n\) with probability \(2^{-(n+1)}\) for \(n \ge 0\). Starting from bankroll \(B\), we want the probability of being able to keep playing forever, which is the same as never letting the bankroll fall below the amount needed for the next ticket. The bankroll can in principle grow without bound, so the state space is infinite. The implementation therefore replaces the raw infinite recursion by a finite residue model together with a self-similar matrix equation.

Mathematical Approach

The numerical method follows the same probability law in all three implementations. Its key idea is that after shifting the bankroll by \(t-1\), the process becomes translation-invariant in blocks of size \(t-1\).

Step 1: Shift the State and Truncate the Tiny Tail

Let

$d=t-1,\qquad c=B-d.$

Here \(c\) is the bankroll measured above the minimum playable threshold. One round changes \(c\) by

$c' = c + 2^n - t.$

If \(c' \le 0\), the player can no longer buy another ticket, so that path is ruined.

The exact distribution has infinitely many outcomes \(n=0,1,2,\dots\), but the implementation keeps only

$0 \le n \lt H,\qquad H=50.$

The omitted probability mass is

$\sum_{n=H}^{\infty} 2^{-(n+1)} = 2^{-50},$

which is far below the requested \(7\)-decimal accuracy.

Step 2: Decompose the State into Layer and Residue

Every playable shifted bankroll \(c \ge 1\) can be written uniquely as

$c = k d + r + 1,\qquad k \ge 0,\qquad 0 \le r \lt d.$

The integer \(k\) is the layer, and \(r\) is the residue inside that layer. This decomposition is useful because

$k d + r + 1 + 2^n - t = (k-1)d + (r + 2^n).$

So increasing the starting layer by one simply shifts all future layers by one. The detailed geometry inside each layer stays the same.

Define

$u_k(r)=\Pr(\text{eventual ruin} \mid c = k d + r + 1).$

For fixed \(t\), each \(u_k\) is a vector of length \(d\).

Step 3: Build the Transfer Matrix for Higher Layers

Because the process repeats the same pattern from one layer to the next, there is a single \(d \times d\) matrix \(X\) such that

$u_{k+1} = X u_k,\qquad u_k = X^k u_0.$

To derive \(X\), start from layer \(1\), where \(c = d + r + 1\). After outcome \(n\),

$c' = d + r + 1 + 2^n - t = r + 2^n.$

Write this new state as

$c' = j_n(r)\,d + \rho_n(r) + 1,$

with

$j_n(r)=\left\lfloor\frac{r+2^n-1}{d}\right\rfloor,\qquad \rho_n(r)=(r+2^n-1)\bmod d.$

If \(j_n(r)=0\), the process lands in the boundary layer. If \(j_n(r)\ge 1\), it lands \(j_n(r)\) layers above the boundary. Therefore the \(r\)-th row of \(X\) satisfies

$X_{r,*}=\sum_{j_n(r)=0} 2^{-(n+1)} e_{\rho_n(r)}^{\mathsf T}+\sum_{j_n(r)\ge 1} 2^{-(n+1)}\bigl(X^{j_n(r)}\bigr)_{\rho_n(r),*},$

where \(e_{\rho}^{\mathsf T}\) is the row vector with a \(1\) in column \(\rho\) and zeros elsewhere.

This is a matrix fixed-point equation. The implementations iterate it numerically until successive matrices differ by less than \(10^{-18}\), or until the safety iteration cap is reached.

Step 4: Solve the Boundary Layer Separately

Now consider the lowest playable layer, where

$c = r + 1,\qquad 0 \le r \lt d.$

After outcome \(n\),

$c' = r + 1 + 2^n - t = r + 2^n - d.$

If \(c' \le 0\), ruin happens immediately. Otherwise write

$c' = j'_n(r)\,d + \rho'_n(r) + 1,$

where

$j'_n(r)=\left\lfloor\frac{c'-1}{d}\right\rfloor,\qquad \rho'_n(r)=(c'-1)\bmod d.$

Conditioning on the first round gives

$u_0(r)=q_r+\sum_{j'_n(r)=0}2^{-(n+1)}u_0\!\left(\rho'_n(r)\right)+\sum_{j'_n(r)\ge 1}2^{-(n+1)}\bigl[X^{j'_n(r)}u_0\bigr]_{\rho'_n(r)},$

where \(q_r\) is the total probability of immediate ruin from residue \(r\).

All boundary-layer unknowns can therefore be grouped into a finite linear system

$(I-C)u_0=q,$

which the implementations solve by Gaussian elimination with pivoting.

Step 5: Recover the Required Probability

For the given starting bankroll \(B\), first convert it to shifted capital

$c_0=B-d.$

If \(c_0 \le 0\), the answer is \(0\) because the first ticket cannot even be bought. Otherwise write

$c_0 = k d + r + 1.$

The ruin probability is then

$\bigl[X^k u_0\bigr]_r,$

so the desired survival probability is

$\boxed{P_t(B)=1-\bigl[X^k u_0\bigr]_r.}$

Step 6: Worked Example for \(t=2\)

When \(t=2\), we have \(d=1\), so there is only one residue class. The transfer matrix becomes a single scalar \(x\), and the boundary ruin vector becomes a single scalar \(u\).

From layer \(1\), outcome \(n=0\) sends the process to the boundary layer, while \(n\ge 1\) jumps to layer \(2^n-1\). Therefore

$x=\frac12+\sum_{n=1}^{H-1}2^{-(n+1)}x^{2^n-1}.$

For the boundary layer, \(n=0\) is immediate ruin and \(n\ge 1\) sends the process to layer \(2^n-2\), so

$u=\frac12+\sum_{n=1}^{H-1}2^{-(n+1)}x^{2^n-2}u.$

Thus

$P_2(2)=1-u,\qquad P_2(5)=1-x^3u.$

The numerical solution gives

$P_2(2)\approx 0.2522,\qquad P_2(5)\approx 0.6873,$

which matches the validation checkpoints used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical plan. First they precompute the retained payouts \(2^n\) and their probabilities \(2^{-(n+1)}\) for \(0 \le n \lt 50\). Next they enumerate every transition from one reference layer and record only the resulting residue together with the layer jump that must be applied later.

To evaluate many powers of the same transfer matrix efficiently, the implementation builds \(X\), \(X^2\), \(X^4\), and so on, and reconstructs each required \(X^j\) from the binary expansion of \(j\). After the fixed-point iteration for \(X\) stabilizes, the boundary layer is folded into one finite linear system and solved by Gaussian elimination. Finally the initial bankroll is decomposed into layer and residue, the needed matrix power is applied to the boundary vector, and the result is subtracted from \(1\). The Python entry point is only a thin wrapper around the same compiled numerical core.

Complexity Analysis

Let \(d=t-1\). Suppose the reference-layer and boundary transitions involve \(E\) distinct positive layer jumps, and let \(J_{\max}\) be the largest of them. One fixed-point iteration needs the squaring chain up to

$L=\left\lfloor\log_2 J_{\max}\right\rfloor+1,$

so its dense-matrix cost is \(O((L+E)d^3)\). Solving the boundary linear system costs \(O(d^3)\). Once the transfer matrix is known, evaluating the answer for a starting bankroll in layer \(k\) needs \(O(d^3 \log k)\) time by binary exponentiation. Memory usage is \(O((L+E)d^2)\). For the actual target \(t=15\), this is entirely manageable because \(d=14\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=499
2. St. Petersburg paradox: Wikipedia — St. Petersburg paradox
3. Gambler's ruin: Wikipedia — Gambler's ruin
4. Numerical note: truncating the payout law at \(H=50\) leaves tail probability \(2^{-50}\), which is negligible at \(7\)-decimal precision.
5. Implementation note: the fixed-point iteration targets entrywise change below \(10^{-18}\), and the boundary problem is solved as a finite linear system.