Problem 511: Sequences with Nice Divisibility Properties

Project Euler Problem 511 Solution

EulerSolve provides an optimized solution for Project Euler Problem 511, Sequences with Nice Divisibility Properties, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $A(n,k)$ be the number of length-$n$ sequences $(a_1,\dots,a_n)$ such that every term divides $n$ and $$a_1+\cdots+a_n \equiv -n \pmod{k}.$$ The answer is required modulo $10^9$. Since $n$ can be enormous, the solution does not build sequences directly. Instead, it compresses the divisor set of $n$ by residue modulo $k$, then raises that residue distribution to the $n$-th power under cyclic convolution. Mathematical Approach Write $D(n)=\{d\in\mathbb{Z}_{>0}: d\mid n\}$, and let $M=10^9$. The key observation is that only residues modulo $k$ matter, so the entire counting problem lives inside the cyclic group $\mathbb{Z}_k$. Step 1: Collapse the divisors into residue counts For each residue $r\in\{0,1,\dots,k-1\}$, define $$c_r=\left|\left\{d\in D(n): d\equiv r \pmod{k}\right\}\right|.$$ The vector $\mathbf{c}=(c_0,\dots,c_{k-1})$ records how many one-step choices land in each residue class. This loses no relevant information, because any future sum is reduced modulo $k$. Step 2: Concatenating choices becomes cyclic convolution Suppose $\mathbf{f}$ and $\mathbf{g}$ describe the residue distributions of two independent blocks of choices. If the first block contributes residue $i$ and the second contributes residue $j$, then the combined residue is $i+j \pmod{k}$....

Detailed mathematical approach

Problem Summary

Let $A(n,k)$ be the number of length-$n$ sequences $(a_1,\dots,a_n)$ such that every term divides $n$ and

$$a_1+\cdots+a_n \equiv -n \pmod{k}.$$

The answer is required modulo $10^9$. Since $n$ can be enormous, the solution does not build sequences directly. Instead, it compresses the divisor set of $n$ by residue modulo $k$, then raises that residue distribution to the $n$-th power under cyclic convolution.

Mathematical Approach

Write $D(n)=\{d\in\mathbb{Z}_{>0}: d\mid n\}$, and let $M=10^9$. The key observation is that only residues modulo $k$ matter, so the entire counting problem lives inside the cyclic group $\mathbb{Z}_k$.

Step 1: Collapse the divisors into residue counts

For each residue $r\in\{0,1,\dots,k-1\}$, define

$$c_r=\left|\left\{d\in D(n): d\equiv r \pmod{k}\right\}\right|.$$

The vector $\mathbf{c}=(c_0,\dots,c_{k-1})$ records how many one-step choices land in each residue class. This loses no relevant information, because any future sum is reduced modulo $k$.

Step 2: Concatenating choices becomes cyclic convolution

Suppose $\mathbf{f}$ and $\mathbf{g}$ describe the residue distributions of two independent blocks of choices. If the first block contributes residue $i$ and the second contributes residue $j$, then the combined residue is $i+j \pmod{k}$. Therefore

$$\left(\mathbf{f}\star\mathbf{g}\right)_t=\sum_{i=0}^{k-1} f_i\,g_{(t-i)\bmod k}.$$

This is circular convolution on $\mathbb{Z}_k$. After one pick the distribution is $\mathbf{c}$; after two picks it is $\mathbf{c}\star\mathbf{c}$; after $m$ picks it is $\mathbf{c}^{\star m}$.

Step 3: The target residue is fixed

The required sequences satisfy

$$a_1+\cdots+a_n \equiv -n \pmod{k}.$$

So the residue index to extract is

$$t^\ast=(-n)\bmod k=\bigl(k-(n\bmod k)\bigr)\bmod k.$$

Hence the desired count is

$$A(n,k)=\left(\mathbf{c}^{\star n}\right)_{t^\ast}\pmod{M}.$$

Step 4: Polynomial viewpoint in a quotient ring

