Problem 927: Prime-ary Tree

Project Euler Problem 927 Solution

EulerSolve provides an optimized solution for Project Euler Problem 927, Prime-ary Tree, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementations do not explore the prime-ary tree by testing every integer up to $N$ separately. Instead, they use a structural reduction: first identify which primes $q$ satisfy the orbit condition that governs the problem, and then sum the squarefree products built from those primes. For the full task, the limit is $N=10^7$. The dynamical object attached to a prime modulus $q$ and a prime exponent $r$ is the orbit $$x_0=1,\qquad x_{k+1}\equiv x_k^r+1 \pmod q.$$ A prime $q$ is kept precisely when this orbit reaches $0$ for every distinct prime divisor $r\mid(q-1)$. Once those primes are known, the answer is the sum of all admissible squarefree products that do not exceed $N$. Mathematical Approach Write $\mathcal P(m)$ for the set of distinct prime divisors of $m$. The computation naturally splits into two layers: classify the primes that survive the orbit test, then enumerate the integers obtained by multiplying those primes without repetition. The orbit attached to a pair $(q,r)$ For a prime $q$ and a prime $r\in\mathcal P(q-1)$, define $$x_0^{(q,r)}=1,\qquad x_{k+1}^{(q,r)}\equiv \left(x_k^{(q,r)}\right)^r+1 \pmod q.$$ The implementations use the following acceptance criterion: $$q\in\mathcal G \iff \forall r\in\mathcal P(q-1),\ \exists k\ge 0\text{ such that }x_k^{(q,r)}=0.$$ Only distinct prime divisors of $q-1$ matter....

Detailed mathematical approach

Problem Summary

The implementations do not explore the prime-ary tree by testing every integer up to $N$ separately. Instead, they use a structural reduction: first identify which primes $q$ satisfy the orbit condition that governs the problem, and then sum the squarefree products built from those primes. For the full task, the limit is $N=10^7$.

The dynamical object attached to a prime modulus $q$ and a prime exponent $r$ is the orbit

$$x_0=1,\qquad x_{k+1}\equiv x_k^r+1 \pmod q.$$

A prime $q$ is kept precisely when this orbit reaches $0$ for every distinct prime divisor $r\mid(q-1)$. Once those primes are known, the answer is the sum of all admissible squarefree products that do not exceed $N$.

Mathematical Approach

Write $\mathcal P(m)$ for the set of distinct prime divisors of $m$. The computation naturally splits into two layers: classify the primes that survive the orbit test, then enumerate the integers obtained by multiplying those primes without repetition.

The orbit attached to a pair $(q,r)$

For a prime $q$ and a prime $r\in\mathcal P(q-1)$, define

$$x_0^{(q,r)}=1,\qquad x_{k+1}^{(q,r)}\equiv \left(x_k^{(q,r)}\right)^r+1 \pmod q.$$

The implementations use the following acceptance criterion:

$$q\in\mathcal G \iff \forall r\in\mathcal P(q-1),\ \exists k\ge 0\text{ such that }x_k^{(q,r)}=0.$$

Only distinct prime divisors of $q-1$ matter. If $r^a$ divides $q-1$, the condition attached to $r$ is still checked only once, because the orbit depends on the prime exponent $r$, not on its multiplicity in $q-1$.

Why the test is finite

Modulo $q$ there are only $q$ possible residues, so every orbit is a finite-state process. Therefore exactly one of two things must happen:

$$\text{either some }x_k^{(q,r)}=0,\qquad\text{or a nonzero state repeats before }0\text{ appears.}$$

If the orbit ever repeats a nonzero value, the future evolution is periodic and $0$ will never be reached afterward. This is why cycle detection with constant memory is enough; there is no need to store a full visited set.

There is also a useful invariant at the moment of success. If $x_k^{(q,r)}=0$, then

$$x_{k+1}^{(q,r)}\equiv 0^r+1\equiv 1 \pmod q.$$

So hitting $0$ means that the orbit starting from $1$ has closed into a cycle containing both $1$ and $0$. The code exploits exactly this hit-or-cycle dichotomy.

From qualifying primes to admissible nodes

Let $\mathcal G$ be the set of primes that pass every required orbit test. The implementations then treat the admissible integers up to $N$ as

$$\mathcal S(N)=\left\{\prod_{q\in T} q:\ T\subseteq\mathcal G,\ \prod_{q\in T} q\le N\right\}.$$

The empty subset contributes the empty product $1$, so $1$ is always included. The desired sum is therefore

