Problem 775: Saving Paper

Project Euler Problem 775 Solution

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

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$....

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.

$$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).$$

$$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).$$

$$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),$$

$$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,$$

$$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

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

Problem 775 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <cmath>
#include <functional>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr i64 MOD = 1'000'000'007LL;
constexpr i64 INV2 = (MOD + 1) / 2;

u64 isqrt_u64(const u64 n) {
    if (n == 0) {
        return 0;
    }
    long double d = static_cast<long double>(n);
    u64 x = static_cast<u64>(std::sqrt(d));
    while (static_cast<u128>(x + 1) * (x + 1) <= n) {
        ++x;
    }
    while (static_cast<u128>(x) * x > n) {
        --x;
    }
    return x;
}

u64 icbrt_u64(const u64 n) {
    if (n == 0) {
        return 0;
    }
    long double d = static_cast<long double>(n);
    u64 x = static_cast<u64>(std::cbrt(d));
    while (static_cast<u128>(x + 1) * (x + 1) * (x + 1) <= n) {
        ++x;
    }
    while (static_cast<u128>(x) * x * x > n) {
        --x;
    }
    return x;
}

u64 p2_prefix(const u64 r) {
    // Sum_{t=1..r} 2*ceil(2*sqrt(t)); p2_prefix(0)=0.
    if (r == 0) {
        return 0;
    }

    u64 s = isqrt_u64(r);
    if (static_cast<u128>(s) * s < r) {
        ++s;  // s = ceil(sqrt(r))
    }
    const u64 k = s - 1;

    const u128 base = (static_cast<u128>(8) * k * k * k + static_cast<u128>(3) * k * k + k) /
                      static_cast<u128>(3);
    const u64 consumed = k * k;
    const u64 rem = r - consumed;

    const u64 len1 = (rem < (s - 1)) ? rem : (s - 1);
    const u64 len2 = rem - len1;

    const u128 add1 = static_cast<u128>(len1) * (4 * s - 2);
    const u128 add2 = static_cast<u128>(len2) * (4 * s);

    return static_cast<u64>(base + add1 + add2);
}

i64 mod_mul(const i64 a, const i64 b) {
    return static_cast<i64>((static_cast<__int128>(a) * b) % MOD);
}

i64 mod_sum_range(const u64 l, const u64 r) {
    if (l > r) {
        return 0;
    }
    const i64 len = static_cast<i64>((r - l + 1) % MOD);
    const i64 ends = static_cast<i64>((l % MOD + r % MOD) % MOD);
    return mod_mul(mod_mul(len, ends), INV2);
}

i64 add_phase(i64 acc, const u64 L, const u64 R, const u64 capN, const u64 base, const u64 rBase) {
    if (L > capN) {
        return acc;
    }
    const u64 left = L;
    const u64 right = (R < capN) ? R : capN;
    if (left > right) {
        return acc;
    }

    const u64 len_u = right - left + 1;
    const i64 len = static_cast<i64>(len_u % MOD);
    const i64 sum_n = mod_sum_range(left, right);

    const u64 rr = right - rBase;
    const u64 ll = left - rBase;
    const u64 sum_p_u = p2_prefix(rr) - p2_prefix((ll == 0) ? 0 : (ll - 1));
    const i64 sum_p = static_cast<i64>(sum_p_u % MOD);

    i64 cur = mod_mul(6, sum_n);
    cur = (cur - mod_mul(len, static_cast<i64>(base % MOD)) + MOD) % MOD;
    cur = (cur - sum_p + MOD) % MOD;

    acc += cur;
    if (acc >= MOD) {
        acc -= MOD;
    }
    return acc;
}

i64 G_mod(const u64 N) {
    i64 ans = 0;
    const u64 mMax = icbrt_u64(N);

    for (u64 m = 1; m <= mMax; ++m) {
        const u64 m2 = m * m;
        const u64 mp1 = m + 1;

        const u64 v0 = m * m2;
        const u64 v1 = m2 * mp1;
        const u64 v2 = m * mp1 * mp1;
        const u64 v3 = mp1 * mp1 * mp1;

        const u64 s0 = 6 * m2;
        const u64 s1 = 2 * (m2 + 2 * m * mp1);
        const u64 s2 = 2 * (2 * m * mp1 + mp1 * mp1);

        ans = add_phase(ans, v0, v1, N, s0, v0);
        ans = add_phase(ans, v1 + 1, v2, N, s1, v1);
        ans = add_phase(ans, v2 + 1, v3 - 1, N, s2, v2);
    }

    return ans;
}

}  // namespace

