Problem 412: Gnomon Numbering

Project Euler Problem 412 Solution

EulerSolve provides an optimized solution for Project Euler Problem 412, Gnomon Numbering, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Start with an $m \times m$ square and remove an $n \times n$ square from the lower-right corner, where $0 \le n \lt m$. The remaining gnomon has $m^2-n^2$ cells. We must count the ways to place the integers $1,2,\dots,m^2-n^2$ so that every row increases from left to right and every column increases from top to bottom, then evaluate the count modulo $76543217$. Mathematical Approach Step 1: Write the Gnomon as a Partition Let $p=m-n$. As a Ferrers diagram, the gnomon is the partition $$\lambda=(\underbrace{m,\dots,m}_{p},\underbrace{p,\dots,p}_{n}).$$ The first $p$ rows have length $m$, while the last $n$ rows have length $p$. Therefore the total number of cells is $$|\lambda|=pm+np=(m-n)m+n(m-n)=m^2-n^2.$$ A valid numbering of the gnomon is exactly a standard Young tableau of this shape. Step 2: Apply the Hook-Length Formula If $T(m,n)$ denotes the number of valid numberings, the Frame-Robinson-Thrall formula gives $$T(m,n)=\frac{(|\lambda|)!}{\prod_{u \in \lambda} h(u)},$$ where $h(u)$ is the hook length of cell $u$: the number of cells to its right, plus the number below it, plus the cell itself. So the problem reduces to evaluating the hook product for this particular gnomon shape....

Detailed mathematical approach

Problem Summary

For integers $m,n$ with $0\le n\lt m$, let $L(m,n)$ be an $m\times m$ grid with the top-right $n\times n$ grid removed; it has $m^2-n^2$ cells. We number its cells with the consecutive integers $1,2,3,\dots$ so that the number in every cell is smaller than the number below it and to the left of it, and $\operatorname{LC}(m,n)$ is the number of valid numberings. The statement pictures $L(5,3)$ with two of its valid numberings and records $\operatorname{LC}(3,0)=42$, $\operatorname{LC}(5,3)=250250$, $\operatorname{LC}(6,3)=406029023400$ and $\operatorname{LC}(10,5)\bmod 76543217=61251715$. We must find $\operatorname{LC}(10000,5000)\bmod 76543217$. ^[1]

The whole solution in one sentence: turning the gnomon upside down and replacing every entry by its complement makes each valid numbering a standard Young tableau of one explicit partition, so the hook length formula counts them, and its hook product collapses into a short ratio of factorials evaluated modulo the prime $76543217$.

Mathematical Approach

Nothing in the final method searches or enumerates. We read off what the rule demands, recognise the counted objects as standard Young tableaux, then evaluate the formula in closed form. Write $p=m-n$ and $C=m^2-n^2$ throughout.

1. The rule on a shape with five cells

Write $(i,j)$ for the cell in row $i$ from the top and column $j$ from the left. Take $m=3$ and $n=2$. Removing the top-right $2\times2$ block from the $3\times3$ grid deletes $(1,2)$, $(1,3)$, $(2,2)$ and $(2,3)$, so $L(3,2)$ consists of $(1,1)$, $(2,1)$ and the whole bottom row $(3,1)$, $(3,2)$, $(3,3)$: five cells, and indeed $3^2-2^2=5$.

The rule compares a cell only with its neighbour below and its neighbour to the left. Here that gives two chains and nothing more: $v(1,1)\lt v(2,1)\lt v(3,1)$ down the first column and $v(3,3)\lt v(3,2)\lt v(3,1)$ along the bottom row, where $v$ denotes the number written in a cell. The entry at the corner $(3,1)$ exceeds every other entry of both chains, so it is $5$. The remaining entries come from $1,2,3,4$, and choosing which two of those go into the column fixes the numbering, because each chain is increasing. Hence $\operatorname{LC}(3,2)=\binom{4}{2}=6$.

