Problem 892: Zebra Circles

Project Euler Problem 892 Solution

EulerSolve provides an optimized solution for Project Euler Problem 892, Zebra Circles, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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

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

Problem 892 source code

C++

#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <vector>

namespace {

constexpr std::int64_t kMod = 1234567891LL;

inline std::int64_t mul_mod(std::int64_t a, std::int64_t b) {
    return static_cast<std::int64_t>((__int128)a * b % kMod);
}

std::vector<int> build_inverses(int limit) {
    std::vector<int> inv(limit + 1, 0);
    if (limit >= 1) {
        inv[1] = 1;
    }
    for (int i = 2; i <= limit; ++i) {
        inv[i] = static_cast<int>(kMod - (kMod / i) * 1LL * inv[kMod % i] % kMod);
    }
    return inv;
}

std::int64_t next_d(int n, std::int64_t d_n, std::int64_t d_n1,
                    const std::vector<int>& inv) {
    const std::int64_t n2 = (static_cast<std::int64_t>(n) * n) % kMod;
    const std::int64_t n1 = n + 1;

    std::int64_t term1 = mul_mod(n2, n1);
    term1 = mul_mod(term1, 2LL * n + 3);
    term1 = mul_mod(term1, 16);
    term1 = mul_mod(term1, d_n);

    std::int64_t quad = (2LL * n2 + 4LL * n + 3) % kMod;
    std::int64_t term2 = mul_mod(quad, d_n1);
    term2 = mul_mod(term2, n1);
    term2 = mul_mod(term2, 4);

    std::int64_t numerator = (term1 + term2) % kMod;

    std::int64_t denom_inv = mul_mod(inv[n], inv[n + 2]);
    denom_inv = mul_mod(denom_inv, inv[n + 3]);
    denom_inv = mul_mod(denom_inv, inv[2 * n + 1]);

    return mul_mod(numerator, denom_inv);
}

void validate(std::int64_t d3, std::int64_t d100) {
    if (d3 != 4) {
        std::cerr << "Validation failed for D(3).\n";
        std::exit(1);
    }
    if (d100 != 1172122931LL) {
        std::cerr << "Validation failed for D(100) modulo " << kMod << ".\n";
        std::exit(1);
    }
}

}  // namespace

int main(int argc, char** argv) {
    int target = 10'000'000;
    bool do_validate = true;
    bool target_set = false;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "--no-validate") {
            do_validate = false;
            continue;
        }
        char* end = nullptr;
        long long val = std::strtoll(arg.c_str(), &end, 10);
        if (end && *end == '\0' && !target_set) {
            target = static_cast<int>(std::max(0LL, val));
            target_set = true;
        }
    }

    if (target <= 0) {
        std::cout << 0 << '\n';
        return 0;
    }

    int compute_n = target;
    if (do_validate) {
        compute_n = std::max(compute_n, 100);
    }

    const int inv_limit = 2 * compute_n + 3;
    std::vector<int> inv = build_inverses(inv_limit);

    std::int64_t d_prev = 0;  // D(1)
    std::int64_t d_curr = 2;  // D(2)
    std::int64_t total = 0;
    if (target >= 1) total = d_prev;
    if (target >= 2) total = (total + d_curr) % kMod;

    std::int64_t d3 = -1;
    std::int64_t d100 = -1;

    for (int n = 1; n + 2 <= compute_n; ++n) {
        std::int64_t d_next = next_d(n, d_prev, d_curr, inv);
        const int idx = n + 2;
        if (idx == 3) d3 = d_next;
        if (idx == 100) d100 = d_next;
        if (idx <= target) {
            total += d_next;
            if (total >= kMod) total -= kMod;
        }
        d_prev = d_curr;
        d_curr = d_next;
    }

    if (do_validate) {
        validate(d3, d100);
    }

    std::cout << total % kMod << '\n';
    return 0;
}

Python

def solve():
    MOD = 1234567891
    target = 10000000
    inv_limit = 2 * target + 3
    inv = [0] * (inv_limit + 1)
    inv[1] = 1
    for i in range(2, inv_limit + 1):
        inv[i] = MOD - (MOD // i) * inv[MOD % i] % MOD

    def next_d(n, d_n, d_n1):
        n2 = n * n % MOD; n1 = (n + 1) % MOD
        t1 = n2 * n1 % MOD * ((2*n+3) % MOD) % MOD * 16 % MOD * d_n % MOD
        quad = (2*n2 + 4*n + 3) % MOD
        t2 = quad * d_n1 % MOD * n1 % MOD * 4 % MOD
        num = (t1 + t2) % MOD
        di = inv[n] * inv[n+2] % MOD * inv[n+3] % MOD * inv[2*n+1] % MOD
        return num * di % MOD

    dp, dc = 0, 2  # D(1)=0, D(2)=2
    total = dp + dc  # sum D(1)+D(2)
    for n in range(1, target - 1):
        dn = next_d(n, dp, dc)
        idx = n + 2
        if idx <= target:
            total = (total + dn) % MOD
        dp, dc = dc, dn

    return str(total % MOD)

if __name__ == '__main__':
    print(solve())

Java

public class Euler892 {
    static final long kMod = 1234567891L;

    static int[] buildInverses(int limit) {
        int[] inv = new int[limit + 1];
        if (limit >= 1)
            inv[1] = 1;
        for (int i = 2; i <= limit; ++i) {
            inv[i] = (int) (kMod - (kMod / i) * inv[(int) (kMod % i)] % kMod);
        }
        return inv;
    }

    static long nextD(int n, long d_n, long d_n1, int[] inv) {
        long n2 = ((long) n * n) % kMod;
        long n1 = n + 1;

        long term1 = (n2 * n1) % kMod;
        term1 = (term1 * (2L * n + 3)) % kMod;
        term1 = (term1 * 16) % kMod;
        term1 = (term1 * d_n) % kMod;

        long quad = (2L * n2 + 4L * n + 3) % kMod;
        long term2 = (quad * d_n1) % kMod;
        term2 = (term2 * n1) % kMod;
        term2 = (term2 * 4) % kMod;

        long numerator = (term1 + term2) % kMod;

        long denomInv = ((long) inv[n] * inv[n + 2]) % kMod;
        denomInv = (denomInv * inv[n + 3]) % kMod;
        denomInv = (denomInv * inv[2 * n + 1]) % kMod;

        return (numerator * denomInv) % kMod;
    }

    public static String solve() {
        int target = 10000000;
        int invLimit = 2 * target + 3;
        int[] inv = buildInverses(invLimit);

        long dPrev = 0;
        long dCurr = 2;
        long total = 2;

        for (int n = 1; n + 2 <= target; ++n) {
            long dNext = nextD(n, dPrev, dCurr, inv);
            total += dNext;
            if (total >= kMod)
                total -= kMod;
            dPrev = dCurr;
            dCurr = dNext;
        }

        return Long.toString(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}