Detailed mathematical approach

Problem Summary

The implementations evaluate a function \(f(a,b)\) for the asymmetric random walk that moves \(+b\) or \(-a\) with probability \(1/2\) at each step. For the two concrete pairs that appear in the implementations, \((1,2)\) for validation and \((89,97)\) for the final answer, the numbers \(a\) and \(b\) are coprime. In that setting, every positive return to the origin happens after a multiple of \(a+b\) steps.

The quantity being computed is not a raw infinite sum for its own sake. That sum is the total return kernel at the origin, and its reciprocal is the probability that the walk never returns to its starting point. The whole problem is therefore about converting return probabilities into a rapidly convergent series that can be evaluated with high-precision arithmetic.

Mathematical Approach

The code reveals a very specific structure: it keeps a running series whose \(n\)-th nontrivial term is the probability of being back at 0 after \((a+b)n\) steps. Once that interpretation is made explicit, the rest of the derivation is standard but elegant.

Return times are forced by step balance

Let the walk start at 0. After \(R\) jumps of size \(+b\) and \(L\) jumps of size \(-a\), the position is

$bR-aL.$

A return to the origin therefore requires

$bR=aL.$

For the coprime inputs used here, this forces

$R=an,\qquad L=bn$

for some integer \(n\ge 0\). Hence every return time is of the form

$N=(a+b)n.$

This is the key invariant of the problem: the walk can only revisit 0 when the numbers of right and left steps occur in the fixed ratio \(a:b\).

Closed form for the return probability

Let

$u_n=\Pr(X_{(a+b)n}=0),\qquad u_0=1.$

To be back at 0 after \((a+b)n\) steps, the walk must contain exactly \(an\) right jumps and \(bn\) left jumps. Every sequence of that type has probability \(2^{-(a+b)n}\), and the number of such sequences is \(\binom{(a+b)n}{an}\). Therefore

$u_n=\frac{1}{2^{(a+b)n}}\binom{(a+b)n}{an}.$

The implementations do not recompute this from factorials every time. Instead they use the ratio

$\frac{u_n}{u_{n-1}}=\frac{\prod_{i=1}^{a+b}\bigl((a+b)(n-1)+i\bigr)}{2^{a+b}\left(\prod_{i=1}^{a}\bigl(a(n-1)+i\bigr)\right)\left(\prod_{i=1}^{b}\bigl(b(n-1)+i\bigr)\right)}.$

This is exactly the multiplicative update visible in all three language versions.

From return probabilities to escape probability

Let \(r_n\) be the probability that the first return to 0 occurs at time \((a+b)n\). Then every return can be decomposed by the time of the first return:

$u_n=\sum_{k=1}^{n} r_k\,u_{n-k}\qquad (n\ge 1).$

Introduce generating functions

$U(z)=\sum_{n\ge 0}u_n z^n,\qquad R(z)=\sum_{n\ge 1}r_n z^n.$

The convolution above gives

$U(z)=1+R(z)U(z),\qquad\text{so}\qquad U(z)=\frac{1}{1-R(z)}.$

Setting \(z=1\) yields

$U(1)=\sum_{n\ge 0}u_n,\qquad R(1)=1-\frac{1}{U(1)}.$

Here \(R(1)\) is the total probability of ever coming back to the origin, so the probability of never returning is

$f(a,b)=1-R(1)=\frac{1}{\sum_{n\ge 0}u_n}=\frac{1}{\sum_{n\ge 0}\binom{(a+b)n}{an}\,2^{-(a+b)n}}.$

This reciprocal is exactly what the implementations return.

Why the symmetric case collapses to zero

If \(a=b\), then

$u_n=\frac{1}{2^{2an}}\binom{2an}{an}\sim \frac{1}{\sqrt{\pi an}}.$

The series \(\sum u_n\) diverges, which means the walk returns to the origin with probability 1. Consequently

$f(a,a)=0.$

When \(a\ne b\), Stirling's formula gives

$u_n\sim \frac{C}{\sqrt{n}}\lambda^n,\qquad \lambda=\frac{(a+b)^{a+b}}{2^{a+b}a^a b^b}.$

Weighted AM-GM shows \(\lambda<1\) unless \(a=b\), so the asymmetric case converges exponentially fast. That is why a strict cutoff such as \(10^{-80}\) is more than enough for a 9-decimal final answer.

Worked example: \(a=1,\ b=2\)

The first positive return can only happen after 3 steps, with one \(+2\) move and two \(-1\) moves. Hence

$u_1=\frac{\binom{3}{1}}{2^3}=\frac{3}{8}.$

After 6 steps, the walk must contain two \(+2\) moves and four \(-1\) moves, so

$u_2=\frac{\binom{6}{2}}{2^6}=\frac{15}{64}.$

Therefore

$U(1)=1+\frac{3}{8}+\frac{15}{64}+\frac{\binom{9}{3}}{2^9}+\cdots,$

and the escape probability is

$f(1,2)=\frac{1}{U(1)}\approx 0.427050983124842272.$

This is the same check value used by the implementations before evaluating the target pair \((89,97)\).

How the Code Works

Build the term recurrence instead of giant factorials

The C++, Python, and Java implementations begin with the trivial symmetric case \(a=b\), where the answer is immediately 0. Otherwise they precompute \(2^{a+b}\) once, initialize the current term to \(u_0=1\), and then update each new \(u_n\) multiplicatively from \(u_{n-1}\) using the ratio above. This avoids repeatedly forming enormous factorials or binomial coefficients.

Accumulate the total return series

A running sum starts at 1, corresponding to the walk being at the origin at time 0. Each loop iteration multiplies by the next ratio, adds the new term to the total, and thereby builds

$U(1)=\sum_{n\ge 0}u_n.$

At the end, the implementation returns \(1/U(1)\), which is the escape probability derived above.

Use high precision and stop when the tail is negligible

All three implementations use about 100 decimal digits of working precision. They stop when the newest term drops below \(10^{-80}\), leaving a very comfortable margin before the final 9-digit rounding. The two built-in checks also reflect the mathematics: the symmetric case \((1,1)\) must give 0, and the asymmetric case \((1,2)\) must reproduce the known probability above.

Complexity Analysis

If \(L\) terms are needed before the threshold is reached, each term requires \(O(a+b)\) high-precision multiplications and divisions, because the ratio contains \(a+b\) numerator factors and \(a+b\) denominator factors in total. The running time is therefore

$O(L(a+b)).$

The extra memory usage is \(O(1)\): the implementation stores only the current term, the cumulative sum, the fixed power \(2^{a+b}\), and a few scalar temporaries. For the target input the convergence is rapid, so the practical cost is very small.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=959
2. Random walk: Wikipedia - Random walk
3. Binomial distribution: Wikipedia - Binomial distribution
4. Renewal theory: Wikipedia - Renewal theory
5. Stirling's approximation: Wikipedia - Stirling's approximation
6. Arbitrary-precision arithmetic: Wikipedia - Arbitrary-precision arithmetic