Problem 763: Amoebas in a 3D Grid

Project Euler Problem 763 Solution

EulerSolve provides an optimized solution for Project Euler Problem 763, Amoebas in a 3D Grid, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem defines a sequence $D(n)$ associated with counting amoeba configurations in a three-dimensional grid. The computational goal used by the implementations is to evaluate $D(10000)$ modulo $10^9$, while several smaller values are checked exactly before the large modular run. A direct enumeration of all relevant 3D shapes would grow far too quickly. The implemented method therefore replaces explicit shape generation with a layered dynamic program whose states encode admissible boundary profiles and whose scalar recurrences extract the final sequence value. Mathematical Approach For exposition, denote the two triangular state families by $A_{n,k}(m)$ and $B_{n,k}(m)$, and denote the scalar sequences by $P(m)$ and $Q(m)$. The desired quantity is $$D(n)=Q(n-1).$$ The implementations never materialize the underlying amoebas one by one. Instead, they count how many partial configurations can end in each admissible profile and then close the whole system with two one-dimensional recurrences. Step 1: A triangular threshold controls when a layer can be active The first structural quantity is the triangular threshold $$T(n)=\frac{(n+1)(n+2)}{2}.$$ This is the smallest size parameter at which the level indexed by $n$ can contribute....

Detailed mathematical approach

Problem Summary

The problem defines a sequence $D(n)$ associated with counting amoeba configurations in a three-dimensional grid. The computational goal used by the implementations is to evaluate $D(10000)$ modulo $10^9$, while several smaller values are checked exactly before the large modular run.

A direct enumeration of all relevant 3D shapes would grow far too quickly. The implemented method therefore replaces explicit shape generation with a layered dynamic program whose states encode admissible boundary profiles and whose scalar recurrences extract the final sequence value.

Mathematical Approach

For exposition, denote the two triangular state families by $A_{n,k}(m)$ and $B_{n,k}(m)$, and denote the scalar sequences by $P(m)$ and $Q(m)$. The desired quantity is

$$D(n)=Q(n-1).$$

The implementations never materialize the underlying amoebas one by one. Instead, they count how many partial configurations can end in each admissible profile and then close the whole system with two one-dimensional recurrences.

Step 1: A triangular threshold controls when a layer can be active

The first structural quantity is the triangular threshold

$$T(n)=\frac{(n+1)(n+2)}{2}.$$

This is the smallest size parameter at which the level indexed by $n$ can contribute. Consequently, for every $n\ge 1$,

$$m<T(n)\implies A_{n,k}(m)=B_{n,k}(m)=0.$$

So when the outer loop is processing a fixed $m$, only levels with $T(n)\le m$ are relevant. Because $T(n)\sim n^2/2$, the largest active level is only $O(\sqrt m)$.

Step 2: Boundary conditions collapse the level $n=0$ into a scalar base stream

Negative sizes are impossible, so every quantity evaluated at $m<0$ is defined to be $0$. At level $n=0$, the two state families do not store a full row of profile data. Instead they reduce to the same scalar stream:

$$A_{0,0}(m)=B_{0,0}(m)=P(m),$$

while

$$A_{0,k}(m)=B_{0,k}(m)=0\qquad (k\ne 0).$$

This identification is essential because the higher-level recurrences refer to $n-1$, so the whole system must know what happens when the index reaches zero. The initial scalar condition is

$$Q(0)=1,$$

and all other missing values are supplied automatically by the zero rules for negative shifts.

Step 3: Every active profile obeys a coupled six-term recurrence

For $1\le k\le n$, define the folded index and the edge multiplier by

$$r(n,k)=\begin{cases} k,&k<n,\\ k-1,&k=n, \end{cases} \qquad \lambda(n,k)=\begin{cases} 1,&k<n,\\ 2,&k=n. \end{cases}$$

Then for every active layer $m\ge T(n)$, the two profile families satisfy

$$\begin{aligned} A_{n,k}(m)=&\,A_{n,k}(m-n-2)+B_{n+1,1}(m-n-3)+A_{n+1,k+1}(m-n-3)\\ &+A_{n-1,r(n,k)}(m-n-1)+B_{n,1}(m-n-2)+A_{n,r(n,k)+1}(m-n-2), \end{aligned}$$