Two things are worth keeping. Since the rule only compares neighbouring cells, a valid numbering is a listing of all cells in an order respecting those comparisons: a linear extension of a partial order. And the numbers increase downwards but decrease to the right, the mirror image of the usual convention; Section 2 removes that mismatch.

2. Turning the gnomon into a Young diagram

A Young diagram of a partition $\lambda=(\lambda_1\ge\lambda_2\ge\dots)$ is the set of cells $(i,j)$ with $1\le j\le\lambda_i$, drawn as left-justified rows of non-increasing length. A standard Young tableau of shape $\lambda$ fills it with $1,2,\dots,|\lambda|$, each once, increasing along every row and down every column. ^[3]

Row $i$ of $L(m,n)$ is short, of length $p$, exactly when $i\le n$, and full, of length $m$, when $i\gt n$: the rows grow downwards, the opposite of a Young diagram. So turn the picture upside down and reverse the numbers at the same time. For a numbering $v$ of $L(m,n)$ define

$$w(i,j)=C+1-v(m+1-i,\,j).$$

Claim. This is a bijection between the valid numberings of $L(m,n)$ and the standard Young tableaux of shape

$$\lambda=(\underbrace{m,\dots,m}_{p},\underbrace{p,\dots,p}_{n}).$$

Proof. The shapes match: row $m+1-i$ of $L(m,n)$ has length $m$ precisely when $m+1-i\gt n$, that is $i\le p$, and length $p$ otherwise, so $(i,j)$ lies in $\lambda$ exactly when $(m+1-i,j)$ lies in $L(m,n)$. This $\lambda$ has $p+n=m$ rows, is non-increasing because $p\le m$, and $|\lambda|=pm+np=m^2-n^2=C$. Since $x\mapsto C+1-x$ permutes $\{1,\dots,C\}$, the array $w$ uses each number once exactly when $v$ does. And $w(i,j)\lt w(i,j+1)$ says $v(m+1-i,j+1)\lt v(m+1-i,j)$, the rule "smaller than the number to the left of it" at $(m+1-i,j+1)$, while $w(i,j)\lt w(i+1,j)$ says $v(m-i,j)\lt v(m+1-i,j)$, the rule "smaller than the number below it" at $(m-i,j)$. Each rule appears once, and the map is its own inverse, so it is a bijection. $\Box$

For $m=5$, $n=3$ we get $p=2$ and $\lambda=(5,5,2,2,2)$ with $|\lambda|=16$. The tableau with rows $(1,2,3,4,5)$, $(6,7,8,9,10)$, $(11,12)$, $(13,14)$, $(15,16)$ is standard, and the bijection sends it to the numbering of $L(5,3)$ reading, from the top, $(2,1)$, $(4,3)$, $(6,5)$, $(11,10,9,8,7)$, $(16,15,14,13,12)$: every column increases downwards and every row decreases to the right, as the rule demands. For $n=0$ the shape is the $m\times m$ square, so $\operatorname{LC}(m,0)$ counts square standard Young tableaux. ^[4]

3. The hook length formula

The hook of a cell $(i,j)$ of $\lambda$ consists of that cell, the cells to its right in row $i$ and the cells below it in column $j$; its hook length $h(i,j)$ is the number of those cells. Writing $\lambda'$ for the conjugate partition, whose $j$-th part is the length of column $j$,

