Problem 559: Permuted Matrices

View on Project Euler

Project Euler Problem 559 Solution

EulerSolve provides an optimized solution for Project Euler Problem 559, Permuted Matrices, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 559 asks for the value of \(Q(50000)\) modulo the prime $$M=1000000123.$$ The matrix-counting statement is encoded through auxiliary quantities \(P(k,r,n)\). The C++, Python, and Java implementations never try to enumerate the matrices directly. Instead, they transform the count into a coefficient-extraction problem built from inverse factorial powers, alternating signs, and one outer sum over \(k\). Mathematical Approach Fix \(n\), \(r\), and \(k\). The implementations first compute a normalized core value and multiply by \((n!)^r\) only at the end. That normalization turns the combinatorial problem into a clean reciprocal-series recurrence. Step 1: Normalize by factorial powers Because \(M\) is prime and the computation only needs factorials up to \(n=50000\lt M\), every \(m!\) is invertible modulo \(M\). Define $$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$ These numbers are the basic weights of the method. All divisions by factorials are replaced by modular inverses, and the missing factor \((n!)^r\) is restored only after the normalized coefficient has been computed. Step 2: Split \(n\) into full \(k\)-blocks and a remainder Write $$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$ The recurrence advances in jumps of size \(k\)....

Detailed mathematical approach

Problem Summary

Problem 559 asks for the value of \(Q(50000)\) modulo the prime

$$M=1000000123.$$

The matrix-counting statement is encoded through auxiliary quantities \(P(k,r,n)\). The C++, Python, and Java implementations never try to enumerate the matrices directly. Instead, they transform the count into a coefficient-extraction problem built from inverse factorial powers, alternating signs, and one outer sum over \(k\).

Mathematical Approach

Fix \(n\), \(r\), and \(k\). The implementations first compute a normalized core value and multiply by \((n!)^r\) only at the end. That normalization turns the combinatorial problem into a clean reciprocal-series recurrence.

Step 1: Normalize by factorial powers

Because \(M\) is prime and the computation only needs factorials up to \(n=50000\lt M\), every \(m!\) is invertible modulo \(M\). Define

$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$

These numbers are the basic weights of the method. All divisions by factorials are replaced by modular inverses, and the missing factor \((n!)^r\) is restored only after the normalized coefficient has been computed.

Step 2: Split \(n\) into full \(k\)-blocks and a remainder

Write

$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$

The recurrence advances in jumps of size \(k\). The quotient \(b\) counts how many full \(k\)-steps fit into \(n\), while the remainder \(s\) records the final incomplete part. This is why the formulas separate naturally into a full-block part and a remainder correction.

Step 3: Build the reciprocal power series

Introduce the formal power series

$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$

Now define its reciprocal

$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$

Comparing coefficients in \(A_k(z)F_k(z)=1\) gives

$$a_0=1,$$

$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$

This alternating convolution is the heart of the algorithm. Once the sequence \(a_0,a_1,\dots,a_b\) is known, the rest is just a short finishing sum.

Step 4: Correct the final partial block

If \(s=0\), the normalized core contribution is simply

$$\gamma(n,k,r)=a_b.$$

If \(s\gt 0\), one shifted series is needed:

$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$

The desired core is then the coefficient of \(z^b\) in the product \(A_k(z)G_{k,s}(z)\):

$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$

So the same recurrence handles all complete \(k\)-blocks, and the shifted tail series accounts for the leftover \(s\) entries.

Step 5: Recover \(P(k,r,n)\) and assemble \(Q(n)\)

After the normalized core has been computed, the required count is

$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$

For the Project Euler target, the implementations set \(r=n\) and sum over every \(k\):

$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$

The full problem is therefore reduced to evaluating one reciprocal-series coefficient for each \(k\), summing those normalized cores, and multiplying once at the end.

Worked Example: \(P(1,2,3)=19\)

Take \(k=1\), \(r=2\), and \(n=3\). Then \(b=3\) and \(s=0\). The relevant weights are

$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$

Now apply the recurrence:

$$a_0=1,$$

$$a_1=a_0u_1=1,$$

$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$

$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$

Because \(s=0\), we have \(\gamma(3,1,2)=a_3\). Therefore

$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$

which matches the built-in checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations follow the same mathematics. They first precompute factorials and inverse factorials, then raise each inverse factorial to the exponent \(r\) to obtain the table \(u_m\). For a fixed \(k\), they fill the coefficient sequence \(a_0,a_1,\dots,a_b\) with the alternating convolution above and evaluate the shifted finishing sum when \(s>0\).