$$\begin{aligned} B_{n,k}(m)=&\,B_{n,k}(m-n-2)+B_{n+1,k+1}(m-n-3)+\lambda(n,k)\,A_{n+1,1}(m-n-3)\\ &+B_{n-1,r(n,k)}(m-n-1)+B_{n,r(n,k)+1}(m-n-2)+\lambda(n,k)\,A_{n,1}(m-n-2). \end{aligned}$$

These are exactly the transitions evaluated by the implementations. The special case $k=n$ is the only point where the coefficient $2$ appears, so the edge of the triangular index set has a slightly different combinatorial weight from the interior.

Step 4: Two scalar recurrences close the dynamic system

Once all $A$ and $B$ states at size $m$ have been filled, the implementations update the scalar layer by

$$P(m)=Q(m-1)+4P(m-2)+2A_{1,1}(m-3)+B_{1,1}(m-3),$$

and, for $m\ge 1$,

$$Q(m)=3Q(m-1)+3P(m-2).$$

The second recurrence produces the sequence whose shifted values are the final answers. Therefore

$$D(n)=Q(n-1).$$

From the dynamic-programming point of view, $P(m)$ is the bridge between the boundary level $n=0$ and the genuinely two-parameter profile states with $n\ge 1$.

Step 5: The same recurrence explains the rolling-memory optimization

Every right-hand side only looks backward by offsets such as $m-n-1$, $m-n-2$, $m-n-3$, $m-1$, $m-2$, and $m-3$. No transition needs the full history of all earlier layers.

$$n_{\max}(M)=\max\{n:T(n)\le M\},$$

then all look-backs for the final target size $M$ fit inside a circular buffer of length $n_{\max}(M)+8$. That is why the program can keep the quadratic-time recurrence while avoiding quadratic storage in $m$.

Worked Example: The first nontrivial values

The initial layers can be computed by hand directly from the recurrence.

At $m=0$, no triangular state is active, so

$$P(0)=0,\qquad Q(0)=1.$$

At $m=1$, still no $n\ge 1$ state survives because $T(1)=3$. Hence

$$P(1)=Q(0)=1,\qquad Q(1)=3Q(0)=3.$$

At $m=2$, the same argument gives

$$P(2)=Q(1)+4P(0)=3,\qquad Q(2)=3Q(1)+3P(0)=9.$$

At $m=3$, the first triangular level appears because $T(1)=3$. For the only index pair $(n,k)=(1,1)$, we have $r(1,1)=0$ and $\lambda(1,1)=2$. Using the boundary rule $A_{0,0}(1)=B_{0,0}(1)=P(1)=1$, every invalid or negatively shifted term vanishes, so

$$A_{1,1}(3)=1,\qquad B_{1,1}(3)=1.$$

The scalar update then yields

$$P(3)=Q(2)+4P(1)=9+4=13,$$

$$Q(3)=3Q(2)+3P(1)=27+3=30.$$

Therefore

$$D(1)=1,\qquad D(2)=3,\qquad D(3)=9,\qquad D(4)=30.$$

This small example shows exactly how the first nonzero triangular state enters when the threshold is met.

How the Code Works

The C++, Python, and Java implementations all evaluate the same recurrence. The compiled solvers allocate two triangular state tables for admissible profile pairs $(n,k)$, together with two one-dimensional arrays for the scalar layer. The outer loop advances the size parameter $m$ from $0$ to the target, clears the current column of the circular history buffer, fills every active profile state, and finally updates the two scalar sequences.

Before the dynamic program starts, the code computes the largest active level $n_{\max}$ from the triangular threshold. That value determines both how many profile rows are needed and how large the circular buffer must be. Because the recurrence reaches one row above the current level, the allocation includes a small safety margin beyond $n_{\max}$.

The C++ implementation is used in two arithmetic modes: exact integer evaluation for small verification checkpoints and modular evaluation for the large final query. The Java implementation keeps the modular variant, which is enough for the published answer. The Python implementation is a thin launcher that builds and runs the compiled solver and then extracts the final numeric output, so the three language versions share the same mathematical core.

