Detailed mathematical approach

Problem Summary

Let \(S(n)\) be the minimum exposed surface area of a connected solid made from \(n\) unit cubes in the cubic lattice. If the cubes were wrapped separately, they would use \(6n\) unit squares of paper, so the saved paper is

$g(n)=6n-S(n).$

The task is to compute

$G(N)=\sum_{n=1}^{N} g(n)\pmod{10^9+7}$

for \(N=10^{16}\). A direct search over all \(n\) or over all polycubes is hopeless, so the solution exploits the fact that optimal shapes stay as balanced as possible and differ from a box only by a thin extra layer on one face.

Mathematical Approach

The entire method comes from turning a three-dimensional minimum-surface problem into a sequence of two-dimensional minimum-perimeter problems.

Step 1: The Extra Layer Becomes a 2D Perimeter Problem

Suppose a box already exists and we attach \(t\) more unit cubes to one of its faces. The contact faces disappear inside the solid, the new cubes contribute top faces, and the only net increase along the boundary comes from the perimeter of the footprint drawn on that face.

Therefore the key two-dimensional quantity is the minimum perimeter \(B(t)\) of a connected polyomino with area \(t\). The optimal footprint is as square as possible, so

$B(0)=0,\qquad B(t)=2\left\lceil 2\sqrt{t}\right\rceil \quad (t\ge 1).$

Equivalently, if \(a=\lfloor\sqrt{t}\rfloor\), then the perimeter jumps only when we cross square or almost-square thresholds:

$B(t)=\begin{cases} 4a, & t=a^2,\\ 4a+2, & a^2<t\le a(a+1),\\ 4a+4, & a(a+1)<t\le (a+1)^2. \end{cases}$

Step 2: Split Space into Cubic Shells

Fix \(m=\lfloor \sqrt[3]{n}\rfloor\). Then \(n\) lies in the shell

$m^3\le n<(m+1)^3.$

Inside this shell the best side lengths stay as equal as possible, so the optimal box grows through the sequence

$m\times m\times m,\qquad m\times m\times (m+1),\qquad m\times (m+1)\times (m+1),\qquad (m+1)\times (m+1)\times (m+1).$

That creates three natural phases:

$m^3\le n\le m^2(m+1),$

$m^2(m+1)<n\le m(m+1)^2,$

$m(m+1)^2<n<(m+1)^3.$

In those three ranges the extra layer lies on a face of size \(m^2\), \(m(m+1)\), and \((m+1)^2\), respectively.

Step 3: Exact Surface Formulas in the Three Phases

Write \(n\) as a base box plus a remainder layer.

If

$n=m^3+t,\qquad 0\le t\le m^2,$

then we start from the cube \(m\times m\times m\), whose surface area is \(6m^2\), so

$S(n)=6m^2+B(t).$

If

$n=m^2(m+1)+t,\qquad 1\le t\le m(m+1),$

then the base box is \(m\times m\times(m+1)\), whose surface area is

$2\bigl(m^2+2m(m+1)\bigr),$

hence

$S(n)=2\bigl(m^2+2m(m+1)\bigr)+B(t).$

If

$n=m(m+1)^2+t,\qquad 1\le t\le (m+1)^2-1,$

then the base box is \(m\times(m+1)\times(m+1)\), whose surface area is

$2\bigl(2m(m+1)+(m+1)^2\bigr),$

so

$S(n)=2\bigl(2m(m+1)+(m+1)^2\bigr)+B(t).$

Step 4: Convert Surface Area into Saved Paper

Since \(g(n)=6n-S(n)\), the three phase formulas become

$g(n)=6n-6m^2-B(n-m^3),\qquad m^3\le n\le m^2(m+1),$

$g(n)=6n-2\bigl(m^2+2m(m+1)\bigr)-B\bigl(n-m^2(m+1)\bigr),\qquad m^2(m+1)<n\le m(m+1)^2,$

$g(n)=6n-2\bigl(2m(m+1)+(m+1)^2\bigr)-B\bigl(n-m(m+1)^2\bigr),\qquad m(m+1)^2<n<(m+1)^3.$

The remaining problem is now purely arithmetic: sum these formulas over intervals instead of handling each \(n\) one by one.

Step 5: Summing the 2D Correction in Constant Time

Define the prefix sum

$Q(r)=\sum_{t=1}^{r} B(t).$

Within the block where \(s=\lceil\sqrt{t}\rceil\), the first \(s-1\) values contribute \(4s-2\) and the next \(s\) values contribute \(4s\). Therefore complete square blocks can be summed exactly.

Let

$s=\left\lceil\sqrt{r}\right\rceil,\qquad k=s-1,\qquad u=r-k^2,$

and split the partial block into

$u_1=\min(u,s-1),\qquad u_2=u-u_1.$

Then

$Q(r)=\frac{8k^3+3k^2+k}{3}+u_1(4s-2)+u_2(4s).$

So any correction sum on a phase interval is obtained by one subtraction of two prefix values.

Worked Example

Take \(n=10\). Since

$2^3=8\le 10\le 2^2\cdot 3=12,$

we are in the first phase with \(m=2\) and \(t=10-8=2\). The two-dimensional correction is

$B(2)=2\left\lceil 2\sqrt{2}\right\rceil=6.$

Therefore

$S(10)=6\cdot 2^2+6=30,$

and the saved paper is

$g(10)=6\cdot 10-30=30.$

A second example from the next phase is \(n=14\). Here

$14=2^2\cdot 3+2,$

so

$S(14)=2(4+12)+B(2)=32+6=38,$

which gives

$g(14)=84-38=46.$

Summing the resulting values from \(n=1\) to \(18\) yields \(G(18)=530\), matching the small checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they compute the exact integer cube root of \(N\), because only shells with \(m^3\le N\) need to be processed. Then for each \(m\) they evaluate the three phase intervals described above.

For each interval, the contribution is written as

$6\sum_{n=L}^{R} n - C(R-L+1) - \sum_{t=A}^{B} B(t),$

where \(C\) is the base surface area of the current box and the last term is handled by the prefix formula \(Q\). The arithmetic progression \(\sum n\) is computed in closed form, the correction term is obtained by a difference of two prefix values, and every operation is reduced modulo \(10^9+7\).

No implementation enumerates polycubes, no implementation performs dynamic programming over volumes, and no implementation stores large tables. The whole speedup comes from recognizing the geometric shell pattern and converting it into constant-time arithmetic per shell phase.

Complexity Analysis

The outer loop runs for

$m=1,2,\dots,\left\lfloor N^{1/3}\right\rfloor.$

Each \(m\) contributes exactly three constant-time phase updates, because both the linear sum and the two-dimensional correction sum have closed forms. Hence the total running time is

$O\!\left(N^{1/3}\right),$

and the auxiliary memory usage is

$O(1).$

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=775
2. Polycube: Wikipedia - Polycube
3. Polyomino: Wikipedia - Polyomino
4. Isoperimetric inequality: Wikipedia - Isoperimetric inequality
5. Floor and ceiling functions: Wikipedia - Floor and ceiling functions