$$A(N)=\sum_{n\in\mathcal S(N)} n=\sum_{\substack{T\subseteq\mathcal G\\ \prod_{q\in T}\le N}}\prod_{q\in T} q.$$

This description is squarefree by construction: each accepted prime can appear at most once. The compiled implementations even contain spot-checks against naively extending an accepted prime to $q^2$, which supports the fact that repeated prime powers are not part of the final search space.

Hand-checking the smallest primes

The first few primes already show how the criterion behaves.

For $q=2$, we have $q-1=1$, so $\mathcal P(q-1)=\varnothing$. The condition is vacuous, hence $2\in\mathcal G$.

For $q=3$, the only relevant exponent is $r=2$, and the orbit is

$$1\to 2\to 2\to 2\to\cdots \pmod 3,$$

so $0$ is never reached. Therefore $3\notin\mathcal G$.

For $q=5$, again only $r=2$ matters, but now

$$1\to 2\to 0\to 1\to\cdots \pmod 5,$$

so $5\in\mathcal G$.

For $q=7$, the relevant exponents are $2$ and $3$. The orbit for $r=3$ is

$$1\to 2\to 2\to 2\to\cdots \pmod 7,$$

which already fails, so $7\notin\mathcal G$.

Hence, up to $20$, the accepted primes are exactly $\{2,5\}$. The admissible squarefree products are

$$1,\ 2,\ 5,\ 10,$$

and therefore

$$A(20)=1+2+5+10=18.$$

This is the first built-in checkpoint in the implementations. A second checkpoint is $A(1000)=2089$, which tests the same logic on a larger range.

How the Code Works

All three language paths follow the same mathematical decomposition: build arithmetic infrastructure first, classify primes second, and enumerate squarefree products last.

Screening primes by orbit detection

The first stage is a linear sieve up to $N$. It produces the prime list together with the smallest prime factor of every integer up to the limit, so factoring $q-1$ is fast for each candidate prime $q$.

For every $q$, the implementation extracts the distinct primes in $\mathcal P(q-1)$ and runs the orbit test modulo $q$ for each of them. The transition $x\mapsto x^r+1\pmod q$ is computed directly when $r=2$, and by modular exponentiation for larger $r$. Because the orbit can only hit $0$ or enter a disjoint cycle, constant-memory cycle detection is sufficient.

The C++ implementation parallelizes this prime-screening phase across several workers. The Java implementation performs the same logic serially. The Python entry point reuses the same compiled computation instead of re-implementing the full search in pure Python.

Summing the squarefree products

After $\mathcal G$ has been constructed, the remaining task is a depth-first enumeration over increasing prime indices. At each recursive state, the current product is added to the running sum, and the search continues only with later accepted primes. That ordering rule automatically prevents repeated use of the same prime.

The essential pruning inequality is

$$\text{current} \gt \frac{N}{q_{\text{next}}}.$$

Once it holds, multiplying by the next available prime would already exceed $N$, so the entire suffix of that branch can be skipped. The compiled implementations also keep the checkpoints $A(20)=18$ and $A(1000)=2089$, plus small spot-checks involving $q^2$, before running the full computation at $N=10^7$.

Complexity Analysis

The sieve stage is $O(N)$ time and $O(N)$ memory in the implemented model, because it stores the smallest prime factor of every integer up to $N$. Once that table exists, extracting $\mathcal P(q-1)$ is cheap.

For a fixed test pair $(q,r)$, the orbit phase inspects at most $q$ residues before it either finds $0$ or proves that the orbit cycles elsewhere. Each step costs $O(1)$ when $r=2$ and $O(\log r)$ modular multiplications otherwise. A direct expression for the screening cost is

$$O\!\left(N+\sum_{\substack{q\le N\\ q\text{ prime}}}\ \sum_{r\in\mathcal P(q-1)} q\log r\right).$$

The final depth-first search is output-sensitive in practice rather than powerset-sized, because it explores only branches whose current squarefree product is still at most $N$. Its memory usage is the recursion depth plus the stored list of accepted primes, both much smaller than the sieve table.

Footnotes and References

Problem page: Project Euler 927
Modular arithmetic: Wikipedia - Modular arithmetic
Modular exponentiation: Wikipedia - Modular exponentiation
Cycle detection: Wikipedia - Brent's algorithm
Squarefree integers: Wikipedia - Square-free integer
Integer factorization: Wikipedia - Integer factorization

Problem 927 source code

C++

#include <pthread.h>
#include <algorithm>
#include <atomic>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <unistd.h>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