The embedded checkpoints are

$$D(10)=44499,\qquad D(20)=9204559704,\qquad D(32)=22037102049132222,$$

and

$$D(100)\equiv 780166455\pmod{10^9}.$$

These values verify that the recurrence has been wired correctly before the computation of $D(10000)\bmod 10^9$.

Complexity Analysis

Let $M$ be the largest size parameter reached by the program, so the final query uses $M=n-1$. The active range satisfies $T(n_{\max})\le M<T(n_{\max}+1)$, hence $n_{\max}=\Theta(\sqrt M)$.

For each $m$, the implementation visits every admissible pair $(n,k)$ with $1\le n\le n_{\max}$ and $1\le k\le n$. The number of profile updates per layer is therefore

$$\sum_{n=1}^{n_{\max}} n=\frac{n_{\max}(n_{\max}+1)}{2}=O(n_{\max}^2),$$

which makes the total running time

$$O(M\,n_{\max}^2)=O(M^2).$$

For memory, the circular buffer length is $O(n_{\max})$, and the total number of stored profile entries is

$$O\!\left(n_{\max}\sum_{n=1}^{n_{\max}} n\right)=O(n_{\max}^3)=O(M^{3/2}).$$

The two scalar sequences add only $O(n_{\max})$ more cells, so they do not change the leading asymptotic term.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=763
Dynamic programming: Wikipedia — Dynamic programming
Triangular number: Wikipedia — Triangular number
Recurrence relation: Wikipedia — Recurrence relation
Circular buffer: Wikipedia — Circular buffer

Problem 763 source code

C++

#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

int threshold(const int n) {
    return (n + 1) * (n + 2) / 2;
}

int max_n_for_m(const int m) {
    int n = 0;
    while (threshold(n + 1) <= m) {
        ++n;
    }
    return n;
}

