Problem 606: Gozinta Chains II

Project Euler Problem 606 Solution

EulerSolve provides an optimized solution for Project Euler Problem 606, Gozinta Chains II, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A gozinta chain for $n$ is a sequence $$1=a_0\lt a_1\lt\cdots\lt a_t=n,$$ where each term divides the next one. Let $g(n)$ be the number of such chains, and let $$S(N)=\sum_{\substack{k\le N\\ g(k)=252}} k.$$ The task is to find the last nine digits of $S(10^{36})$. A direct scan up to $10^{36}$ is hopeless, so the solution first characterizes exactly which prime-exponent patterns produce 252 chains, then sums those integers in a much more structured way. Mathematical Approach Write $$n=\prod_{j=1}^{r} p_j^{e_j},\qquad \Omega=e_1+\cdots+e_r.$$ The crucial observation is that the number of gozinta chains depends only on the exponent pattern $(e_1,\dots,e_r)$, not on the actual primes $p_j$. Step 1: Convert Divisors into Exponent Vectors Every divisor of $n$ can be written uniquely as $$d=\prod_{j=1}^{r} p_j^{b_j},\qquad 0\le b_j\le e_j.$$ So divisors correspond to lattice points $b=(b_1,\dots,b_r)$ inside the box $[0,e_1]\times\cdots\times[0,e_r]$. A gozinta chain is therefore a strictly increasing coordinatewise chain $$ (0,\dots,0)=v^{(0)}\lt v^{(1)}\lt\cdots\lt v^{(k)}=(e_1,\dots,e_r), $$ where at least one coordinate increases at every step. Since this description uses only the exponents, the chain count is completely determined by the exponent pattern....

Detailed mathematical approach

Problem Summary

A gozinta chain for $n$ is a sequence

$$1=a_0\lt a_1\lt\cdots\lt a_t=n,$$

where each term divides the next one. Let $g(n)$ be the number of such chains, and let

$$S(N)=\sum_{\substack{k\le N\\ g(k)=252}} k.$$

The task is to find the last nine digits of $S(10^{36})$. A direct scan up to $10^{36}$ is hopeless, so the solution first characterizes exactly which prime-exponent patterns produce 252 chains, then sums those integers in a much more structured way.

Mathematical Approach

Write

$$n=\prod_{j=1}^{r} p_j^{e_j},\qquad \Omega=e_1+\cdots+e_r.$$

The crucial observation is that the number of gozinta chains depends only on the exponent pattern $(e_1,\dots,e_r)$, not on the actual primes $p_j$.

Step 1: Convert Divisors into Exponent Vectors

Every divisor of $n$ can be written uniquely as

$$d=\prod_{j=1}^{r} p_j^{b_j},\qquad 0\le b_j\le e_j.$$

So divisors correspond to lattice points $b=(b_1,\dots,b_r)$ inside the box $[0,e_1]\times\cdots\times[0,e_r]$. A gozinta chain is therefore a strictly increasing coordinatewise chain

$$ (0,\dots,0)=v^{(0)}\lt v^{(1)}\lt\cdots\lt v^{(k)}=(e_1,\dots,e_r), $$

where at least one coordinate increases at every step. Since this description uses only the exponents, the chain count is completely determined by the exponent pattern.

Step 2: Count Chains with Exactly $k$ Steps

Fix a chain length $k$, meaning $k$ strict divisibility moves from $1$ to $n$. Let the increment in coordinate $j$ during step $t$ be $\delta_{j,t}$. Then

$$\delta_{j,t}\ge 0,\qquad \sum_{t=1}^{k}\delta_{j,t}=e_j.$$

For one fixed prime $p_j$, the number of weak compositions of $e_j$ into $k$ parts is

$$\binom{e_j+k-1}{k-1}.$$

Multiplying over all prime coordinates gives the number of increment tables if zero steps are temporarily allowed:

$$T_k(e_1,\dots,e_r)=\prod_{j=1}^{r}\binom{e_j+k-1}{k-1}.$$

However, $T_k$ still counts forbidden cases where some step has $\delta_{1,t}=\cdots=\delta_{r,t}=0$, which would repeat a divisor instead of moving to a strictly larger one.

Step 3: Enforce Strictness by Inclusion-Exclusion

To obtain genuine gozinta chains, every step must increase at least one coordinate. If exactly $i$ of the $k$ step positions are declared active, the same stars-and-bars argument gives