$$h(i,j)=(\lambda_i-j)+(\lambda'_j-i)+1.$$

The hook length formula, found in 1953 by Frame, Robinson and Thrall, states that the number of standard Young tableaux of shape $\lambda$ is ^[2]

$$f^{\lambda}=\frac{|\lambda|!}{\prod_{(i,j)\in\lambda}h(i,j)}.$$

Two checks against facts we already have. For the $3\times3$ square every hook length is $7-i-j$, so the hooks read $5,4,3$, then $4,3,2$, then $3,2,1$; their product is $60\cdot24\cdot6=8640$ and $9!/8640=42$, the statement's $\operatorname{LC}(3,0)$. For $n=m-1$ we have $p=1$ and $\lambda=(m,1,1,\dots,1)$, one row and one column crossing at a corner: the corner hook is $2m-1$, and the rest of the row and the rest of the column each contribute $m-1,m-2,\dots,1$, giving

$$f^{\lambda}=\frac{(2m-1)!}{(2m-1)\,\bigl((m-1)!\bigr)^2}=\binom{2m-2}{m-1}.$$

For $m=3$ this is $\binom{4}{2}=6$, the number counted by hand in Section 1; for $m=4,5$ it gives $20$ and $70$.

4. The hook lengths of the gnomon, in three blocks

Our $\lambda$ has $\lambda_i=m$ for $i\le p$ and $\lambda_i=p$ for $i\gt p$; counting columns gives $\lambda'_j=m$ for $j\le p$ and $\lambda'_j=p$ for $j\gt p$, so $\lambda'=\lambda$ and the shape is self-conjugate. Substituting into Section 3 splits the diagram into three blocks.

The corner square $i,j\le p$, row and column both full: $h(i,j)=(m-j)+(m-i)+1$${}=2m+1-i-j$.
The arm to its right, $i\le p$ and $j=p+b$ with $1\le b\le n$, row full and column short: $h=(m-p-b)+(p-i)+1$${}=m-i-b+1$.
The arm below it, $i=p+a$ with $1\le a\le n$ and $j\le p$, row short and column full: $h=(p-j)+(m-p-a)+1$${}=m-a-j+1$.

The last two expressions are one function of two indices, taken at $(i,b)$ and at $(j,a)$, and both pairs run over the rectangle $\{1,\dots,p\}\times\{1,\dots,n\}$. The arms therefore carry equal multisets of hook lengths — as self-conjugacy predicts — so

$$H=\prod_{(i,j)\in\lambda}h(i,j)=P_A\,P_B^{2},$$

with $P_A$ the product over the corner square and $P_B$ the product over one arm.

5. Each row collapses into a ratio of factorials

Fix a row $i\le p$ of the corner square. As $j$ runs from $1$ to $p$, the value $2m+1-i-j$ runs through the consecutive integers from $2m-i$ down to $2m+1-i-p=m+n+1-i$, so that row contributes

$$\prod_{j=1}^{p}(2m+1-i-j)=\frac{(2m-i)!}{(m+n-i)!}.$$

Multiplying over $i=1,\dots,p$, the numerators are the factorials of $2m-1,\dots,2m-p=m+n$ and the denominators those of $m+n-1,\dots,m+n-p=2n$:

$$\begin{aligned}P_A&=\prod_{i=1}^{p}\frac{(2m-i)!}{(m+n-i)!}\\&=\frac{\prod_{t=m+n}^{2m-1}t!}{\prod_{s=2n}^{m+n-1}s!}.\end{aligned}$$

One row of an arm behaves the same way: for fixed $i\le p$ the values $m-i-b+1$ with $1\le b\le n$ are the consecutive integers from $m-i$ down to $m-i-n+1=p-i+1$, all positive, so

$$\begin{aligned}P_B&=\prod_{i=1}^{p}\frac{(m-i)!}{(p-i)!}\\&=\frac{\prod_{t=n}^{m-1}t!}{\prod_{u=0}^{p-1}u!}.\end{aligned}$$

Here $0!=1$ and an empty product is $1$. For $n=0$ both arms are empty, $p=m$, and the two products in $P_B$ coincide, so $P_B=1$ and only the corner square remains. With Section 3 this is the closed form the programs evaluate:

$$\operatorname{LC}(m,n)=\frac{(m^2-n^2)!}{P_A\,P_B^{2}}.$$

The gain is dramatic: the diagram has $m^2-n^2$ cells, but $P_A$ and $P_B$ need only $4p$ factorials of arguments below $2m$ — $20000$ of them here — each read from one table of $2m$ entries.

6. The two exact checkpoints worked out

For $m=5$, $n=3$ we have $p=2$, $m+n=8$, $2m-1=9$ and $2n=6$, so

$$P_A=\frac{8!\,9!}{6!\,7!}=4032,$$

$$P_B=\frac{3!\,4!}{0!\,1!}=144.$$

Then $H=4032\cdot144^2=83607552$, and with $16!=20922789888000$ the formula gives $\operatorname{LC}(5,3)=250250$, the value in the statement. For $m=6$, $n=3$ the same computation gives $P_A=359251200$ and $P_B=8640$, hence $H=26817958379520000$, and dividing $27!$${}=10888869450418352160768000000$ by it leaves $406029023400$, again as stated.

A third fact detects wrong formulas cheaply. The largest entry of a square tableau must sit in its bottom-right corner, since every other cell has a neighbour to the right or below; for $m\ge2$, deleting that corner leaves $\lambda$ for $n=1$, and adjoining it again with the entry $m^2$ restores the square. Hence $\operatorname{LC}(m,1)=\operatorname{LC}(m,0)$ for every $m\ge2$, for instance $\operatorname{LC}(4,1)=\operatorname{LC}(4,0)=24024$.

7. Evaluating the quotient modulo 76543217

The modulus $M=76543217$ is prime, as trial division by every integer from $2$ to $\lfloor\sqrt M\rfloor=8748$ confirms. Dividing by a nonzero residue $x$ means multiplying by a residue $y$ with $xy\equiv1\pmod M$, which undoes the multiplication by $x$; modulo $7$, for example, $y=5$ inverts $x=3$, since $3\cdot5=15\equiv1$. Fermat's little theorem gives $x^{M-1}\equiv1\pmod M$ whenever $M\nmid x$, ^[5] so $x\cdot x^{M-2}\equiv1$ and

$$x^{-1}\equiv x^{M-2}\pmod M,$$

so an inverse costs one modular exponentiation. ^[6] The square-and-multiply loops of these programs consume one binary digit of the exponent per iteration, squaring the running base every time, so with $27$ binary digits in $M-2$, $18$ of them $1$, each call costs $27$ squarings and $18$ multiplications. ^[7] The computation is therefore

$$\operatorname{LC}(m,n)\equiv C!\cdot H^{M-2}\pmod M.$$

This needs an argument, because reducing numerator and denominator separately is not always allowed. Every hook length lies between $1$ and $h(1,1)=2m-1$, here $19999$ and far below $M$, so no hook length is a multiple of the prime $M$ and $H$ is invertible; the same bound covers the factorials inside $P_A$ and $P_B$, whose arguments stay below $2m$. And $C=10000^2-5000^2=75000000$ is smaller than $M$, by $1543217$, so $C!$ carries no factor $M$ either and the reduced value is not $0$. The checkpoint $\operatorname{LC}(10,5)\bmod M=61251715$ tests this pipeline on $75$ cells.

8. The cost that remains, and a shortcut not taken

After Section 5 the hook side costs $O(m)$ operations, so the only expensive object left is $C!\bmod M$, which all three programs obtain from the obvious loop of $C$ modular multiplications. A shortcut does exist. Wilson's theorem says that $(M-1)!\equiv-1\pmod M$ for a prime $M$, since every nonzero residue pairs off with its inverse and only $1$ and $M-1$ are their own inverses. ^[10] Splitting that product at $C$ gives

$$C!\cdot\prod_{k=C+1}^{M-1}k\equiv-1\pmod M,$$

so $C!$ is $-1$ times the inverse of the tail product, and the tail has only $M-1-C=1543216$ factors, about $48$ times fewer than $C$. The published programs do not take this route; a separate implementation of it returns the same value as they do.

How the Code Works

The C++ program ^[8] keeps the modulus in MOD and the parameters in an Options struct defaulting to $m=10000$, $n=5000$ with the checkpoints enabled. parse_int_after_prefix reads the integer after a prefix such as --m= and reports failure instead of throwing, and parse_arguments uses it for --m= and --n=, accepts --skip-checkpoints and enforces $m\ge1$ and $0\le n\lt m$. It does not also require $2m-1\lt M$, the bound Section 7 needs for the division; the published parameters satisfy it with room to spare, but for accepted inputs beyond it LC_mod is not generally correct. mod_pow is the square-and-multiply loop of Section 7.

LC_mod is Sections 5 and 7 in order. It builds fact_small with the factorials of $0,1,\dots,2m-1$ modulo MOD and accumulates fact_cells in one loop over all cells; it then forms pa_num, pa_den, pb_num and pb_den as the four products of Section 5, divides through mod_pow(..., MOD - 2), multiplies hooks together as $P_A P_B^2$, and returns fact_cells times the inverse of hooks. Three exponentiations are used where one would do; at $27$ squarings and $18$ multiplications each, that is invisible beside the $75000000$ multiplications of the factorial loop. The factor (i % MOD) in that loop changes nothing here, every $i$ being below $M$; it keeps the loop correct for larger $m$, where a multiple of $M$ among the cells makes fact_cells zero — and zero is then the right answer, because $H$ stays coprime to $M$ as long as $2m-1\lt M$.

LC_exact_small is the independent checker. It builds row_len and col_len explicitly, reads each hook length off them as (row_len[i] - j) + (col_len[j] - i) - 1 in $0$-based indices — the $-1$ because the cell itself is counted twice — and divides the exact factorial by the exact hook product in Boost's arbitrary-precision cpp_int. ^[9] run_checkpoints calls it for $\operatorname{LC}(3,0)$, $\operatorname{LC}(5,3)$ and $\operatorname{LC}(6,3)$, then checks $\operatorname{LC}(10,5)\bmod M$ through LC_mod, so both paths run before main prints the result.

The Python and Java versions implement LC_mod only, with the parameters fixed in the source: solve in Python mirrors the C++ function statement by statement with its own mod_pow, and Euler412.solve in Java does the same over an int[] factSmall. Neither port carries the exact checks, a choice rather than a limitation: Python's integers are unbounded and would need no library, while Java would want BigInteger. C++ and Java both store residues in $32$-bit fields, and since $M\lt2^{27}$ a product of two residues stays below $2^{54}$, so widening one operand — 1LL * in C++, 1L * in Java — holds every intermediate inside the signed $64$-bit range, whose limit exceeds $(M-1)^2=5858863915622656$. Python needs no such care, its list holding references to integer objects rather than packed words.

Complexity and Verification

The factorial loop performs $C=m^2-n^2$ modular multiplications, $75000000$ for the target, and dominates everything else: $2m-1=19999$ for the factorial table, $4p=20000$ for the four products, and three exponentiations of $27$ squarings and $18$ multiplications each. The method is $O(m^2)$ in time and $O(m)$ in memory, that memory being a table of $2m=20000$ entries — $80000$ bytes of $32$-bit fields in C++ and Java, more in Python, where the list holds object references. The exact checker costs $O(C)$ big-integer multiplications on numbers of $O(C\log C)$ bits, so it is kept for the small checkpoints; at $27$ cells it is instant. All three programs print the same value, the C++ and Java ones reaching it in a fraction of a second and Python in a few seconds on the machine used here.

Beyond the packaged checkpoints, these were verified independently for this article. A direct count of linear extensions of the gnomon order, built straight from the words of the statement, agrees with the hook length formula applied to $\lambda$ for every $m\le5$ and $0\le n\lt m$, including $\operatorname{LC}(5,3)=250250$ and the hand count $\operatorname{LC}(3,2)=6$; reading the rule as "below and to the right" instead yields $784$ for $L(5,3)$, so the orientation of Section 2 is the one meant. A formula-free count of all paths from the empty diagram to $\lambda$ in Young's lattice matches the formula for every $m\le12$, and the modular pipeline, reimplemented separately, reproduces the exact value modulo $M$ for every $m\le45$ and $0\le n\lt m$. The block formulas of Section 4 were checked cell by cell against the hooks of $\lambda$ and $\lambda'$, together with $\lambda'=\lambda$, equal arm multisets and $\max h=2m-1$, for every $m\lt30$. The identity $\operatorname{LC}(m,1)=\operatorname{LC}(m,0)$ was checked for $2\le m\le6$, the square counts $\operatorname{LC}(m,0)$ for $m=1,\dots,6$ being 1, 2, 42, 24024, 701149020, 1671643033734960. Finally the Wilson route of Section 8 was implemented separately and returns, at the target parameters, exactly what the three published programs print. ^[4]

Footnotes and References

The references supply the background facts used above. The bijection of Section 2, the three-block decomposition of the hook product, its collapse into factorials and the invertibility argument of Section 7 are derived in this article.

Project Euler 412 — Gnomon Numbering. The official statement defines L(m, n) as an m × m grid with the top-right n × n grid removed, pictures L(5, 3) with two valid numberings, asks for numberings by consecutive integers in which every cell holds a number smaller than the number below it and to the left of it, writes LC(m, n) for their count, supplies the four checkpoints quoted in the summary, and asks for LC(10000, 5000) modulo 76543217.
Wikipedia — Hook length formula. The hook of a cell consists of the cell itself, the cells to its right in its row and the cells below it in its column, and the number of standard Young tableaux of shape λ is n! divided by the product of all hook lengths. The page credits the discovery to Frame, Robinson and Thrall in 1953.
Wikipedia — Young tableau. A Young, or Ferrers, diagram is a collection of boxes in left-justified rows of non-increasing length; counting the boxes in each column gives the conjugate partition, whose diagram is the reflection of the original in its main diagonal; and a filling by 1 to n is standard when every row and every column increases.
OEIS A039622 — number of n × n Young tableaux. For n ≥ 1 the sequence counts arrangements of 1, 2, …, n² in an n × n array with every row and every column increasing, which Section 2 identifies with LC(n, 0); its values for n = 1, …, 6 are 1, 2, 42, 24024, 701149020, 1671643033734960, listed there after an initial term for n = 0.
Wikipedia — Fermat's little theorem. For a prime p and any integer a, the number a^p − a is a multiple of p; when a is not divisible by p this is equivalent to a^(p−1) ≡ 1 (mod p), the form used in Section 7.
Wikipedia — Modular multiplicative inverse. For a prime modulus m the page gives the inverse explicitly as a^(m−2) mod m, and notes that this is an alternative to the extended Euclidean algorithm, used when modular exponentiation is already available — the trade-off these programs accept.
Wikipedia — Exponentiation by squaring. The iterative method consumes the exponent from its least significant bit, using ⌊log₂ n⌋ squarings and at most as many further multiplications, and notes its use for large exponents modulo a number. The published mod_pow and modPow are this loop, squaring the base on every iteration, the last one included.
C++ — Euler412.cpp. The linked immutable revision contains parse_arguments, mod_pow, LC_mod, LC_exact_small and run_checkpoints, and includes Boost's cpp_int header. It is the implementation record for this article; the independent counts above are not part of it, nor is the Wilson route of Section 8.
Boost.Multiprecision — cpp_int. The page documents cpp_int as a signed integer type of arbitrary precision that expands as needed, as LC_exact_small needs for 27! and the exact hook product.
Wikipedia — Wilson's theorem. For a natural number n > 1 the product of all positive integers below n is one less than a multiple of n exactly when n is prime, that is (p − 1)! ≡ −1 (mod p); the elementary proof pairs each nonzero residue with its inverse, only 1 and p − 1 inverting themselves. Section 8 splits that product.

Problem 412 source code

C++

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

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using boost::multiprecision::cpp_int;

constexpr int MOD = 76543217;

struct Options {
    int m = 10000;
    int n = 5000;
    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;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    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 >= 0 && options.n < options.m;
}

int mod_pow(i64 base, i64 exp) {
    i64 result = 1;
    i64 cur = base % MOD;
    i64 e = exp;
    while (e > 0) {
        if (e & 1LL) {
            result = result * cur % MOD;
        }
        cur = cur * cur % MOD;
        e >>= 1LL;
    }
    return static_cast<int>(result);
}

int LC_mod(const int m, const int n) {
    const int p = m - n;
    const i64 cells = 1LL * m * m - 1LL * n * n;

    const int max_small = 2 * m - 1;
    std::vector<int> fact_small(static_cast<std::size_t>(max_small + 1), 1);
    for (int i = 1; i <= max_small; ++i) {
        fact_small[static_cast<std::size_t>(i)] =
            static_cast<int>(1LL * fact_small[static_cast<std::size_t>(i - 1)] * i % MOD);
    }

    int fact_cells = 1;
    for (i64 i = 1; i <= cells; ++i) {
        fact_cells = static_cast<int>(1LL * fact_cells * (i % MOD) % MOD);
    }

    // Hook-product decomposition for shape λ = (m repeated p times, p repeated n times):
    // P = P_A * P_B^2, where
    // P_A = prod_{t=m+n}^{2m-1} t! / prod_{s=2n}^{m+n-1} s!
    // P_B = prod_{t=n}^{m-1} t! / prod_{u=0}^{p-1} u!
    int pa_num = 1;
    for (int t = m + n; t <= 2 * m - 1; ++t) {
        pa_num = static_cast<int>(1LL * pa_num * fact_small[static_cast<std::size_t>(t)] % MOD);
    }
    int pa_den = 1;
    for (int s = 2 * n; s <= m + n - 1; ++s) {
        pa_den = static_cast<int>(1LL * pa_den * fact_small[static_cast<std::size_t>(s)] % MOD);
    }
    const int pa = static_cast<int>(1LL * pa_num * mod_pow(pa_den, MOD - 2) % MOD);

    int pb_num = 1;
    for (int t = n; t <= m - 1; ++t) {
        pb_num = static_cast<int>(1LL * pb_num * fact_small[static_cast<std::size_t>(t)] % MOD);
    }
    int pb_den = 1;
    for (int u = 0; u <= p - 1; ++u) {
        pb_den = static_cast<int>(1LL * pb_den * fact_small[static_cast<std::size_t>(u)] % MOD);
    }
    const int pb = static_cast<int>(1LL * pb_num * mod_pow(pb_den, MOD - 2) % MOD);

    int hooks = pa;
    hooks = static_cast<int>(1LL * hooks * pb % MOD);
    hooks = static_cast<int>(1LL * hooks * pb % MOD);

    return static_cast<int>(1LL * fact_cells * mod_pow(hooks, MOD - 2) % MOD);
}

cpp_int LC_exact_small(const int m, const int n) {
    const int p = m - n;
    const int rows = m;
    std::vector<int> row_len(static_cast<std::size_t>(rows), p);
    for (int i = 0; i < p; ++i) {
        row_len[static_cast<std::size_t>(i)] = m;
    }

    std::vector<int> col_len(static_cast<std::size_t>(m), p);
    for (int j = 0; j < p; ++j) {
        col_len[static_cast<std::size_t>(j)] = m;
    }

    int cells = 0;
    for (int i = 0; i < rows; ++i) {
        cells += row_len[static_cast<std::size_t>(i)];
    }

    cpp_int numer = 1;
    for (int i = 1; i <= cells; ++i) {
        numer *= i;
    }

    cpp_int denom = 1;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < row_len[static_cast<std::size_t>(i)]; ++j) {
            const int hook =
                (row_len[static_cast<std::size_t>(i)] - j) + (col_len[static_cast<std::size_t>(j)] - i) - 1;
            denom *= hook;
        }
    }

    return numer / denom;
}

