Problem 380: Amazing Mazes!

View on Project Euler

Project Euler Problem 380 Solution

EulerSolve provides an optimized solution for Project Euler Problem 380, Amazing Mazes!, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(C(m,n)\) denote the number of spanning trees of the \(m \times n\) rectangular grid graph. The Project Euler instance asks for the value at \((m,n)=(100,500)\), but the exact integer has more than twenty-five thousand decimal digits, so the implementations return the answer in scientific notation instead of constructing the full integer explicitly. The local solution files all exploit the same fact: for this particular graph family, the Matrix-Tree theorem converts the counting problem into a product over explicit Laplacian eigenvalues. Mathematical Approach Step 1: Use the Matrix-Tree Theorem For any connected graph \(G\) with Laplacian eigenvalues $$0=\lambda_0\lt\lambda_1\le \lambda_2 \le \dots \le \lambda_{|V|-1},$$ Kirchhoff's Matrix-Tree theorem gives the spanning-tree count as $$\tau(G)=\frac{1}{|V|}\prod_{k=1}^{|V|-1}\lambda_k.$$ In our case, the grid graph has \(mn\) vertices, so once the nonzero Laplacian eigenvalues are known, we immediately obtain \(C(m,n)\). Step 2: Laplacian Eigenvalues of the Rectangular Grid The \(m \times n\) rectangular grid is the Cartesian product of two path graphs, one of length \(m\) and one of length \(n\)....

Detailed mathematical approach

Problem Summary

Let \(C(m,n)\) denote the number of spanning trees of the \(m \times n\) rectangular grid graph. The Project Euler instance asks for the value at \((m,n)=(100,500)\), but the exact integer has more than twenty-five thousand decimal digits, so the implementations return the answer in scientific notation instead of constructing the full integer explicitly.

The local solution files all exploit the same fact: for this particular graph family, the Matrix-Tree theorem converts the counting problem into a product over explicit Laplacian eigenvalues.

Mathematical Approach

Step 1: Use the Matrix-Tree Theorem

For any connected graph \(G\) with Laplacian eigenvalues

$$0=\lambda_0\lt\lambda_1\le \lambda_2 \le \dots \le \lambda_{|V|-1},$$

Kirchhoff's Matrix-Tree theorem gives the spanning-tree count as

$$\tau(G)=\frac{1}{|V|}\prod_{k=1}^{|V|-1}\lambda_k.$$

In our case, the grid graph has \(mn\) vertices, so once the nonzero Laplacian eigenvalues are known, we immediately obtain \(C(m,n)\).

Step 2: Laplacian Eigenvalues of the Rectangular Grid

The \(m \times n\) rectangular grid is the Cartesian product of two path graphs, one of length \(m\) and one of length \(n\). For the path graph with \(m\) vertices, the Laplacian eigenvalues are

$$\alpha_i=2-2\cos\frac{\pi i}{m},\qquad i=0,1,\dots,m-1,$$

and similarly for the path graph with \(n\) vertices,

$$\beta_j=2-2\cos\frac{\pi j}{n},\qquad j=0,1,\dots,n-1.$$

For a Cartesian product, Laplacian eigenvalues add, so the grid eigenvalues are

$$\lambda_{i,j}=\alpha_i+\beta_j=4-2\cos\frac{\pi i}{m}-2\cos\frac{\pi j}{n}.$$

The pair \((i,j)=(0,0)\) gives the single zero eigenvalue corresponding to the constant vector. Every other pair gives a positive eigenvalue because the grid graph is connected.

Step 3: Closed Product Formula for \(C(m,n)\)

Substituting the grid eigenvalues into Kirchhoff's formula yields

$$\boxed{C(m,n)=\frac{1}{mn}\prod_{\substack{0 \le i \lt m\\0 \le j \lt n\\ (i,j)\ne(0,0)}}\left(4-2\cos\frac{\pi i}{m}-2\cos\frac{\pi j}{n}\right).}$$

This is the exact formula implemented by all three solution files. The main difficulty is no longer graph theory but numerical scale: the product contains \(mn-1\) positive factors and becomes astronomically large even for moderate \(m\) and \(n\).

Step 4: Why the Code Works in Base-10 Logarithms

Instead of multiplying the eigenvalues directly, the programs sum their base-10 logarithms:

$$L=\log_{10} C(m,n)=\sum_{\substack{0 \le i \lt m\\0 \le j \lt n\\ (i,j)\ne(0,0)}}\log_{10}\left(4-2\cos\frac{\pi i}{m}-2\cos\frac{\pi j}{n}\right)-\log_{10}(mn).$$