To compute \(P(k,r,n)\), the implementation multiplies the normalized core by \((n!)^r\). To compute \(Q(n)\), it uses \(r=n\), repeats the same per-\(k\) core computation for every \(1\le k\le n\), adds the cores, and only then multiplies by \((n!)^n\). The C++ version parallelizes the outer loop over \(k\), the Java version runs the same logic serially, and the Python version acts as a thin bridge that invokes the compiled C++ solver and returns its numeric result.

Complexity Analysis

Precomputing factorials and inverse factorials costs \(O(n)\). Building the table \(u_m=(m!)^{-r}\) with repeated fast exponentiation costs \(O(n\log r)\). For one fixed \(k\), the recurrence length is \(b=\lfloor n/k\rfloor\), and the nested alternating sum costs \(O(b^2)\) time. Therefore

$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$

so \(Q(n)\) is evaluated in overall \(O(n^2)\) time. The main tables use \(O(n)\) memory; the multithreaded C++ version adds one reusable work array per worker thread.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=559
  2. Permutation: Wikipedia — Permutation
  3. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
  4. Formal power series: Wikipedia — Formal power series
  5. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse

Problem 559 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <vector>

using int64 = long long;

constexpr int64 MOD = 1000000123LL;

static int64 mod_pow(int64 base, int64 exp) {
    int64 result = 1 % MOD;
    base %= MOD;
    while (exp > 0) {
        if (exp & 1LL) result = (result * base) % MOD;
        base = (base * base) % MOD;
        exp >>= 1LL;
    }
    return result;
}

static void build_factorials(int n, std::vector<int64>& fact, std::vector<int64>& inv_fact) {
    fact.assign(n + 1, 0);
    inv_fact.assign(n + 1, 0);
    fact[0] = 1;
    for (int i = 1; i <= n; ++i) {
        fact[i] = (fact[i - 1] * i) % MOD;
    }
    inv_fact[n] = mod_pow(fact[n], MOD - 2);
    for (int i = n; i >= 1; --i) {
        inv_fact[i - 1] = (inv_fact[i] * i) % MOD;
    }
}

static std::vector<int64> build_pow_invfact(const std::vector<int64>& inv_fact, int n, int r) {
    std::vector<int64> pow_invfact(n + 1, 0);
    for (int i = 0; i <= n; ++i) {
        pow_invfact[i] = mod_pow(inv_fact[i], r);
    }
    return pow_invfact;
}

static int64 compute_core_with_workspace(int n, int k, const int64* pow_invfact, std::vector<int64>& dp) {
    const int blocks = n / k;
    const int rem = n - blocks * k;

    dp[0] = 1;
    for (int i = 1; i <= blocks; ++i) {
        int64 sum = 0;
        int sign = 1;
        int idx = k;
        for (int j = 1; j <= i; ++j) {
            int64 term = (dp[i - j] * pow_invfact[idx]) % MOD;
            sum += sign * term;
            sign = -sign;
            idx += k;
        }
        sum %= MOD;
        if (sum < 0) sum += MOD;
        dp[i] = sum;
    }

    if (rem == 0) {
        return dp[blocks];
    }

    int64 total = 0;
    int sign = 1;
    int idx = rem;
    for (int j = 0; j <= blocks; ++j) {
        int64 term = (dp[blocks - j] * pow_invfact[idx]) % MOD;
        total += sign * term;
        sign = -sign;
        idx += k;
    }
    total %= MOD;
    if (total < 0) total += MOD;
    return total;
}

static int64 compute_P(int n, int r, int k, const std::vector<int64>& fact,
                       const std::vector<int64>& inv_fact) {
    std::vector<int64> pow_invfact = build_pow_invfact(inv_fact, n, r);
    std::vector<int64> dp(n + 1, 0);
    int64 core = compute_core_with_workspace(n, k, pow_invfact.data(), dp);
    int64 factor = mod_pow(fact[n], r);
    return (core * factor) % MOD;
}

static int64 compute_Q(int n, const std::vector<int64>& fact, const std::vector<int64>& inv_fact,
                       unsigned threads) {
    std::vector<int64> pow_invfact = build_pow_invfact(inv_fact, n, n);
    const int64* w = pow_invfact.data();

    if (threads == 0) threads = 1;
    threads = std::min<unsigned>(threads, static_cast<unsigned>(n));

    std::atomic<int> next_k(1);
    std::vector<int64> partial(threads, 0);

    auto worker = [&](unsigned tid) {
        std::vector<int64> dp(n + 1, 0);
        int64 local = 0;
        while (true) {
            int k = next_k.fetch_add(1, std::memory_order_relaxed);
            if (k > n) break;

            int64 core = compute_core_with_workspace(n, k, w, dp);
            local += core;
        }
        partial[tid] = local % MOD;
    };

    std::vector<std::thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    for (auto& th : pool) th.join();

    int64 total = 0;
    for (int64 v : partial) {
        total += v;
    }
    total %= MOD;
    int64 factor = mod_pow(fact[n], n);
    return (total * factor) % MOD;
}

