Problem 1016: Squareful Matrix

Detailed mathematical approach

Problem Summary

An \(n\times n\) matrix is called squareful if every number from \(1\) to \(n^2\) appears exactly once, every two consecutive numbers stand in cells that share a side, and, for every \(k\lt n\), the numbers \(1\) to \(k^2\) fill a contiguous \(k\times k\) square inside the matrix. The matrix

$$\begin{pmatrix}13&14&15&16\\12&1&4&5\\11&2&3&6\\10&9&8&7\end{pmatrix}$$

is squareful: \(1\) alone is a \(1\times1\) square, \(1,2,3,4\) fill a \(2\times2\) square, \(1\) to \(9\) fill a \(3\times3\) square, and each number sits next to the one before it. Its trace, the sum along the main diagonal from the top left to the bottom right, is \(13+1+3+7=24\). ^[2] Let \(T(n)\) be the sum of the traces of all squareful matrices of size at most \(n\times n\). We are given \(T(2)=41\), \(T(12)=9063089\) and \(T(123)\equiv1139320302\)\({}\ (\mathrm{mod}\ 1234567891)\), and we must find \(T(123456789)\) modulo \(1234567891\). ^[1]

Here is the plan. A squareful matrix is nothing but a chain of nested squares, each grown out of the previous one by adding an L-shaped ring at one of its four corners, and for \(n\ge2\) there are only \(2^{n+1}\) of them. The trace cannot see where the diagonal cells sit inside a ring, only which ring they belong to. That splits the matrices into two families of equal size: in one of them every matrix has the same trace, while in the other the traces are governed by a race between two counters. Adding everything up leaves a formula built from \(2^n\), one central binomial coefficient and the running total of central binomial coefficients, so \(T(n)\) costs a single pass of about \(n/2\) steps. The sections below build this up from the smallest cases, so that every step can be checked by hand. Size \(1\), which has a single matrix instead of four, is treated separately throughout.

Mathematical Approach

We first describe every squareful matrix exactly, then find what the trace depends on, then count.

1. Nested squares

For \(1\le k\le n\) let \(B_k\) be the set of cells holding the numbers \(1\) to \(k^2\). The third condition says that \(B_k\) is a contiguous \(k\times k\) square for every \(k\lt n\), and \(B_n\) is the whole matrix, so this holds for every \(k\). Since the numbers \(1\) to \(k^2\) are among the numbers \(1\) to \((k+1)^2\), we have \(B_k\subseteq B_{k+1}\). A \(k\times k\) square inside a \((k+1)\times(k+1)\) square must touch two of its sides, so:

\(B_k\) occupies one of the four corners of \(B_{k+1}\).

What is left over, \(R_k=B_{k+1}\setminus B_k\), is an L-shaped ring of \(2k+1\) cells: one full row and one full column of \(B_{k+1}\), meeting at the corner of \(B_{k+1}\) opposite to \(B_k\). In the example above, \(B_2\) is the square holding \(1,2,3,4\), and \(R_2\) is the ring \(5,6,7,8,9\).

2. Every ring is a bent path

A ring is one cell wide, so if we join cells that share a side, \(R_k\) is a path: a chain of \(2k+1\) cells with two ends and no branches. A numbering in which consecutive numbers share a side is a Hamiltonian path of the grid of cells, and on a chain such a numbering has nowhere to start but one of the two ends. ^[3] The numbers \(k^2+1,\dots,(k+1)^2\) fill the ring and each stands next to the one before, so they march along this chain from one of its ends to the other. If \(B_k\) sits at the top left of \(B_{k+1}\), the ring is the rightmost column plus the bottom row, and its two ends are the top-right and the bottom-left corner of \(B_{k+1}\), with the bend at the bottom-right corner. In general the two ends of \(R_k\) are the two corners of \(B_{k+1}\) next to the corner occupied by \(B_k\), and the bend is the corner opposite to it. So:

\(k^2+1\) sits at one end of the ring and \((k+1)^2\) at the other end.

3. Where the largest number sits

Each end of the ring has exactly one neighbour inside \(B_k\). With \(B_k\) at the top left, the end at the top-right corner of \(B_{k+1}\) touches only the top-right corner of \(B_k\), and the end at the bottom-left corner touches only the bottom-left corner of \(B_k\). Since \(k^2+1\) must stand next to \(k^2\), and \(k^2\) lies in \(B_k\), we learn two things at once. First, \(k^2\) sits in a corner of \(B_k\) — one of the two corners of \(B_k\) that touch the corner shared with \(B_{k+1}\). Second, the march across the ring ends at the other end. To say what that means, name each corner of a square by a pair of bits \((u,v)\), where \(u=0\) stands for the top edge and \(u=1\) for the bottom edge, \(v=0\) for the left edge and \(v=1\) for the right edge. Then:

