Problem 338: Cutting Rectangular Grid Paper

View on Project Euler

Project Euler Problem 338 Solution

EulerSolve provides an optimized solution for Project Euler Problem 338, Cutting Rectangular Grid Paper, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer rectangle \(w \times h\) with \(0 \lt h \le w \le N\), let \(F(w,h)\) be the number of distinct rectangles that can be obtained by cutting the sheet along grid lines into two pieces and rearranging those two pieces into a new rectangle. Rectangles congruent to the original one are excluded, and \(w \times h\) is identified with \(h \times w\). We must compute $$G(N)=\sum_{1\le h\le w\le N} F(w,h).$$ Since the problem asks for \(N=10^{12}\), neither iterating over all pairs \((w,h)\) nor testing all cuts is remotely feasible. The solution therefore converts the geometry into an arithmetic counting problem. Mathematical Approach Arithmetic Reformulation of One Cut For a fixed rectangle \(w \ge h\), let \(x\) be the smaller side of a candidate new rectangle. Area is preserved, so \(x\mid wh\) and \(x\le\sqrt{wh}\). A useful equivalent criterion is $$F(w,h)=\#\Bigl\{x:\ x\mid wh,\ x\le\sqrt{wh},\ x\ne h,\ (w-x)\mid w\ \text{or}\ (x-h)\mid x\Bigr\}.$$ This already shows that the problem is really about divisibility, not about simulating cuts square by square....

Detailed mathematical approach

Problem Summary

A sheet of grid paper has integer dimensions \(w \times h\) and grid spacing \(1\). We cut it along the grid lines into two pieces and rearrange those two pieces, without overlap, into a rectangle of different dimensions. Then \(F(w,h)\) is the number of distinct rectangles that can be made this way: a rectangle congruent to the sheet does not count, and \(x \times y\) is the same rectangle as \(y \times x\). The statement gives \(F(2,1)=0\), \(F(2,2)=1\), \(F(9,4)=3\) — from the \(9 \times 4\) sheet one can make \(18 \times 2\), \(12 \times 3\) and \(6 \times 6\) — and \(F(9,8)=2\). We must evaluate

$$G(N)=\sum_{0 \lt h \le w \le N} F(w,h)$$

for \(N=10^{12}\) and report the result modulo \(10^8\). The checkpoints \(G(10)=55\), \(G(10^3)=971745\) and \(G(10^5)=9992617687\) are given. [1]

The whole solution in one sentence: the rectangles a sheet can be made into are exactly those a staircase of \(k\) steps produces, turning sides \((k+1)a\) and \(kb\) into sides \(ka\) and \((k+1)b\), so we count those staircase moves, remove the ones that return the sheet and the ones that repeat a result, and evaluate the three resulting divisor sums in about \(N^{2/3}\) operations.

Mathematical Approach

1. What the area already tells us

Cutting and rearranging preserves area, so every rectangle made from a \(w \times h\) sheet has sides \(x,y\) with \(xy=wh\). For the smallest sheets that settles everything: the area \(2\) factors only as \(2 \times 1\), which is the sheet again, so \(F(2,1)=0\); the sheet \(2 \times 2\), cut into two \(2 \times 1\) pieces laid end to end, gives \(1 \times 4\), the only other factorisation of \(4\), so \(F(2,2)=1\); for \(3 \times 3\) only \(9 \times 1\) is available and no two-piece cut reaches it, so \(F(3,3)=0\). The real question is which factorisations of \(wh\) are reachable.

2. The staircase cut

Fix integers \(k \ge 1\) and \(c,d \ge 1\) and take the sheet of width \((k+1)c\) and height \(kd\). Divide it into a \((k+1) \times k\) array of \(c \times d\) blocks, block \((i,j)\) covering the columns \(ic+1,\dots,(i+1)c\) and the rows \(jd+1,\dots,(j+1)d\) for \(0 \le i \le k\) and \(0 \le j \le k-1\). Let piece \(A\) consist of the blocks with \(i+j \le k-1\) and piece \(B\) of those with \(i+j \ge k\). Their boundary is a staircase of \(k\) steps, each \(c\) wide and \(d\) high, from the top edge of the sheet to the bottom edge; both pieces are connected, and each holds \(k(k+1)/2\) blocks.

Now translate \(B\) by \(c\) to the left and \(d\) upward. In block coordinates that is the map \((i,j)\mapsto(i-1,j+1)\), which sends the blocks of \(B\) to

$$\{(i,j):i+j\ge k,\ i\le k-1,\ j\le k\},$$

while \(A\) stays where it is, with \(i \le k-1\) and \(j \le k-1\) for all of its blocks. The two images are disjoint and together form the whole \(k \times (k+1)\) array of blocks, so the pieces now tile the rectangle of width \(kc\) and height \((k+1)d\). Every staircase therefore realises

$$\bigl((k+1)c\bigr)\times(kd)\ \to\ (kc)\times\bigl((k+1)d\bigr).$$

