Problem 519: Tricoloured Coin Fountains

Project Euler Problem 519 Solution

EulerSolve provides an optimized solution for Project Euler Problem 519, Tricoloured Coin Fountains, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for the tricoloured coin-fountain count $T(20000)$ modulo $10^9$. Exhaustive generation is completely impractical at this size, so the solution works with truncated generating functions and formal power-series algebra instead of exploring fountain states directly. Mathematical Approach The computation is organized around one auxiliary series $F_2(x)$. After that series is known to the required degree, a second rational transformation produces the generating function whose coefficients are the desired values $T(n)$. Step 1: Start from the continued fraction The auxiliary series is $$F_2(x)=\frac{1}{1-\frac{x^2}{1-\frac{x^3}{1-\frac{x^4}{\ddots}}}}.$$ This continued fraction is the compact object behind the counting problem. Expanding it directly to degree $20000$ would be awkward, so the implementations switch to an equivalent quotient of two explicit $q$-series. Step 2: Replace it by a quotient of two $q$-series Define the finite $q$-Pochhammer product $$(x;x)_m=\prod_{k=1}^{m}(1-x^k).$$ The key identity used by the implementations is $$F_2(x)=\frac{\sum_{m\ge 0}(-1)^m\dfrac{x^{m(m+2)}}{(x;x)_m}}{\sum_{m\ge 0}(-1)^m\dfrac{x^{m(m+1)}}{(x;x)_m}}.$$ So the continued fraction is converted into a numerator series and a denominator series....

Detailed mathematical approach

Problem Summary

The problem asks for the tricoloured coin-fountain count $T(20000)$ modulo $10^9$. Exhaustive generation is completely impractical at this size, so the solution works with truncated generating functions and formal power-series algebra instead of exploring fountain states directly.

Mathematical Approach

The computation is organized around one auxiliary series $F_2(x)$. After that series is known to the required degree, a second rational transformation produces the generating function whose coefficients are the desired values $T(n)$.

Step 1: Start from the continued fraction

The auxiliary series is

$$F_2(x)=\frac{1}{1-\frac{x^2}{1-\frac{x^3}{1-\frac{x^4}{\ddots}}}}.$$

This continued fraction is the compact object behind the counting problem. Expanding it directly to degree $20000$ would be awkward, so the implementations switch to an equivalent quotient of two explicit $q$-series.

Step 2: Replace it by a quotient of two $q$-series

Define the finite $q$-Pochhammer product

$$(x;x)_m=\prod_{k=1}^{m}(1-x^k).$$

The key identity used by the implementations is

$$F_2(x)=\frac{\sum_{m\ge 0}(-1)^m\dfrac{x^{m(m+2)}}{(x;x)_m}}{\sum_{m\ge 0}(-1)^m\dfrac{x^{m(m+1)}}{(x;x)_m}}.$$

So the continued fraction is converted into a numerator series and a denominator series. For a truncation to degree $N$, only finitely many indices $m$ contribute: once both $m(m+1)$ and $m(m+2)$ exceed $N$, that index can no longer affect any coefficient up to $x^N$.

Step 3: Expand $1/(x;x)_m$ incrementally

The reciprocal products are updated one step at a time through

$$\frac{1}{(x;x)_{m+1}}=\frac{1}{(x;x)_m}\cdot\frac{1}{1-x^{m+1}}=\frac{1}{(x;x)_m}\left(1+x^{m+1}+x^{2(m+1)}+\cdots\right).$$

This means that once the truncated coefficients of $1/(x;x)_m$ are known, the next one is obtained by adding shifted copies spaced by $m+1$. The code therefore carries a single evolving truncated series instead of rebuilding each reciprocal product from scratch.

Step 4: Recover coefficients by formal series division

Suppose

$$A(x)=\frac{N(x)}{D(x)},\qquad D(x)=1+\sum_{i\ge 1} d_i x^i,\qquad N(x)=\sum_{n\ge 0} n_n x^n.$$

$$A(x)=\sum_{n\ge 0} a_n x^n,$$

then comparing coefficients in $A(x)D(x)=N(x)$ gives

$$a_0=n_0,\qquad a_n=n_n-\sum_{i=1}^{n} d_i a_{n-i}\pmod{10^9}\quad(n\ge 1).$$