Encode the residue vector as

$$P(x)=\sum_{r=0}^{k-1} c_r x^r$$

inside the ring $(\mathbb{Z}/M\mathbb{Z})[x]/(x^k-1)$. There, $x^a x^b=x^{(a+b)\bmod k}$, so polynomial multiplication is exactly cyclic convolution. Therefore the coefficient of $x^t$ in $P(x)^n$ counts length-$n$ divisor sequences whose sum is congruent to $t \pmod{k}$. The answer is the coefficient of $x^{t^\ast}$.

Step 5: Binary exponentiation removes the linear dependence on $n$

A direct dynamic program would perform one convolution per position and cost $O(nk^2)$, which is impossible for the real input. Because convolution is associative, repeated squaring works: compute $\mathbf{c}, \mathbf{c}^{\star 2}, \mathbf{c}^{\star 4}, \dots$, and combine only the powers corresponding to the set bits of $n$. This reduces the main cost to $O(k^2\log n)$.

Worked Example: $A(3,4)=4$

The divisors of $3$ are $1$ and $3$, so modulo $4$ the one-step vector is

$$\mathbf{c}=(0,1,0,1).$$

One more convolution gives

$$\mathbf{c}^{\star 2}=(2,0,2,0),$$

and the third choice gives

$$\mathbf{c}^{\star 3}=(0,4,0,4).$$

Since

$$t^\ast=(-3)\bmod 4=1,$$

the required component is $4$. Equivalently, every term is either $1$ or $3$, and the total is $1 \pmod 4$ exactly when an odd number of terms are $3$, giving $\binom31+\binom33=4$.

How the Code Works

The C++ and Java implementations first factor $n$ by trial division and then enumerate every divisor recursively from the prime factorization. Those divisors are collapsed into a length-$k$ count vector according to their residues modulo $k$.

Next, the implementation performs circular convolution between two length-$k$ vectors, always reducing arithmetic modulo $10^9$. The compiled implementations use wider temporary accumulators and reduce periodically so that intermediate products remain safe.

Repeated squaring raises the one-step residue vector to the $n$-th convolution power, and the final answer is the component at residue $(-n)\bmod k$. The Python implementation returns the same value by delegating to the compiled solver, so the C++, Python, and Java implementations all follow the same mathematical method.

Complexity Analysis

Factoring $n$ by trial division costs $O(\sqrt{n})$ time in the worst case. If $\tau(n)$ is the number of divisors of $n$, then enumerating all divisors and filling the residue-count vector costs $O(\tau(n))$ time and $O(\tau(n))$ storage.

Each circular convolution of two vectors of length $k$ costs $O(k^2)$ time and $O(k)$ additional working memory. Binary exponentiation uses $O(\log n)$ such convolutions, so the full method runs in

$$O\!\left(\sqrt{n}+\tau(n)+k^2\log n\right)$$

time and uses

$$O\!\left(k+\tau(n)\right)$$

memory. For the actual input size, the convolution phase is the dominant cost.

Footnotes and References

Problem page: https://projecteuler.net/problem=511
Circular convolution: Wikipedia — Circular convolution
Exponentiation by squaring: Wikipedia — Exponentiation by squaring
Divisor: Wikipedia — Divisor
Modular arithmetic: Wikipedia — Modular arithmetic

Problem 511 source code

C++

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