The example of the statement is \(k=2\), \(c=3\), \(d=2\). Piece \(A\) is the three leftmost columns at full height together with the bottom two rows of the next three columns; the cut leaves the top edge at \(x=3\), drops \(2\), runs \(3\) to the right and drops \(2\) to the bottom edge, and moving \(B\) three left and two up produces the \(6 \times 6\) square. Reading the same sheet as \(4\) wide and \(9\) high gives the other two rectangles of the statement: \(k=1\), \(c=2\), \(d=9\) is a straight cut into two \(2 \times 9\) halves stacked into \(2 \times 18\), and \(k=3\), \(c=1\), \(d=3\) gives \(3 \times 12\).

The construction is old: Pacioli described around 1500 how to cut a \(4 \times 24\) rectangle into two pieces forming \(3 \times 32\), which is \(k=3\), \(c=1\), \(d=8\), Leonardo da Vinci drew the matching zigzag, and a Japanese book of 1721 illustrates the dissection of \(3a \times 4b\) into \(4a \times 3b\), again \(k=3\). [2]

3. Why nothing else works

Equal area is not enough: the Wallace–Bolyai–Gerwien theorem guarantees a dissection into finitely many pieces, not into two. [3] We prove that the dimensions reachable from a sheet are exactly those of Section 2. Let the two rectangles have sides \(w \ge h\) and \(W \ge H\) with \(WH=wh\); rearranging is reversible, so take the second to be the more elongated one, \(W \gt w \ge h \gt H\). As in the figure of the statement, the pieces are rearranged on the grid: each placement is a quarter turn, possibly a reflection, and a translation by whole cells.

A piece whose bounding box has sides \(p \ge q\) fits into a rectangle with sides \(a \ge b\) only if \(p \le a\) and \(q \le b\), since the turned placement would need the stronger \(p \le b\); as each piece lies in both rectangles, \(p \le w\) and \(q \le H\). A piece is connected, so the columns it meets in the new rectangle form an interval of at most \(p \le w \lt W\) columns; if both pieces met the left edge, neither could reach the right edge. Hence one piece contains the whole leftmost column and the other the whole rightmost column, and both boxes are \(H\) high there.

In the sheet each box therefore has height \(H\) or width \(H\). If one has height \(H\), a row of the sheet outside it is covered by the other piece alone, which gives that piece the full width \(w\); as \(w \gt H\) its box is \(w \times H\), and a row outside that box gives the first one the width \(w\) too. If instead one box has width \(H\), the same argument with columns makes both boxes \(H\) wide and \(h\) high, and turning the sheet a quarter turn — the same sheet — returns us to the first case. So the sheet can be read as \(s \times t\) with \(\{s,t\}=\{w,h\}\) and both boxes measuring \(s \times H\); for the move taking \(9 \times 8\) to \(6 \times 12\) it is \(s=8\), not \(9\). As \(s \ge h \gt H\), no piece turns a quarter turn between that reading and the new rectangle.

Put \(\delta=W-s\) and \(\varepsilon=t-H\), both positive. Two boxes of height \(H\) cover the \(t\) rows only if one takes the bottom \(H\) and the other the top \(H\), so each piece contains a full-width band of the \(\varepsilon\) rows the other box misses; likewise one box takes the left \(s\) columns of the new rectangle and the other the right \(s\), each piece holding a full-height band of \(\delta\) columns. Covering forces \(\delta \le s\) and \(\varepsilon \le H\), and \(s(H+\varepsilon)=(s+\delta)H\) gives \(s\varepsilon=H\delta\). If \(\delta=s\) then \(\varepsilon=H\): the pieces are the two halves \(s \times H\) and the cut is straight, \(k=1\).

Otherwise \(\delta \lt s\) and \(\varepsilon \lt H\), so the boxes share \(s-\delta\) columns of the new rectangle and \(H-\varepsilon\) rows of the sheet. Reflected differently left to right, the two \(\delta\)-bands would fall in the same columns of the sheet, where one piece fills every row of the bottom box and the other every row of the top box, and those boxes share \(H-\varepsilon \gt 0\) rows: impossible. Reflected differently top to bottom, the two \(\varepsilon\)-bands would fall in the same rows of the new rectangle, where the pieces would overlap in the shared columns. Both pieces therefore undergo the same reflection, which we absorb into the orientation of the new rectangle: the placements become translations, by \(0\) and by \(\delta\) columns.

Number the columns from \(0\) and group them by their remainder modulo \(\delta\); a translation by \(0\) or \(\delta\) columns keeps every cell in its group, so each group holds as many cells in the sheet as in the new rectangle. If \(c_r\) of the \(s\) columns of the sheet lie in group \(r\), then \(c_r+1\) of the \(W=s+\delta\) columns of the new rectangle do, and counting that group's cells in both rectangles gives

$$t\,c_r=H(c_r+1),\qquad \varepsilon\,c_r=H.$$

so every \(c_r\) is the same integer \(k=H/\varepsilon \ge 1\). Summing over the \(\delta\) groups gives \(s=k\delta\), and with it \(H=k\varepsilon\), \(t=(k+1)\varepsilon\) and \(W=(k+1)\delta\). Read with its other side first, the sheet is \(\bigl((k+1)\varepsilon\bigr)\times(k\delta)\) and the new rectangle \((k\varepsilon)\times\bigl((k+1)\delta\bigr)\): the staircase of Section 2 with \(c=\varepsilon\), \(d=\delta\). Reversing a staircase exchanges \(c\) and \(d\), so taking the more elongated rectangle second costs nothing, and the classification is complete.