Because the constant term of the denominator is $1$, every new coefficient depends only on earlier ones. The implementations use this recurrence first to recover $F_2(x)$, and later to recover the final target series.

Step 5: Transform $F_2(x)$ into the target sequence

After $F_2(x)$ has been computed to the required degree, the target series is

$$t(x)=\frac{x(2F_2(x)-1)}{1-2x(2F_2(x)-1)}.$$

The desired counting sequence is then extracted through

$$T(n)=3\,[x^n]\,t(x),$$

where $[x^n]$ denotes the coefficient of $x^n$. So the whole computation is: build the two $q$-series for $F_2$, divide them, form the numerator and denominator of $t(x)$, divide again, read the coefficient of $x^n$, and multiply by $3$.

Worked Example: Recover $T(4)$

To degree $4$, the continued fraction contributes only enough to give

$$F_2(x)=1+x^2+x^4+O(x^5).$$

Hence

$$A(x)=2F_2(x)-1=1+2x^2+2x^4+O(x^5),$$

$$N(x)=xA(x)=x+2x^3+O(x^5),$$

$$D(x)=1-2xA(x)=1-2x-4x^3+O(x^5).$$

Writing $t(x)=\sum_{n\ge 0} t_n x^n$, the division recurrence yields

$$t_0=0,\qquad t_1=1,\qquad t_2=2,\qquad t_3=6,\qquad t_4=16.$$

Therefore

$$T(4)=3\,t_4=48,$$

which matches the low-degree checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. The truncated numerator and denominator of $F_2(x)$ are accumulated term by term, with alternating signs, while a rolling truncated series stores the current reciprocal product $1/(x;x)_m$. After the first formal division, the code builds the numerator and denominator of $t(x)$ from the coefficients of $F_2(x)$ and applies the same division recurrence again.

The final answer is the coefficient of $x^{20000}$ in $t(x)$, multiplied by $3$ and reduced modulo $10^9$. The implementations also perform small internal checks: they compare the $q$-series coefficients of $F_2$ with a direct low-degree continued-fraction expansion, and they verify checkpoint values such as $T(4)=48$ and $T(10)=17760$. The Java implementation delegates the heavy computation to the compiled implementation of the same formulas.

Complexity Analysis

Let $N$ be the truncation degree. Building the two $q$-series for $F_2(x)$ takes $O(N\sqrt{N})$ coefficient updates because only $m\le O(\sqrt{N})$ contribute, and each contributing index touches $O(N)$ coefficients. The dominant cost is the two formal series divisions, each of which is $O(N^2)$. Therefore the overall running time is $O(N^2)$, and the memory usage is $O(N)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=519
Formal power series: Wikipedia — Formal power series
$q$-Pochhammer symbol: Wikipedia — q-Pochhammer symbol
Continued fraction: Wikipedia — Continued fraction

Problem 519 source code

C++

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

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;
using u32 = std::uint32_t;

constexpr u64 kMod = 1'000'000'000ULL;  // last 9 digits

u64 mod_norm(i128 v) {
    const i128 m = static_cast<i128>(kMod);
    v %= m;
    if (v < 0) {
        v += m;
    }
    return static_cast<u64>(v);
}

std::vector<u32> series_divide(const std::vector<u32>& num, const std::vector<u32>& den,
                               const int degree) {
    // den[0] must be 1. Returns num/den modulo kMod up to x^degree.
    std::vector<u32> out(static_cast<std::size_t>(degree + 1), 0U);
    out[0] = num[0];
    for (int n = 1; n <= degree; ++n) {
        i128 acc = static_cast<i128>(num[static_cast<std::size_t>(n)]);
        for (int i = 1; i <= n; ++i) {
            acc -= static_cast<i128>(static_cast<u64>(den[static_cast<std::size_t>(i)]) *
                                     out[static_cast<std::size_t>(n - i)]);
        }
        out[static_cast<std::size_t>(n)] = static_cast<u32>(mod_norm(acc));
    }
    return out;
}