if \(k^2\) sits in the corner \((u,v)\) of \(B_k\), with \(k\ge2\), then \((k+1)^2\) sits in the corner \((1-u,1-v)\) of \(B_{k+1}\),

the corner opposite to it; for \(k=1\) the rule needs no bits, because the single cell of \(B_1\) is all four of its corners at once. The two corners of \(B_k\) that touch the corner shared with \(B_{k+1}\) are the two that differ from the shared one in exactly one bit.

In the example, \(4\) sits at the top-right corner of \(B_2\), that is \((0,1)\); then \(9\) sits at the bottom-left corner of \(B_3\), which is \((1,0)\); then \(16\) sits at the top-right corner of \(B_4\), which is \((0,1)\) again. Both bits flip at every step.

4. How many squareful matrices there are

Write \(c_k\) for the corner of \(B_k\) holding \(k^2\) and \(X_k\) for the corner of \(B_{k+1}\) occupied by \(B_k\). Section 3 says \(c_{k+1}\) is the opposite of \(c_k\), and that \(c_k\) must be one of the two corners of \(B_k\) that touch \(X_k\), which means \(X_k\) differs from \(c_k\) in exactly one bit. So once \(c_k\) is known, there are exactly two choices for \(X_k\), and either choice fixes everything about the ring: which cells it uses, from which end the march starts, and therefore all of its numbers. Going from \(B_1\) to \(B_n\) there are \(n-1\) steps. The first step is special because \(B_1\) is a single cell, which is all four of its own corners at once: all four placements of \(B_1\) in \(B_2\) are allowed, and both ends of \(R_1\) touch the single cell, so the first step has \(8\) outcomes instead of \(2\). Hence there are \(8\cdot2^{n-2}=2^{n+1}\) squareful matrices of size \(n\ge2\), and exactly one of size \(1\). A direct search over all self-avoiding walks confirms the counts \(8,16,32,64,128\) for \(n=2,\dots,6\).

5. Two growth patterns

Everything in this section and the next five assumes \(n\ge2\); size \(1\) is handled on its own. Give a corner \((u,v)\) the parity of \(u+v\): the top-left and bottom-right corners, \((0,0)\) and \((1,1)\), are even, the top-right and bottom-left corners are odd. Passing from \(c_k\) to \(c_{k+1}\) flips both bits, which keeps the parity; passing from \(c_k\) to \(X_k\) flips exactly one bit, which changes it. The chain starts at \(k=1\): \(c_1\) is meaningless, since the single cell of \(B_1\) is all four of its own corners, but \(1\) and \(4\) always sit in neighbouring corners of \(B_2\) — a march along the four cells of \(B_2\) ends next to where it began — so \(X_1\) differs from \(c_2\) in exactly one bit as well. Hence all the \(c_k\) with \(k\ge2\) share one parity, all the \(X_k\) share the other one, and there are only two kinds of squareful matrix:

• diagonal growth: every \(B_k\) sits at the top-left or the bottom-right corner of \(B_{k+1}\), and \(n^2\) sits at the top-right or the bottom-left corner of the matrix;
• anti-diagonal growth: every \(B_k\) sits at the top-right or the bottom-left corner of \(B_{k+1}\), and \(n^2\) sits at the top-left or the bottom-right corner.

Each kind holds \(2^n\) of the \(2^{n+1}\) matrices, because the kind is decided by which of the four corners holds \(n^2\). The example has \(16\) at the top right, so it has diagonal growth, and indeed \(B_1,B_2\) sit at top-left corners and \(B_3\) at a bottom-right corner.

6. What the trace sees

Put \(B_0=\varnothing\), so that \(R_0=B_1\) is the single cell holding \(1\) and the rings \(R_0,R_1,\dots,R_{n-1}\) cover the matrix. Let \(f_k\) be how many cells of the main diagonal of the whole matrix lie inside \(B_k\), with \(f_0=0\) and \(f_n=n\). The ring \(R_k\) therefore contains \(f_{k+1}-f_k\) diagonal cells. Write

$$g(k)=k^2+k+1.$$

Numbers on a ring fall off by one with every step away from the cell of \((k+1)^2\), and a step away always changes the row or the column by one, so the number in a ring cell is \((k+1)^2\) minus its row-and-column distance from the cell of \((k+1)^2\). With rows and columns of \(B_{k+1}\) numbered \(0\) to \(k\), this gives two cases. Under diagonal growth \(B_k\) and \(B_{k+1}\) share a corner that lies on the main diagonal of the matrix, so the main diagonals of all the nested squares lie on one and the same line, the main diagonal of the matrix. The ring then meets that line only in its bend, which is at distance \(k\) from either end, so that cell holds \((k+1)^2-k=g(k)\); and \((k+1)^2\) sits at a top-right or bottom-left corner, by Section 5. Under anti-diagonal growth \((k+1)^2\) sits at a top-left or bottom-right corner, so the number in the cell with row \(y\) and column \(x\) is \((k+1)^2-(y+x)\) or \(k^2+1+(y+x)\); the main diagonal of the matrix meets \(B_{k+1}\) in a line of cells with \(x-y\) constant, and that line meets the ring in cells with \(y+x\) equal to \(\lvert x-y\rvert\) or \(2k-\lvert x-y\rvert\), one from each arm. Those two numbers add up to \(2g(k)\) in both cases. A ring may also miss the line completely, and when the line only grazes the corner of \(B_{k+1}\) the single cell it meets has \(y+x=k\) and holds \(g(k)\). So, however the rings lie:

