Detailed mathematical approach

Problem Summary

Problem 892, Zebra Circles, asks for the cumulative total of a sequence \(D(n)\) modulo

$p=1{,}234{,}567{,}891,$

with initial values \(D(1)=0\) and \(D(2)=2\). The required output is

$S(N)=\sum_{k=1}^{N} D(k)\pmod p,\qquad N=10^7.$

The difficulty is that \(N\) is very large and the defining update is a rational recurrence. The successful approach is to keep everything in modular arithmetic, turn each division into multiplication by a precomputed modular inverse, and stream the sequence in one forward pass.

Mathematical Approach

The sequence is generated from the three-term recurrence

$D_{n+2}=\frac{16n^2(n+1)(2n+3)D_n+4(n+1)(2n^2+4n+3)D_{n+1}}{n(n+2)(n+3)(2n+1)}\qquad (n\ge 1).$

So the problem is not to discover \(D(n)\) term by term with exact rational arithmetic, but to evaluate this recurrence efficiently modulo \(p\) and accumulate the total sum.

Step 1: Isolate the algebraic pieces of the recurrence

Write

$A_n=16n^2(n+1)(2n+3),\qquad B_n=4(n+1)(2n^2+4n+3),$

$C_n=n(n+2)(n+3)(2n+1).$

Then the recurrence becomes

$D_{n+2}=\frac{A_nD_n+B_nD_{n+1}}{C_n}.$

This rewriting makes the structure clear: each new term depends only on the previous two terms and on explicit polynomials in \(n\). That immediately suggests a constant-state iteration rather than a large dynamic-programming table.

Step 2: Replace division by modular inversion

Over arithmetic modulo \(p\), division by an invertible value \(x\) is multiplication by its inverse:

$\frac{u}{x}\equiv u\,x^{-1}\pmod p.$

Therefore the recurrence can be evaluated as

$D_{n+2}\equiv \left(A_nD_n+B_nD_{n+1}\right)\cdot C_n^{-1}\pmod p.$

Because the denominator factors as \(n\), \(n+2\), \(n+3\), and \(2n+1\), the implementation only needs inverses for ordinary integers appearing in that range. A safe upper bound is \(2N+3\), which covers every denominator factor used during the loop.

Step 3: Precompute all inverses in linear time

Let \(I(i)\) denote the modular inverse of \(i\) modulo \(p\). The implementation uses

$I(1)=1,$

$I(i)\equiv \left(p-\left\lfloor\frac{p}{i}\right\rfloor\right)I(p\bmod i)\pmod p,\qquad i\ge 2.$

This comes from the Euclidean identity

$p=\left\lfloor\frac{p}{i}\right\rfloor i+(p\bmod i),$

which implies

$(p\bmod i)\equiv -\left\lfloor\frac{p}{i}\right\rfloor i\pmod p.$

After multiplying by the relevant inverses, one gets a recurrence for \(I(i)\) in terms of a smaller argument \(p\bmod i\). This avoids expensive repeated inverse computations inside the main loop.

Step 4: Stream both the sequence and the total

Once the inverse table exists, each new term requires only a fixed number of modular multiplications and additions. The running state consists of the current pair \((D_n,D_{n+1})\) and the partial sum

$S_m=\sum_{k=1}^{m} D(k)\pmod p.$

After computing \(D_{n+2}\), the implementation immediately updates

$S_{n+2}=S_{n+1}+D_{n+2}\pmod p.$

No old term is needed again after it has been shifted out of the two-term window, so the recurrence itself uses constant live state.

Worked Example

The first few terms show how the rational formula collapses to ordinary integers before reduction modulo \(p\).

For \(n=1\), using \(D_1=0\) and \(D_2=2\),

$D_3=\frac{16\cdot 1^2\cdot 2\cdot 5\cdot D_1+4\cdot 2\cdot 9\cdot D_2}{1\cdot 3\cdot 4\cdot 3} =\frac{0+144}{36}=4.$

For \(n=2\),

$D_4=\frac{16\cdot 2^2\cdot 3\cdot 7\cdot D_2+4\cdot 3\cdot 19\cdot D_3}{2\cdot 4\cdot 5\cdot 5} =\frac{2688+912}{200}=18.$

For \(n=3\),

$D_5=\frac{16\cdot 3^2\cdot 4\cdot 9\cdot D_3+4\cdot 4\cdot 33\cdot D_4}{3\cdot 5\cdot 6\cdot 7} =\frac{20736+9504}{630}=48.$

So the sequence begins

$0,\ 2,\ 4,\ 18,\ 48,\dots$

and the partial sums begin

$0,\ 2,\ 6,\ 24,\ 72,\dots$

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They fix the modulus \(p\), the target \(N\), and an inverse-table limit of \(2N+3\). Then they build a table of modular inverses for every integer in that range. After that they initialize the recurrence with the two known starting terms and the initial running sum.

Inside the loop, the implementation evaluates the two polynomial coefficients in the numerator, combines them with the previous two sequence values, multiplies by the product of the four denominator inverses, and reduces modulo \(p\). The new term is added to the running total immediately, and the two-term state window is shifted forward. This design avoids floating-point arithmetic, avoids per-step inverse recomputation, and keeps the recurrence state independent of \(N\).

Complexity Analysis

The inverse table has size \(2N+3\), so building it costs \(O(N)\) time and \(O(N)\) memory. The recurrence loop then performs \(N-2\) iterations, each doing only a constant amount of modular arithmetic, so the loop is also \(O(N)\). The full method is therefore linear in \(N\), with \(O(N)\) memory dominated by the inverse table and only \(O(1)\) additional state for the recurrence and the running sum.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=892
2. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
3. Modular arithmetic: Wikipedia — Modular arithmetic
4. Recurrence relation: Wikipedia — Recurrence relation
5. Euclidean algorithm: Wikipedia — Euclidean algorithm