Problem 935: Rolling Square

Project Euler Problem 935 Solution

EulerSolve provides an optimized solution for Project Euler Problem 935, Rolling Square, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementations do not simulate a square rolling step by step. Instead, they use a number-theoretic parametrization in which each admissible contribution is represented by an integer $q \ge 1$, one of four companion values $k \in \{4q+1,4q+2,4q+3,4q+4\}$, and a second parameter $m$ with $1 \le m \le q$. The key invariant is primitiveness: a contribution is valid exactly when $\gcd(m,k)=1$. The bound in the problem is imposed through the associated size parameter $n$: the families $k=4q+1$, $4q+2$, and $4q+4$ contribute at $n=4q$, while $k=4q+3$ contributes at $n=4q+2$. So the total answer up to $N$ is a cumulative sum over even $n$, with no odd case at all. Mathematical Approach Once the rolling-square geometry has been reduced to this parametrization, the problem becomes a clean coprime-counting question. The whole job is to evaluate the right coprime counts quickly and add them in the right four families....

Detailed mathematical approach

Problem Summary

The implementations do not simulate a square rolling step by step. Instead, they use a number-theoretic parametrization in which each admissible contribution is represented by an integer $q \ge 1$, one of four companion values $k \in \{4q+1,4q+2,4q+3,4q+4\}$, and a second parameter $m$ with $1 \le m \le q$.

The key invariant is primitiveness: a contribution is valid exactly when $\gcd(m,k)=1$. The bound in the problem is imposed through the associated size parameter $n$: the families $k=4q+1$, $4q+2$, and $4q+4$ contribute at $n=4q$, while $k=4q+3$ contributes at $n=4q+2$. So the total answer up to $N$ is a cumulative sum over even $n$, with no odd case at all.

Mathematical Approach

Once the rolling-square geometry has been reduced to this parametrization, the problem becomes a clean coprime-counting question. The whole job is to evaluate the right coprime counts quickly and add them in the right four families.

The Four Admissible Families

Define

$$R(q,k)=\#\{1 \le m \le q : \gcd(m,k)=1\}.$$

The arithmetic structure visible in the implementations is then

$$S(N)=\sum_{q=1}^{\lfloor N/4 \rfloor}\Bigl(R(q,4q+1)+R(q,4q+2)+R(q,4q+4)\Bigr)+\sum_{q=1}^{\lfloor (N-2)/4 \rfloor}R(q,4q+3).$$

Equivalently, if $F(n)$ denotes the contribution attached to one value of $n$, then

$$F(4q)=R(q,4q+1)+R(q,4q+2)+R(q,4q+4), \qquad F(4q+2)=R(q,4q+3),$$

and $F(n)=0$ for odd $n$. That is the first important invariant: only the residue classes $0$ and $2$ modulo 4 can occur.

Primitiveness Becomes a Coprime Count

For fixed $q$ and $k$, the only thing that matters about $m$ is whether it shares a prime factor with $k$. Prime powers do not create new exclusions: if a prime $p$ divides $k$, then every multiple of $p$ must be removed, regardless of whether $p$ occurs once or many times in $k$.

That is why the count depends only on the distinct prime divisors of $k$. Writing

$$\operatorname{rad}(k)=\prod_{p \mid k} p,$$

the forbidden values of $m$ are exactly the multiples of the primes dividing $\operatorname{rad}(k)$.

Inclusion-Exclusion Over the Distinct Prime Factors

Counting the complement of those forbidden multiples gives the exact formula

$$R(q,k)=\sum_{d \mid \operatorname{rad}(k)} \mu(d)\left\lfloor \frac{q}{d}\right\rfloor,$$

where $\mu$ is the Möbius function. This is simply inclusion-exclusion written compactly: start with all $q$ candidates, subtract multiples of each prime divisor of $k$, add back multiples of each product of two distinct prime divisors, and continue alternating.

The implementations do not build a global Möbius table. They factor $k$, list its distinct prime divisors, enumerate all subsets of that list, multiply the chosen primes to obtain $d$, and add or subtract $\left\lfloor q/d \right\rfloor$ according to the subset parity. That direct subset view matches the formula exactly.

Worked Example

Take $q=2$. Then the three $n=4q=8$ branches are $k=9$, $10$, and $12$.

For $k=12$, the distinct prime divisors are $2$ and $3$, so

$$R(2,12)=\left\lfloor\frac{2}{1}\right\rfloor-\left\lfloor\frac{2}{2}\right\rfloor-\left\lfloor\frac{2}{3}\right\rfloor+\left\lfloor\frac{2}{6}\right\rfloor=2-1-0+0=1.$$

Also, $R(2,9)=2$ because both $1$ and $2$ are coprime to $9$, and $R(2,10)=1$ because only $1$ is coprime to $10$. Hence

$$F(8)=R(2,9)+R(2,10)+R(2,12)=2+1+1=4.$$

The $4q+2$ branch gives $n=10$ with $k=11$, so

$$F(10)=R(2,11)=2.$$