the numbers of a ring that fall on the main diagonal add up to \(g(k)\) times how many there are.

Writing \(\tau\) for the trace, it is therefore

$$\tau=\sum_{k=0}^{n-1}(f_{k+1}-f_k)\,g(k),$$

and since \(g(k)-g(k-1)=2k\), regrouping the sum gives the form we shall use:

$$\tau=n(n^2-n+1)-2\sum_{k=1}^{n-1}k\,f_k.$$

Nothing but the counts \(f_1,\dots,f_{n-1}\) survives. Where a diagonal cell sits inside its ring, and which of the two arms it is on, do not matter at all.

7. Diagonal growth: one trace for all

Under diagonal growth \(B_k\) and \(B_{k+1}\) share a corner on the main diagonal, so the main diagonal of \(B_k\) is part of the main diagonal of the whole matrix and \(f_k=k\) for every \(k\). Every ring then contains exactly one diagonal cell, and

$$\tau=\sum_{k=0}^{n-1}g(k)=\frac{n(n^2+2)}{3}.$$

For the example this is \(1+3+7+13=24\), which is \(4\cdot18/3\), the trace we were given. All \(2^n\) matrices with diagonal growth have this one trace, so together they contribute \(S_1(n)=2^n n(n^2+2)/3\) to the sum of all traces.

8. Anti-diagonal growth: peeling with two counters

For the other family it is easier to run the chain backwards and peel the matrix down to its centre. Under anti-diagonal growth \(B_k\) sits at the top-right or the bottom-left corner of \(B_{k+1}\), so each peeling step removes either the top row together with the rightmost column, or the leftmost column together with the bottom row. Let \(r\) count the steps of the first kind so far and \(s\) those of the second kind. Rows and columns of the matrix numbered \(0\) to \(n-1\), after \(r+s\) steps what is left is

$$B_k=\{(y,x):r\le y\le n-1-s,\ s\le x\le n-1-r\},\qquad k=n-r-s,$$

where \(y\) is the row and \(x\) the column of a cell.

A diagonal cell in row and column \(i\) has survived exactly when \(r\le i\le n-1-s\) and \(s\le i\le n-1-r\), that is when

$$\max(r,s)\le\min(i,\,n-1-i).$$

So the surviving diagonal cells always form a block in the middle of the diagonal, of length \(n-2\max(r,s)\), which gives \(f_k=\max(0,\,n-2\max(r,s))\). A step shaves one cell off each end of that block when it raises the larger of the two counters, and shaves nothing when it raises the smaller one. Writing \(d=\min(i,n-1-i)\) for the distance of a diagonal cell from the nearer end of the diagonal:

a cell at distance \(d\) leaves at the first step at which the larger counter passes \(d\).

Apart from the middle cell of an odd size, the diagonal cells come in symmetric pairs at the same distance, and the two cells of a pair always leave together. A step that raises the larger counter removes the two cells at the ends of the surviving block, or the lone middle cell if that is all there is; a step that raises the smaller counter removes nothing, and once the block is empty no step removes anything.

9. When a diagonal cell leaves

If a cell leaves at step number \(i\), it was inside \(B_{n-i+1}\) and is outside \(B_{n-i}\), so it belongs to the ring \(R_{n-i}\). Section 6 does not say that this one cell holds \(g(n-i)\): it says that a departing pair holds \(2g(n-i)\) between them, and a departing lone cell holds \(g(n-i)\). The squareful matrix

$$\begin{pmatrix}4&3\\1&2\end{pmatrix}$$

has anti-diagonal growth, and its two diagonal cells hold \(4\) and \(2\), not \(3\) and \(3\), although \(4+2=2g(1)\). Since we only ever add traces, we may give each departing cell the accounting weight \(g(n-i)\). The distances \(d\) that occur are read off the diagonal: for \(n=2M\) each of \(d=0,1,\dots,M-1\) belongs to two cells; for \(n=2M+1\) the same is true for \(d=0,\dots,M-1\) while the middle cell alone has \(d=M\).

10. Counting the peeling orders