$$\prod_{j=1}^{r}\binom{e_j+i-1}{i-1}$$

possibilities. Choosing the active positions and applying inclusion-exclusion yields

$$ F_k(e_1,\dots,e_r)=\sum_{i=1}^{k}(-1)^{k-i}\binom{k}{i}\prod_{j=1}^{r}\binom{e_j+i-1}{i-1}. $$

This is the number of chains with exactly $k$ strict moves. Because each move raises the total exponent sum by at least $1$, we must have $1\le k\le \Omega$. Therefore

$$ G(e_1,\dots,e_r)=\sum_{k=1}^{\Omega} F_k(e_1,\dots,e_r) $$

is the total number of gozinta chains for any integer having exponent pattern $(e_1,\dots,e_r)$.

Step 4: Identify the Unique Pattern Giving 252

The implementations enumerate exponent partitions and evaluate the exact formula above. The outcome is strikingly simple: the only pattern with

$$G(e_1,\dots,e_r)=252$$

$$ (e_1,\dots,e_r)=(3,3). $$

Hence every qualifying integer has the form

$$ n=p^3q^3=(pq)^3,\qquad p\lt q \text{ prime}. $$

The primes must be distinct: $p^6$ has pattern $(6)$, not $(3,3)$, so it does not belong to the target set.

Step 5: Worked Example for the Pattern $(3,3)$

Take $n=2^3\cdot 3^3=216$. Here $\Omega=6$, and for a fixed step count $k$ we get

$$ T_k=\binom{3+k-1}{k-1}^2=\binom{k+2}{2}^2. $$

Applying inclusion-exclusion gives

$$ F_k=\sum_{i=1}^{k}(-1)^{k-i}\binom{k}{i}\binom{i+2}{2}^2. $$

The values are

$$ \begin{aligned} F_1&=1,\\ F_2&=14,\\ F_3&=55,\\ F_4&=92,\\ F_5&=70,\\ F_6&=20. \end{aligned} $$

Summing them shows

$$ 1+14+55+92+70+20=252. $$

So $216$ really has exactly 252 gozinta chains, and by the previous step every other qualifying number is obtained by replacing $2$ and $3$ with any distinct prime pair.

Step 6: Reduce $S(N)$ to a Sum over Prime Pairs

Let

$$ X=\left\lfloor N^{1/3}\right\rfloor. $$

Then $(pq)^3\le N$ is equivalent to $pq\le X$, so

$$ S(N)=\sum_{\substack{p\lt q\\ pq\le X}} (pq)^3. $$

For the actual problem, $N=10^{36}$, hence $X=10^{12}$.

Now define the prime-cube prefix sum

$$ P_3(y)=\sum_{\substack{q\le y\\ q\text{ prime}}} q^3. $$

Grouping terms by the smaller prime $p$ gives

$$ S(N)=\sum_{p\le \sqrt{X}} p^3\left(P_3\!\left(\left\lfloor\frac{X}{p}\right\rfloor\right)-P_3(p)\right). $$

The subtraction by $P_3(p)$ enforces $q>p$, so every unordered prime pair is counted exactly once.

Step 7: Compute $P_3$ with a Quotient-Based Prime Summatory Sieve

A full sieve up to $10^{12}$ would be too large, so the implementation works only with the distinct quotient values

$$ v=\left\lfloor\frac{X}{i}\right\rfloor, $$

of which there are only $O(\sqrt{X})$. For each such $v$, it starts from the ordinary cube sum

$$ \sum_{m=2}^{v} m^3=\left(\frac{v(v+1)}{2}\right)^2-1. $$

Then it performs Min_25-style prime corrections. When sieving by a prime $p$, every state with $v\ge p^2$ is updated by subtracting the contribution of numbers whose smallest remaining prime factor is $p$:

$$ H(v)\leftarrow H(v)-p^3\left(H\!\left(\left\lfloor\frac{v}{p}\right\rfloor\right)-\sum_{q<p} q^3\right). $$

After all primes up to $\sqrt{X}$ have been processed, $H(v)$ equals $P_3(v)$. Because only the last nine digits are needed, every arithmetic operation is reduced modulo $10^9$.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical reduction. They first verify the structural claim by evaluating the exact chain-count formula on exponent partitions and confirming that $(3,3)$ is the only pattern with 252 chains.

