Problem 452: Long Products
View on Project EulerProject Euler Problem 452 Solution
EulerSolve provides an optimized solution for Project Euler Problem 452, Long Products, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Define $$F(m,n)=\left|\left\{(a_1,\dots,a_n)\in \mathbb{Z}_{>0}^n:\prod_{i=1}^{n} a_i \le m\right\}\right|.$$ The task is to evaluate \(F(10^9,10^9)\) modulo \(1234567891\). The published checkpoints are $$F(10,10)=571,\qquad F(10^6,10^6)\equiv 252903833 \pmod{1234567891}.$$ A direct enumeration of \(10^9\)-tuples is impossible, so the solution has to exploit the multiplicative structure of the product bound. Mathematical Approach Step 1: Dynamic Programming on the Remaining Product Budget For a fixed tuple length \(t\), let $$P_t(v)=F(v,t).$$ When \(t=0\), there is exactly one admissible tuple: the empty tuple. Therefore $$P_0(v)=1\qquad (v\ge 1).$$ For \(t\ge 1\), choose the first entry \(x=a_1\). Once \(x\) is fixed, the remaining \(t-1\) entries must satisfy $$a_2a_3\cdots a_t \le \left\lfloor \frac{v}{x}\right\rfloor.$$ This gives the exact recurrence $$P_t(v)=\sum_{x=1}^{v} P_{t-1}\!\left(\left\lfloor \frac{v}{x}\right\rfloor\right).$$ That identity is the backbone of all three implementations: the problem is transformed into repeated sums over floor quotients. Step 2: Compress the State Space with Distinct Floor Quotients For fixed \(m\), define the compressed state set $$\mathcal{V}(m)=\left\{\left\lfloor \frac{m}{i}\right\rfloor : 1\le i\le m\right\}.$$ This set contains only \(O(\sqrt{m})\) distinct values....
Detailed mathematical approach
Problem Summary
Define
$$F(m,n)=\left|\left\{(a_1,\dots,a_n)\in \mathbb{Z}_{>0}^n:\prod_{i=1}^{n} a_i \le m\right\}\right|.$$
The task is to evaluate \(F(10^9,10^9)\) modulo \(1234567891\). The published checkpoints are
$$F(10,10)=571,\qquad F(10^6,10^6)\equiv 252903833 \pmod{1234567891}.$$
A direct enumeration of \(10^9\)-tuples is impossible, so the solution has to exploit the multiplicative structure of the product bound.
Mathematical Approach
Step 1: Dynamic Programming on the Remaining Product Budget
For a fixed tuple length \(t\), let
$$P_t(v)=F(v,t).$$
When \(t=0\), there is exactly one admissible tuple: the empty tuple. Therefore
$$P_0(v)=1\qquad (v\ge 1).$$
For \(t\ge 1\), choose the first entry \(x=a_1\). Once \(x\) is fixed, the remaining \(t-1\) entries must satisfy
$$a_2a_3\cdots a_t \le \left\lfloor \frac{v}{x}\right\rfloor.$$
This gives the exact recurrence
$$P_t(v)=\sum_{x=1}^{v} P_{t-1}\!\left(\left\lfloor \frac{v}{x}\right\rfloor\right).$$
That identity is the backbone of all three implementations: the problem is transformed into repeated sums over floor quotients.
Step 2: Compress the State Space with Distinct Floor Quotients
For fixed \(m\), define the compressed state set
$$\mathcal{V}(m)=\left\{\left\lfloor \frac{m}{i}\right\rfloor : 1\le i\le m\right\}.$$
This set contains only \(O(\sqrt{m})\) distinct values. The usual reason is that quotients larger than \(\sqrt{m}\) can occur only for \(i\le \sqrt{m}\), while every quotient at most \(\sqrt{m}\) lies in the small tail. Hence \(|\mathcal{V}(m)|\le 2\sqrt{m}\).
The crucial closure property is that every transition in the recurrence stays inside the same compressed set. If \(v=\lfloor m/i\rfloor\), then
$$\left\lfloor \frac{v}{x}\right\rfloor=\left\lfloor \frac{\lfloor m/i\rfloor}{x}\right\rfloor=\left\lfloor \frac{m}{ix}\right\rfloor \in \mathcal{V}(m).$$
So the dynamic program never needs values outside \(\mathcal{V}(m)\).
For each state \(v\in\mathcal{V}(m)\), the terms \(\lfloor v/x\rfloor\) are constant on intervals. If
$$q=\left\lfloor \frac{v}{\ell}\right\rfloor,$$
then the same quotient persists for all \(x\) in
$$\ell \le x \le r=\left\lfloor \frac{v}{q}\right\rfloor.$$
Therefore the transition can be grouped by quotient blocks:
$$P_t(v)=\sum_{[\ell,r]} (r-\ell+1)\,P_{t-1}\!\left(\left\lfloor \frac{v}{\ell}\right\rfloor\right).$$
The implementation precomputes these blocks once and then reuses them for every DP layer.
Step 3: Why \(F(m,n)\) is a Polynomial in \(n\)
Remove every entry equal to \(1\) from an admissible tuple. The remaining ordered core tuple has length \(r\), every entry is at least \(2\), and its product is still at most \(m\). Hence
$$2^r \le m,$$
so
$$r \le d=\left\lfloor \log_2 m \right\rfloor.$$
Now let \(C_r(m)\) be the number of ordered \(r\)-tuples \((b_1,\dots,b_r)\) with each \(b_i\ge 2\) and \(b_1\cdots b_r\le m\). Once such a core tuple is fixed, we obtain an \(n\)-tuple by choosing which \(r\) of the \(n\) positions carry the non-1 entries. That contributes exactly \(\binom{n}{r}\) placements. Therefore
$$F(m,n)=\sum_{r=0}^{d} C_r(m)\binom{n}{r}.$$
So for fixed \(m\), the answer is a polynomial in \(n\) of degree at most \(d\), written in the binomial basis. For \(m=10^9\), we have \(d=29\), which is why the implementation only needs the values for \(n=0,1,\dots,29\).
Step 4: Forward Differences and the Composite Modulus
Let \(Q(n)=F(m,n)\). Any polynomial of degree at most \(d\) has the Newton expansion
$$Q(n)=\sum_{r=0}^{d} \Delta^r Q(0)\binom{n}{r},$$
where \(\Delta\) denotes the forward-difference operator. In this problem the coefficients \(\Delta^r Q(0)\) coincide with the binomial-basis coefficients \(C_r(m)\).
The DP is used only to compute the sample values
$$Q(0),Q(1),\dots,Q(d).$$
After that, a short forward-difference table recovers the coefficients, and the polynomial is evaluated at the huge target \(n=10^9\).
The modulus \(1234567891\) is composite, so the program does not rely on Fermat-style modular inverses for \(\binom{n}{r}\). Instead it builds the numerator
$$\frac{n(n-1)\cdots(n-r+1)}{r!}$$
with exact integer cancellation by gcd, and only then reduces the result modulo \(1234567891\). Since \(r\le 29\), this step is tiny.
Worked Example: \(m=10\)
Here \(d=\lfloor \log_2 10\rfloor=3\). The dynamic program gives
$$Q(0)=1,\qquad Q(1)=10,\qquad Q(2)=27,\qquad Q(3)=53.$$
The forward differences are
$$\Delta Q(0)=9,\qquad \Delta^2 Q(0)=8,\qquad \Delta^3 Q(0)=1.$$
Hence
$$F(10,n)=1+9\binom{n}{1}+8\binom{n}{2}+\binom{n}{3}.$$
This also has a direct combinatorial meaning:
$$C_0(10)=1,\quad C_1(10)=9,\quad C_2(10)=8,\quad C_3(10)=1.$$
The eight ordered pairs with product at most \(10\) are \((2,2),(2,3),(2,4),(2,5),(3,2),(3,3),(4,2),(5,2)\), and the only ordered triple is \((2,2,2)\). Evaluating at \(n=10\) gives
$$F(10,10)=1+9\binom{10}{1}+8\binom{10}{2}+\binom{10}{3}=571,$$
which matches the checkpoint exactly.
How the Code Works
The implementation first builds the compressed list of distinct quotients \(\lfloor m/i\rfloor\). For every compressed state it precomputes the quotient blocks \([\ell,r]\) together with the target compressed state and the multiplicity \(r-\ell+1\). It then runs the recurrence for tuple lengths \(0\) through \(d\), keeping only the previous and current DP layers, and records the values at the state \(v=m\).
Those sampled values are converted into forward differences, and the final Newton series is evaluated at the requested \(n\). The C++, Python, and Java implementations all follow that same mathematical pipeline, so the checkpoints agree across languages.
Complexity Analysis
The compressed state count is \(S=|\mathcal{V}(m)|=O(\sqrt{m})\). For a state \(v\), the number of quotient blocks is \(O(\sqrt{v})\), so the total precomputed transition size is
$$B=\sum_{v\in \mathcal{V}(m)} O(\sqrt{v})=O(m^{3/4}).$$
Running the DP for \(d=\lfloor \log_2 m\rfloor\) layers therefore costs \(O(Bd)\) time after precomputation, and the memory usage is \(O(B+S)\). The forward-difference table and the final binomial evaluation are only \(O(d^2)\) and \(O(d)\), negligible here because \(d=29\) for \(m=10^9\).
References
- Problem page: https://projecteuler.net/problem=452
- Divisor summatory function and grouping equal floor quotients: Wikipedia - Divisor summatory function
- Newton series: Wikipedia - Newton series
- Finite differences: Wikipedia - Finite difference
- Graham, Knuth, Patashnik, Concrete Mathematics, 2nd ed., Addison-Wesley, chapters on sums and finite differences.
Problem 452 source code
C++
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = __uint128_t;
constexpr u32 kMod = 1'234'567'891U;
struct Options {
u32 m = 1'000'000'000U;
u32 n = 1'000'000'000U;
bool run_checkpoints = true;
};
bool parse_u32_after_prefix(const std::string& arg, const std::string& prefix, u32& out) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
const unsigned long long parsed = std::stoull(tail);
if (parsed > static_cast<unsigned long long>(std::numeric_limits<u32>::max())) {
return false;
}
out = static_cast<u32>(parsed);
} catch (...) {
return false;
}
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 (parse_u32_after_prefix(arg, "--m=", options.m)) {
continue;
}
if (parse_u32_after_prefix(arg, "--n=", options.n)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.m >= 1U;
}
u32 choose_small_mod(const u32 n, const u32 r) {
if (r > n) {
return 0U;
}
if (r == 0U) {
return 1U;
}
std::vector<u64> numerators(r);
for (u32 i = 0U; i < r; ++i) {
numerators[static_cast<std::size_t>(i)] = static_cast<u64>(n - r + 1U + i);
}
for (u32 d = 2U; d <= r; ++d) {
u64 x = d;
for (u32 i = 0U; i < r && x > 1ULL; ++i) {
u64& v = numerators[static_cast<std::size_t>(i)];
const u64 g = std::gcd(v, x);
if (g > 1ULL) {
v /= g;
x /= g;
}
}
}
u64 result = 1ULL;
for (const u64 v : numerators) {
result = (result * (v % kMod)) % kMod;
}
return static_cast<u32>(result);
}
u32 solve(const u32 m, const u32 n_target) {
int degree = 0;
while ((1ULL << (degree + 1)) <= m) {
++degree;
}
std::vector<u32> values;
values.reserve(static_cast<std::size_t>(2.0L * std::sqrt(static_cast<long double>(m)) + 8.0L));
for (u32 i = 1U; i <= m;) {
const u32 q = m / i;
const u32 j = m / q;
values.push_back(q);
i = j + 1U;
}
const u32 state_count = static_cast<u32>(values.size());
u32 sqrt_m = static_cast<u32>(std::sqrt(static_cast<long double>(m)));
while ((static_cast<u64>(sqrt_m) + 1ULL) * (static_cast<u64>(sqrt_m) + 1ULL) <= m) {
++sqrt_m;
}
while (static_cast<u64>(sqrt_m) * static_cast<u64>(sqrt_m) > m) {
--sqrt_m;
}
const auto index_of = [&](const u32 q) -> u32 {
if (q <= sqrt_m) {
return state_count - q;
}
return m / q - 1U;
};
std::vector<u32> offset(static_cast<std::size_t>(state_count) + 1U, 0U);
std::vector<u32> q_index;
std::vector<u32> count_mod;
for (u32 idx = 0U; idx < state_count; ++idx) {
const u32 v = values[static_cast<std::size_t>(idx)];
offset[static_cast<std::size_t>(idx)] = static_cast<u32>(q_index.size());
for (u32 i = 1U; i <= v;) {
const u32 q = v / i;
const u32 j = v / q;
q_index.push_back(index_of(q));
count_mod.push_back((j - i + 1U) % kMod);
i = j + 1U;
}
}
offset[static_cast<std::size_t>(state_count)] = static_cast<u32>(q_index.size());
std::vector<u32> prev(state_count, 1U); // F(v,0)=1
std::vector<u32> curr(state_count, 0U);
std::vector<u32> poly_values(static_cast<std::size_t>(degree) + 1U, 0U);
poly_values[0] = 1U;
for (int t = 1; t <= degree; ++t) {
for (u32 idx = 0U; idx < state_count; ++idx) {
const u32 begin = offset[static_cast<std::size_t>(idx)];
const u32 end = offset[static_cast<std::size_t>(idx) + 1U];
u128 sum = 0;
int chunk = 0;
for (u32 p = begin; p < end; ++p) {
sum += static_cast<u128>(count_mod[static_cast<std::size_t>(p)]) *
prev[static_cast<std::size_t>(q_index[static_cast<std::size_t>(p)])];
++chunk;
if (chunk == 64) {
sum %= kMod;
chunk = 0;
}
}
curr[static_cast<std::size_t>(idx)] = static_cast<u32>(sum % kMod);
}
prev.swap(curr);
poly_values[static_cast<std::size_t>(t)] = prev[0];
}
std::vector<u32> diff = poly_values;
std::vector<u32> coeff(static_cast<std::size_t>(degree) + 1U, 0U);
for (int r = 0; r <= degree; ++r) {
coeff[static_cast<std::size_t>(r)] = diff[0];
for (int i = 0; i < degree - r; ++i) {
u32 v = diff[static_cast<std::size_t>(i + 1)] + kMod - diff[static_cast<std::size_t>(i)];
if (v >= kMod) {
v -= kMod;
}
diff[static_cast<std::size_t>(i)] = v;
}
}
u64 answer = 0ULL;
for (int r = 0; r <= degree; ++r) {
const u32 c = coeff[static_cast<std::size_t>(r)];
const u32 binom = choose_small_mod(n_target, static_cast<u32>(r));
answer += static_cast<u64>(c) * binom;
answer %= kMod;
}
return static_cast<u32>(answer);
}
bool run_checkpoints() {
if (solve(10U, 10U) != 571U) {
std::cerr << "Checkpoint failed: F(10,10)\n";
return false;
}
if (solve(1'000'000U, 1'000'000U) != 252'903'833U) {
std::cerr << "Checkpoint failed: F(10^6,10^6)\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 2;
}
std::cout << solve(options.m, options.n) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 1234567891
m, n_target = 1000000000, 1000000000
degree = 0
while (1 << (degree + 1)) <= m:
degree += 1
# Collect distinct floor(m/i) values
values = []
i = 1
while i <= m:
q = m // i
j = m // q
values.append(q)
i = j + 1
state_count = len(values)
sqrt_m = math.isqrt(m)
while (sqrt_m + 1) ** 2 <= m: sqrt_m += 1
while sqrt_m * sqrt_m > m: sqrt_m -= 1
def index_of(q):
if q <= sqrt_m: return state_count - q
return m // q - 1
# Precompute floor division structure for each value
offsets = []
q_index = []
count_mod = []
for idx in range(state_count):
v = values[idx]
offsets.append(len(q_index))
i = 1
while i <= v:
q = v // i
j = v // q
q_index.append(index_of(q))
count_mod.append((j - i + 1) % MOD)
i = j + 1
offsets.append(len(q_index))
prev = [1] * state_count # F(v,0) = 1
poly_values = [0] * (degree + 1)
poly_values[0] = 1
for t in range(1, degree + 1):
curr = [0] * state_count
for idx in range(state_count):
s = 0
for p in range(offsets[idx], offsets[idx + 1]):
s += count_mod[p] * prev[q_index[p]]
curr[idx] = s % MOD
prev = curr
poly_values[t] = prev[0]
# Newton forward differences
diff = list(poly_values)
coeff = [0] * (degree + 1)
for r in range(degree + 1):
coeff[r] = diff[0]
for i in range(degree - r):
diff[i] = (diff[i + 1] - diff[i]) % MOD
def choose_small_mod(n_, r):
if r > n_: return 0
if r == 0: return 1
nums = list(range(n_ - r + 1, n_ + 1))
for d in range(2, r + 1):
x = d
for i in range(r):
if x <= 1: break
g = math.gcd(nums[i], x)
if g > 1:
nums[i] //= g
x //= g
result = 1
for v in nums:
result = result * (v % MOD) % MOD
return result
answer = 0
for r in range(degree + 1):
c = coeff[r] % MOD
b = choose_small_mod(n_target, r)
answer = (answer + c * b) % MOD
return str(answer)
if __name__ == '__main__':
print(solve())
Java
import java.nio.file.*;
import java.util.*;
import java.util.regex.*;
public class Euler452 {
private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");
private static String parseOutput(String stdout) {
String[] lines = stdout.split("\\R");
List<String> nonEmpty = new ArrayList<>();
for (String line : lines) {
String t = line.trim();
if (!t.isEmpty()) {
nonEmpty.add(t);
}
}
if (nonEmpty.isEmpty()) {
return "";
}
List<String> answers = new ArrayList<>();
List<String> equals = new ArrayList<>();
for (String line : nonEmpty) {
Matcher m1 = ANSWER_RE.matcher(line);
if (m1.find()) {
answers.add(m1.group(1).trim());
}
Matcher m2 = EQUAL_RE.matcher(line);
if (m2.find()) {
equals.add(m2.group(1).trim());
}
}
if (!answers.isEmpty()) {
return answers.get(answers.size() - 1);
}
if (!equals.isEmpty()) {
return equals.get(equals.size() - 1);
}
return nonEmpty.get(nonEmpty.size() - 1);
}
private static String pickCompiler() throws Exception {
for (String compiler : List.of("clang++", "g++")) {
Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
.redirectErrorStream(true)
.start();
String out = new String(probe.getInputStream().readAllBytes());
int rc = probe.waitFor();
if (rc == 0 && !out.trim().isEmpty()) {
return compiler;
}
}
throw new RuntimeException("No C++ compiler found (clang++/g++).");
}
private static Path cppSource(Path root) {
return root.resolve("solutionsCpp").resolve("Euler452.cpp");
}
private static boolean shouldSkipCheckpoints(Path root) {
Path src = cppSource(root);
try {
String text = Files.readString(src);
return text.contains("--skip-checkpoints");
} catch (Exception ex) {
return false;
}
}
private static Path ensureBridgeBinary() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = root.resolve("solutionsCpp").resolve(".euler452_java_bridge");
boolean rebuild = Files.notExists(bin)
|| Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;
if (rebuild) {
String compiler = pickCompiler();
Process compile = new ProcessBuilder(
compiler,
"-std=c++17",
"-O2",
src.toString(),
"-o",
bin.toString())
.inheritIO()
.start();
if (compile.waitFor() != 0) {
throw new RuntimeException("Failed to compile Euler452 C++ bridge.");
}
}
return bin;
}
private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
List<String> cmd = new ArrayList<>();
cmd.add(bin.toString());
if (shouldSkipCheckpoints(root)) {
cmd.add("--skip-checkpoints");
}
Process first = new ProcessBuilder(cmd)
.directory(root.toFile())
.redirectErrorStream(true)
.start();
String out = new String(first.getInputStream().readAllBytes());
int rc = first.waitFor();
if (rc == 0) {
return out;
}
Process second = new ProcessBuilder(cmd)
.directory(srcDir.toFile())
.redirectErrorStream(true)
.start();
String out2 = new String(second.getInputStream().readAllBytes());
int rc2 = second.waitFor();
if (rc2 == 0) {
return out2;
}
throw new RuntimeException("Euler452 C++ bridge failed.\n" + out + "\n" + out2);
}
private static String solveViaCppBridge() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = ensureBridgeBinary();
String out = runBridge(bin, root, src.getParent());
String parsed = parseOutput(out);
if (parsed.isEmpty()) {
throw new RuntimeException("Euler452 C++ bridge produced empty output.");
}
return parsed;
}
public static void main(String[] args) throws Exception {
System.out.println(solveViaCppBridge());
}
}