u64 mod_pow(u64 base, int exp, int mod) {
    if (mod == 1) {
        return 0;
    }
    u64 result = 1 % mod;
    base %= static_cast<u64>(mod);
    while (exp > 0) {
        if (exp & 1) {
            result = (result * base) % static_cast<u64>(mod);
        }
        base = (base * base) % static_cast<u64>(mod);
        exp >>= 1;
    }
    return result;
}

u64 next_term(u64 x, int p, int mod) {
    if (p == 2) {
        return (x * x + 1ULL) % static_cast<u64>(mod);
    }
    return (mod_pow(x, p, mod) + 1ULL) % static_cast<u64>(mod);
}

bool hits_zero(int mod, int p) {
    if (mod == 1) {
        return true;
    }

    const u64 mod_u = static_cast<u64>(mod);
    const u64 start = 1ULL % mod_u;
    if (start == 0ULL) {
        return true;
    }

    u64 tortoise = start;
    u64 hare = next_term(start, p, mod);
    if (hare == 0ULL) {
        return true;
    }

    u64 power = 1ULL;
    u64 lam = 1ULL;
    while (tortoise != hare) {
        if (power == lam) {
            tortoise = hare;
            power <<= 1U;
            lam = 0ULL;
        }
        hare = next_term(hare, p, mod);
        if (hare == 0ULL) {
            return true;
        }
        ++lam;
    }

    return false;
}

struct Sieve {
    std::vector<int> primes;
    std::vector<int> spf;
};

Sieve build_sieve(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);

    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const i64 v = static_cast<i64>(p) * i;
            if (v > n || p > spf[i]) {
                break;
            }
            spf[static_cast<int>(v)] = p;
        }
    }

    return {std::move(primes), std::move(spf)};
}

int factor_list(int x, const std::vector<int>& spf, int* factors, const int limit) {
    int count = 0;
    while (x > 1) {
        const int p = spf[x];
        if (count < limit) {
            factors[count] = p;
        }
        ++count;

        while (x % p == 0) {
            x /= p;
        }
    }
    return count;
}

bool is_good_prime(int q, const std::vector<int>& spf) {
    if (q == 2) {
        return true;
    }

    int factor_buffer[32];
    const int divisor_count = factor_list(q - 1, spf, factor_buffer,
                                          static_cast<int>(std::size(factor_buffer)));
    std::sort(factor_buffer, factor_buffer + divisor_count, std::greater<int>());

    for (int i = 0; i < divisor_count; ++i) {
        const int p = factor_buffer[i];
        if (!hits_zero(q, p)) {
            return false;
        }
    }
    return true;
}

int detect_thread_count(std::size_t work_items) {
    long cores = ::sysconf(_SC_NPROCESSORS_ONLN);
    int threads = (cores > 0) ? static_cast<int>(cores) : 4;
    if (threads < 1) threads = 1;
    if (work_items == 0) return 1;
    if (static_cast<std::size_t>(threads) > work_items) {
        threads = static_cast<int>(work_items);
    }
    return threads;
}

struct PrimeWorkerTask {
    const std::vector<int>* primes = nullptr;
    const std::vector<int>* spf = nullptr;
    std::atomic<std::size_t>* next_idx = nullptr;
    std::vector<int> local_good;
};

void* prime_worker_entry(void* raw) {
    auto* task = static_cast<PrimeWorkerTask*>(raw);
    while (true) {
        const std::size_t idx = task->next_idx->fetch_add(1, std::memory_order_relaxed);
        if (idx >= task->primes->size()) break;
        const int q = (*task->primes)[idx];
        if (is_good_prime(q, *task->spf)) {
            task->local_good.push_back(q);
        }
    }
    return nullptr;
}

std::vector<int> good_primes_upto(int limit) {
    const Sieve sieve = build_sieve(limit);

    const int threads = detect_thread_count(sieve.primes.size());
    if (threads <= 1 || sieve.primes.size() < 5000U) {
        std::vector<int> good;
        good.reserve(sieve.primes.size() / 8);
        for (int q : sieve.primes) {
            if (is_good_prime(q, sieve.spf)) {
                good.push_back(q);
            }
        }
        return good;
    }

    std::atomic<std::size_t> next_idx{0};
    std::vector<pthread_t> handles(static_cast<std::size_t>(threads));
    std::vector<PrimeWorkerTask> tasks(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        auto& task = tasks[static_cast<std::size_t>(t)];
        task.primes = &sieve.primes;
        task.spf = &sieve.spf;
        task.next_idx = &next_idx;
        task.local_good.clear();
        task.local_good.reserve(sieve.primes.size() / (8U * static_cast<std::size_t>(threads)) + 8U);
        pthread_create(&handles[static_cast<std::size_t>(t)], nullptr, prime_worker_entry, &task);
    }

    std::vector<int> good;
    good.reserve(sieve.primes.size() / 8);
    for (int t = 0; t < threads; ++t) {
        pthread_join(handles[static_cast<std::size_t>(t)], nullptr);
        auto& vec = tasks[static_cast<std::size_t>(t)].local_good;
        good.insert(good.end(), vec.begin(), vec.end());
    }
    std::sort(good.begin(), good.end());
    return good;
}