Next they set $X=10^{12}$, generate all primes up to $\sqrt{X}$, and build the table of distinct quotient values $\lfloor X/i\rfloor$. Each table entry is initialized with the closed form for $\sum m^3$, then corrected prime by prime until only prime cubes remain in the summatory table.

Finally they loop over the smaller prime $p$, query the prime-cube prefix sum at $\lfloor X/p\rfloor$ and at $p$, and accumulate

$$ p^3\left(P_3\!\left(\left\lfloor\frac{X}{p}\right\rfloor\right)-P_3(p)\right) $$

modulo $10^9$. Before attacking $10^{36}$, the implementation also checks smaller validated cases such as $S(10^6)=8462952$ and $S(10^{12})=623291998881978$.

Complexity Analysis

The exponent-pattern verification is tiny compared with the main computation. For the main solver, sieving primes up to $\sqrt{X}$ costs $O(\sqrt{X}\log\log X)$ time and $O(\sqrt{X})$ memory. The quotient table contains only $O(\sqrt{X})$ states. During the prime-correction phase, a prime $p$ touches only those states with $v\ge p^2$, so the number of table updates is

$$ \sum_{p\le \sqrt{X}} O\!\left(\min\!\left(\sqrt{X},\frac{X}{p^2}\right)\right), $$

which is sublinear in $X$ and follows the usual complexity profile of Min_25-style prime-summatory methods. The final accumulation over primes up to $\sqrt{X}$ is linear in $\pi(\sqrt{X})$. Overall, the implementation uses $O(\sqrt{X})$ memory and practical sublinear time for $X=10^{12}$.

Footnotes and References

Problem page: https://projecteuler.net/problem=606
Partially ordered sets and chains: Wikipedia - Partially ordered set
Stars and bars for weak compositions: Wikipedia - Stars and bars
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Prime sieving background: Wikipedia - Sieve of Eratosthenes
Prime-counting and related summatory methods: Wikipedia - Prime-counting function

Problem 606 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>
#include <cmath>

// Project Euler 606: only numbers of the form (p*q)^3 (p<q primes) have 252 gozinta chains,
// so S(N) reduces to sum_{p<q, pq <= cbrt(N)} (pq)^3. For cbrt(10^36)=10^12 we compute this
// modulo 10^9 using a Min_25 style prime power-sum sieve for sum_{p<=x} p^3.

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128_t;

static constexpr u64 MOD = 1000000000ULL;  // last 9 digits

static inline u64 add_mod(u64 a, u64 b) {
    a += b;
    if (a >= MOD) a -= MOD;
    return a;
}
static inline u64 sub_mod(u64 a, u64 b) { return (a >= b) ? (a - b) : (a + MOD - b); }
static inline u64 mul_mod(u64 a, u64 b) { return (u64)((u128)a * b % MOD); }

static u64 sum_cubes_mod(u64 n) {
    // sum_{i=1..n} i^3 = (n(n+1)/2)^2
    const u128 t = (u128)n * (n + 1) / 2;
    const u64 tm = (u64)(t % MOD);
    return mul_mod(tm, tm);
}