Once \(L\) is known, write

$$e=\lfloor L \rfloor,\qquad \mu=10^{L-e},\qquad C(m,n)=\mu \times 10^e,$$

where \(1 \le \mu \lt 10\). The code rounds \(\mu\) to five significant digits. If rounding pushes the mantissa to \(10.0000\), it divides by \(10\) and increments the exponent, exactly as the local implementations do.

The C++ version uses Kahan summation for the logarithm sum, which is a small but worthwhile refinement when tens of thousands of floating-point terms are added. The Python and Java versions use the same formula without that extra compensation step.

Worked Checkpoints

The C++ file contains several checkpoints that are useful both mathematically and as implementation tests.

For \(m=n=1\), there is only one vertex and therefore exactly one spanning tree, so \(C(1,1)=1\).

For \(m=n=2\), the nonzero eigenvalues are \(2\), \(2\), and \(4\), hence

$$C(2,2)=\frac{2 \cdot 2 \cdot 4}{4}=4.$$

The next exact checkpoint used in code is

$$C(3,4)=2415,$$

and symmetry of the rectangular grid requires \(C(3,4)=C(4,3)\), which the C++ program checks numerically.

For a larger validation point, the code formats

$$C(9,12)\approx 2.5720 \times 10^{46},$$

and for the actual Project Euler target it prints

$$C(100,500)\approx 6.3202 \times 10^{25093}.$$

How the Code Works

The C++ solution exposes optional arguments --m=, --n=, and --skip-checkpoints. It precomputes two cosine tables, loops over all \((i,j)\ne(0,0)\), forms the eigenvalue \(4-2c_i-2c_j\), accumulates \(\log_{10}\) with Kahan compensation, subtracts \(\log_{10}(mn)\), and finally reconstructs the scientific notation string. The Python solution follows the same mathematical pipeline and also caches cosine tables, but keeps the code minimal and does not include checkpoints. The Java solution uses the same formula and the same scientific-format rounding rule, but computes the cosine values directly inside the loops instead of storing arrays.

Complexity Analysis

The dominant work is the double loop over the \(mn-1\) nonzero eigenvalue positions, so the running time is \(O(mn)\). The C++ and Python versions store cosine tables of lengths \(m\) and \(n\), so they use \(O(m+n)\) extra memory. The Java translation recomputes cosines and therefore only needs \(O(1)\) extra memory beyond a few scalars. In every version, this is vastly cheaper than forming a giant determinant or performing exact big-integer linear algebra.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=380
  2. Matrix-Tree theorem: Wikipedia — Kirchhoff's theorem
  3. Laplacian matrix of a graph: Wikipedia — Laplacian matrix
  4. Cartesian product of graphs: Wikipedia — Cartesian product of graphs
  5. Path graph: Wikipedia — Path graph

Problem 380 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

