Problem 362: Squarefree Factors
View on Project EulerProject Euler Problem 362 Solution
EulerSolve provides an optimized solution for Project Euler Problem 362, Squarefree Factors, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(n \ge 2\), let \(\operatorname{fsf}(n)\) denote the number of unordered factorizations of \(n\) into squarefree factors greater than 1. Repeated factors are allowed as long as each factor itself is squarefree. The Project Euler quantity is $$S(N)=\sum_{n=2}^{N}\operatorname{fsf}(n),\qquad N=10^{10}.$$ As a small example, \(12=2^2\cdot 3\) has exactly two valid factorizations: $$12=2\cdot 2\cdot 3=2\cdot 6,$$ so \(\operatorname{fsf}(12)=2\). A direct scan up to \(10^{10}\) is impossible, so the implementation groups integers by their exponent pattern and counts entire families at once. Mathematical Approach Write the prime factorization of \(n\) as $$n=\prod_{i=1}^{r} p_i^{e_i},\qquad e_i\ge 1.$$ The multiset of exponents \((e_1,\dots,e_r)\) is the only part that matters for the squarefree-factor count. Step 1: Replace Unordered Factorizations by a System on Subsets Every squarefree factor of \(n\) is obtained by choosing a nonempty subset \(J\subseteq\{1,\dots,r\}\) and taking $$f_J=\prod_{i\in J} p_i.$$ Because the factorization is unordered, all we need is the multiplicity \(x_J\in\mathbb{Z}_{\ge 0}\) of each subset type....
Detailed mathematical approach
Problem Summary
For each integer \(n \ge 2\), let \(\operatorname{fsf}(n)\) denote the number of unordered factorizations of \(n\) into squarefree factors greater than 1. Repeated factors are allowed as long as each factor itself is squarefree. The Project Euler quantity is
$$S(N)=\sum_{n=2}^{N}\operatorname{fsf}(n),\qquad N=10^{10}.$$
As a small example, \(12=2^2\cdot 3\) has exactly two valid factorizations:
$$12=2\cdot 2\cdot 3=2\cdot 6,$$
so \(\operatorname{fsf}(12)=2\). A direct scan up to \(10^{10}\) is impossible, so the implementation groups integers by their exponent pattern and counts entire families at once.
Mathematical Approach
Write the prime factorization of \(n\) as
$$n=\prod_{i=1}^{r} p_i^{e_i},\qquad e_i\ge 1.$$
The multiset of exponents \((e_1,\dots,e_r)\) is the only part that matters for the squarefree-factor count.
Step 1: Replace Unordered Factorizations by a System on Subsets
Every squarefree factor of \(n\) is obtained by choosing a nonempty subset \(J\subseteq\{1,\dots,r\}\) and taking
$$f_J=\prod_{i\in J} p_i.$$
Because the factorization is unordered, all we need is the multiplicity \(x_J\in\mathbb{Z}_{\ge 0}\) of each subset type. Prime \(p_i\) must appear exactly \(e_i\) times across all chosen factors, hence
$$\sum_{J\ni i} x_J=e_i\qquad (1\le i\le r).$$
So \(\operatorname{fsf}(n)\) is precisely the number of nonnegative integer solutions of this linear system. In particular, it depends only on the exponent signature, not on the actual prime values. If \(\sigma=(e_1,\dots,e_r)\), we may therefore write \(F(\sigma)\) for this common value.
Step 2: Count \(F(\sigma)\) with Memoized DFS
The C++ and Java implementations build all nonempty bitmasks \(J\), sort them by decreasing size, and recurse on a vector of remaining exponents. At one mask \(J\), the multiplicity \(x_J\) can be any value from \(0\) up to
$$\min_{i\in J} \operatorname{remaining}_i.$$
For each choice, the code subtracts that amount from every coordinate in \(J\) and continues with the next mask. The base case returns 1 if all remaining exponents are zero and 0 otherwise. Memoization uses the pair
$$\bigl(\text{mask position},\text{remaining exponent vector}\bigr)$$
as the state key. This is exactly what fsf_from_signature and its DFS helper implement.
Step 3: Worked Signature Example
For \(12=2^2\cdot 3\), the signature is \((2,1)\). There are three subset types: \(\{1\}\), \(\{2\}\), and \(\{1,2\}\). Writing their multiplicities as \(x_1,x_2,x_{12}\), the constraints become
$$x_1+x_{12}=2,\qquad x_2+x_{12}=1,\qquad x_1,x_2,x_{12}\ge 0.$$
The two solutions are
$$ (x_1,x_2,x_{12})=(2,1,0)\quad\text{and}\quad (1,0,1), $$
corresponding to \(2\cdot 2\cdot 3\) and \(2\cdot 6\). Hence \(F(2,1)=2\). Every integer with exponent multiset \(\{2,1\}\) has the same squarefree-factor count.
Step 4: Enumerate Only Feasible Exponent Signatures
To sum over all \(n\le N\), the program does not iterate over integers. Instead it enumerates exponent signatures \(\sigma=(e_1,\dots,e_r)\) in nonincreasing order
$$e_1\ge e_2\ge \dots\ge e_r\ge 1.$$
For such a signature, the smallest integer realizing it is obtained by pairing the largest exponent with the smallest prime:
$$m_{\min}(\sigma)=2^{e_1}3^{e_2}5^{e_3}\cdots p_r^{e_r}.$$
If \(m_{\min}(\sigma)>N\), then no integer up to \(N\) can have that signature. If \(m_{\min}(\sigma)\le N\), then the signature is feasible. This is exactly why generate_signatures_dfs grows signatures along the seed primes \(2,3,5,\dots\) while forcing exponents to stay nonincreasing.
Step 5: Count Integers for One Ordered Exponent Tuple
A signature is an unordered multiset of exponents, but actual integers assign those exponents to increasing primes. For one ordered tuple \(a=(a_1,\dots,a_r)\), define
$$C(a,N)=\left|\left\{(q_1,\dots,q_r): q_1\lt q_2\lt \dots\lt q_r,\ \prod_{i=1}^{r} q_i^{a_i}\le N\right\}\right|.$$
The recursion in OrderedPrimeTupleCounter chooses the primes left to right. If
$$T_k=a_k+a_{k+1}+\dots+a_r$$
is the remaining exponent sum at depth \(k\), then the current prime must satisfy
$$q_k\le \left\lfloor R^{1/T_k}\right\rfloor,$$
where \(R\) is the remaining limit after previous prime choices. This bound is sharp because every later prime is at least \(q_k\). On the last position there is no deeper recursion; the number of valid choices is obtained directly from the prime-counting function:
$$\pi\!\left(\left\lfloor R^{1/a_r}\right\rfloor\right)-\text{startIndex}.$$
If some exponents repeat, the same ordered tuple would arise many times, so the code generates only unique permutations of the signature.
Step 6: A Concrete Family Below 100
The signature \((2,1)\) illustrates the split between \(F(\sigma)\) and \(C(a,N)\). We already know \(F(2,1)=2\). For \(N=100\), the ordered tuple \((2,1)\) counts integers of the form \(p^2q\) with \(p\lt q\), while \((1,2)\) counts integers of the form \(pq^2\).
For \((2,1)\): \(2^2q\le 100\) gives 8 possibilities for \(q\), and \(3^2q\le 100\) gives 3 more, so
$$C((2,1),100)=11.$$
For \((1,2)\): \(2q^2\le 100\) gives 3 possibilities and \(3q^2\le 100\) gives 1 more, so
$$C((1,2),100)=4.$$
Therefore this entire signature family contributes
$$F(2,1)\bigl(C((2,1),100)+C((1,2),100)\bigr)=2(11+4)=30$$
to \(S(100)\). This is exactly the logic used for every signature.
Step 7: Prime Counting and the Final Sum
The last remaining ingredient is a fast \(\pi(x)\). The implementation builds a sieve up to \(\sqrt{N}\), stores the prime list and the small values of \(\pi(x)\), and for larger arguments uses a Lehmer-style prime counter with cached
$$\phi(x,s)=\left|\left\{1\le m\le x:\ p_1,\dots,p_s\nmid m\right\}\right|.$$
This makes the many root-bound queries inside the tuple recursion fast enough. If \(\Sigma(N)\) is the set of feasible signatures and \(\mathcal{U}(\sigma)\) is the set of unique permutations of \(\sigma\), then the program computes
$$\boxed{S(N)=\sum_{\sigma\in\Sigma(N)} F(\sigma)\sum_{a\in\mathcal{U}(\sigma)} C(a,N).}$$
How the Code Works
The C++ file contains the full optimized solver. PrimeTable builds the sieve and small \(\pi(x)\) table, PrimeCounter provides the Lehmer-style prime-counting routine, fsf_from_signature counts squarefree factorizations for one signature, generate_signatures enumerates feasible exponent patterns, generate_unique_permutations expands repeated exponents into distinct orderings, and OrderedPrimeTupleCounter::count_for_exponents counts the prime tuples for one ordered exponent vector.
The Java solution mirrors the same mathematics and almost the same function decomposition in a single-threaded implementation. The Python file is intentionally thin: it recompiles and runs the C++ source when necessary, then parses the final numeric output. In other words, the mathematical source of truth is the shared C++/Java algorithm, while Python is just a bridge.
The C++ program also validates the logic with internal checkpoints: it verifies the known value \(S(100)=193\), compares fast and brute-force computations on a small limit, and checks that single-threaded and multi-threaded runs agree.
Complexity Analysis
The sieve and small prime table up to \(\sqrt{N}\) cost \(O(\sqrt{N}\log\log \sqrt{N})\) time and \(O(\sqrt{N})\) memory. For \(N=10^{10}\), that preprocessing size is only \(10^5\).
The rest of the runtime is dominated by the signature tasks. There is no simple single closed form, because the DFS state counts depend on the exponent pattern and on the distribution of root-bound prime queries. However, the key compression is dramatic: the algorithm works on signatures, permutations, and memoized prime-tuple states rather than on all integers up to \(N\).
In this problem, the number of distinct prime factors is automatically small: the product of the first ten primes is \(6469693230\), while including the eleventh prime already exceeds \(10^{10}\). Hence every signature has length at most 10, which keeps the subset DFS, the permutation generation, and the tuple recursion practical. Memory is dominated by the memo tables for the factorization DFS and the prime-counting cache.
Footnotes and References
- Problem page: https://projecteuler.net/problem=362
- Squarefree integers: Wikipedia — Square-free integer
- Prime-counting function: Wikipedia — Prime-counting function
- Meissel-Lehmer prime counting: Wikipedia — Meissel-Lehmer algorithm
Problem 362 source code
C++
#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u16 = std::uint16_t;
using u8 = std::uint8_t;
using u128 = unsigned __int128;
constexpr u64 kDefaultLimit = 10'000'000'000ULL;
constexpr u64 kCheckpointKnownLimit = 100ULL;
constexpr u64 kCheckpointKnownExpected = 193ULL;
constexpr u64 kCheckpointBruteLimit = 500ULL;
constexpr u64 kThreadConsistencyLimit = 200'000ULL;
struct Options {
u64 limit = kDefaultLimit;
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0U;
};
bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const u64 digit = static_cast<u64>(c - '0');
if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
return false;
}
parsed = parsed * 10ULL + digit;
}
value = parsed;
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u64 parsed = 0ULL;
if (!parse_u64_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
u64 parsed_u64 = 0ULL;
if (parse_u64_after_prefix(arg, "--limit=", parsed_u64)) {
options.limit = parsed_u64;
continue;
}
unsigned parsed_unsigned = 0U;
if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
options.requested_threads = parsed_unsigned;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.limit < 2ULL) {
std::cerr << "--limit must be at least 2.\n";
return false;
}
return true;
}
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload) {
if (!allow_multithreading || workload < 2ULL) {
return 1U;
}
unsigned threads = requested_threads;
if (threads == 0U) {
threads = std::thread::hardware_concurrency();
if (threads == 0U) {
threads = 1U;
}
}
return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}
u64 isqrt_u64(u64 x) {
u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
while ((r + 1ULL) <= x / (r + 1ULL)) {
++r;
}
while (r > x / r) {
--r;
}
return r;
}
u64 pow_u64_capped(u64 base, int exp, u64 cap) {
u64 result = 1ULL;
for (int i = 0; i < exp; ++i) {
if (base != 0ULL && result > cap / base) {
return cap + 1ULL;
}
result *= base;
}
return result;
}
u64 kth_root_floor(u64 n, int k) {
if (k <= 1 || n <= 1ULL) {
return n;
}
const long double approx =
std::pow(static_cast<long double>(n), 1.0L / static_cast<long double>(k));
u64 r = static_cast<u64>(approx);
if (r == 0ULL) {
r = 1ULL;
}
while (pow_u64_capped(r + 1ULL, k, n) <= n) {
++r;
}
while (pow_u64_capped(r, k, n) > n) {
--r;
}
return r;
}
std::string to_string_u128(u128 value) {
if (value == 0) {
return "0";
}
std::string s;
while (value > 0) {
const u32 digit = static_cast<u32>(value % 10U);
s.push_back(static_cast<char>('0' + digit));
value /= 10U;
}
std::reverse(s.begin(), s.end());
return s;
}
struct PrimeTable {
int limit = 0;
std::vector<int> primes;
std::vector<int> pi_small;
explicit PrimeTable(int sieve_limit) { build(sieve_limit); }
void build(int sieve_limit) {
limit = sieve_limit;
std::vector<int> lp(static_cast<std::size_t>(limit) + 1ULL, 0);
primes.reserve(static_cast<std::size_t>(limit) / 10ULL + 32ULL);
for (int i = 2; i <= limit; ++i) {
if (lp[static_cast<std::size_t>(i)] == 0) {
lp[static_cast<std::size_t>(i)] = i;
primes.push_back(i);
}
for (const int p : primes) {
const long long v = 1LL * p * i;
if (p > lp[static_cast<std::size_t>(i)] || v > limit) {
break;
}
lp[static_cast<std::size_t>(v)] = p;
}
}
pi_small.assign(static_cast<std::size_t>(limit) + 1ULL, 0);
int ptr = 0;
for (int i = 1; i <= limit; ++i) {
pi_small[static_cast<std::size_t>(i)] =
pi_small[static_cast<std::size_t>(i - 1)];
if (ptr < static_cast<int>(primes.size()) && primes[static_cast<std::size_t>(ptr)] == i) {
++pi_small[static_cast<std::size_t>(i)];
++ptr;
}
}
}
};
u64 icbrt_u64(u64 x) {
u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
while (static_cast<u128>(r + 1ULL) * (r + 1ULL) * (r + 1ULL) <= x) {
++r;
}
while (static_cast<u128>(r) * r * r > x) {
--r;
}
return r;
}
u64 iroot4_u64(u64 x) {
u64 r = static_cast<u64>(std::sqrt(std::sqrt(static_cast<long double>(x))));
while (static_cast<u128>(r + 1ULL) * (r + 1ULL) * (r + 1ULL) * (r + 1ULL) <= x) {
++r;
}
while (static_cast<u128>(r) * r * r * r > x) {
--r;
}
return r;
}
class PrimeCounter {
public:
explicit PrimeCounter(const PrimeTable& table) : table_(table) { build_phi_table(); }
u64 pi(u64 x) const {
if (x <= static_cast<u64>(table_.limit)) {
return static_cast<u64>(table_.pi_small[static_cast<std::size_t>(x)]);
}
auto& cache = pi_cache();
const auto it = cache.find(x);
if (it != cache.end()) {
return it->second;
}
const u64 a = pi(iroot4_u64(x));
const u64 b = pi(isqrt_u64(x));
const u64 c = pi(icbrt_u64(x));
i64 sum = static_cast<i64>(phi(x, static_cast<int>(a)));
sum += static_cast<i64>((b + a - 2ULL) * (b - a + 1ULL) / 2ULL);
for (u64 i = a + 1ULL; i <= b; ++i) {
const u64 w = x / static_cast<u64>(table_.primes[static_cast<std::size_t>(i - 1ULL)]);
sum -= static_cast<i64>(pi(w));
if (i <= c) {
const u64 bi = pi(isqrt_u64(w));
for (u64 j = i; j <= bi; ++j) {
const u64 q =
w / static_cast<u64>(table_.primes[static_cast<std::size_t>(j - 1ULL)]);
sum -= static_cast<i64>(pi(q)) - static_cast<i64>(j - 1ULL);
}
}
}
const u64 ans = static_cast<u64>(sum);
cache.emplace(x, ans);
return ans;
}
private:
using i64 = std::int64_t;
static constexpr int kSmallPhiPrimes = 7; // 2,3,5,7,11,13,17
static constexpr int kSmallPhiMod = 510510;
static constexpr int kPhiCacheSMax = 128;
static constexpr u64 kPhiCacheXMax = 50'000'000'000ULL;
const PrimeTable& table_;
std::vector<std::array<u32, kSmallPhiPrimes + 1>> phi_table_;
void build_phi_table() {
phi_table_.assign(static_cast<std::size_t>(kSmallPhiMod) + 1ULL, {});
for (int n = 0; n <= kSmallPhiMod; ++n) {
phi_table_[static_cast<std::size_t>(n)][0] = static_cast<u32>(n);
}
for (int s = 1; s <= kSmallPhiPrimes; ++s) {
const int p = table_.primes[static_cast<std::size_t>(s - 1)];
for (int n = 0; n <= kSmallPhiMod; ++n) {
phi_table_[static_cast<std::size_t>(n)][static_cast<std::size_t>(s)] =
phi_table_[static_cast<std::size_t>(n)][static_cast<std::size_t>(s - 1)] -
phi_table_[static_cast<std::size_t>(n / p)][static_cast<std::size_t>(s - 1)];
}
}
}
static std::unordered_map<u64, u64>& pi_cache() {
static thread_local std::unordered_map<u64, u64> cache;
return cache;
}
static std::unordered_map<u64, u64>& phi_cache() {
static thread_local std::unordered_map<u64, u64> cache;
return cache;
}
u64 phi(u64 x, int s) const {
if (s == 0) {
return x;
}
if (s <= kSmallPhiPrimes) {
const u64 block = x / static_cast<u64>(kSmallPhiMod);
const u64 rem = x % static_cast<u64>(kSmallPhiMod);
return block * static_cast<u64>(phi_table_[kSmallPhiMod][static_cast<std::size_t>(s)]) +
static_cast<u64>(
phi_table_[static_cast<std::size_t>(rem)][static_cast<std::size_t>(s)]);
}
if (x <= static_cast<u64>(table_.limit) &&
static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]) *
static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]) >
x) {
return static_cast<u64>(table_.pi_small[static_cast<std::size_t>(x)] - s + 1);
}
const u64 ps = static_cast<u64>(table_.primes[static_cast<std::size_t>(s - 1)]);
if (ps * ps > x) {
return pi(x) - static_cast<u64>(s) + 1ULL;
}
if (s <= kPhiCacheSMax && x <= kPhiCacheXMax) {
const u64 key = (x << 8U) ^ static_cast<u64>(s);
auto& cache = phi_cache();
const auto it = cache.find(key);
if (it != cache.end()) {
return it->second;
}
const u64 ans = phi(x, s - 1) - phi(x / ps, s - 1);
cache.emplace(key, ans);
return ans;
}
return phi(x, s - 1) - phi(x / ps, s - 1);
}
};
struct FsfMemoKey {
u16 mask_pos = 0U;
u64 packed_state = 0ULL;
bool operator==(const FsfMemoKey& other) const {
return mask_pos == other.mask_pos && packed_state == other.packed_state;
}
};
struct FsfMemoKeyHash {
std::size_t operator()(const FsfMemoKey& key) const {
const u64 mixed = key.packed_state ^ (static_cast<u64>(key.mask_pos) * 0x9E3779B97F4A7C15ULL);
return static_cast<std::size_t>(mixed ^ (mixed >> 33U));
}
};
u64 pack_small_vector(const std::array<u8, 10>& values, int size) {
u64 packed = 0ULL;
for (int i = 0; i < size; ++i) {
packed |= static_cast<u64>(values[static_cast<std::size_t>(i)]) << (6U * i);
}
return packed;
}
u64 fsf_from_signature(const std::vector<u8>& signature) {
const int r = static_cast<int>(signature.size());
if (r == 0) {
return 0ULL;
}
std::vector<int> masks;
masks.reserve((1 << r) - 1);
for (int mask = 1; mask < (1 << r); ++mask) {
masks.push_back(mask);
}
std::sort(masks.begin(), masks.end(), [](int a, int b) {
const int pa = __builtin_popcount(static_cast<unsigned>(a));
const int pb = __builtin_popcount(static_cast<unsigned>(b));
if (pa != pb) {
return pa > pb;
}
return a < b;
});
std::vector<std::vector<u8>> mask_bits;
mask_bits.reserve(masks.size());
for (const int mask : masks) {
std::vector<u8> bits;
for (int i = 0; i < r; ++i) {
if ((mask >> i) & 1) {
bits.push_back(static_cast<u8>(i));
}
}
mask_bits.push_back(std::move(bits));
}
std::array<u8, 10> remaining{};
for (int i = 0; i < r; ++i) {
remaining[static_cast<std::size_t>(i)] = signature[static_cast<std::size_t>(i)];
}
std::unordered_map<FsfMemoKey, u64, FsfMemoKeyHash> memo;
memo.reserve(4096);
std::function<u64(u16)> dfs = [&](u16 pos) -> u64 {
if (pos == static_cast<u16>(mask_bits.size())) {
for (int i = 0; i < r; ++i) {
if (remaining[static_cast<std::size_t>(i)] != 0U) {
return 0ULL;
}
}
return 1ULL;
}
const FsfMemoKey key{pos, pack_small_vector(remaining, r)};
const auto it = memo.find(key);
if (it != memo.end()) {
return it->second;
}
const std::vector<u8>& bits = mask_bits[static_cast<std::size_t>(pos)];
u8 max_take = std::numeric_limits<u8>::max();
for (const u8 bit : bits) {
max_take = std::min(max_take, remaining[static_cast<std::size_t>(bit)]);
}
u64 total = 0ULL;
for (u32 take = 0U;; ++take) {
total += dfs(static_cast<u16>(pos + 1U));
if (take == max_take) {
break;
}
for (const u8 bit : bits) {
--remaining[static_cast<std::size_t>(bit)];
}
}
for (const u8 bit : bits) {
remaining[static_cast<std::size_t>(bit)] =
static_cast<u8>(remaining[static_cast<std::size_t>(bit)] + max_take);
}
memo.emplace(key, total);
return total;
};
return dfs(0U);
}
void generate_signatures_dfs(const std::vector<int>& seed_primes,
int prime_index,
int last_exp,
u64 current,
u64 limit,
std::vector<u8>& cur,
std::vector<std::vector<u8>>& out) {
if (prime_index >= static_cast<int>(seed_primes.size())) {
return;
}
const u64 p = static_cast<u64>(seed_primes[static_cast<std::size_t>(prime_index)]);
u64 p_pow = 1ULL;
for (int e = 1; e <= last_exp; ++e) {
if (p_pow > limit / p) {
break;
}
p_pow *= p;
if (current > limit / p_pow) {
break;
}
const u64 next_value = current * p_pow;
cur.push_back(static_cast<u8>(e));
out.push_back(cur);
generate_signatures_dfs(seed_primes, prime_index + 1, e, next_value, limit, cur, out);
cur.pop_back();
}
}
std::vector<std::vector<u8>> generate_signatures(u64 limit, const PrimeTable& table) {
std::vector<std::vector<u8>> signatures;
std::vector<u8> cur;
generate_signatures_dfs(table.primes, 0, 63, 1ULL, limit, cur, signatures);
return signatures;
}
void generate_unique_permutations_rec(std::array<int, 64>& freq,
std::vector<u8>& cur,
std::size_t pos,
std::vector<std::vector<u8>>& out) {
if (pos == cur.size()) {
out.push_back(cur);
return;
}
for (int value = 1; value < static_cast<int>(freq.size()); ++value) {
if (freq[static_cast<std::size_t>(value)] == 0) {
continue;
}
--freq[static_cast<std::size_t>(value)];
cur[pos] = static_cast<u8>(value);
generate_unique_permutations_rec(freq, cur, pos + 1ULL, out);
++freq[static_cast<std::size_t>(value)];
}
}
std::vector<std::vector<u8>> generate_unique_permutations(const std::vector<u8>& signature) {
std::array<int, 64> freq{};
for (const u8 e : signature) {
++freq[static_cast<std::size_t>(e)];
}
std::vector<std::vector<u8>> out;
std::vector<u8> cur(signature.size(), 0U);
generate_unique_permutations_rec(freq, cur, 0ULL, out);
return out;
}
struct CountMemoKey {
u64 limit = 0ULL;
u32 start_index = 0U;
u8 pos = 0U;
bool operator==(const CountMemoKey& other) const {
return limit == other.limit && start_index == other.start_index && pos == other.pos;
}
};
struct CountMemoKeyHash {
std::size_t operator()(const CountMemoKey& key) const {
u64 h = key.limit;
h ^= static_cast<u64>(key.start_index) * 0x9E3779B97F4A7C15ULL;
h ^= static_cast<u64>(key.pos) * 0xC2B2AE3D27D4EB4FULL;
h ^= (h >> 31U);
return static_cast<std::size_t>(h);
}
};
class OrderedPrimeTupleCounter {
public:
OrderedPrimeTupleCounter(const PrimeTable& table, const PrimeCounter& prime_counter)
: table_(table), prime_counter_(prime_counter) {}
u64 count_for_exponents(const std::vector<u8>& exponents, u64 limit) const {
const int r = static_cast<int>(exponents.size());
if (r == 0) {
return 0ULL;
}
if (r == 1) {
const u64 root = kth_root_floor(limit, static_cast<int>(exponents[0]));
return prime_counter_.pi(root);
}
std::array<int, 10> tail_sum{};
int running = 0;
for (int i = r - 1; i >= 0; --i) {
running += static_cast<int>(exponents[static_cast<std::size_t>(i)]);
tail_sum[static_cast<std::size_t>(i)] = running;
}
std::unordered_map<CountMemoKey, u64, CountMemoKeyHash> memo;
memo.reserve(4096);
std::function<u64(int, u32, u64)> dfs = [&](int pos, u32 start_index, u64 rem_limit) -> u64 {
if (pos == r) {
return 1ULL;
}
const CountMemoKey key{rem_limit, start_index, static_cast<u8>(pos)};
const auto it = memo.find(key);
if (it != memo.end()) {
return it->second;
}
u64 result = 0ULL;
if (pos == r - 1) {
const int exp = static_cast<int>(exponents[static_cast<std::size_t>(pos)]);
const u64 root = kth_root_floor(rem_limit, exp);
const u64 upto = prime_counter_.pi(root);
result = (upto > static_cast<u64>(start_index))
? (upto - static_cast<u64>(start_index))
: 0ULL;
memo.emplace(key, result);
return result;
}
const int tail = tail_sum[static_cast<std::size_t>(pos)];
const u64 max_p = kth_root_floor(rem_limit, tail);
u64 upto = prime_counter_.pi(max_p);
if (upto > table_.primes.size()) {
upto = static_cast<u64>(table_.primes.size());
}
const int exp = static_cast<int>(exponents[static_cast<std::size_t>(pos)]);
for (u64 idx = static_cast<u64>(start_index); idx < upto; ++idx) {
const u64 p = static_cast<u64>(table_.primes[static_cast<std::size_t>(idx)]);
const u64 p_pow = pow_u64_capped(p, exp, rem_limit);
if (p_pow > rem_limit) {
break;
}
result += dfs(pos + 1, static_cast<u32>(idx + 1ULL), rem_limit / p_pow);
}
memo.emplace(key, result);
return result;
};
return dfs(0, 0U, limit);
}
private:
const PrimeTable& table_;
const PrimeCounter& prime_counter_;
};
struct Task {
std::vector<u8> exponents;
u64 fsf_value = 0ULL;
};
u128 solve_fast(u64 limit, bool allow_multithreading, unsigned requested_threads) {
const u64 root_limit = std::max<u64>(isqrt_u64(limit), 100ULL);
if (root_limit > static_cast<u64>(std::numeric_limits<int>::max())) {
throw std::runtime_error("Sieve limit too large for this implementation.");
}
const PrimeTable table(static_cast<int>(root_limit));
const PrimeCounter prime_counter(table);
const OrderedPrimeTupleCounter tuple_counter(table, prime_counter);
const std::vector<std::vector<u8>> signatures = generate_signatures(limit, table);
std::vector<Task> tasks;
tasks.reserve(16'000ULL);
for (const auto& signature : signatures) {
const u64 fsf = fsf_from_signature(signature);
const std::vector<std::vector<u8>> permutations = generate_unique_permutations(signature);
for (const auto& p : permutations) {
tasks.push_back(Task{p, fsf});
}
}
const unsigned threads =
choose_thread_count(allow_multithreading, requested_threads, tasks.size());
if (threads == 1U) {
u128 total = 0;
for (const Task& task : tasks) {
const u64 count = tuple_counter.count_for_exponents(task.exponents, limit);
total += static_cast<u128>(task.fsf_value) * static_cast<u128>(count);
}
return total;
}
std::atomic<std::size_t> next_task{0ULL};
std::vector<u128> partials(threads, 0);
std::vector<std::thread> workers;
workers.reserve(threads);
for (unsigned tid = 0; tid < threads; ++tid) {
workers.emplace_back([&, tid]() {
const OrderedPrimeTupleCounter local_counter(table, prime_counter);
u128 local = 0;
while (true) {
const std::size_t idx = next_task.fetch_add(1ULL, std::memory_order_relaxed);
if (idx >= tasks.size()) {
break;
}
const Task& task = tasks[idx];
const u64 count = local_counter.count_for_exponents(task.exponents, limit);
local += static_cast<u128>(task.fsf_value) * static_cast<u128>(count);
}
partials[static_cast<std::size_t>(tid)] = local;
});
}
for (auto& th : workers) {
th.join();
}
u128 total = 0;
for (const u128 value : partials) {
total += value;
}
return total;
}
std::vector<int> squarefree_numbers(int limit) {
std::vector<int> sqf;
std::vector<int> square_div(static_cast<std::size_t>(limit) + 1ULL, 0);
for (int p = 2; p * p <= limit; ++p) {
const int pp = p * p;
for (int n = pp; n <= limit; n += pp) {
square_div[static_cast<std::size_t>(n)] = 1;
}
}
for (int n = 2; n <= limit; ++n) {
if (square_div[static_cast<std::size_t>(n)] == 0) {
sqf.push_back(n);
}
}
return sqf;
}
u64 brute_fsf_single(u64 n,
std::size_t start_idx,
const std::vector<int>& sqf,
std::unordered_map<u64, u64>& memo) {
const u64 key = (n << 16U) ^ static_cast<u64>(start_idx);
const auto it = memo.find(key);
if (it != memo.end()) {
return it->second;
}
u64 total = 0ULL;
for (std::size_t i = start_idx; i < sqf.size(); ++i) {
const u64 d = static_cast<u64>(sqf[i]);
if (d > n) {
break;
}
if (n % d != 0ULL) {
continue;
}
if (d == n) {
++total;
} else {
total += brute_fsf_single(n / d, i, sqf, memo);
}
}
memo.emplace(key, total);
return total;
}
u64 brute_s(u64 limit) {
if (limit < 2ULL) {
return 0ULL;
}
if (limit > 100'000ULL) {
throw std::runtime_error("brute_s is intended only for small limits.");
}
const std::vector<int> sqf = squarefree_numbers(static_cast<int>(limit));
u64 total = 0ULL;
for (u64 n = 2ULL; n <= limit; ++n) {
std::unordered_map<u64, u64> memo;
total += brute_fsf_single(n, 0ULL, sqf, memo);
}
return total;
}
void run_checkpoints() {
const u64 brute_100 = brute_s(kCheckpointKnownLimit);
if (brute_100 != kCheckpointKnownExpected) {
throw std::runtime_error("Brute-force checkpoint S(100)=193 failed.");
}
const u64 fast_100 = static_cast<u64>(solve_fast(kCheckpointKnownLimit, false, 1U));
if (fast_100 != kCheckpointKnownExpected) {
throw std::runtime_error("Fast checkpoint S(100)=193 failed.");
}
const u64 brute_small = brute_s(kCheckpointBruteLimit);
const u64 fast_small = static_cast<u64>(solve_fast(kCheckpointBruteLimit, false, 1U));
if (fast_small != brute_small) {
throw std::runtime_error("Brute-force consistency checkpoint failed.");
}
const u64 single = static_cast<u64>(solve_fast(kThreadConsistencyLimit, false, 1U));
const u64 threaded = static_cast<u64>(solve_fast(kThreadConsistencyLimit, true, 2U));
if (single != threaded) {
throw std::runtime_error("Thread consistency checkpoint failed.");
}
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
try {
if (options.run_checkpoints) {
run_checkpoints();
}
const auto start = std::chrono::steady_clock::now();
const u128 answer = solve_fast(
options.limit, options.allow_multithreading, options.requested_threads);
const auto end = std::chrono::steady_clock::now();
const std::chrono::duration<double> elapsed = end - start;
std::cout << "S(" << options.limit << ") = " << to_string_u128(answer) << '\n';
std::cout << "Elapsed: " << elapsed.count() << " s\n";
} catch (const std::exception& ex) {
std::cerr << "Error: " << ex.what() << '\n';
return 1;
}
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 Euler362 {
static long isqrt(long n) {
if (n < 0)
return 0;
long r = (long) Math.sqrt(n);
while ((r + 1) <= n / (r + 1))
r++;
while (r > n / r)
r--;
return r;
}
static long icbrt(long n) {
if (n < 0)
return 0;
long r = (long) Math.cbrt(n);
while ((r + 1) * (r + 1) * (r + 1) <= n)
r++;
while (r * r * r > n)
r--;
return r;
}
static long iroot4(long n) {
if (n < 0)
return 0;
long r = (long) Math.sqrt(Math.sqrt(n));
while ((r + 1) * (r + 1) * (r + 1) * (r + 1) <= n)
r++;
while (r * r * r * r > n)
r--;
return r;
}
static long powCapped(long base, int exp, long cap) {
long res = 1;
for (int i = 0; i < exp; i++) {
if (base != 0 && res > cap / base)
return cap + 1;
res *= base;
}
return res;
}
static long kthRootFloor(long n, int k) {
if (k <= 1 || n <= 1)
return n;
long r = (long) Math.pow(n, 1.0 / k);
if (r == 0)
r = 1;
while (powCapped(r + 1, k, n) <= n)
r++;
while (powCapped(r, k, n) > n)
r--;
return r;
}
static class PrimeTable {
int limit;
int[] primes;
int[] piSmall;
PrimeTable(int limit) {
this.limit = limit;
int[] lp = new int[limit + 1];
List<Integer> prs = new ArrayList<>();
for (int i = 2; i <= limit; i++) {
if (lp[i] == 0) {
lp[i] = i;
prs.add(i);
}
for (int p : prs) {
long v = (long) p * i;
if (p > lp[i] || v > limit)
break;
lp[(int) v] = p;
}
}
primes = new int[prs.size()];
for (int i = 0; i < prs.size(); i++)
primes[i] = prs.get(i);
piSmall = new int[limit + 1];
int ptr = 0;
for (int i = 1; i <= limit; i++) {
piSmall[i] = piSmall[i - 1];
if (ptr < primes.length && primes[ptr] == i) {
piSmall[i]++;
ptr++;
}
}
}
}
static class PrimeCounter {
PrimeTable table;
Map<Long, Long> piCache = new HashMap<>();
Map<Long, Long> phiCache = new HashMap<>();
static final int K_SMALL_PHI_PRIMES = 7;
static final int K_SMALL_PHI_MOD = 510510;
static final int K_PHI_CACHE_S_MAX = 128;
static final long K_PHI_CACHE_X_MAX = 50000000000L;
int[][] phiTable;
PrimeCounter(PrimeTable table) {
this.table = table;
phiTable = new int[K_SMALL_PHI_MOD + 1][K_SMALL_PHI_PRIMES + 1];
for (int n = 0; n <= K_SMALL_PHI_MOD; n++)
phiTable[n][0] = n;
for (int s = 1; s <= K_SMALL_PHI_PRIMES; s++) {
int p = table.primes[s - 1];
for (int n = 0; n <= K_SMALL_PHI_MOD; n++) {
phiTable[n][s] = phiTable[n][s - 1] - phiTable[n / p][s - 1];
}
}
}
long phi(long x, int s) {
if (s == 0)
return x;
if (s <= K_SMALL_PHI_PRIMES) {
long block = x / K_SMALL_PHI_MOD;
long rem = x % K_SMALL_PHI_MOD;
return block * phiTable[K_SMALL_PHI_MOD][s] + phiTable[(int) rem][s];
}
if (x <= table.limit && (long) table.primes[s - 1] * table.primes[s - 1] > x) {
return table.piSmall[(int) x] - s + 1;
}
long ps = table.primes[s - 1];
if (ps * ps > x)
return pi(x) - s + 1;
if (s <= K_PHI_CACHE_S_MAX && x <= K_PHI_CACHE_X_MAX) {
long key = (x << 8) ^ s;
if (phiCache.containsKey(key))
return phiCache.get(key);
long ans = phi(x, s - 1) - phi(x / ps, s - 1);
phiCache.put(key, ans);
return ans;
}
return phi(x, s - 1) - phi(x / ps, s - 1);
}
long pi(long x) {
if (x <= table.limit)
return table.piSmall[(int) x];
if (piCache.containsKey(x))
return piCache.get(x);
long a = pi(iroot4(x));
long b = pi(isqrt(x));
long c = pi(icbrt(x));
long ans = phi(x, (int) a) + (b + a - 2) * (b - a + 1) / 2;
for (long i = a + 1; i <= b; i++) {
long w = x / table.primes[(int) (i - 1)];
ans -= pi(w);
if (i <= c) {
long bi = pi(isqrt(w));
for (long j = i; j <= bi; j++) {
long q = w / table.primes[(int) (j - 1)];
ans -= pi(q) - (j - 1);
}
}
}
piCache.put(x, ans);
return ans;
}
}
static long fsfFromSignature(List<Integer> signature) {
int r = signature.size();
if (r == 0)
return 0;
List<Integer> masks = new ArrayList<>();
for (int m = 1; m < (1 << r); m++)
masks.add(m);
masks.sort((a, b) -> {
int pa = Integer.bitCount(a);
int pb = Integer.bitCount(b);
if (pa != pb)
return Integer.compare(pb, pa);
return Integer.compare(a, b);
});
int[][] maskBits = new int[masks.size()][];
for (int i = 0; i < masks.size(); i++) {
int mask = masks.get(i);
int[] bits = new int[Integer.bitCount(mask)];
int idx = 0;
for (int k = 0; k < r; k++) {
if (((mask >> k) & 1) == 1)
bits[idx++] = k;
}
maskBits[i] = bits;
}
byte[] remaining = new byte[r];
for (int i = 0; i < r; i++)
remaining[i] = signature.get(i).byteValue();
Map<Long, Long> memo = new HashMap<>();
return fsfDfs(0, maskBits, remaining, r, memo);
}
static long packRemaining(byte[] remaining, int r) {
long packed = 0;
for (int i = 0; i < r; i++)
packed |= (long) remaining[i] << (6 * i);
return packed;
}
static long fsfDfs(int pos, int[][] maskBits, byte[] remaining, int r, Map<Long, Long> memo) {
if (pos == maskBits.length) {
for (int i = 0; i < r; i++)
if (remaining[i] != 0)
return 0;
return 1;
}
long key = packRemaining(remaining, r) ^ ((long) pos * 0x9E3779B97F4A7C15L);
if (memo.containsKey(key))
return memo.get(key);
int[] bits = maskBits[pos];
int maxTake = Integer.MAX_VALUE;
for (int bit : bits)
maxTake = Math.min(maxTake, remaining[bit]);
long total = 0;
for (int take = 0;; take++) {
total += fsfDfs(pos + 1, maskBits, remaining, r, memo);
if (take == maxTake)
break;
for (int bit : bits)
remaining[bit]--;
}
for (int bit : bits)
remaining[bit] += maxTake;
memo.put(key, total);
return total;
}
static void generateSignaturesDfs(int[] primes, int pIdx, int lastExp, long current, long limit, List<Integer> cur,
List<List<Integer>> out) {
if (pIdx >= primes.length)
return;
long p = primes[pIdx];
long pPow = 1;
for (int e = 1; e <= lastExp; e++) {
if (pPow > limit / p)
break;
pPow *= p;
if (current > limit / pPow)
break;
long nextVal = current * pPow;
cur.add(e);
out.add(new ArrayList<>(cur));
generateSignaturesDfs(primes, pIdx + 1, e, nextVal, limit, cur, out);
cur.remove(cur.size() - 1);
}
}
static List<List<Integer>> generateSignatures(long limit, PrimeTable table) {
List<List<Integer>> signatures = new ArrayList<>();
generateSignaturesDfs(table.primes, 0, 63, 1, limit, new ArrayList<>(), signatures);
return signatures;
}
static void generateUniquePermutationsRec(Map<Integer, Integer> freq, int[] cur, int pos, List<List<Integer>> out) {
if (pos == cur.length) {
List<Integer> list = new ArrayList<>();
for (int v : cur)
list.add(v);
out.add(list);
return;
}
for (int val : new ArrayList<>(freq.keySet())) {
int f = freq.get(val);
if (f > 0) {
freq.put(val, f - 1);
cur[pos] = val;
generateUniquePermutationsRec(freq, cur, pos + 1, out);
freq.put(val, f);
}
}
}
static List<List<Integer>> generateUniquePermutations(List<Integer> signature) {
Map<Integer, Integer> freq = new HashMap<>();
for (int e : signature)
freq.put(e, freq.getOrDefault(e, 0) + 1);
List<List<Integer>> out = new ArrayList<>();
generateUniquePermutationsRec(freq, new int[signature.size()], 0, out);
return out;
}
static class OrderedPrimeTupleCounter {
PrimeTable table;
PrimeCounter primeCounter;
OrderedPrimeTupleCounter(PrimeTable table, PrimeCounter primeCounter) {
this.table = table;
this.primeCounter = primeCounter;
}
long countForExponents(List<Integer> exponents, long limit) {
int r = exponents.size();
if (r == 0)
return 0;
if (r == 1) {
long root = kthRootFloor(limit, exponents.get(0));
return primeCounter.pi(root);
}
int[] tailSum = new int[r];
int running = 0;
for (int i = r - 1; i >= 0; i--) {
running += exponents.get(i);
tailSum[i] = running;
}
Map<Long, Long> memo = new HashMap<>();
return countDfs(0, 0, limit, r, exponents, tailSum, memo);
}
long countDfs(int pos, int startIndex, long remLimit, int r, List<Integer> exponents, int[] tailSum,
Map<Long, Long> memo) {
if (pos == r)
return 1;
long key = remLimit ^ ((long) startIndex * 0x9E3779B97F4A7C15L) ^ ((long) pos * 0xC2B2AE3D27D4EB4FL);
if (memo.containsKey(key))
return memo.get(key);
long result = 0;
if (pos == r - 1) {
int exp = exponents.get(pos);
long root = kthRootFloor(remLimit, exp);
long upto = primeCounter.pi(root);
result = (upto > startIndex) ? (upto - startIndex) : 0;
memo.put(key, result);
return result;
}
int tail = tailSum[pos];
long maxP = kthRootFloor(remLimit, tail);
long upto = primeCounter.pi(maxP);
if (upto > table.primes.length)
upto = table.primes.length;
int exp = exponents.get(pos);
for (long idx = startIndex; idx < upto; idx++) {
long p = table.primes[(int) idx];
long pPow = powCapped(p, exp, remLimit);
if (pPow > remLimit)
break;
result += countDfs(pos + 1, (int) (idx + 1), remLimit / pPow, r, exponents, tailSum, memo);
}
memo.put(key, result);
return result;
}
}
static long solveFast(long limit) {
long rootLimit = Math.max(isqrt(limit), 100);
PrimeTable table = new PrimeTable((int) rootLimit);
PrimeCounter primeCounter = new PrimeCounter(table);
OrderedPrimeTupleCounter tupleCounter = new OrderedPrimeTupleCounter(table, primeCounter);
List<List<Integer>> signatures = generateSignatures(limit, table);
long total = 0;
for (List<Integer> signature : signatures) {
long fsf = fsfFromSignature(signature);
List<List<Integer>> perms = generateUniquePermutations(signature);
for (List<Integer> p : perms) {
long cnt = tupleCounter.countForExponents(p, limit);
total += fsf * cnt;
}
}
return total;
}
static String solve() {
return Long.toString(solveFast(10000000000L));
}
public static void main(String[] args) {
System.out.println(solve());
}
}