At the smallest nontrivial level, $q=1$ gives $F(4)=3$ and $F(6)=1$, so $S(6)=4$, exactly matching the checkpoint embedded in the implementations.

How the Code Works

The C++, Python, and Java implementations all evaluate the same arithmetic decomposition. They differ only in execution strategy, not in the mathematics.

Prime-Factor Preprocessing

A smallest-prime-factor table is precomputed up to $N+4$. This is a linear sieve, so every later factorization of a candidate $k$ can peel off its distinct prime divisors in near-constant amortized time per divisor instead of testing divisibility up to $\sqrt{k}$.

Counting One Branch

For each required value of $k$, the implementation extracts the distinct prime divisors, enumerates all subset products, and evaluates the inclusion-exclusion sum for $R(q,k)$. Because the bound is $m \le q$, each subset contributes exactly one floor term $\left\lfloor q/d \right\rfloor$.

No symbolic simplification is required. The same routine works uniformly for $4q+1$, $4q+2$, $4q+3$, and $4q+4$.

Outer Accumulation

The outer loop runs over $q$ and evaluates the four admissible families, adding a branch only when its associated $n$ is at most the requested bound. The native implementations split the $q$-interval into chunks, let each worker accumulate an exact local sum, and then reduce the partial sums at the end. The Python entry point delegates to the same compiled arithmetic, so it produces the same count rather than a different algorithm.

Complexity Analysis

Let $K=N+4$ and $Q=\lfloor K/4 \rfloor$. The sieve phase is $O(K)$ time and $O(K)$ memory. That memory cost is the dominant storage use, because the smallest-prime-factor table has one entry per integer up to the limit.

The accumulation phase evaluates up to four candidate values of $k$ for each $q$. A single coprime count costs $O(2^{\omega(k)})$, where $\omega(k)$ is the number of distinct prime divisors of $k$. For the target limit $N=10^8$, every $k \le 100000004$ has at most 8 distinct prime divisors, so each inclusion-exclusion pass uses at most $2^8=256$ subset terms. In practice this makes the arithmetic phase very close to linear in the number of processed branches, and parallel chunking improves wall-clock time without changing the exact count.

Footnotes and References

Problem page: Project Euler 935
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Möbius function: Wikipedia - Möbius function
Coprime integers: Wikipedia - Coprime integers
Sieve of Eratosthenes and prime preprocessing: Wikipedia - Sieve of Eratosthenes

Problem 935 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

using u32 = std::uint32_t;
using u64 = std::uint64_t;

static std::vector<u32> build_spf(u32 n) {
    std::vector<u32> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(n / 10));
    if (n >= 1) spf[1] = 1;
    for (u32 i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (u32 p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > n) break;
            spf[static_cast<std::size_t>(v)] = p;
            if (p == spf[static_cast<std::size_t>(i)]) break;
        }
    }
    return spf;
}

static int factor_distinct(u32 x, const std::vector<u32>& spf, u32 out[10]) {
    int cnt = 0;
    while (x > 1) {
        const u32 p = spf[static_cast<std::size_t>(x)];
        out[cnt++] = p;
        do {
            x /= p;
        } while (x > 1 && x % p == 0);
    }
    return cnt;
}

static u64 count_relprimes_pie(u32 pmax, u32 k, const std::vector<u32>& spf) {
    u32 fac[10];
    const int m = factor_distinct(k, spf, fac);
    const int lim = 1 << m;
    u64 res = 0;
    for (int mask = 0; mask < lim; ++mask) {
        u64 prod = 1;
        int bits = 0;
        int mm = mask;
        bool ok = true;
        while (mm) {
            const int b = __builtin_ctz(mm);
            mm &= (mm - 1);
            ++bits;
            const u32 p = fac[b];
            if (prod > pmax / p) {
                ok = false;
                break;
            }
            prod *= p;
        }
        if (!ok) continue;
        const u64 term = pmax / prod;
        if (bits & 1) {
            res -= term;
        } else {
            res += term;
        }
    }
    return res;
}