namespace {

struct Options {
    int m = 100;
    int n = 500;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--m=", options.m) ||
            parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.m >= 1 && options.n >= 1;
}

long double compute_log10_spanning_trees(const int m, const int n) {
    const long double pi = acosl(-1.0L);

    std::vector<long double> cos_m(static_cast<std::size_t>(m));
    std::vector<long double> cos_n(static_cast<std::size_t>(n));
    for (int i = 0; i < m; ++i) {
        cos_m[static_cast<std::size_t>(i)] = cosl(pi * static_cast<long double>(i) / m);
    }
    for (int j = 0; j < n; ++j) {
        cos_n[static_cast<std::size_t>(j)] = cosl(pi * static_cast<long double>(j) / n);
    }

    // Kahan summation keeps rounding error tiny across many logarithms.
    long double sum_log10 = 0.0L;
    long double compensation = 0.0L;
    for (int i = 0; i < m; ++i) {
        const long double ci = cos_m[static_cast<std::size_t>(i)];
        for (int j = 0; j < n; ++j) {
            if (i == 0 && j == 0) {
                continue;
            }
            const long double cj = cos_n[static_cast<std::size_t>(j)];
            const long double eigenvalue = 4.0L - 2.0L * ci - 2.0L * cj;
            const long double y = log10l(eigenvalue) - compensation;
            const long double t = sum_log10 + y;
            compensation = (t - sum_log10) - y;
            sum_log10 = t;
        }
    }

    sum_log10 -= log10l(static_cast<long double>(m) * static_cast<long double>(n));
    return sum_log10;
}

std::string format_scientific_from_log10(const long double log10_value, const int significant_digits) {
    const long double exponent_ld = floorl(log10_value);
    long long exponent = static_cast<long long>(exponent_ld);
    long double mantissa = powl(10.0L, log10_value - exponent_ld);

    long double scale = 1.0L;
    for (int i = 1; i < significant_digits; ++i) {
        scale *= 10.0L;
    }

    long double rounded = roundl(mantissa * scale);
    if (rounded >= 10.0L * scale) {
        rounded /= 10.0L;
        ++exponent;
    }

    std::ostringstream oss;
    oss << std::fixed << std::setprecision(significant_digits - 1)
        << static_cast<double>(rounded / scale) << 'e' << exponent;
    return oss.str();
}

bool is_nearly_equal(const long double a, const long double b, const long double rel_tol) {
    const long double diff = fabsl(a - b);
    const long double scale = std::max(fabsl(a), fabsl(b));
    if (scale == 0.0L) {
        return diff <= rel_tol;
    }
    return diff <= rel_tol * scale;
}

bool run_checkpoints() {
    if (!is_nearly_equal(compute_log10_spanning_trees(1, 1), 0.0L, 1e-18L)) {
        std::cerr << "Checkpoint failed: C(1,1) != 1\n";
        return false;
    }

    const long double c22 = powl(10.0L, compute_log10_spanning_trees(2, 2));
    if (llround(c22) != 4LL) {
        std::cerr << "Checkpoint failed: C(2,2)\n";
        return false;
    }

    const long double c34 = powl(10.0L, compute_log10_spanning_trees(3, 4));
    if (llround(c34) != 2415LL) {
        std::cerr << "Checkpoint failed: C(3,4)\n";
        return false;
    }

    const std::string c912 = format_scientific_from_log10(compute_log10_spanning_trees(9, 12), 5);
    if (c912 != "2.5720e46") {
        std::cerr << "Checkpoint failed: C(9,12), got " << c912 << '\n';
        return false;
    }

    const long double symmetry_lhs = compute_log10_spanning_trees(3, 4);
    const long double symmetry_rhs = compute_log10_spanning_trees(4, 3);
    if (!is_nearly_equal(symmetry_lhs, symmetry_rhs, 1e-16L)) {
        std::cerr << "Checkpoint failed: symmetry C(m,n)=C(n,m)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const long double answer_log10 = compute_log10_spanning_trees(options.m, options.n);
    std::cout << format_scientific_from_log10(answer_log10, 5) << '\n';
    return 0;
}

Python

import math

def solve():
    m, n = 100, 500
    pi = math.pi

    cos_m = [math.cos(pi * i / m) for i in range(m)]
    cos_n = [math.cos(pi * j / n) for j in range(n)]

    sum_log10 = 0.0
    for i in range(m):
        ci = cos_m[i]
        for j in range(n):
            if i == 0 and j == 0:
                continue
            cj = cos_n[j]
            eigenvalue = 4.0 - 2.0 * ci - 2.0 * cj
            sum_log10 += math.log10(eigenvalue)

    sum_log10 -= math.log10(m * n)

    # Format as scientific with 5 significant digits
    exponent = int(math.floor(sum_log10))
    mantissa = 10.0 ** (sum_log10 - exponent)
    scale = 10000.0
    rounded = round(mantissa * scale)
    if rounded >= 10.0 * scale:
        rounded /= 10.0
        exponent += 1
    result = rounded / scale
    return f"{result:.4f}e{exponent}"

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

Java

public class Euler380 {
    public static String solve() {
        int m = 100, n = 500;
        double pi = Math.PI;
        double sumLog10 = 0.0;
        for (int i = 0; i < m; i++) {
            double ci = Math.cos(pi * i / m);
            for (int j = 0; j < n; j++) {
                if (i == 0 && j == 0)
                    continue;
                double cj = Math.cos(pi * j / n);
                double eigenvalue = 4.0 - 2.0 * ci - 2.0 * cj;
                sumLog10 += Math.log10(eigenvalue);
            }
        }
        sumLog10 -= Math.log10((double) m * n);
        int exponent = (int) Math.floor(sumLog10);
        double mantissa = Math.pow(10.0, sumLog10 - exponent);
        double scale = 10000.0;
        long rounded = Math.round(mantissa * scale);
        if (rounded >= (long) (10.0 * scale)) {
            rounded /= 10;
            exponent++;
        }
        double result = rounded / scale;
        return String.format("%.4fe%d", result, exponent);
    }

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