static bool run_validations(const std::vector<int64>& fact, const std::vector<int64>& inv_fact) {
    struct CheckP {
        int k;
        int r;
        int n;
        int64 expected;
    };
    const std::vector<CheckP> checks = {
        {1, 2, 3, 19},
        {2, 4, 6, 65508751},
        {7, 5, 30, 161858102},
    };

    for (const auto& check : checks) {
        int64 got = compute_P(check.n, check.r, check.k, fact, inv_fact);
        if (got != check.expected) {
            std::cerr << "Validation failed for P(" << check.k << ", " << check.r << ", "
                      << check.n << "): got=" << got << " expected=" << check.expected << "\n";
            return false;
        }
    }

    const int64 expected_q5 = 879391168;
    int64 q5 = compute_Q(5, fact, inv_fact, 1);
    if (q5 != expected_q5) {
        std::cerr << "Validation failed for Q(5): got=" << q5 << " expected=" << expected_q5
                  << "\n";
        return false;
    }

    const int64 expected_q50 = 819573537;
    int64 q50 = compute_Q(50, fact, inv_fact, 1);
    if (q50 != expected_q50) {
        std::cerr << "Validation failed for Q(50): got=" << q50 << " expected=" << expected_q50
                  << "\n";
        return false;
    }

    return true;
}

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n = 50000;
    if (argc > 1) n = std::atoi(argv[1]);
    if (n <= 0) {
        std::cerr << "n must be positive.\n";
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 4;
    if (argc > 2) threads = static_cast<unsigned>(std::max(1, std::atoi(argv[2])));

    int max_n = std::max(n, 50);
    std::vector<int64> fact;
    std::vector<int64> inv_fact;
    build_factorials(max_n, fact, inv_fact);

    if (!run_validations(fact, inv_fact)) return 1;

    int64 answer = compute_Q(n, fact, inv_fact, threads);
    std::cout << answer << "\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 Euler559 {
    static final long MOD = 1000000123L;

    static long modPow(long base, long exp) {
        long result = 1 % MOD;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) != 0)
                result = (result * base) % MOD;
            base = (base * base) % MOD;
            exp >>= 1;
        }
        return result;
    }

    static class Factorials {
        long[] fact;
        long[] invFact;

        Factorials(int n) {
            fact = new long[n + 1];
            invFact = new long[n + 1];
            fact[0] = 1;
            for (int i = 1; i <= n; i++) {
                fact[i] = (fact[i - 1] * i) % MOD;
            }
            invFact[n] = modPow(fact[n], MOD - 2);
            for (int i = n; i >= 1; i--) {
                invFact[i - 1] = (invFact[i] * i) % MOD;
            }
        }
    }

    static long[] buildPowInvFact(long[] invFact, int n, int r) {
        long[] powInvFact = new long[n + 1];
        for (int i = 0; i <= n; i++) {
            powInvFact[i] = modPow(invFact[i], r);
        }
        return powInvFact;
    }

    static long computeCore(int n, int k, long[] powInvFact) {
        int blocks = n / k;
        int rem = n % k;

        long[] dp = new long[blocks + 1];
        dp[0] = 1;

        for (int i = 1; i <= blocks; i++) {
            long sum = 0;
            int idx = k;
            for (int j = 1; j <= i;) {
                sum += dp[i - j] * powInvFact[idx];
                j++;
                idx += k;
                if (j > i)
                    break;

                sum -= dp[i - j] * powInvFact[idx];
                j++;
                idx += k;

                if ((j & 7) == 1) {
                    sum %= MOD;
                }
            }
            sum %= MOD;
            if (sum < 0)
                sum += MOD;
            dp[i] = sum;
        }

        if (rem == 0)
            return dp[blocks];

        long total = 0;
        int idx = rem;
        for (int j = 0; j <= blocks;) {
            total += dp[blocks - j] * powInvFact[idx];
            j++;
            idx += k;
            if (j > blocks)
                break;

            total -= dp[blocks - j] * powInvFact[idx];
            j++;
            idx += k;

            if ((j & 7) == 0) {
                total %= MOD;
            }
        }
        total %= MOD;
        if (total < 0)
            total += MOD;
        return total;
    }

    static long computeQ(int n) {
        Factorials f = new Factorials(n);
        long[] powInvFact = buildPowInvFact(f.invFact, n, n);

        long total = 0;
        for (int k = 1; k <= n; k++) {
            total = (total + computeCore(n, k, powInvFact)) % MOD;
        }

        long factor = modPow(f.fact[n], n);
        return (total * factor) % MOD;
    }

    public static String solve() {
        return Long.toString(computeQ(50000));
    }

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