A peeling order is a sequence of \(n-1\) choices, so one corner of the matrix carries \(2^{n-1}\) of them, and the family has \(2\cdot2^{n-1}=2^n\) members. For a fixed \(d\), in how many orders does the larger counter first pass \(d\) exactly at step \(i\)? The counter that gets there has \(d+1\) steps, the last of them being step \(i\); of the earlier \(i-1\) steps it owns \(d\), in any positions; and it can be either of the two counters. That is \(2\binom{i-1}{d}\) prefixes of length \(i\), and they exist for \(d+1\le i\le 2d+1\), because the other counter owns \(i-1-d\) steps and must not have passed \(d\) itself. ^[5] The remaining \(n-1-i\) choices are free, which multiplies by \(2^{n-1-i}\) and turns each prefix into that many whole orders. For even \(n=2M\) every distance satisfies \(2d+1\le n-1\), so every diagonal cell does leave, and the whole family contributes

$$S_2(n)=8\sum_{d=0}^{M-1}\ \sum_{i=d+1}^{2d+1}\binom{i-1}{d}2^{\,n-1-i}g(n-i),$$

where the factor \(8\) is two corners of the matrix, times two cells at each distance, times the two counters. For an odd size \(n=2M+1\) the distances \(d\le M-1\) behave in exactly the same way, but the middle cell, which is alone at its distance \(d=M\), needs two terms of its own: since \(2M+1\) exceeds the \(2M\) steps, the middle cell either leaves at some step \(i\) with \(M+1\le i\le2M\) or survives to the very end, and it survives in \(\binom{2M}{M}\) of the orders of each corner, contributing \(g(0)=1\). Altogether, for \(M\ge1\),

$$\begin{aligned}S_2(2M+1)={}&8\sum_{d=0}^{M-1}\ \sum_{i=d+1}^{2d+1}\binom{i-1}{d}2^{\,2M-i}g(2M+1-i)\\&+4\sum_{i=M+1}^{2M}\binom{i-1}{M}2^{\,2M-i}g(2M+1-i)+2\binom{2M}{M},\end{aligned}$$

where the middle terms carry a \(4\) and a \(2\) instead of an \(8\) because the middle cell has no partner.

11. Collapsing the double sum

Exchanging the two sums turns the inner one into a partial row of Pascal's triangle: ^[6] the condition \(i\le2d+1\) is \(d\ge(i-1)/2\), so for a fixed \(i\) the sum runs over \(\binom{i-1}{d}\) for \(\lceil(i-1)/2\rceil\le d\le\min(i-1,M-1)\). Writing \(L=i-1\) and

$$A_L=\sum_{d=\lceil L/2\rceil}^{\min(L,\,M-1)}\binom{L}{d},$$

the double sum for an even size \(n=2M\) becomes the single sum

$$S_2(n)=8\sum_{L=0}^{n-2}2^{\,n-2-L}g(n-1-L)\,A_L.$$

Three identities then finish the job, each of them a short computation.

First, a full upper half of a row is known exactly:

$$\sum_{d=\lceil L/2\rceil}^{L}\binom{L}{d}=\frac{2^L+c_L}{2}\qquad(L\ge0),$$

where \(c_L=\binom{L}{L/2}\) for even \(L\) and \(c_L=0\) for odd \(L\); it is the symmetry \(\binom{L}{d}=\binom{L}{L-d}\) applied to the whole row. Second, the part of the row above \(M-1\), which the cap cuts off, is a tail of a row, and a tail of a row is the same thing as a sum over first arrivals:

$$2^{-L}\sum_{a=M}^{L}\binom{L}{a}=\sum_{l=M}^{L}\binom{l-1}{M-1}2^{-l}\qquad(1\le M\le L),$$

because both sides give the probability, for \(L\) equally likely choices, that the \(M\)-th step of one kind has arrived by step \(L\): the left side counts by how many such steps there are, the right side by when the \(M\)-th one came. Third, sums of central binomial coefficients against a power of \(4\) telescope; ^[4] the plainest case is

$$\sum_{j=0}^{M}\binom{2j}{j}4^{-j}=(2M+1)\binom{2M}{M}4^{-M}\qquad(M\ge0),$$

which has a one-line certificate: with \(a_j=\binom{2j}{j}4^{-j}\) the recurrence \(2j\,a_j=(2j-1)a_{j-1}\) gives \((2j+1)a_j-(2j-1)a_{j-1}=a_j\), and summing that from \(1\) to \(M\) leaves exactly the identity. The same device handles a polynomial weight: for \(r(j)=j\,w(j)\) one has \(a_{j+1}r(j+1)-a_jr(j)\)\({}=a_j\left[\tfrac12(2j+1)w(j+1)-j\,w(j)\right]\), whose leading terms cancel, so every polynomial weight is matched by a \(w\) of the same degree and the sum closes as a polynomial in \(M\) times \(\binom{2M}{M}4^{-M}\).