std::vector<u32> compute_F2_qseries(const int degree) {
    // F2(x) = 1/(1 - x^2/(1 - x^3/(1 - x^4/(...))))
    // Using the known identity (specializing the A047998 bivariate g.f. at y=x):
    //   F2(x) = (Sum_{n>=0} (-1)^n x^{n(n+2)} / (x;x)_n) / (Sum_{n>=0} (-1)^n x^{n(n+1)} / (x;x)_n),
    // where (x;x)_n = Prod_{k=1..n} (1 - x^k).
    std::vector<u32> den(static_cast<std::size_t>(degree + 1), 0U);
    std::vector<u32> num(static_cast<std::size_t>(degree + 1), 0U);

    std::vector<u32> invprod(static_cast<std::size_t>(degree + 1), 0U);
    invprod[0] = 1U;  // 1/(x;x)_0

    for (int n = 0;; ++n) {
        const int eD = n * (n + 1);
        const int eN = n * (n + 2);
        if (eD > degree && eN > degree) {
            break;
        }

        const bool neg = (n & 1) != 0;
        if (eD <= degree) {
            for (int i = 0; i + eD <= degree; ++i) {
                const u32 add = invprod[static_cast<std::size_t>(i)];
                u32& cell = den[static_cast<std::size_t>(i + eD)];
                if (!neg) {
                    const u64 s = static_cast<u64>(cell) + add;
                    cell = static_cast<u32>(s >= kMod ? s - kMod : s);
                } else {
                    cell = (cell >= add) ? static_cast<u32>(cell - add)
                                         : static_cast<u32>(cell + kMod - add);
                }
            }
        }
        if (eN <= degree) {
            for (int i = 0; i + eN <= degree; ++i) {
                const u32 add = invprod[static_cast<std::size_t>(i)];
                u32& cell = num[static_cast<std::size_t>(i + eN)];
                if (!neg) {
                    const u64 s = static_cast<u64>(cell) + add;
                    cell = static_cast<u32>(s >= kMod ? s - kMod : s);
                } else {
                    cell = (cell >= add) ? static_cast<u32>(cell - add)
                                         : static_cast<u32>(cell + kMod - add);
                }
            }
        }

        const int next = n + 1;
        if (next == 0) {
            continue;
        }
        for (int t = next; t <= degree; ++t) {
            const u64 s = static_cast<u64>(invprod[static_cast<std::size_t>(t)]) +
                          invprod[static_cast<std::size_t>(t - next)];
            invprod[static_cast<std::size_t>(t)] = static_cast<u32>(s % kMod);
        }
    }

    // den[0] == 1 so formal division is well-defined mod kMod.
    return series_divide(num, den, degree);
}

std::vector<u32> compute_F2_contfrac_small(const int degree) {
    // Direct series expansion of the continued fraction, for checkpointing only.
    const int k_max = degree + 5;
    std::vector<std::vector<u32>> F(static_cast<std::size_t>(k_max + 2));
    F[static_cast<std::size_t>(k_max + 1)] =
        std::vector<u32>(static_cast<std::size_t>(degree + 1), 0U);
    F[static_cast<std::size_t>(k_max + 1)][0] = 1U;

    for (int k = k_max; k >= 2; --k) {
        std::vector<u32> den(static_cast<std::size_t>(degree + 1), 0U);
        den[0] = 1U;
        const auto& next = F[static_cast<std::size_t>(k + 1)];
        for (int i = k; i <= degree; ++i) {
            const u32 v = next[static_cast<std::size_t>(i - k)];
            den[static_cast<std::size_t>(i)] = (v == 0U) ? 0U : static_cast<u32>(kMod - v);
        }

        // Invert den: inv[0]=1, inv[n] = -sum_{i=1..n} den[i]*inv[n-i]
        std::vector<u32> inv(static_cast<std::size_t>(degree + 1), 0U);
        inv[0] = 1U;
        for (int n = 1; n <= degree; ++n) {
            i128 acc = 0;
            for (int i = 1; i <= n; ++i) {
                acc -= static_cast<i128>(static_cast<u64>(den[static_cast<std::size_t>(i)]) *
                                         inv[static_cast<std::size_t>(n - i)]);
            }
            inv[static_cast<std::size_t>(n)] = static_cast<u32>(mod_norm(acc));
        }
        F[static_cast<std::size_t>(k)] = std::move(inv);
    }
    return F[2];
}