namespace {

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

constexpr u64 kMod = 1'000'000'000ULL;

std::vector<std::pair<u64, int>> factorize(u64 n) {
    std::vector<std::pair<u64, int>> factors;
    u64 x = n;
    for (u64 p = 2ULL; p * p <= x; p += (p == 2ULL ? 1ULL : 2ULL)) {
        if (x % p != 0ULL) {
            continue;
        }
        int e = 0;
        while (x % p == 0ULL) {
            x /= p;
            ++e;
        }
        factors.push_back({p, e});
    }
    if (x > 1ULL) {
        factors.push_back({x, 1});
    }
    return factors;
}

void build_divisors_rec(const std::vector<std::pair<u64, int>>& factors, const int idx,
                        const u64 cur, std::vector<u64>& out) {
    if (idx == static_cast<int>(factors.size())) {
        out.push_back(cur);
        return;
    }
    const auto [p, e] = factors[static_cast<std::size_t>(idx)];
    u64 value = 1ULL;
    for (int i = 0; i <= e; ++i) {
        build_divisors_rec(factors, idx + 1, cur * value, out);
        value *= p;
    }
}

std::vector<u64> divisors_of(u64 n) {
    const auto factors = factorize(n);
    std::vector<u64> divisors;
    build_divisors_rec(factors, 0, 1ULL, divisors);
    return divisors;
}

std::vector<u32> cyclic_convolution(const std::vector<u32>& a, const std::vector<u32>& b,
                                    const u64 mod) {
    const int k = static_cast<int>(a.size());
    std::vector<u64> accum(static_cast<std::size_t>(k), 0ULL);
    constexpr u64 kReduceThreshold = (1ULL << 62);

    for (int i = 0; i < k; ++i) {
        const u64 ai = a[static_cast<std::size_t>(i)];
        if (ai == 0ULL) {
            continue;
        }
        const int split = k - i;
        for (int j = 0; j < split; ++j) {
            u64& cell = accum[static_cast<std::size_t>(i + j)];
            cell += ai * static_cast<u64>(b[static_cast<std::size_t>(j)]);
            if (cell >= kReduceThreshold) {
                cell %= mod;
            }
        }
        for (int j = split; j < k; ++j) {
            u64& cell = accum[static_cast<std::size_t>(i + j - k)];
            cell += ai * static_cast<u64>(b[static_cast<std::size_t>(j)]);
            if (cell >= kReduceThreshold) {
                cell %= mod;
            }
        }
    }

    std::vector<u32> out(static_cast<std::size_t>(k), 0U);
    for (int i = 0; i < k; ++i) {
        out[static_cast<std::size_t>(i)] =
            static_cast<u32>(accum[static_cast<std::size_t>(i)] % mod);
    }
    return out;
}

std::vector<u32> cyclic_power(std::vector<u32> base, u64 exp, const u64 mod) {
    const int k = static_cast<int>(base.size());
    std::vector<u32> result(static_cast<std::size_t>(k), 0U);
    result[0] = 1U;

    while (exp > 0ULL) {
        if (exp & 1ULL) {
            result = cyclic_convolution(result, base, mod);
        }
        exp >>= 1ULL;
        if (exp > 0ULL) {
            base = cyclic_convolution(base, base, mod);
        }
    }
    return result;
}

u64 seq_mod(const u64 n, const int k) {
    const std::vector<u64> divisors = divisors_of(n);
    std::vector<u32> base(static_cast<std::size_t>(k), 0U);
    for (const u64 d : divisors) {
        const int residue = static_cast<int>(d % static_cast<u64>(k));
        base[static_cast<std::size_t>(residue)] += 1U;
    }

    const std::vector<u32> ways = cyclic_power(base, n, kMod);
    const int target = static_cast<int>((static_cast<u64>(k) - (n % static_cast<u64>(k))) %
                                        static_cast<u64>(k));
    return static_cast<u64>(ways[static_cast<std::size_t>(target)]);
}

u64 seq_mod_small_dp(const int n, const int k) {
    const std::vector<u64> divisors = divisors_of(static_cast<u64>(n));
    std::vector<u64> dp(static_cast<std::size_t>(k), 0ULL);
    dp[0] = 1ULL;

    for (int step = 0; step < n; ++step) {
        std::vector<u64> next(static_cast<std::size_t>(k), 0ULL);
        for (int r = 0; r < k; ++r) {
            if (dp[static_cast<std::size_t>(r)] == 0ULL) {
                continue;
            }
            for (const u64 d : divisors) {
                const int nr = (r + static_cast<int>(d % static_cast<u64>(k))) % k;
                next[static_cast<std::size_t>(nr)] += dp[static_cast<std::size_t>(r)];
            }
        }
        dp.swap(next);
    }
    const int target = (k - (n % k)) % k;
    return dp[static_cast<std::size_t>(target)] % kMod;
}

bool run_checkpoints() {
    if (seq_mod(3ULL, 4) != 4ULL) {
        std::cerr << "Checkpoint failed: Seq(3,4)\n";
        return false;
    }
    if (seq_mod(4ULL, 11) != 8ULL) {
        std::cerr << "Checkpoint failed: Seq(4,11)\n";
        return false;
    }
    if (seq_mod(1'111ULL, 24) != 840'643'584ULL) {
        std::cerr << "Checkpoint failed: Seq(1111,24) mod 1e9\n";
        return false;
    }
    if (seq_mod(8ULL, 15) != seq_mod_small_dp(8, 15)) {
        std::cerr << "Checkpoint failed: fast/DP mismatch\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 n = 1'234'567'898'765ULL;
    constexpr int k = 4'321;
    const u64 answer = seq_mod(n, k);
    std::cout << std::setw(9) << std::setfill('0') << answer << '\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 Euler511 {
    static class Factor {
        long p;
        int e;

        Factor(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Long> getDivisors(long n) {
        List<Factor> factors = new ArrayList<>();
        long x = n;
        for (long p = 2; p * p <= x; p += (p == 2 ? 1 : 2)) {
            if (x % p == 0) {
                int e = 0;
                while (x % p == 0) {
                    e++;
                    x /= p;
                }
                factors.add(new Factor(p, e));
            }
        }
        if (x > 1) {
            factors.add(new Factor(x, 1));
        }

        List<Long> divs = new ArrayList<>();
        buildDivisors(factors, 0, 1L, divs);
        return divs;
    }

    static void buildDivisors(List<Factor> factors, int idx, long current, List<Long> divs) {
        if (idx == factors.size()) {
            divs.add(current);
            return;
        }
        Factor f = factors.get(idx);
        long val = 1;
        for (int i = 0; i <= f.e; i++) {
            buildDivisors(factors, idx + 1, current * val, divs);
            val *= f.p;
        }
    }

    static int[] cyclicConvolution(int[] a, int[] b, int mod) {
        int k = a.length;
        long[] out = new long[k];
        for (int i = 0; i < k; i++) {
            long ai = a[i];
            if (ai == 0)
                continue;
            int split = k - i;
            for (int j = 0; j < split; j++) {
                out[i + j] += ai * b[j];
                if (out[i + j] >= 4000000000000000000L) { // Prevent overflow
                    out[i + j] %= mod;
                }
            }
            for (int j = split; j < k; j++) {
                out[i + j - k] += ai * b[j];
                if (out[i + j - k] >= 4000000000000000000L) {
                    out[i + j - k] %= mod;
                }
            }
        }
        int[] res = new int[k];
        for (int i = 0; i < k; i++) {
            res[i] = (int) (out[i] % mod);
        }
        return res;
    }

    static int[] cyclicPower(int[] base, long exp, int mod) {
        int k = base.length;
        int[] res = new int[k];
        res[0] = 1;

        while (exp > 0) {
            if ((exp & 1) == 1) {
                res = cyclicConvolution(res, base, mod);
            }
            exp >>= 1;
            if (exp > 0) {
                base = cyclicConvolution(base, base, mod);
            }
        }
        return res;
    }

    public static void main(String[] args) {
        long n = 1234567898765L;
        int k = 4321;
        int mod = 1000000000;

        List<Long> divisors = getDivisors(n);
        int[] base = new int[k];
        for (long d : divisors) {
            base[(int) (d % k)]++;
        }

        int[] ways = cyclicPower(base, n, mod);
        int target = (int) ((k - (n % k)) % k);

        System.out.printf("%09d\n", ways[target]);
    }
}