Now put them to work. The first identity gives \(A_L=\tfrac12(2^L+c_L)-T_L\) with \(T_L=\sum_{a=M}^{L}\binom{L}{a}\), a tail that is empty for \(L\lt M\), so \(S_2(n)\) splits into three sums. The halves of the rows give

$$2^n\sum_{L=0}^{n-2}g(n-1-L),$$

a polynomial in \(n\) times \(2^n\). The middle terms \(c_L\) survive only for even \(L=2j\), where \(2^{-L}c_L=\binom{2j}{j}4^{-j}\), and give

$$2^n\sum_{j=0}^{M-1}\binom{2j}{j}4^{-j}g(n-1-2j),$$

in which \(g(n-1-2j)\) is a quadratic in \(j\), so the third identity closes it into a polynomial in \(M\) times \(\binom{2M}{M}\). The tails give the third sum: the second identity turns \(T_L\) into \(2^L\sum_{l=M}^{L}\binom{l-1}{M-1}2^{-l}\), and exchanging the two sums leaves \(\sum_{L=l}^{n-2}g(n-1-L)=P(n-1-l)\) inside, so that part is

$$-2^{\,n+1}\sum_{l=M}^{n-2}\binom{l-1}{M-1}2^{-l}P(n-1-l),\qquad P(X)=\frac{X(X^2+3X+5)}{3}.$$

It closes with a certificate of the same kind. Writing \(b_l=\binom{l-1}{M-1}2^{-l}\), the recurrence \(2(l-M+1)b_{l+1}=l\,b_l\) makes \(r(l)=(l-M)v(l)\) an antidifference for the weight \(\tfrac12l\,v(l+1)-(l-M)v(l)\), which has one degree more than \(v\); a cubic weight is therefore matched up to a single constant, and that constant multiplies \(\sum_lb_l\), which the second identity and then the first turn into the tail of the row \(n-2\) above its middle. The odd sizes are treated in the same way, with the two middle-cell sums of Section 10 added. The rest is finite bookkeeping, and the outcome, with the \(2^n n(n^2+2)/3\) of Section 7 added back, is the total trace \(S(n)\) of all squareful matrices of one size \(n\ge2\):

$$S(n)=\frac{2^n}{3}n(2n^2+7)+\frac{8}{15}\binom{n}{n/2}n(n+1)(n-2)$$

for even \(n\), and

$$S(n)=\frac{2^n}{3}(2n^3+7n+3)+\frac{4}{15}\binom{n-1}{\frac{n-1}{2}}n(4n^2-6n-13)$$

for odd \(n\), while \(S(1)=1\) because size \(1\) has a single matrix instead of four. The first values are \(1,40,216,960,3528,11904,36752\), which agree with a direct search over all squareful matrices up to size \(7\).

12. Adding up the sizes

What the problem asks for is \(T(n)=S(1)+S(2)+\dots+S(n)\). The parts with a power of two add up in the usual way: the polynomial in front of \(2^m\) is \((2m^3+7m)/3\) for even \(m\) and \((2m^3+7m+3)/3\) for odd \(m\), and

$$\sum_{m=1}^{n}2^m p(m)=\frac{2^n}{3}\left(4n^3-12n^2+50n-63-(-1)^n\right)+\frac{64}{3},$$

which can be checked by subtracting two consecutive values. The parts with a binomial coefficient all carry \(\binom{2j}{j}\) with \(j=\lfloor m/2\rfloor\): the size \(m=2j\) brings \(\frac{32}{15}j(j-1)(2j+1)\) and the size \(m=2j+1\) brings \(\frac{4}{15}(2j+1)(16j^2+4j-15)\), so together

$$Q(j)=\frac{4}{15}(2j+1)(24j^2-4j-15).$$

Now comes the only step with a surprise in it. Because \(\binom{2j+2}{j+1}=\frac{2(2j+1)}{j+1}\binom{2j}{j}\), for any polynomial \(s\) the difference of \(\binom{2j}{j}\,j\,s(j)\) between \(j+1\) and \(j\) equals \(\binom{2j}{j}\left[2(2j+1)s(j+1)-j\,s(j)\right]\). With

$$s(j)=\frac{8}{15}(8j^2-24j+21)$$

one has the polynomial identity

$$Q(j)=2(2j+1)s(j+1)-j\,s(j)-\frac{28}{3},$$

which the reader can verify by multiplying out, and therefore

$$\sum_{j=0}^{M}Q(j)\binom{2j}{j}=2(2M+1)s(M+1)\binom{2M}{M}-\frac{28}{3}\sum_{j=0}^{M}\binom{2j}{j}.$$