bool run_checkpoints() {
    if (LC_exact_small(3, 0) != 42) {
        std::cerr << "Checkpoint failed: LC(3,0)\n";
        return false;
    }
    if (LC_exact_small(5, 3) != cpp_int("250250")) {
        std::cerr << "Checkpoint failed: LC(5,3)\n";
        return false;
    }
    if (LC_exact_small(6, 3) != cpp_int("406029023400")) {
        std::cerr << "Checkpoint failed: LC(6,3)\n";
        return false;
    }
    if (LC_mod(10, 5) != 61251715) {
        std::cerr << "Checkpoint failed: LC(10,5) mod 76543217\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;
    }

    std::cout << LC_mod(options.m, options.n) << '\n';
    return 0;
}

Python

def solve():
    MOD = 76543217
    m, n = 10000, 5000
    p = m - n
    cells = m * m - n * n

    def mod_pow(base, exp, mod=MOD):
        r = 1; base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod
            exp >>= 1
        return r

    max_small = 2 * m - 1
    fact_small = [1] * (max_small + 1)
    for i in range(1, max_small + 1):
        fact_small[i] = fact_small[i-1] * i % MOD

    fact_cells = 1
    for i in range(1, cells + 1):
        fact_cells = fact_cells * (i % MOD) % MOD

    pa_num = 1
    for t in range(m + n, 2 * m):
        pa_num = pa_num * fact_small[t] % MOD
    pa_den = 1
    for s in range(2 * n, m + n):
        pa_den = pa_den * fact_small[s] % MOD
    pa = pa_num * mod_pow(pa_den, MOD - 2) % MOD

    pb_num = 1
    for t in range(n, m):
        pb_num = pb_num * fact_small[t] % MOD
    pb_den = 1
    for u in range(p):
        pb_den = pb_den * fact_small[u] % MOD
    pb = pb_num * mod_pow(pb_den, MOD - 2) % MOD

    hooks = pa * pb % MOD * pb % MOD
    return str(fact_cells * mod_pow(hooks, MOD - 2) % MOD)

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