int main() {
    assert(G_mod(18) == 530);
    assert(G_mod(1'000'000) == 951'640'919LL);

    std::cout << G_mod(10'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

import math

MOD = 1000000007
INV2 = (MOD + 1) // 2

def p2_prefix(r):
    if r == 0: return 0
    s = math.isqrt(r)
    if s * s < r:
        s += 1
    k = s - 1
    
    base = (8 * k**3 + 3 * k**2 + k) // 3
    consumed = k * k
    rem = r - consumed
    
    len1 = rem if rem < (s - 1) else (s - 1)
    len2 = rem - len1
    
    add1 = len1 * (4 * s - 2)
    add2 = len2 * (4 * s)
    
    return base + add1 + add2

def mod_sum_range(l, r):
    if l > r: return 0
    length = (r - l + 1) % MOD
    ends = (l + r) % MOD
    return (length * ends * INV2) % MOD

def add_phase(acc, L, R, capN, base, rBase):
    if L > capN: return acc
    left = L
    right = R if R < capN else capN
    if left > right: return acc
    
    len_u = right - left + 1
    length = len_u % MOD
    sum_n = mod_sum_range(left, right)
    
    rr = right - rBase
    ll = left - rBase
    sum_p_u = p2_prefix(rr) - p2_prefix(ll - 1 if ll > 0 else 0)
    sum_p = sum_p_u % MOD
    
    cur = (6 * sum_n) % MOD
    cur = (cur - length * (base % MOD) + MOD) % MOD
    cur = (cur - sum_p + MOD) % MOD
    
    return (acc + cur) % MOD

def G_mod(N):
    ans = 0
    mMax = int(math.floor(math.pow(N, 1.0/3.0)))
    while (mMax + 1)**3 <= N: mMax += 1
    while mMax**3 > N: mMax -= 1
    
    for m in range(1, mMax + 1):
        m2 = m * m
        mp1 = m + 1
        
        v0 = m * m2
        v1 = m2 * mp1
        v2 = m * mp1 * mp1
        v3 = mp1 * mp1 * mp1
        
        s0 = 6 * m2
        s1 = 2 * (m2 + 2 * m * mp1)
        s2 = 2 * (2 * m * mp1 + mp1 * mp1)
        
        ans = add_phase(ans, v0, v1, N, s0, v0)
        ans = add_phase(ans, v1 + 1, v2, N, s1, v1)
        ans = add_phase(ans, v2 + 1, v3 - 1, N, s2, v2)
        
    return ans

def solve():
    return str(G_mod(10**16))

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

Java

public class Euler775 {
    static final long MOD = 1000000007L;
    static final long INV2 = (MOD + 1) / 2;

    static long isqrt(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.sqrt((double) n);
        while ((x + 1) * (x + 1) <= n)
            x++;
        while (x * x > n && x > 0)
            x--;
        return x;
    }

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

    static long p2Prefix(long r) {
        if (r == 0)
            return 0;
        long s = isqrt(r);
        if (s * s < r)
            s++;
        long k = s - 1;

        long base = (8L * k * k * k + 3L * k * k + k) / 3L;
        long consumed = k * k;
        long rem = r - consumed;

        long len1 = (rem < (s - 1)) ? rem : (s - 1);
        long len2 = rem - len1;

        long add1 = len1 * (4L * s - 2L);
        long add2 = len2 * (4L * s);

        return base + add1 + add2;
    }

    static long modMul(long a, long b) {
        return (a * b) % MOD;
    }

    static long modSumRange(long l, long r) {
        if (l > r)
            return 0;
        long len = (r - l + 1) % MOD;
        long ends = ((l % MOD) + (r % MOD)) % MOD;
        return modMul(modMul(len, ends), INV2);
    }

    static long addPhase(long acc, long L, long R, long capN, long base, long rBase) {
        if (L > capN)
            return acc;
        long left = L;
        long right = (R < capN) ? R : capN;
        if (left > right)
            return acc;

        long lenU = right - left + 1;
        long len = lenU % MOD;
        long sumN = modSumRange(left, right);

        long rr = right - rBase;
        long ll = left - rBase;
        long sumPU = p2Prefix(rr) - p2Prefix((ll == 0) ? 0 : (ll - 1));
        long sumP = sumPU % MOD;

        long cur = modMul(6, sumN);
        cur = (cur - modMul(len, base % MOD) + MOD) % MOD;
        cur = (cur - sumP + MOD) % MOD;

        acc += cur;
        if (acc >= MOD)
            acc -= MOD;
        return acc;
    }

    static long GMod(long N) {
        long ans = 0;
        long mMax = icbrt(N);

        for (long m = 1; m <= mMax; m++) {
            long m2 = m * m;
            long mp1 = m + 1;

            long v0 = m * m2;
            long v1 = m2 * mp1;
            long v2 = m * mp1 * mp1;
            long v3 = mp1 * mp1 * mp1;

            long s0 = 6 * m2;
            long s1 = 2 * (m2 + 2 * m * mp1);
            long s2 = 2 * (2 * m * mp1 + mp1 * mp1);

            ans = addPhase(ans, v0, v1, N, s0, v0);
            ans = addPhase(ans, v1 + 1, v2, N, s1, v1);
            ans = addPhase(ans, v2 + 1, v3 - 1, N, s2, v2);
        }

        return ans;
    }

    public static String solve() {
        return Long.toString(GMod(10000000000000000L));
    }

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