static u64 solve(u32 n_max) {
    const u32 k_max = n_max + 4;
    const std::vector<u32> spf = build_spf(k_max);
    const u32 q_max = (k_max - 1) / 4;

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 4;
    if (threads > 12) threads = 12;
    if (q_max < 200'000) threads = 1;

    std::vector<std::thread> workers;
    std::vector<u64> partial(static_cast<std::size_t>(threads), 0);
    workers.reserve(threads);

    const u32 chunk = (q_max + threads - 1) / threads;
    for (unsigned t = 0; t < threads; ++t) {
        const u32 lo = static_cast<u32>(t) * chunk + 1;
        const u32 hi = std::min<u32>(q_max, static_cast<u32>(t + 1) * chunk);
        if (lo > hi) continue;
        workers.emplace_back([&, t, lo, hi]() {
            u64 local = 0;
            for (u32 q = lo; q <= hi; ++q) {
                u32 k = (q << 2) + 1;
                u32 n = k - 1;
                if (n <= n_max) local += count_relprimes_pie(q, k, spf);

                k = (q << 2) + 2;
                n = k - 2;
                if (n <= n_max) local += count_relprimes_pie(q, k, spf);

                k = (q << 2) + 3;
                n = k - 1;
                if (n <= n_max) local += count_relprimes_pie(q, k, spf);

                k = (q << 2) + 4;
                n = k - 4;
                if (n <= n_max) local += count_relprimes_pie(q, k, spf);
            }
            partial[static_cast<std::size_t>(t)] = local;
        });
    }

    for (auto& th : workers) th.join();
    u64 ans = 0;
    for (u64 v : partial) ans += v;
    return ans;
}

int main() {
    assert(solve(6) == 4);
    assert(solve(100) == 805);
    std::cout << solve(100'000'000) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

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

public class Euler935 {

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        if (n >= 1)
            spf[1] = 1;
        int[] primes = new int[n / 10 + 10]; // Rough upper bound
        int primeCount = 0;

        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                if (primeCount < primes.length)
                    primes[primeCount++] = i;
            }
            for (int j = 0; j < primeCount; ++j) {
                int p = primes[j];
                long v = (long) i * p;
                if (v > n)
                    break;
                spf[(int) v] = p;
                if (p == spf[i])
                    break;
            }
        }
        return spf;
    }

    static int factorDistinct(int x, int[] spf, int[] out) {
        int cnt = 0;
        while (x > 1) {
            int p = spf[x];
            out[cnt++] = p;
            do {
                x /= p;
            } while (x > 1 && x % p == 0);
        }
        return cnt;
    }

    static long countRelprimesPie(int pmax, int k, int[] spf) {
        int[] fac = new int[10];
        int m = factorDistinct(k, spf, fac);
        int lim = 1 << m;
        long res = 0;

        for (int mask = 0; mask < lim; ++mask) {
            long prod = 1;
            int bits = 0;
            int mm = mask;
            boolean ok = true;

            while (mm != 0) {
                int b = Integer.numberOfTrailingZeros(mm);
                mm &= (mm - 1);
                bits++;
                int p = fac[b];

                if (prod > pmax / p) {
                    ok = false;
                    break;
                }
                prod *= p;
            }

            if (!ok)
                continue;
            long term = pmax / prod;
            if ((bits & 1) != 0) {
                res -= term;
            } else {
                res += term;
            }
        }
        return res;
    }

    public static String solve(int nMax) {
        int kMax = nMax + 4;
        int[] spf = buildSpf(kMax);
        int qMax = (kMax - 1) / 4;

        int threads = Runtime.getRuntime().availableProcessors();
        if (threads < 1)
            threads = 1;
        if (threads > 12)
            threads = 12;
        if (qMax < 200000)
            threads = 1;

        int chunk = (qMax + threads - 1) / threads;
        long[] partial = new long[threads];

        if (threads == 1) {
            long local = 0;
            for (int q = 1; q <= qMax; ++q) {
                int k = (q << 2) + 1;
                int n = k - 1;
                if (n <= nMax)
                    local += countRelprimesPie(q, k, spf);

                k = (q << 2) + 2;
                n = k - 2;
                if (n <= nMax)
                    local += countRelprimesPie(q, k, spf);

                k = (q << 2) + 3;
                n = k - 1;
                if (n <= nMax)
                    local += countRelprimesPie(q, k, spf);

                k = (q << 2) + 4;
                n = k - 4;
                if (n <= nMax)
                    local += countRelprimesPie(q, k, spf);
            }
            return Long.toString(local);
        }

        try {
            ExecutorService executor = Executors.newFixedThreadPool(threads);
            List<Callable<Void>> tasks = new ArrayList<>();

            for (int t = 0; t < threads; ++t) {
                final int threadIdx = t;
                final int lo = t * chunk + 1;
                final int hi = Math.min(qMax, (t + 1) * chunk);
                if (lo > hi)
                    continue;

                tasks.add(() -> {
                    long local = 0;
                    for (int q = lo; q <= hi; ++q) {
                        int k = (q << 2) + 1;
                        int n = k - 1;
                        if (n <= nMax)
                            local += countRelprimesPie(q, k, spf);

                        k = (q << 2) + 2;
                        n = k - 2;
                        if (n <= nMax)
                            local += countRelprimesPie(q, k, spf);

                        k = (q << 2) + 3;
                        n = k - 1;
                        if (n <= nMax)
                            local += countRelprimesPie(q, k, spf);

                        k = (q << 2) + 4;
                        n = k - 4;
                        if (n <= nMax)
                            local += countRelprimesPie(q, k, spf);
                    }
                    partial[threadIdx] = local;
                    return null;
                });
            }

            List<Future<Void>> results = executor.invokeAll(tasks);
            for (Future<Void> res : results) {
                res.get(); // wait for completion
            }
            executor.shutdown();

        } catch (Exception e) {
            e.printStackTrace();
        }

        long ans = 0;
        for (long val : partial) {
            ans += val;
        }
        return Long.toString(ans);
    }

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