4. Moves

Call a move a triple \((k,a,b)\) of positive integers. It describes the sheet with sides \((k+1)a\) and \(kb\) and the rectangle with sides \(ka\) and \((k+1)b\): the first side shrinks by the factor \(k/(k+1)\), the second grows by \((k+1)/k\). By Sections 2 and 3, the rectangles counted by \(F(w,h)\) are exactly the results of the moves whose sheet is \(\{w,h\}\) in either of its two readings, with the result congruent to the sheet excluded.

Once the sheet and the result are read in a fixed order, the move is determined. Let the sheet be read as \((A,B)\) and the result as \((X,Y)\) with \(X \lt A\). Equal areas give \(AB=XY\) and hence \(A/X=Y/B\); write \(A/X=p/q\) in lowest terms. A move needs \(A/X=(k+1)/k\), a ratio of consecutive integers, so \(p=q+1\) is forced, and with it \(k=q\), \(a=A/p\) and \(b=B/q\). A sheet–result pair therefore has at most one move per reading, and at most two of the four readings can have the first side shrink and the second grow.

5. Worked example: the sheets 9×4 and 9×8

For the sheet \(9 \times 4\) the reading \((9,4)\) needs \(k+1 \mid 9\) and \(k \mid 4\), leaving only \(k=2\) with \(a=3\), \(b=2\) and the result \(6 \times 6\); the reading \((4,9)\) needs \(k+1 \mid 4\) and \(k \mid 9\), leaving \(k=1\) with \(2 \times 18\) and \(k=3\) with \(3 \times 12\). Three moves, three different results, so \(F(9,4)=3\).

For the sheet \(9 \times 8\) the reading \((9,8)\) allows \(k=2\), with \(a=3\), \(b=4\) and the result \(6 \times 12\), and \(k=8\), with \(a=b=1\) and the result \(8 \times 9\) — the sheet itself, turned. The reading \((8,9)\) allows \(k=1\), giving \(4 \times 18\), and \(k=3\), with \(a=2\), \(b=3\), giving \(6 \times 12\) again. Four moves, but one gives the sheet back and two describe the same pair, so \(F(9,8)=2\). The two corrections counted below are exactly these two phenomena.

6. Counting all moves

A sheet is admissible when both of its sides are at most \(N\). For a fixed \(k\) the move \((k,a,b)\) is admissible precisely when \((k+1)a \le N\) and \(kb \le N\), and these two conditions are independent, so with \(d=k+1\) the number of admissible moves is

$$P(N)=\sum_{d=2}^{N}\Bigl\lfloor\frac{N}{d}\Bigr\rfloor\Bigl\lfloor\frac{N}{d-1}\Bigr\rfloor.$$

Only the sheet is bounded; the result may stick out beyond \(N\), as the sheet \(9 \times 4\) contributes \(18 \times 2\) to \(G(10)\) even though \(18 \gt 10\). We computed \(P(10)=81\), \(P(10^3)=994101\) and \(P(10^5)=9998969787\).

7. The moves that give the sheet back

The move \((k,a,b)\) returns a rectangle congruent to its sheet exactly when \(a=b\). Since \(ka \ne (k+1)a\), the unordered pairs \(\{ka,(k+1)b\}\) and \(\{(k+1)a,kb\}\) can only agree through \(ka=kb\) together with \((k+1)b=(k+1)a\), that is \(a=b\); and conversely \(a=b\) makes the result the sheet read backwards. The move \(k=8\), \(a=b=1\) of Section 5 is such a case: eight unit steps carry the \(9 \times 8\) sheet into the \(8 \times 9\) rectangle. Each pair \((k,a)\) with \((k+1)a \le N\) gives one such move, and \(ka \le N\) then holds automatically, so their number is

$$\begin{aligned}\sum_{d\ge2}\Bigl\lfloor\frac{N}{d}\Bigr\rfloor&=D(N)-N,\\ D(n)&=\sum_{k=1}^{n}\Bigl\lfloor\frac{n}{k}\Bigr\rfloor,\end{aligned}$$

where \(D(n)\) counts the pairs \((u,v)\) of positive integers with \(uv \le n\), the lattice points under a hyperbola. [5] For \(N=10\) this removes \(27-10=17\) of the \(81\) moves.

8. The pairs that two moves describe