std::vector<u32> compute_t_from_F2(const std::vector<u32>& F2, const int degree) {
    // t(x) = x*(2*F2(x) - 1) / (1 - 2x(2*F2(x) - 1)) = T(x)/3
    std::vector<u32> A(static_cast<std::size_t>(degree + 1), 0U);     // A = 2F2 - 1
    std::vector<u32> den(static_cast<std::size_t>(degree + 1), 0U);   // 1 - 2xA
    std::vector<u32> num(static_cast<std::size_t>(degree + 1), 0U);   // xA

    for (int i = 0; i <= degree; ++i) {
        const u64 v = 2ULL * F2[static_cast<std::size_t>(i)];
        A[static_cast<std::size_t>(i)] = static_cast<u32>(v % kMod);
    }
    A[0] = (A[0] + static_cast<u32>(kMod - 1ULL)) % static_cast<u32>(kMod);

    den[0] = 1U;
    for (int i = 1; i <= degree; ++i) {
        const u64 v = 2ULL * A[static_cast<std::size_t>(i - 1)];
        den[static_cast<std::size_t>(i)] = static_cast<u32>(v == 0ULL ? 0ULL : (kMod - (v % kMod)) % kMod);
        num[static_cast<std::size_t>(i)] = A[static_cast<std::size_t>(i - 1)];
    }

    return series_divide(num, den, degree);
}

bool run_checkpoints() {
    constexpr int small_deg = 80;
    const auto f2_q = compute_F2_qseries(small_deg);
    const auto f2_cf = compute_F2_contfrac_small(small_deg);
    if (f2_q != f2_cf) {
        std::cerr << "Checkpoint failed: F2 q-series / contfrac mismatch\n";
        return false;
    }

    const auto t = compute_t_from_F2(f2_q, 20);
    const u64 T4 = (3ULL * t[4]) % kMod;
    const u64 T10 = (3ULL * t[10]) % kMod;
    if (T4 != 48ULL) {
        std::cerr << "Checkpoint failed: T(4)\n";
        return false;
    }
    if (T10 != 17'760ULL) {
        std::cerr << "Checkpoint failed: T(10)\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr int n = 20'000;
    const auto f2 = compute_F2_qseries(n);
    const auto t = compute_t_from_F2(f2, n);
    const u64 answer = (3ULL * t[static_cast<std::size_t>(n)]) % kMod;
    std::cout << std::setw(9) << std::setfill('0') << answer << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000000000
    n = 20000

    def mod_norm(v):
        return v % MOD

    def series_divide(num, den, degree):
        out = [0] * (degree + 1)
        out[0] = num[0] % MOD
        for i in range(1, degree + 1):
            acc = num[i] if i < len(num) else 0
            for j in range(1, i + 1):
                if j < len(den): acc -= den[j] * out[i-j]
            out[i] = mod_norm(acc)
        return out

    def compute_F2(degree):
        den = [0] * (degree + 1)
        num = [0] * (degree + 1)
        invprod = [0] * (degree + 1)
        invprod[0] = 1
        for nn in range(10000):
            eD = nn * (nn + 1); eN = nn * (nn + 2)
            if eD > degree and eN > degree: break
            neg = nn & 1
            if eD <= degree:
                for i in range(degree - eD + 1):
                    v = invprod[i]
                    if not neg: den[i+eD] = (den[i+eD] + v) % MOD
                    else: den[i+eD] = (den[i+eD] - v) % MOD
            if eN <= degree:
                for i in range(degree - eN + 1):
                    v = invprod[i]
                    if not neg: num[i+eN] = (num[i+eN] + v) % MOD
                    else: num[i+eN] = (num[i+eN] - v) % MOD
            nxt = nn + 1
            for t in range(nxt, degree + 1):
                invprod[t] = (invprod[t] + invprod[t-nxt]) % MOD
        return series_divide(num, den, degree)

    def compute_t(f2, degree):
        A = [(2 * f2[i]) % MOD for i in range(degree + 1)]
        A[0] = (A[0] - 1) % MOD
        den = [0] * (degree + 1); num = [0] * (degree + 1)
        den[0] = 1
        for i in range(1, degree + 1):
            den[i] = (-2 * A[i-1]) % MOD
            num[i] = A[i-1]
        return series_divide(num, den, degree)

    f2 = compute_F2(n)
    t = compute_t(f2, n)
    return str(3 * t[n] % MOD).zfill(9)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler519 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler519.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler519_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler519 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler519 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler519 C++ bridge produced empty output.");
        }
        return parsed;
    }

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