The constant \(-28/3\) cannot be absorbed. If \(s\) has degree \(e\), then \(2(2j+1)s(j+1)-j\,s(j)\) has degree \(e+1\) with three times the leading coefficient of \(s\), so matching \(Q\), of degree \(3\), forces \(s\) to be quadratic; its three coefficients are then fixed by the three leading coefficients of \(Q\), and the displayed constant is what is left over. The running total of central binomial coefficients it leaves behind is not a hypergeometric term either: Gosper's algorithm, which finds a hypergeometric antidifference whenever one exists, returns none for \(\binom{2j}{j}\), as a symbolic run confirms. ^[7] We therefore keep that finite running total and compute it, which is what the program spends its single pass on. Putting the two parts together, for \(n\ge1\) and with \(M=\lfloor n/2\rfloor\),

$$15\,T(n)=5\cdot2^n\left(4n^3-12n^2+50n-63-(-1)^n\right)+275+\gamma(n)\binom{2M}{M}-140\sum_{j=0}^{M}\binom{2j}{j},$$

$$\gamma(n)=\begin{cases}16n(2n^2-8n+11),&n\equiv1\ (\mathrm{mod}\ 2)\\[2pt]4(n+1)(4n^2-18n+35),&n\equiv0\ (\mathrm{mod}\ 2)\end{cases}$$

The two shapes of \(\gamma\) differ only because for even \(n\) the size \(2M+1\) is not included and its term must be taken back out; the constant \(275\) is \(5\cdot64\) from the sum of powers of two, less \(45\) for the three that the formula for odd sizes over-counts at size \(1\). As a check, \(n=2\) gives \(M=1\), \(\binom21=2\) and \(1+2=3\), so \(15\,T(2)=5\cdot4\cdot20\)\({}+275+180\cdot2-140\cdot3\)\({}=400+275+360-420=615\), and \(T(2)=41\). The identity is meant for \(n\ge1\): at \(n=0\) its right-hand side gives \(-45\), that is \(T(0)=-3\), so the empty range is set to \(T(0)=0\) separately, exactly as the program does.

13. Working modulo 1234567891

The modulus \(p=1234567891\) is prime, which the self-test verifies by trial division, so every non-zero remainder has an inverse and \(15^{-1}=15^{\,p-2}\) by Fermat's little theorem. ^[8] The programs accept \(0\le n\lt p\), and that restriction is what makes the division below safe: every factor of the accumulated denominator is one of \(1,2,\dots,M\) with \(M=\lfloor n/2\rfloor\lt p\), so none of them is \(0\) modulo \(p\) and neither is their product. The power \(2^n\) costs about \(\log n\) multiplications by repeated squaring, and the polynomials cost a handful. Everything therefore hinges on the two binomial quantities at \(M=\lfloor n/2\rfloor\), which is about \(6.17\cdot10^7\) for the size we are asked about: the central binomial coefficient itself and the running total of all central binomial coefficients up to it.

Evaluating \(\binom{2k}{k}=\frac{4k-2}{k}\binom{2k-2}{k-1}\) step by step with a fresh inverse each time would cost one modular power per term. Instead the program keeps three numbers for a block \(a\le k\le b\) of the range. Each block starts from the relative value \(1\), and with

$$b_{a,k}=\prod_{t=a}^{k}\frac{4t-2}{t}$$

it keeps

$$N=\prod_{t=a}^{b}(4t-2),\qquad D=\prod_{t=a}^{b}t,\qquad \Sigma=D\sum_{k=a}^{b}b_{a,k},$$

so that \(N/D=b_{a,b}\) is the block's last relative value and \(\Sigma/D\) its relative sum. For the block that starts at \(a=1\) these are the central binomial coefficient and its sum, because \(b_{1,k}=\binom{2k}{k}\); a later block carries relative values only, and the merge below supplies the missing scale from the blocks in front of it. One step of a block updates the three numbers by

$$N\leftarrow N(4k-2),\qquad D\leftarrow Dk,\qquad \Sigma\leftarrow\Sigma k+N,$$

where \(\Sigma k\) rescales the old total to the new common denominator and the \(N\) that is added is the new one. No division happens inside the loop; a single inverse of \(D\) at the end turns both \(N\) and \(\Sigma\) into the numbers we want, and the term \(\binom{0}{0}=1\) is added to the total separately. Two neighbouring blocks merge by