template <bool UseMod>
u64 solve_a2(const int target_m, const u64 mod_value = 0ULL) {
    const int nmax = max_n_for_m(target_m);
    const int ring = nmax + 8;
    const int nalloc = nmax + 2;

    std::vector<std::vector<u64>> u(nalloc + 1);
    std::vector<std::vector<u64>> v(nalloc + 1);
    for (int n = 0; n <= nalloc; ++n) {
        u[n].assign(static_cast<std::size_t>(n + 2) * static_cast<std::size_t>(ring), 0ULL);
        v[n].assign(static_cast<std::size_t>(n + 2) * static_cast<std::size_t>(ring), 0ULL);
    }

    std::vector<u64> a(static_cast<std::size_t>(ring), 0ULL);
    std::vector<u64> f0(static_cast<std::size_t>(ring), 0ULL);
    a[0] = 1ULL;

    auto norm = [&](const u128 x) -> u64 {
        if constexpr (UseMod) {
            return static_cast<u64>(x % static_cast<u128>(mod_value));
        } else {
            return static_cast<u64>(x);
        }
    };

    auto get_a = [&](const int m) -> u64 {
        if (m < 0) {
            return 0ULL;
        }
        return a[static_cast<std::size_t>(m % ring)];
    };

    auto get_f0 = [&](const int m) -> u64 {
        if (m < 0) {
            return 0ULL;
        }
        return f0[static_cast<std::size_t>(m % ring)];
    };

    auto get_u = [&](const int n, const int k, const int m) -> u64 {
        if (m < 0) {
            return 0ULL;
        }
        if (n == 0) {
            return (k == 0) ? get_f0(m) : 0ULL;
        }
        if (n < 0 || n > nalloc || k < 1 || k > n || m < threshold(n)) {
            return 0ULL;
        }
        const std::size_t idx = static_cast<std::size_t>(k) * static_cast<std::size_t>(ring)
                              + static_cast<std::size_t>(m % ring);
        return u[static_cast<std::size_t>(n)][idx];
    };

    auto get_v = [&](const int n, const int k, const int m) -> u64 {
        if (m < 0) {
            return 0ULL;
        }
        if (n == 0) {
            return (k == 0) ? get_f0(m) : 0ULL;
        }
        if (n < 0 || n > nalloc || k < 1 || k > n || m < threshold(n)) {
            return 0ULL;
        }
        const std::size_t idx = static_cast<std::size_t>(k) * static_cast<std::size_t>(ring)
                              + static_cast<std::size_t>(m % ring);
        return v[static_cast<std::size_t>(n)][idx];
    };

    for (int m = 0; m <= target_m; ++m) {
        const int cur = m % ring;

        for (int n = 1; n <= nmax; ++n) {
            for (int k = 1; k <= n; ++k) {
                const std::size_t idx = static_cast<std::size_t>(k) * static_cast<std::size_t>(ring)
                                      + static_cast<std::size_t>(cur);
                u[static_cast<std::size_t>(n)][idx] = 0ULL;
                v[static_cast<std::size_t>(n)][idx] = 0ULL;
            }
        }

        for (int n = 1; n <= nmax; ++n) {
            if (m < threshold(n)) {
                continue;
            }

            for (int k = 1; k <= n; ++k) {
                const int h = (k < n) ? k : (k - 1);

                const u128 u_sum = static_cast<u128>(get_u(n, k, m - n - 2))
                                 + static_cast<u128>(get_v(n + 1, 1, m - n - 3))
                                 + static_cast<u128>(get_u(n + 1, k + 1, m - n - 3))
                                 + static_cast<u128>(get_u(n - 1, h, m - n - 1))
                                 + static_cast<u128>(get_v(n, 1, m - n - 2))
                                 + static_cast<u128>(get_u(n, h + 1, m - n - 2));

                const u128 v_sum = static_cast<u128>(get_v(n, k, m - n - 2))
                                 + static_cast<u128>(get_v(n + 1, k + 1, m - n - 3))
                                 + static_cast<u128>((k == n) ? 2ULL : 1ULL)
                                       * static_cast<u128>(get_u(n + 1, 1, m - n - 3))
                                 + static_cast<u128>(get_v(n - 1, h, m - n - 1))
                                 + static_cast<u128>(get_v(n, h + 1, m - n - 2))
                                 + static_cast<u128>((k == n) ? 2ULL : 1ULL)
                                       * static_cast<u128>(get_u(n, 1, m - n - 2));

                const std::size_t idx = static_cast<std::size_t>(k) * static_cast<std::size_t>(ring)
                                      + static_cast<std::size_t>(cur);
                u[static_cast<std::size_t>(n)][idx] = norm(u_sum);
                v[static_cast<std::size_t>(n)][idx] = norm(v_sum);
            }
        }

        const u128 f0_sum = static_cast<u128>(get_a(m - 1))
                          + static_cast<u128>(4ULL) * static_cast<u128>(get_f0(m - 2))
                          + static_cast<u128>(2ULL) * static_cast<u128>(get_u(1, 1, m - 3))
                          + static_cast<u128>(get_v(1, 1, m - 3));
        f0[static_cast<std::size_t>(cur)] = norm(f0_sum);

        if (m >= 1) {
            const u128 a_sum = static_cast<u128>(3ULL) * static_cast<u128>(get_a(m - 1))
                             + static_cast<u128>(3ULL) * static_cast<u128>(get_f0(m - 2));
            a[static_cast<std::size_t>(cur)] = norm(a_sum);
        }
    }

    return a[static_cast<std::size_t>(target_m % ring)];
}

u64 D_exact(const int n) {
    assert(n >= 1);
    return solve_a2<false>(n - 1);
}

u64 D_mod(const int n, const u64 mod) {
    assert(n >= 1);
    return solve_a2<true>(n - 1, mod);
}

}  // namespace