Suppose the move \((k,a,b)\) has sheet \((A,B)=((k+1)a,kb)\) and result \((X,Y)=(ka,(k+1)b)\), and a second move describes the same sheet and result. By Section 4 it must pair the sides the other way round, matching \(A\) with \(Y\) and \(B\) with \(X\). If \(a \gt b\), that reading does have \(A \gt Y\) and \(B \lt X\), and the second move exists exactly when \(A/Y=a/b\) is a ratio of consecutive integers; if \(a \lt b\) the same holds for \(b/a\). With \(g=\gcd(a,b)\) the condition is \(|a-b|=g\), and then \(\min(a,b)=k'g\) and \(\max(a,b)=(k'+1)g\) for some \(k' \ge 1\).

A doubly described pair is thus given by two step counts \(k \ne k'\) and a scale \(g\), and consists of the rectangles

$$\begin{aligned}R_1&=(k+1)(k'+1)g\times kk'g,\\ R_2&=(k+1)k'g\times k(k'+1)g,\end{aligned}$$

described by the four moves \((k,(k'+1)g,k'g)\), \((k',(k+1)g,kg)\), \((k,k'g,(k'+1)g)\) and \((k',kg,(k+1)g)\): the first two carry \(R_1\) to \(R_2\), the last two carry \(R_2\) to \(R_1\). For \(k=1\), \(k'=2\), \(g=1\) this is the pair \(6 \times 2\) and \(4 \times 3\), reached from the \(6 \times 2\) sheet by \((1,3,2)\) and \((2,2,1)\); for \(k=3\), \(k'=2\), \(g=1\) it is the pair \(12 \times 6\) and \(8 \times 9\) of Section 5. The requirement \(k \ne k'\) matters: for \(k=k'\) one of the two rectangles is a square and the two readings collapse into a single move, as for the move \((2,3,2)\) of the \(9 \times 4\) sheet, which has \(|a-b|=1=\gcd(3,2)\) with \(k'=2=k\) and whose result \(6 \times 6\) is counted once.

Now count these pairs with \(k \gt k'\), so that each is met once. The largest side of \(R_1\) is \((k+1)(k'+1)g\), the largest side of \(R_2\) is \(k(k'+1)g\), and the pair must be subtracted once for each of its two rectangles that is an admissible sheet. Substituting \((u,v)=(k,k'+1)\) for \(R_2\) and \((u,v)=(k'+1,k+1)\) for \(R_1\) turns both conditions into \(uvg \le N\), and the first substitution produces every pair with \(u \ge v \ge 2\) exactly once, the second every pair with \(2 \le u \lt v\) exactly once. The number of sheet–result pairs with two descriptions is therefore

$$M(N)=\#\{(u,v,g):u,v\ge2,\ uvg\le N\}.$$

Let \(T(n)\) count all triples with \(uvg \le n\). The triples with \(u=1\) and those with \(v=1\) number \(D(N)\) each, and the \(N\) triples with \(u=v=1\) are then removed twice, so inclusion–exclusion [7] gives

$$M(N)=T(N)-2D(N)+N.$$

For \(N=10\) that is \(53-54+10=9\).

9. The identity

Every admissible move is counted once in \(P(N)\). Discarding the moves that return the sheet and the second description of each doubly described pair leaves one count per sheet and distinct result, which is what \(G(N)\) sums:

$$\begin{aligned}G(N)=P(N)&-\bigl(D(N)-N\bigr)\\&-\bigl(T(N)-2D(N)+N\bigr),\end{aligned}$$

and the two terms \(N\) cancel, leaving the identity the programs evaluate,

$$G(N)=P(N)-T(N)+D(N).$$

At \(N=10\) the three steps read \(81-17-9=55\), and the closed form \(81-53+27=55\). At \(N=10^3\) it gives \(994101-29425+7069=971745\) and at \(N=10^5\) it gives \(9998969787-7518850+1166750\)\({}=9992617687\), both as stated. Note that \(T(n)\) is the summatory function of the three-fold divisor function \(d_3\), which counts the ordered factorisations \(n=rst\). [6]

10. Evaluating D(n) by quotient blocks

The quotient \(\lfloor n/k\rfloor\) takes only about \(2\sqrt n\) distinct values, and the indices sharing one value \(q=\lfloor n/k\rfloor\) form the block \(k \le j \le r\) with \(r=\lfloor n/q\rfloor\). Adding \(q(r-k+1)\) and jumping to \(r+1\) evaluates \(D(n)\) in \(O(\sqrt n)\) steps; the classical hyperbola form \(D(n)=2\sum_{k\le u}\lfloor n/k\rfloor-u^2\) with \(u=\lfloor\sqrt n\rfloor\) costs the same. [4] For \(n=10^{12}\) the block loop takes \(1999999\) steps.

11. Evaluating T(n)

Let \(L=\lfloor n^{1/3}\rfloor\). Three coordinates all exceeding \(L\) would give a product of at least \((L+1)^3 \gt n\), so every triple with \(uvg \le n\) has a coordinate at most \(L\). Let \(S_1,S_2,S_3\) be the sets of triples whose first, second or third coordinate is at most \(L\). Then

$$\begin{aligned}|S_1|&=\sum_{a=1}^{L}D\Bigl(\Bigl\lfloor\frac{n}{a}\Bigr\rfloor\Bigr),\\ |S_1\cap S_2|&=\sum_{a=1}^{L}\sum_{b=1}^{L}\Bigl\lfloor\frac{n}{ab}\Bigr\rfloor,\end{aligned}$$

and \(|S_1\cap S_2\cap S_3|=L^3\), because \(abc \le L^3 \le n\) for all \(a,b,c \le L\). The three single sets have equal sizes, as do the three pairwise intersections, so inclusion–exclusion gives

$$\begin{aligned}T(n)=\ &3\sum_{a=1}^{L}D\Bigl(\Bigl\lfloor\frac{n}{a}\Bigr\rfloor\Bigr)\\ &-3\sum_{a=1}^{L}\sum_{b=1}^{L}\Bigl\lfloor\frac{n}{ab}\Bigr\rfloor+L^3.\end{aligned}$$

For \(n=10^{12}\) the cube root is exact: \(L=10^4\) and \(L^3=n\).

12. Evaluating P(N)

In \(P(N)\) two quotients move at once. Both \(\lfloor N/i\rfloor\) and \(\lfloor N/(i-1)\rfloor\) are constant on the intersection of their blocks, so the loop multiplies the product of the two quotients by the length of that intersection and jumps to the end of whichever block finishes first. The number of iterations stays \(O(\sqrt N)\): \(86\) for \(N=10^3\) and \(2828424\) for \(N=10^{12}\).

How the Code Works

The C++ program [10] evaluates the three counts of Sections 6 to 8 by the methods of Sections 10 to 12. iroot3 computes \(L=\lfloor n^{1/3}\rfloor\) from a cbrt estimate and corrects it upward and downward with exact cubes held in 128-bit integers, so the floating-point value cannot leave an off-by-one. divisor_summatory is the quotient-block loop of Section 10 behind a memo table (std::unordered_map); inside triple_count the arguments \(\lfloor n/a\rfloor\) for \(a \le L\) are pairwise distinct, so that table is a safety net rather than a saving.

sum_floor_products is Section 12: for the current \(i\) it reads \(a=\lfloor n/i\rfloor\) and \(b=\lfloor n/(i-1)\rfloor\), takes the first indices after the current blocks as n / a + 1 and n / b + 2 — the shift by two accounts for the argument \(i-1\) — clamps the jump at \(n+1\) and adds cnt * a * b. sum_double_floor evaluates the double sum of Section 11, splitting \(1 \le a \le L\) into consecutive blocks whose partial sums are added at the end, and triple_count assembles \(3\,|S_1|-3\,|S_1\cap S_2|+L^3\). compute_G returns \(P(n)-T(n)+D(n)\) as an exact integer: the accumulator widths hold every intermediate value, and the cast of cnt to the wide type widens the product before the multiplication. Only mod_norm reduces modulo \(10^8\), shifting a negative residue back into range.

main runs compute_G on \(10\), \(10^3\) and \(10^5\), compares the residues with \(55\), \(971745\) and \(9992617687\) modulo \(10^8\) and stops with a validation message if any differs, before computing the value for \(N=10^{12}\). The Python and Java ports repeat the same six functions, Java in camel case, but compute the target directly, without those checks.

The integer widths differ because the terms of \(P(N)\) sum to roughly \(N^2\), about \(10^{24}\) at the target, far beyond the \(9.2 \cdot 10^{18}\) of a signed 64-bit integer. C++ accumulates in __int128_t, a GCC and Clang extension. [8] Java uses BigInteger for \(P\), the double sum, \(T\) and \(G\), but keeps the loop variables, the denominators, the cube-root correction and divisorSummatory with its Map cache in long, which is safe because \(D(n)\le n(1+\ln n)\) stays far below the 64-bit limit at the target. The double sum would fit a long too, so its \(10^8\) BigInteger additions are the price of uniformity. Python relies on its unlimited integers. [9] One detail of the Java port deserves notice before it is reused: while collecting the partial sums of the double sum it catches and discards every exception, so a failure in one part would silently leave the sum too small. All three ports keep \(ab \le L^2=10^8\), so the inner divisions are safe.

Complexity and Verification

Two parts dominate. The double sum of Section 11 has \(L^2=\lfloor n^{1/3}\rfloor^2 \le n^{2/3}\) terms, which is \(10^8\) at \(n=10^{12}\). The divisor-summatory evaluations cost \(O(\sqrt{n/a})\) each, and

$$\sum_{a=1}^{L}\sqrt{\frac{n}{a}}\le2\sqrt{nL}\le2n^{2/3},$$

so they are of the same order. Counted exactly at \(n=10^{12}\): the triple count asks for \(10^4\) distinct divisor-summatory arguments and runs \(397074327\) quotient blocks, the double sum runs \(10^8\) terms, and \(P\) and \(D(n)\) add \(2828424\) and \(1999999\) blocks — \(501902750\) loop iterations, each quotient block costing two integer divisions and each product block four. The method is thus \(O(n^{2/3})\) in arithmetic operations. Compiled with optimisation the C++ program took about \(2\) seconds on the machine used here, the Java port about the same and the Python port about \(31\) seconds; such times depend on the machine, not on the method. Memory is \(O(n^{1/3})\): the memo table holds one entry per distinct \(\lfloor n/a\rfloor\) with \(a \le L\), at most \(10^4\), and nothing else grows with \(n\).

The three checkpoints built into the C++ program pass. Beyond them, the following was verified independently for this article.

  • A separate program summed \(F(w,h)\) over all \(500500\) sheets with \(w \le 1000\), listing for each sheet the distinct results of its moves; the total is \(971745\), as stated.
  • The identity was checked against a complete enumeration of moves for every \(N \le 150\): the admissible moves number \(P(N)\), those with \(a=b\) number \(D(N)-N\), the doubly described pairs number \(T(N)-2D(N)+N\), no pair is ever described more than twice, and the difference is \(G(N)\).
  • The staircase of Section 2 was built cell by cell for all \(54\) triples with \(k \le 6\) and \(c,d \le 3\): both pieces connected, equal in area, and the translated union exactly the target rectangle.
  • The classification of Section 3 was tested on all \(62\) rectangles of area at most \(32\): a program enumerated every partition of the cells into two edge-connected pieces — at least as permissive as a single cut — identified the pieces up to the eight symmetries of the square, and compared the shape pairs with those of every rectangle of equal area. The reachable rectangles were exactly the results of the moves, the largest count being \(F(12,2)=3\).
  • The block loops agree with direct summation, \(P(n)\) for \(n \le 200\) and the inclusion–exclusion form of \(T(n)\) for \(n \le 300\); \(D(10)=27\) and \(T(10)=53\) match the published summatory functions. [5] [6]

Footnotes and References

The references below supply the background facts used above; the staircase construction, the argument of Section 3, the move calculus and the three counts of the identity are derived in the article itself.

  1. Project Euler 338 — Cutting Rectangular Grid Paper. The official statement defines F(w, h) for a w × h sheet cut along grid lines into two pieces and rearranged, excluding rectangles congruent to the sheet and identifying w × h with h × w, and G(N) as the sum over 0 < h ≤ w ≤ N. It gives F(2,1) = 0, F(2,2) = 1, F(9,4) = 3 with the rearrangements 18 × 2, 12 × 3, 6 × 6, and F(9,8) = 2, the checkpoints G(10) = 55, G(10³) = 971745, G(10⁵) = 9992617687, and asks for G(10¹²) mod 10⁸.
  2. G. N. Frederickson — updates to chapter 7 of Dissections: Plane & Fancy. The page records the history of the step dissection: Pacioli described around 1500, in De Viribus Quantitatis, how to cut a 4 × 24 rectangle into two pieces forming 3 × 32; Leonardo da Vinci outlined the matching zigzag in the Codex Atlanticus; and Kan Chu Sen's Wakoku Chiyekurabe of 1721 illustrates the dissection of a (3a × 4b)-rectangle into a (4a × 3b)-rectangle — all cases k = 3 of Section 2.
  3. Wikipedia — Wallace–Bolyai–Gerwien theorem. Two polygons can be cut into finitely many pieces that reassemble into one another, by translations and rotations, exactly when their areas are equal; the article adds that the standard construction uses far more pieces than the minimum, so equal area gives no two-piece dissection.
  4. Wikipedia — Divisor summatory function. The article describes D(x) as the number of lattice points under the hyperbola jk = x, gives D(x) = 2·Σ(k ≤ u) ⌊x/k⌋ − u² with u = ⌊√x⌋ and states that this allows D(x) to be computed in O(√x) time; its Piltz divisor problem section treats the same sums for d_k, of which T(n) is the case 3.
  5. OEIS A006218 — Σ ⌊n/k⌋. The entry gives a(n) = Σ(k = 1..n) ⌊n/k⌋, equal to the number of solutions of x·y = z with x, y, z ≤ n, and lists a(10) = 27, the value of D(10) used in Sections 7 and 9.
  6. OEIS A061201 — partial sums of the three-fold divisor function. The entry sums d₃, the number of ordered factorisations n = r·s·t, and states that a(n) counts the ordered triples with a·b·c ≤ n. Its terms begin 1, 4, 7, 13, 16, 25, 28, 38, 44, 53, so a(10) = 53 = T(10).
  7. Wikipedia — Inclusion–exclusion principle. The article states the alternating formula for the size of a union and, in complementary form, the count of elements having none of several properties: Section 8 uses the latter to remove the triples with a coordinate 1, Section 11 the former over the three sets of triples with a small coordinate.
  8. GCC manual — 128-bit integers. The manual documents __int128 and unsigned __int128 as an extension available on targets with a wide enough integer mode; the C++ solution uses the signed form for the accumulators of Sections 6 and 11.
  9. Python — numeric types states that integers have unlimited precision and that // is the floored quotient, on which the Python port relies. Java — BigInteger documents the arbitrary-precision integers used by the Java port.
  10. C++ — Euler338.cpp. The linked immutable revision contains the functions discussed above and the three checkpoint validations in main; it is the implementation record for this article, and the brute-force comparisons listed above are separate programs.

Problem 338 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <mutex>
#include <thread>
#include <unordered_map>
#include <vector>
#include <cmath>
#include <functional>

using int64 = long long;
using i128 = __int128_t;

namespace {

constexpr int64 MOD = 100000000LL;

int64 iroot3(int64 n) {
    long double approx = cbrt(static_cast<long double>(n));
    int64 x = static_cast<int64>(approx);
    auto cube = [](int64 v) -> i128 {
        return static_cast<i128>(v) * v * v;
    };
    while (cube(x + 1) <= n) ++x;
    while (cube(x) > n) --x;
    return x;
}

int64 divisor_summatory(int64 n, std::unordered_map<int64, int64>& cache) {
    if (n <= 0) return 0;
    auto it = cache.find(n);
    if (it != cache.end()) return it->second;

    int64 res = 0;
    int64 k = 1;
    while (k <= n) {
        int64 q = n / k;
        int64 r = n / q;
        res += q * (r - k + 1);
        k = r + 1;
    }

    cache[n] = res;
    return res;
}

i128 sum_floor_products(int64 n) {
    if (n < 2) return 0;
    i128 sum = 0;
    int64 i = 2;
    while (i <= n) {
        int64 a = n / i;
        int64 b = n / (i - 1);

        int64 next_a = n / a + 1;
        int64 next_b = n / b + 2;
        int64 nxt = std::min(next_a, next_b);
        if (nxt > n + 1) nxt = n + 1;

        int64 cnt = nxt - i;
        sum += static_cast<i128>(cnt) * a * b;
        i = nxt;
    }
    return sum;
}

i128 sum_double_floor(int64 n, int64 limit) {
    // sum_{a=1..limit} sum_{b=1..limit} floor(n/(a*b))
    int thread_count = static_cast<int>(std::thread::hardware_concurrency());
    if (thread_count <= 0) thread_count = 1;
    if (limit < thread_count) thread_count = static_cast<int>(limit);

    std::vector<i128> partial(thread_count, 0);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    for (int t = 0; t < thread_count; ++t) {
        int64 start = 1 + (limit * t) / thread_count;
        int64 end = (limit * (t + 1)) / thread_count;
        workers.emplace_back([&, t, start, end]() {
            i128 local = 0;
            for (int64 a = start; a <= end; ++a) {
                for (int64 b = 1; b <= limit; ++b) {
                    int64 denom = a * b;
                    local += n / denom;
                }
            }
            partial[t] = local;
        });
    }
    for (auto& th : workers) th.join();

    i128 total = 0;
    for (auto v : partial) total += v;
    return total;
}

i128 triple_count(int64 n) {
    // T(n) = #{(a,b,c) >= 1 : a*b*c <= n}
    // Use inclusion-exclusion with L = floor(n^(1/3)).
    int64 L = iroot3(n);

    std::unordered_map<int64, int64> cache;
    cache.reserve(static_cast<size_t>(L * 2 + 100));

    i128 sum_a = 0;
    int64 a = 1;
    while (a <= L) {
        int64 q = n / a;
        int64 r = n / q;
        if (r > L) r = L;
        int64 cnt = r - a + 1;
        sum_a += static_cast<i128>(cnt) * divisor_summatory(q, cache);
        a = r + 1;
    }

    i128 sum_ab = sum_double_floor(n, L);
    i128 l3 = static_cast<i128>(L) * L * L;

    return 3 * sum_a - 3 * sum_ab + l3;
}

i128 compute_G(int64 n) {
    if (n < 2) return 0;

    std::unordered_map<int64, int64> cache;
    cache.reserve(1 << 16);

    i128 s_prod = sum_floor_products(n);
    int64 d_n = divisor_summatory(n, cache);
    i128 t_n = triple_count(n);

    return s_prod - t_n + d_n;
}

int64 mod_norm(i128 v) {
    v %= MOD;
    if (v < 0) v += MOD;
    return static_cast<int64>(v);
}

}  // namespace

int main() {
    const int64 N = 1000000000000LL;

    struct Check { int64 n; int64 expected; };
    const Check checks[] = {
        {10, 55},
        {1000, 971745},
        {100000, 9992617687LL},
    };

    for (const auto& chk : checks) {
        int64 got = mod_norm(compute_G(chk.n));
        if (got != (chk.expected % MOD)) {
            std::cerr << "Validation failure: G(" << chk.n << ") mod 1e8 = "
                      << got << ", expected " << (chk.expected % MOD) << '\n';
            return 1;
        }
    }

    int64 answer = mod_norm(compute_G(N));
    std::cout << answer << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 100_000_000
    N = 10**12

    def iroot3(n):
        x = round(n ** (1/3))
        while (x+1)**3 <= n: x += 1
        while x**3 > n: x -= 1
        return x

    def divisor_summatory(n, cache):
        if n <= 0: return 0
        if n in cache: return cache[n]
        res = 0
        k = 1
        while k <= n:
            q = n // k
            r = n // q
            res += q * (r - k + 1)
            k = r + 1
        cache[n] = res
        return res

    def sum_floor_products(n):
        if n < 2: return 0
        s = 0
        i = 2
        while i <= n:
            a = n // i
            b = n // (i - 1)
            next_a = n // a + 1
            next_b = n // b + 2
            nxt = min(next_a, next_b, n + 1)
            cnt = nxt - i
            s += cnt * a * b
            i = nxt
        return s

    def sum_double_floor(n, limit):
        total = 0
        for a in range(1, limit + 1):
            for b in range(1, limit + 1):
                total += n // (a * b)
        return total

    def triple_count(n):
        L = iroot3(n)
        cache = {}
        sum_a = 0
        a = 1
        while a <= L:
            q = n // a
            r = n // q
            if r > L: r = L
            cnt = r - a + 1
            sum_a += cnt * divisor_summatory(q, cache)
            a = r + 1
        sum_ab = sum_double_floor(n, L)
        return 3 * sum_a - 3 * sum_ab + L**3

    def compute_G(n):
        if n < 2: return 0
        cache = {}
        s_prod = sum_floor_products(n)
        d_n = divisor_summatory(n, cache)
        t_n = triple_count(n)
        return s_prod - t_n + d_n

    answer = compute_G(N) % MOD
    return str(answer)

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

Java

import java.util.*;
import java.util.concurrent.*;
import java.util.concurrent.atomic.*;

public class Euler338 {

    static final long MOD = 100000000L;

    static long iroot3(long n) {
        long approx = (long) Math.cbrt((double) n);
        while ((approx + 1) * (approx + 1) * (approx + 1) <= n && (approx + 1) > 0)
            approx++;
        while (approx * approx * approx > n)
            approx--;
        return approx;
    }

    static long divisorSummatory(long n, Map<Long, Long> cache) {
        if (n <= 0)
            return 0;
        if (cache.containsKey(n))
            return cache.get(n);

        long res = 0;
        long k = 1;
        while (k <= n) {
            long q = n / k;
            long r = n / q;
            res += q * (r - k + 1);
            k = r + 1;
        }

        cache.put(n, res);
        return res;
    }

    static java.math.BigInteger sumFloorProducts(long n) {
        if (n < 2)
            return java.math.BigInteger.ZERO;
        java.math.BigInteger sum = java.math.BigInteger.ZERO;
        long i = 2;
        while (i <= n) {
            long a = n / i;
            long b = n / (i - 1);

            long next_a = n / a + 1;
            long next_b = n / b + 2;
            long nxt = Math.min(next_a, next_b);
            if (nxt > n + 1)
                nxt = n + 1;

            long cnt = nxt - i;
            sum = sum.add(java.math.BigInteger.valueOf(cnt).multiply(java.math.BigInteger.valueOf(a))
                    .multiply(java.math.BigInteger.valueOf(b)));
            i = nxt;
        }
        return sum;
    }

    static java.math.BigInteger sumDoubleFloor(long n, long limit) {
        int threads = Runtime.getRuntime().availableProcessors();
        if (limit < threads)
            threads = (int) limit;

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<java.math.BigInteger>> futures = new ArrayList<>();

        for (int t = 0; t < threads; t++) {
            final long start = 1 + (limit * t) / threads;
            final long end = (limit * (t + 1)) / threads;

            futures.add(executor.submit(() -> {
                java.math.BigInteger local = java.math.BigInteger.ZERO;
                for (long a = start; a <= end; a++) {
                    for (long b = 1; b <= limit; b++) {
                        local = local.add(java.math.BigInteger.valueOf(n / (a * b)));
                    }
                }
                return local;
            }));
        }

        java.math.BigInteger total = java.math.BigInteger.ZERO;
        for (Future<java.math.BigInteger> f : futures) {
            try {
                total = total.add(f.get());
            } catch (Exception e) {
            }
        }
        executor.shutdown();
        return total;
    }

    static java.math.BigInteger tripleCount(long n) {
        long L = iroot3(n);
        Map<Long, Long> cache = new HashMap<>();

        java.math.BigInteger sum_a = java.math.BigInteger.ZERO;
        long a = 1;
        while (a <= L) {
            long q = n / a;
            long r = n / q;
            if (r > L)
                r = L;
            long cnt = r - a + 1;
            sum_a = sum_a.add(java.math.BigInteger.valueOf(cnt)
                    .multiply(java.math.BigInteger.valueOf(divisorSummatory(q, cache))));
            a = r + 1;
        }

        java.math.BigInteger sum_ab = sumDoubleFloor(n, L);
        java.math.BigInteger l3 = java.math.BigInteger.valueOf(L).multiply(java.math.BigInteger.valueOf(L))
                .multiply(java.math.BigInteger.valueOf(L));

        return sum_a.multiply(java.math.BigInteger.valueOf(3))
                .subtract(sum_ab.multiply(java.math.BigInteger.valueOf(3)))
                .add(l3);
    }

    static java.math.BigInteger computeG(long n) {
        if (n < 2)
            return java.math.BigInteger.ZERO;
        Map<Long, Long> cache = new HashMap<>();

        java.math.BigInteger sProd = sumFloorProducts(n);
        long dN = divisorSummatory(n, cache);
        java.math.BigInteger tN = tripleCount(n);

        return sProd.subtract(tN).add(java.math.BigInteger.valueOf(dN));
    }

    public static String solve() {
        long N = 1000000000000L;
        java.math.BigInteger ans = computeG(N).remainder(java.math.BigInteger.valueOf(MOD));
        long val = ans.longValue();
        if (val < 0)
            val += MOD;
        return String.valueOf(val);
    }

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