$$(N,D,\Sigma)\ast(N',D',\Sigma')=(NN',\,DD',\,\Sigma D'+N\Sigma'),$$

because the values of the second block are its own relative values scaled by the last value of the first. The merge is associative but not commutative, so the range \(1\le k\le M\) can be cut into as many blocks as there are processors and each block run on its own thread, as long as the blocks are merged in increasing order; merging the whole range this way starts at \(a=1\) and therefore does give \(\binom{2M}{M}\) and its sum, with the term at \(0\) added separately. The whole computation is \(\Theta(n)\) multiplications, a couple of modular powers, and one triple of numbers per thread.

How the Code Works

The C++ program ^[9] is a direct transcription of Sections 12 and 13. power is the usual square-and-multiply modulo \(p\), and INVERSE_15 is \(15^{\,p-2}\), computed at compile time and checked by a static_assert. The struct Partial holds the triple \((N,D,\Sigma)\) of Section 13, combine is the merge rule, and binomial_worker is the three-line loop that runs one block of the range. central_binomials splits \(1\le k\le M\) into equal blocks, one per thread, starts the blocks after the first with pthread_create, runs the first block on the calling thread, runs any block whose thread could not be created, joins the rest, merges them in order, inverts \(D\) once and returns the pair \(\left(\binom{2M}{M},\sum_{j\le M}\binom{2j}{j}\right)\). It asks for one thread per million values of \(k\), up to the number of processors, so small inputs stay serial. solve evaluates the formula of Section 12: the polynomial in \(2^n\), the two shapes of \(\gamma(n)\) by the parity of \(n\), then one multiplication by INVERSE_15. Every intermediate value stays inside \(64\) bits, and the bound that covers the additions as well is

$$2(p-1)^2\lt2^{63}\lt2^{64}.$$

Both summands inside combine are at most \((p-1)^2\), the worker's numerator product is below \(2(p-1)^2\) because \(k\le(p-1)/2\), its sum update is at most \((p-1)^2/2+(p-1)\), and the polynomial evaluation keeps every subtraction preceded by a multiple of \(p\), so nothing is ever negative. Signed and unsigned \(64\)-bit arithmetic therefore agree, and Java needs no unsigned types.

The rest of the file exists to check the mathematics. Enumeration walks over every squareful matrix of one size directly: it extends a self-avoiding walk cell by cell, keeps the smallest and largest row and column used so far, and whenever the next number is a perfect square \(j^2\) it demands that those bounds span exactly \(j\) rows and \(j\) columns, which is the third condition of the problem. It adds up the numbers it places on the main diagonal, and it also counts the matrices, so the count \(2^{n+1}\) of Section 4 is tested as well. direct_trace_sums is an independent method with no formula in it: for each size and each of the four corners it stores, for every cell, the total of the numbers in that cell over all squareful matrices of that size whose largest number sits at that corner. Growing the size by one, the cells inside the inner square inherit the table of the opposite corner, which is Section 3, and a ring cell at row-and-column distance \(\delta\) from that corner gets \(q(n^2-\delta)\), where \(q\) is the number of inner matrices, which is Section 6; the two placements of the inner square are added. Reading the main diagonal of the four tables and adding the previous total gives \(T(n)\) in \(O(n^3)\) steps. run_tests checks the primality of \(p\), the three given values, agreement with the direct search through size \(6\), agreement with the \(O(n^3)\) method through size \(128\), the binomial pair against a Pascal triangle built row by row for every \(M\le512\), and that splitting the range into \(2\), \(3\) or \(7\) blocks gives the same pair as the serial run. Without arguments the program prints \(T(123456789)\) modulo \(p\); given \(n\), and optionally a number of threads, it prints \(T(n)\).

The Python and Java versions follow the C++ program function by function in their arithmetic, their block split, their merge order and their self-tests; the thread handling differs in the details of each language. Python runs the blocks on its own threads, falls back to the calling thread for a block whose thread will not start, carries a failed worker's exception back to the caller and refuses to merge a block that produced no result, since Thread.join does not pass an exception on by itself. Java submits every block to an ExecutorService, even when there is only one, and Future.get propagates a failure on its own; a check at class initialisation records that every intermediate value fits into a signed long. Python uses the process affinity where the operating system offers it and the system processor count otherwise, so its automatic number of blocks can differ from the C++ one on the same machine. None of this changes the result of a successful run.

Complexity and Verification

The main loop runs \(\lfloor n/2\rfloor\) times and does three multiplications and three divisions with remainder per step, so the time grows like \(n\); besides the loop there are two modular powers in C++, and three in Python and Java, which work out the inverse of \(15\) when they start instead of at compile time; the memory is one triple of numbers per thread. For \(n=123456789\) that is about \(6.17\cdot10^7\) steps spread over the blocks. Measured as the wall-clock time of the whole program, with the default number of threads, on an Apple M4 Max with \(16\) cores, the C++ program built by clang++ -std=c++20 -O3 -pthread takes about \(0.03\) seconds, the Java version on OpenJDK \(17\) about \(0.05\) seconds including the start of the virtual machine, and CPython \(3.11\), where the interpreter dominates and the global lock keeps the threads from helping, about \(10\) seconds. The checking methods inside the programs are far slower by design: the direct search over all squareful matrices visits every self-avoiding walk of the grid and is only used up to size \(6\), and the \(O(n^3)\) aggregation up to size \(128\).

All three programs print the same value and pass their self-tests. Further checks, which are not part of the published solutions, confirm the mathematics rather than the code. A separate search enumerates every squareful matrix up to size \(7\) and finds the counts \(2^{n+1}\) for \(2\le n\le7\), the single matrix of size \(1\), and the totals \(T(1),\dots,T(7)=1,41,257,\)\(1217,4745,16649,53401\), which the formula reproduces; the same search confirms, matrix by matrix, that the trace equals \(\sum(f_{k+1}-f_k)g(k)\) of Section 6, that the two growth patterns of Section 5 are the only ones, and that diagonal growth always gives \(n(n^2+2)/3\). A separate program that only counts peeling orders, as in Sections 9 and 10, reproduces the sum of all traces for every size up to \(25\), including the middle-cell terms at odd sizes, and its cumulative totals, computed as exact whole numbers with no remainders, agree with the closed formula for every size up to \(399\). A third route avoids the first-passage count altogether: it adds up the survivors directly. With \(F_{n,k}\) the total of \(f_k\) over the whole anti-diagonal family,

$$F_{n,k}=2^k\sum_{r=0}^{n-k}\binom{n-k}{r}\max\left(0,\,n-2\max(r,n-k-r)\right),$$

where the sum over \(r\) runs over the first \(n-k\) steps and \(2^k\) is the two corners times the \(2^{k-1}\) free choices after them, and feeding this into the second form of the trace in Section 6 it gives the same two closed formulas. The identities of Sections 11 and 12 were each checked term by term over the ranges in which they are used, including the certificate \(Q(j)=2(2j+1)s(j+1)\)\({}-j\,s(j)-28/3\) and the telescoped sum it implies, and a symbolic run of Gosper's algorithm reproduced the failure quoted in Section 12. Finally, the three programs and an independent exact-arithmetic computation agree modulo \(p\) at \(n=200\), \(1000\) and \(5000\), and the three programs agree at \(n=123456789\).

Footnotes and References

The references below supply the background facts used above. The description of all squareful matrices, the count \(2^{n+1}\) for \(n\ge2\), the two growth patterns, the ring lemma for the trace, the peeling picture and the two closed formulas are derived explicitly in this article.

1. Project Euler 1016 — Squareful Matrix. The official statement: an n × n matrix is squareful if each number from 1 to n² appears exactly once, each pair of consecutive numbers are orthogonally adjacent, and for each k less than n the numbers 1 to k² form a contiguous square submatrix. The 4 × 4 example given there has trace 24, T(n) is the sum of the traces of all squareful matrices whose size does not exceed n × n, T(2) = 41, T(12) = 9063089 and T(123) ≡ 1139320302 mod 1234567891, and the task asks for T(123456789) mod 1234567891.
2. Wikipedia — Trace (linear algebra). The trace of a square matrix is the sum of the elements on its main diagonal, and it is defined only for square matrices.
3. Wikipedia — Hamiltonian path. A Hamiltonian path is a path in a graph that visits every vertex exactly once. Numbering the cells of the matrix so that consecutive numbers share a side is exactly a Hamiltonian path of the grid graph, and inside a ring, whose graph is a path, such a numbering must start at one of the two ends.
4. Wikipedia — Central binomial coefficient. The central binomial coefficients are the numbers C(2n, n) = (2n)!/(n!)², namely 1, 2, 6, 20, 70, 252, …, and their generating function is 1/√(1 − 4x). The page is background for the definition; the finite identity of Section 11 is proved there by its own certificate.
5. Wikipedia — Negative binomial distribution. Also called the Pascal distribution, it counts the failures before a fixed number of successes in independent trials with two outcomes; the number of sequences in which the (d + 1)-st success of one kind falls exactly on trial i is C(i − 1, d).
6. Wikipedia — Pascal's rule. The identity C(n − 1, k) + C(n − 1, k − 1) = C(n, k) for the numbers of Pascal's triangle; it is what turns a partial row into a recurrence in the length of the row.
7. Wikipedia — Gosper's algorithm. Gosper's algorithm decides, for a hypergeometric term a(j), whether it has a hypergeometric antidifference, that is whether a(1) + … + a(n) can be written as S(n) − S(0) for a hypergeometric term S, and finds S when it exists; the page describes the general algorithm and says nothing about this particular sum. Running it on a(j) = C(2j, j) returns no antidifference, a calculation carried out for this article, which is why Section 12 keeps the finite running total instead of a closed form for it.
8. Wikipedia — Fermat's little theorem. For a prime p and an integer a, the number a^p − a is a multiple of p; for a not divisible by p this gives a^(p−1) ≡ 1 and hence the inverse a^(p−2) mod p.
9. C++ — Euler1016.cpp. The linked immutable C++ revision contains central_binomials, solve, Enumeration, direct_trace_sums and the checks run by --self-test. It is the implementation record for this article.