Detailed mathematical approach

Problem Summary

For each cycle length \(m\), two walkers move independently on an \(m\)-cycle. In every round each walker chooses one of three actions with equal probability: one step clockwise, no move, or one step counterclockwise. If the chase lasts \(T\) rounds before the walkers meet, the payment for that game is \(T^2\). Let \(E_m\) be the expected payment when the starting positions are chosen uniformly at random, and define

$G(n)=\sum_{m=2}^{n} E_m.$

The task is to evaluate this quantity up to \(n=500\). The implementations do not simulate random games. Instead, they solve the exact expectation equations for the relative-separation walk on the cycle.

Mathematical Approach

The key simplification is that absolute positions do not matter. Only the relative separation on the cycle matters, so the original two-walker process becomes a single Markov chain on \(\mathbb{Z}_m\).

Step 1: Reduce the Chase to Relative Separation

Let \(D_t\in\mathbb{Z}_m\) be the directed separation after \(t\) rounds. The chase ends exactly when

$D_t=0.$

Each walker contributes an increment from \(\{-1,0,+1\}\) with probability \(1/3\), so the relative increment \(X_t=D_{t+1}-D_t\) has distribution

$\mathbb{P}(X_t=0)=\frac13,\qquad \mathbb{P}(X_t=\pm1)=\frac29,\qquad \mathbb{P}(X_t=\pm2)=\frac19.$

Therefore the transition operator acting on a function \(f:\mathbb{Z}_m\to\mathbb{R}\) is

$ (Pf)(d)=\frac13 f(d)+\frac29\bigl(f(d+1)+f(d-1)\bigr)+\frac19\bigl(f(d+2)+f(d-2)\bigr), $

where all indices are understood modulo \(m\).

Step 2: Diagonalize the Transition Operator with Fourier Modes

Because the operator is circulant, the discrete Fourier modes

$u_k(d)=e^{2\pi i k d/m},\qquad k=0,1,\dots,m-1,$

are eigenvectors of \(P\). Writing

$\theta_k=\frac{2\pi k}{m},$

one gets

$ \lambda_k=\frac13+\frac29\left(e^{i\theta_k}+e^{-i\theta_k}\right)+\frac19\left(e^{2i\theta_k}+e^{-2i\theta_k}\right) =\frac13+\frac49\cos\theta_k+\frac29\cos(2\theta_k). $

These eigenvalues are exactly the denominators used by the implementation. The mode \(k=0\) has \(\lambda_0=1\), so constants lie in the kernel of \(I-P\); that is why every spectral solve must be written with zero total source.

Step 3: Solve the First Moment of the Meeting Time

Let \(T_d\) be the number of rounds until the first meeting when the initial separation is \(d\), and define

$a_d=\mathbb{E}[T_d],\qquad a_0=0.$

For every nonzero separation, one round passes and the process restarts from the next state, so

$a_d=1+(Pa)_d\qquad (d\ne0).$

To make the equation compatible with Fourier inversion on the whole cycle, rewrite it as the zero-sum Poisson equation

$ (I-P)a=\mathbf{1}-m e_0, $

where \(e_0\) is \(1\) at separation \(0\) and \(0\) elsewhere. After removing the zero mode and taking the real part, the first moment becomes

$ a_d=\sum_{k=1}^{m-1}\frac{1-\cos\left(2\pi k d/m\right)}{1-\lambda_k}. $

This is the first large spectral sum carried out by the implementations.

Step 4: Turn the Second Moment into Another Poisson Equation

Now define the second moment

$b_d=\mathbb{E}[T_d^2],\qquad b_0=0.$

Using \(T_d=1+T_{D_1}\) after one round, we get for \(d\ne0\)

$ b_d=1+2(Pa)_d+(Pb)_d. $

Since \(a_d=1+(Pa)_d\), this simplifies to

$ (I-P)b=2a-1 $

away from the absorbing state. Again we need a zero-sum source, so define

$ c_d=2a_d-1\quad(d\ne0),\qquad c_0=-\sum_{d=1}^{m-1} c_d. $

Let \(z\) be the unique zero-mean solution of

$ (I-P)z=c,\qquad \sum_{d=0}^{m-1} z_d=0. $

Then \(b_d=z_d-z_0\), which enforces \(b_0=0\). If

$ \widehat{c}_k=\sum_{d=0}^{m-1} c_d\,e^{-2\pi i k d/m}, $

then the average payment for a fixed cycle size is

$ E_m=\frac1m\sum_{d=0}^{m-1} b_d =-z_0 =-\frac1m\sum_{k=1}^{m-1}\frac{\widehat{c}_k}{1-\lambda_k}. $

That is the second spectral solve implemented in all three languages.

Step 5: Sum Over All Cycle Sizes

Once \(E_m\) is available for one cycle length, the global function is just

$G(n)=\sum_{m=2}^{n} E_m.$

The implementations evaluate this directly for \(m=2,3,\dots,500\) and accumulate the running total.

Worked Example: \(m=2\)

For \(m=2\), there is only one nontrivial Fourier mode, \(k=1\), with

$ \lambda_1=\frac13+\frac49\cos\pi+\frac29\cos(2\pi)=\frac19, \qquad 1-\lambda_1=\frac89. $

The first moment from separation \(1\) is therefore

$ a_1=\frac{1-\cos\pi}{8/9}=\frac{2}{8/9}=\frac94. $

Next,

$ c_1=2a_1-1=\frac72, \qquad c_0=-\frac72, $

so the nonzero Fourier coefficient is

$ \widehat{c}_1=c_0-c_1=-7. $

Hence

$ E_2=-\frac12\cdot\frac{-7}{8/9}=\frac{63}{16}=3.9375. $

This is the average of \(T^2\) over the two possible starting separations \(0\) and \(1\). The separation \(0\) contributes \(0\), so the conditional second moment from the nontrivial start is twice as large.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical plan. For each cycle length \(m\), they first compute the spectral factors \(1-\lambda_k\) for all nonzero modes. Then they evaluate the first-moment formula by fixing one mode at a time and generating the powers of \(e^{2\pi i k/m}\) incrementally across all separations, which avoids a fresh trigonometric call inside the inner loop.

After that, the implementation builds the balanced source for the second-moment equation, performs a direct discrete Fourier transform of that source, divides each nonzero mode by the same spectral factor, and reads off the average payment from the zero-coordinate of the zero-mean solution. The imaginary parts cancel by conjugate symmetry, so only the real part is kept numerically. Finally, the values for \(m=2\) through \(m=500\) are accumulated and printed in scientific notation.

Complexity Analysis

For one fixed cycle length \(m\), the implementation performs two explicit mode-by-distance sweeps, so the time cost is \(O(m^2)\). The working memory is \(O(m)\), because only a few arrays of length \(m\) are stored at once. Summing over all cycle sizes from \(2\) to \(n\) gives

$ \sum_{m=2}^{n} O(m^2)=O(n^3), $

with peak memory \(O(n)\). For the required bound \(n=500\), this direct spectral method is easily practical.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=683
2. Circulant matrices: Wikipedia - Circulant matrix
3. Discrete Fourier transform: Wikipedia - Discrete Fourier transform
4. Markov chains: Wikipedia - Markov chain
5. Hitting times: Wikipedia - Hitting time