int main() {
    constexpr u64 MOD = 1'000'000'000ULL;

    assert(D_exact(10) == 44'499ULL);
    assert(D_exact(20) == 9'204'559'704ULL);
    assert(D_exact(32) == 22'037'102'049'132'222ULL);
    assert(D_mod(100, MOD) == 780'166'455ULL);

    const u64 ans = D_mod(10'000, MOD);
    std::cout << std::setw(9) << std::setfill('0') << ans << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler763 {
    static int threshold(int n) {
        return (n + 1) * (n + 2) / 2;
    }

    static int maxNForM(int m) {
        int n = 0;
        while (threshold(n + 1) <= m) {
            ++n;
        }
        return n;
    }

    static long solveA2(int targetM, long modValue) {
        int nmax = maxNForM(targetM);
        int ring = nmax + 8;
        int nalloc = nmax + 2;

        long[][] u = new long[nalloc + 1][(nalloc + 2) * ring];
        long[][] v = new long[nalloc + 1][(nalloc + 2) * ring];
        long[] a = new long[ring];
        long[] f0 = new long[ring];
        a[0] = 1L;

        for (int m = 0; m <= targetM; ++m) {
            int cur = m % ring;

            for (int n = 1; n <= nmax; ++n) {
                for (int k = 1; k <= n; ++k) {
                    int idx = k * ring + cur;
                    u[n][idx] = 0L;
                    v[n][idx] = 0L;
                }
            }

            for (int n = 1; n <= nmax; ++n) {
                if (m < threshold(n)) {
                    continue;
                }

                for (int k = 1; k <= n; ++k) {
                    int h = (k < n) ? k : (k - 1);

                    long uSum = getU(u, f0, nalloc, ring, n, k, m - n - 2)
                            + getV(v, f0, nalloc, ring, n + 1, 1, m - n - 3)
                            + getU(u, f0, nalloc, ring, n + 1, k + 1, m - n - 3)
                            + getU(u, f0, nalloc, ring, n - 1, h, m - n - 1)
                            + getV(v, f0, nalloc, ring, n, 1, m - n - 2)
                            + getU(u, f0, nalloc, ring, n, h + 1, m - n - 2);

                    long vSum = getV(v, f0, nalloc, ring, n, k, m - n - 2)
                            + getV(v, f0, nalloc, ring, n + 1, k + 1, m - n - 3)
                            + ((k == n) ? 2L : 1L) * getU(u, f0, nalloc, ring, n + 1, 1, m - n - 3)
                            + getV(v, f0, nalloc, ring, n - 1, h, m - n - 1)
                            + getV(v, f0, nalloc, ring, n, h + 1, m - n - 2)
                            + ((k == n) ? 2L : 1L) * getU(u, f0, nalloc, ring, n, 1, m - n - 2);

                    int idx = k * ring + cur;
                    u[n][idx] = uSum % modValue;
                    v[n][idx] = vSum % modValue;
                }
            }

            long f0Sum = getA(a, ring, m - 1)
                    + 4L * getF0(f0, ring, m - 2)
                    + 2L * getU(u, f0, nalloc, ring, 1, 1, m - 3)
                    + getV(v, f0, nalloc, ring, 1, 1, m - 3);
            f0[cur] = f0Sum % modValue;

            if (m >= 1) {
                long aSum = 3L * getA(a, ring, m - 1) + 3L * getF0(f0, ring, m - 2);
                a[cur] = aSum % modValue;
            }
        }

        return a[targetM % ring];
    }

    static long getA(long[] a, int ring, int m) {
        if (m < 0)
            return 0L;
        return a[m % ring];
    }

    static long getF0(long[] f0, int ring, int m) {
        if (m < 0)
            return 0L;
        return f0[m % ring];
    }

    static long getU(long[][] u, long[] f0, int nalloc, int ring, int n, int k, int m) {
        if (m < 0)
            return 0L;
        if (n == 0)
            return (k == 0) ? getF0(f0, ring, m) : 0L;
        if (n < 0 || n > nalloc || k < 1 || k > n || m < threshold(n))
            return 0L;
        int idx = k * ring + (m % ring);
        return u[n][idx];
    }

    static long getV(long[][] v, long[] f0, int nalloc, int ring, int n, int k, int m) {
        if (m < 0)
            return 0L;
        if (n == 0)
            return (k == 0) ? getF0(f0, ring, m) : 0L;
        if (n < 0 || n > nalloc || k < 1 || k > n || m < threshold(n))
            return 0L;
        int idx = k * ring + (m % ring);
        return v[n][idx];
    }

    static long DMod(int n, long mod) {
        return solveA2(n - 1, mod);
    }

    public static String solve() {
        long ans = DMod(10000, 1000000000L);
        return String.format("%09d", ans);
    }

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