static std::vector<u64> sieve_primes(u64 n) {
    std::vector<bool> is_prime(n + 1, true);
    if (n >= 0) is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;
    for (u64 p = 2; p * p <= n; ++p) {
        if (!is_prime[p]) continue;
        for (u64 j = p * p; j <= n; j += p) is_prime[j] = false;
    }
    std::vector<u64> primes;
    primes.reserve(n / 10);
    for (u64 i = 2; i <= n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

struct PrimeCubeSummatory {
    u64 n = 0;
    u64 sq = 0;
    std::vector<u64> vals;   // distinct floor(n/i)
    std::vector<u64> g;      // sum_{p<=vals[j]} p^3 (mod MOD) after sieve
    std::vector<int> id1;    // vals index for x<=sq
    std::vector<int> id2;    // vals index for x>sq, keyed by n/x
    std::vector<u64> primes; // primes up to sq
    std::vector<u64> pref_p3;

    explicit PrimeCubeSummatory(u64 n_) : n(n_) {
        sq = (u64)std::sqrt((long double)n);
        primes = sieve_primes(sq);

        for (u64 l = 1; l <= n;) {
            const u64 w = n / l;
            vals.push_back(w);
            l = n / w + 1;
        }

        const int m = (int)vals.size();
        g.resize(m);
        id1.assign((size_t)sq + 1, -1);
        id2.assign((size_t)sq + 1, -1);
        for (int i = 0; i < m; ++i) {
            const u64 w = vals[i];
            g[i] = sub_mod(sum_cubes_mod(w), 1); // exclude 1^3
            if (w <= sq) id1[w] = i;
            else id2[n / w] = i;
        }

        pref_p3.assign(primes.size() + 1, 0);
        for (std::size_t i = 0; i < primes.size(); ++i) {
            const u64 p = primes[i];
            const u64 p3 = (u64)((u128)p * p % MOD * p % MOD);
            pref_p3[i + 1] = add_mod(pref_p3[i], p3);
        }

        for (std::size_t i = 0; i < primes.size(); ++i) {
            const u64 p = primes[i];
            const u64 p2 = p * p;
            if (p2 > n) break;
            const u64 p3 = (u64)((u128)p * p % MOD * p % MOD);

            for (int j = 0; j < m && vals[j] >= p2; ++j) {
                const u64 w = vals[j];
                const int idx = id(w / p);
                const u64 term = sub_mod(g[idx], pref_p3[i]);
                g[j] = sub_mod(g[j], mul_mod(p3, term));
            }
        }
    }

    inline int id(u64 x) const { return (x <= sq) ? id1[x] : id2[n / x]; }
    inline u64 sum_p3(u64 x) const { return g[id(x)]; }
};

static u64 g_from_pattern(const std::vector<int>& exps, u64 limit) {
    if (exps.empty()) return 1;
    int omega = 0;
    for (int e : exps) omega += e;

    std::vector<std::vector<u64>> C((size_t)omega + 1, std::vector<u64>((size_t)omega + 1, 0));
    C[0][0] = 1;
    for (int n = 1; n <= omega; ++n) {
        C[n][0] = C[n][n] = 1;
        for (int k = 1; k < n; ++k) C[n][k] = C[n - 1][k - 1] + C[n - 1][k];
    }

    auto binom_u64 = [](int n, int k) -> u64 {
        if (k < 0 || k > n) return 0;
        k = std::min(k, n - k);
        u128 r = 1;
        for (int i = 1; i <= k; ++i) r = (r * (u128)(n - k + i)) / (u128)i;
        return (u64)r;
    };

    std::vector<u128> total((size_t)omega + 1, 0);
    for (int k = 1; k <= omega; ++k) {
        u128 prod = 1;
        for (int a : exps) prod *= (u128)binom_u64(a + k - 1, k - 1);
        total[(size_t)k] = prod;
    }

    u128 g = 0;
    for (int k = 1; k <= omega; ++k) {
        i128 fk = 0;
        for (int i = 1; i <= k; ++i) {
            const u128 term = (u128)C[k][i] * total[(size_t)i];
            if (((k - i) & 1) != 0) fk -= (i128)term;
            else fk += (i128)term;
        }
        g += (u128)fk;
        if (g > (u128)limit) return limit + 1;
    }
    return (u64)g;
}

static void gen_partitions(int rem, int max_part, std::vector<int>& cur, std::vector<std::vector<int>>& out) {
    if (rem == 0) {
        out.push_back(cur);
        return;
    }
    for (int x = std::min(rem, max_part); x >= 1; --x) {
        cur.push_back(x);
        gen_partitions(rem - x, x, cur, out);
        cur.pop_back();
    }
}

static u128 brute_sum_pairs_cubed(u64 X) {
    std::vector<u64> primes = sieve_primes(X);
    u128 sum = 0;
    for (std::size_t i = 0; i < primes.size(); ++i) {
        const u64 p = primes[i];
        for (std::size_t j = i + 1; j < primes.size(); ++j) {
            const u64 q = primes[j];
            if ((u128)p * q > X) break;
            const u128 pq = (u128)p * q;
            sum += pq * pq * pq;
        }
    }
    return sum;
}

int main() {
    // Verify the exponent-pattern characterization of 252.
    std::vector<std::vector<int>> patterns;
    for (int omega = 1; omega <= 20; ++omega) {
        std::vector<std::vector<int>> parts;
        std::vector<int> cur;
        gen_partitions(omega, omega, cur, parts);
        for (auto& p : parts) {
            std::sort(p.begin(), p.end(), std::greater<int>());
            if (g_from_pattern(p, 252) == 252) patterns.push_back(p);
        }
    }
    std::sort(patterns.begin(), patterns.end());
    patterns.erase(std::unique(patterns.begin(), patterns.end()), patterns.end());
    assert(patterns.size() == 1 && patterns[0].size() == 2 && patterns[0][0] == 3 && patterns[0][1] == 3);

    // Statement validations.
    assert((u64)brute_sum_pairs_cubed(100) == 8462952ULL); // S(10^6)
    assert((u64)brute_sum_pairs_cubed(10000) == 623291998881978ULL); // S(10^12)

    const u64 X = 1000000000000ULL; // cbrt(10^36)
    PrimeCubeSummatory pc(X);

    u64 ans = 0;
    for (u64 p : pc.primes) { // primes up to sqrt(X)
        const u64 lim = X / p;
        if (lim <= p) break;
        const u64 sum_q = sub_mod(pc.sum_p3(lim), pc.sum_p3(p));
        const u64 p3 = (u64)((u128)p * p % MOD * p % MOD);
        ans = add_mod(ans, mul_mod(p3, sum_q));
    }

    std::cout << std::setw(9) << std::setfill('0') << ans << "\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.ArrayList;
import java.util.List;

public class Euler606 {
    static final long MOD = 1000000000L;

    static long addMod(long a, long b) {
        long s = a + b;
        if (s >= MOD)
            s -= MOD;
        return s;
    }

    static long subMod(long a, long b) {
        return (a >= b) ? (a - b) : (a + MOD - b);
    }

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

    static long sumCubesMod(long n) {
        long tm = 0;
        if (n % 2 == 0) {
            tm = (n / 2 % MOD) * ((n + 1) % MOD) % MOD;
        } else {
            tm = (n % MOD) * ((n + 1) / 2 % MOD) % MOD;
        }
        return mulMod(tm, tm);
    }

    static List<Long> sievePrimes(long n) {
        int limit = (int) n;
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrime[i] = true;
        for (int p = 2; p * p <= limit; p++) {
            if (isPrime[p]) {
                for (int j = p * p; j <= limit; j += p) {
                    isPrime[j] = false;
                }
            }
        }
        List<Long> primes = new ArrayList<>();
        for (int i = 2; i <= limit; i++) {
            if (isPrime[i])
                primes.add((long) i);
        }
        return primes;
    }

    public static String solve() {
        long X = 1000000000000L;
        long sq = (long) Math.sqrt(X);
        List<Long> primes = sievePrimes(sq);

        List<Long> vals = new ArrayList<>();
        for (long l = 1; l <= X;) {
            long w = X / l;
            vals.add(w);
            l = X / w + 1;
        }

        int m = vals.size();
        long[] g = new long[m];
        int[] id1 = new int[(int) sq + 1];
        int[] id2 = new int[(int) sq + 1];

        for (int i = 0; i < m; i++) {
            long w = vals.get(i);
            g[i] = subMod(sumCubesMod(w), 1);
            if (w <= sq) {
                id1[(int) w] = i;
            } else {
                id2[(int) (X / w)] = i;
            }
        }

        long[] prefP3 = new long[primes.size() + 1];
        for (int i = 0; i < primes.size(); i++) {
            long p = primes.get(i);
            long p3 = p * p % MOD * p % MOD;
            prefP3[i + 1] = addMod(prefP3[i], p3);
        }

        for (int i = 0; i < primes.size(); i++) {
            long p = primes.get(i);
            long p2 = p * p;
            if (p2 > X)
                break;
            long p3 = p * p % MOD * p % MOD;

            for (int j = 0; j < m; j++) {
                long w = vals.get(j);
                if (w < p2)
                    break;
                long x = w / p;
                int idx = (x <= sq) ? id1[(int) x] : id2[(int) (X / x)];
                long term = subMod(g[idx], prefP3[i]);
                g[j] = subMod(g[j], mulMod(p3, term));
            }
        }

        long ans = 0;
        for (int i = 0; i < primes.size(); i++) {
            long p = primes.get(i);
            long lim = X / p;
            if (lim <= p)
                break;

            int idxLim = (lim <= sq) ? id1[(int) lim] : id2[(int) (X / lim)];
            int idxP = (p <= sq) ? id1[(int) p] : id2[(int) (X / p)];

            long sumQ = subMod(g[idxLim], g[idxP]);
            long p3 = p * p % MOD * p % MOD;
            ans = addMod(ans, mulMod(p3, sumQ));
        }

        return String.format("%09d", ans);
    }

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