Java

public class Euler412 {
    private static final int MOD = 76543217;

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

    public static String solve() {
        int m = 10000;
        int n = 5000;
        int p = m - n;
        long cells = 1L * m * m - 1L * n * n;

        int maxSmall = 2 * m - 1;
        int[] factSmall = new int[maxSmall + 1];
        factSmall[0] = 1;
        for (int i = 1; i <= maxSmall; ++i) {
            factSmall[i] = (int) (1L * factSmall[i - 1] * i % MOD);
        }

        int factCells = 1;
        for (long i = 1; i <= cells; ++i) {
            factCells = (int) (1L * factCells * (i % MOD) % MOD);
        }

        int paNum = 1;
        for (int t = m + n; t <= 2 * m - 1; ++t) {
            paNum = (int) (1L * paNum * factSmall[t] % MOD);
        }
        int paDen = 1;
        for (int s = 2 * n; s <= m + n - 1; ++s) {
            paDen = (int) (1L * paDen * factSmall[s] % MOD);
        }
        int pa = (int) (1L * paNum * modPow(paDen, MOD - 2) % MOD);

        int pbNum = 1;
        for (int t = n; t <= m - 1; ++t) {
            pbNum = (int) (1L * pbNum * factSmall[t] % MOD);
        }
        int pbDen = 1;
        for (int u = 0; u <= p - 1; ++u) {
            pbDen = (int) (1L * pbDen * factSmall[u] % MOD);
        }
        int pb = (int) (1L * pbNum * modPow(pbDen, MOD - 2) % MOD);

        int hooks = pa;
        hooks = (int) (1L * hooks * pb % MOD);
        hooks = (int) (1L * hooks * pb % MOD);

        long ans = 1L * factCells * modPow(hooks, MOD - 2) % MOD;
        return String.valueOf(ans);
    }

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