Problem 338: Cutting Rectangular Grid Paper

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

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.

Writing the first divisibility condition as $w=ad$ and $x=a(d-1)$ gives the standard parameterization

$$\bigl(ad,\ b(d-1)\bigr)\longrightarrow \bigl(a(d-1),\ bd\bigr),\qquad d\ge 2,$$

and, after swapping the roles of the two original sides, every valid rearrangement is captured by the same adjacent-factor exchange $d \leftrightarrow d-1$.

This matches the examples from the statement. For instance, $9\times 4=(3\cdot 3)\times(2\cdot 2)$ with $d=3$ gives $6\times 6$, and the same framework also produces $12\times 3$ and $18\times 2$ after choosing the appropriate orientation.

Raw Global Counting

Once the geometry is parameterized by $(a,b,d)$, the global sum becomes a floor-sum. For a fixed $d\ge 2$, we may choose

$$1\le a\le \left\lfloor\frac{N}{d}\right\rfloor,\qquad 1\le b\le \left\lfloor\frac{N}{d-1}\right\rfloor,$$

so the raw number of candidate constructions is

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

This is exactly the quantity computed by sum_floor_products in the code. It counts valid staircase descriptions, but not yet distinct rectangles: the same rectangle can arise from more than one divisor-chain description, and self-congruent boundary cases also need correction.

Pair and Triple Counting Corrections

Two classical summatory functions clean up those overlaps:

$$D(n)=\sum_{k=1}^{n}\left\lfloor\frac{n}{k}\right\rfloor =\#\{(u,v)\in\mathbb N^2:\ uv\le n\},$$

$$T(n)=\#\{(a,b,c)\in\mathbb N^3:\ abc\le n\}.$$

The exact identity implemented by the solver is

$$\boxed{G(N)=P(N)-T(N)+D(N).}$$

Intuitively, $P(N)$ counts all adjacent-factor exchanges, $T(N)$ removes the resulting overcount coming from divisor chains of length three, and $D(N)$ restores the pair-level boundary cases that would otherwise be subtracted once too many. This identity agrees with the given checks $G(10)=55$, $G(10^3)=971745$, and $G(10^5)=9992617687$.

Fast Evaluation of $D(n)$

The divisor summatory function is not computed term by term. If

$$q=\left\lfloor\frac{n}{k}\right\rfloor,$$

then the same quotient persists for all $k$ in the interval

$$k\le j\le r,\qquad r=\left\lfloor\frac{n}{q}\right\rfloor.$$

Hence one block contributes $q(r-k+1)$, and the loop jumps directly from one quotient block to the next. This is the standard harmonic-grouping or quotient-block trick, reducing the cost from linear to about $O(\sqrt n)$ distinct blocks.

Fast Evaluation of $T(n)$

A naive triple count is far too slow. Let

$$L=\left\lfloor n^{1/3}\right\rfloor.$$

Every triple $(a,b,c)$ with $abc\le n$ has at least one coordinate at most $L$. Counting triples by how many coordinates are $\le L$ gives the inclusion-exclusion identity

$$T(n)=3\sum_{a=1}^{L} D\!\left(\left\lfloor\frac{n}{a}\right\rfloor\right) -3\sum_{a=1}^{L}\sum_{b=1}^{L}\left\lfloor\frac{n}{ab}\right\rfloor +L^3.$$

The first term counts triples with at least one small coordinate, the second removes those with at least two small coordinates counted twice, and the final $L^3$ adds back the triples where all three coordinates are small.

How the Code Works

The helper iroot3 computes the exact integer cube root $L=\lfloor n^{1/3}\rfloor$. The function divisor_summatory evaluates $D(n)$ by quotient blocks and memoizes repeated calls. The function sum_floor_products evaluates $P(n)$ by jumping over intervals where both $\lfloor n/i\rfloor$ and $\lfloor n/(i-1)\rfloor$ stay constant. The function triple_count implements the formula for $T(n)$, with sum_double_floor handling the double floor sum. Finally, compute_G returns $P(n)-T(n)+D(n)$, and only the last step reduces modulo $10^8$.

Complexity Analysis

The dominant cost is the double sum up to $L=\lfloor n^{1/3}\rfloor$, so the overall running time is about $O(L^2)=O(n^{2/3})$. The remaining divisor-summatory pieces are sublinear thanks to quotient blocks and caching. Memory usage is about $O(n^{1/3})$ plus cache overhead. This is why $N=10^{12}$ is practical here, whereas direct enumeration over all rectangles would be hopeless.

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());
    }
}