u64 sum_squarefree_products_leq(const std::vector<int>& primes, i64 limit) {
    u64 total = 0;

    auto dfs = [&](auto&& self, std::size_t idx, i64 current) -> void {
        total += static_cast<u64>(current);
        for (std::size_t i = idx; i < primes.size(); ++i) {
            const i64 p = primes[i];
            if (current > limit / p) {
                break;
            }
            self(self, i + 1, current * p);
        }
    };

    dfs(dfs, 0, 1);
    return total;
}

u64 solve(i64 n) {
    const std::vector<int> good = good_primes_upto(static_cast<int>(n));

    for (int q : good) {
        if (static_cast<i64>(q) * q > n) {
            continue;
        }
        assert(!hits_zero(q * q, 2));
    }

    return sum_squarefree_products_leq(good, n);
}

void run_validations() {
    assert(solve(20) == 18);
    assert(solve(1000) == 2089);
}

}  // namespace

int main() {
    run_validations();
    constexpr i64 kN = 10'000'000;
    std::cout << solve(kN) << '\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.*;

public class Euler927 {

    static long modPow(long base, int exp, int mod) {
        if (mod == 1)
            return 0;
        long res = 1 % mod;
        long b = base % mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * b) % mod;
            b = (b * b) % mod;
            exp >>= 1;
        }
        return res;
    }

    static long nextTerm(long x, int p, int mod) {
        if (p == 2) {
            return (x * x + 1L) % mod;
        }
        return (modPow(x, p, mod) + 1L) % mod;
    }

    static boolean hitsZero(int mod, int p) {
        if (mod == 1)
            return true;
        long start = 1L % mod;
        if (start == 0)
            return true;

        long tortoise = start;
        long hare = nextTerm(start, p, mod);
        if (hare == 0)
            return true;

        long power = 1;
        long lam = 1;
        while (tortoise != hare) {
            if (power == lam) {
                tortoise = hare;
                power <<= 1;
                lam = 0;
            }
            hare = nextTerm(hare, p, mod);
            if (hare == 0)
                return true;
            lam++;
        }
        return false;
    }

    static class Sieve {
        List<Integer> primes = new ArrayList<>();
        int[] spf;
    }

    static Sieve buildSieve(int n) {
        Sieve sieve = new Sieve();
        sieve.spf = new int[n + 1];
        for (int i = 2; i <= n; i++) {
            if (sieve.spf[i] == 0) {
                sieve.spf[i] = i;
                sieve.primes.add(i);
            }
            for (int p : sieve.primes) {
                long v = (long) p * i;
                if (v > n || p > sieve.spf[i])
                    break;
                sieve.spf[(int) v] = p;
            }
        }
        return sieve;
    }

    static boolean isGoodPrime(int q, int[] spf) {
        if (q == 2)
            return true;

        int x = q - 1;
        List<Integer> factors = new ArrayList<>();
        while (x > 1) {
            int p = spf[x];
            factors.add(p);
            while (x % p == 0)
                x /= p;
        }

        Collections.sort(factors, Collections.reverseOrder());
        // Deduplicate
        List<Integer> uniqueFactors = new ArrayList<>();
        int last = -1;
        for (int p : factors) {
            if (p != last) {
                uniqueFactors.add(p);
                last = p;
            }
        }

        for (int p : uniqueFactors) {
            if (!hitsZero(q, p))
                return false;
        }
        return true;
    }

    static long total = 0;

    static void dfs(int idx, long current, List<Integer> good, long limit) {
        total += current;
        for (int i = idx; i < good.size(); ++i) {
            long p = good.get(i);
            if (current > limit / p)
                break;
            dfs(i + 1, current * p, good, limit);
        }
    }

    public static String solve(long n) {
        Sieve sieve = buildSieve((int) n);
        List<Integer> good = new ArrayList<>();
        for (int q : sieve.primes) {
            if (isGoodPrime(q, sieve.spf)) {
                good.add(q);
            }
        }

        total = 0;
        dfs(0, 1, good, n);
        return Long.toString(total);
    }

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