Project Euler

Long Products

Define F(m, n) as the number of n -tuples (a_1, a_2,..., a_n) of positive integers for which the product does not exceed m: F(m, n) = |bigl{(a_1,..., a_n) in Z_(>0)^n: a_1 a_2... a_n <= mbigr}|. Gi...

Source sync Apr 19, 2026

Problem #0452

Level Level 19

Solved By 701

Languages C++, Python

Answer 345558983

Length 257 words

modular_arithmeticnumber_theorylinear_algebra

Define $F(m,n)$ as the number of $n$-tuples of positive integers for which the product of the elements doesn’t exceed $m$.

You are given that $F(10, 10) = 571$ and $F(10^6, 10^6) \bmod 1\,234\,567\,891 = 252903833$.

Find $F(10^9, 10^9) \bmod 1\,234\,567\,891$.

Problem 452: Long Products

Mathematical Foundation

Theorem 1 (Reduction to Piltz divisor sum). Grouping tuples by product value,

F(m, n) = \sum_{k=1}^{m} d_n(k),

where $d_n(k)$ counts the number of ordered factorizations of $k$ into exactly $n$ positive integer factors (the $n$ -th Piltz divisor function).

Proof. Each $n$ -tuple $(a_1, \ldots, a_n)$ with product $k$ is an ordered factorization of $k$ into $n$ parts. Summing over all $k \le m$ counts every valid tuple exactly once. $\square$

Lemma 1 (Multiplicativity of $d_n$ ). The function $d_n$ is multiplicative: for $\gcd(a, b) = 1$ , $d_n(ab) = d_n(a) d_n(b)$ . At a prime power,

d_n(p^a) = \binom{n + a - 1}{a}.

Proof. Multiplicativity follows from the Chinese Remainder Theorem applied to ordered factorizations. For $p^a$ , distributing $a$ copies of prime $p$ among $n$ positions is a stars-and-bars count: $\binom{n + a - 1}{a}$ . $\square$

Lemma 2 (Exponent bound). Since $m = 10^9$ and $2^{30} > 10^9$ , every prime power $p^a \le m$ satisfies $a \le 29$ . Hence $\binom{n + a - 1}{a}$ is a polynomial of degree at most $29$ in $n$ , evaluable modulo the prime $1\,234\,567\,891$ .

Proof. If $p^a \le 10^9$ then $a \le \lfloor \log_2(10^9) \rfloor = 29$ . The binomial coefficient $\binom{n + a - 1}{a} = \frac{(n+a-1)(n+a-2)\cdots n}{a!}$ is a polynomial of degree $a$ in $n$ . $\square$

Theorem 2 (Bucket decomposition). The set $\{\lfloor m/k \rfloor : k = 1, \ldots, m\}$ has at most $2\lfloor\sqrt{m}\rfloor$ distinct values. Storing partial sums $S[v]$ for each such value $v$ suffices to compute $F(m, n)$ with $O(\sqrt{m})$ storage.

Proof. For $k \le \sqrt{m}$ , the value $\lfloor m/k \rfloor$ can be any of at most $\sqrt{m}$ values. For $k > \sqrt{m}$ , $\lfloor m/k \rfloor < \sqrt{m}$ , giving at most $\sqrt{m}$ values. The union has at most $2\sqrt{m}$ elements. $\square$

Editorial

We initialize bucket arrays. We then precompute C[a] = binomial(n+a-1, a) mod P for a = 0, …, 29. Finally, process small primes (p <= B).

Pseudocode

Input: m, n, modulus P
Output: F(m, n) mod P
Let B = floor(sqrt(m))
Initialize bucket arrays:
Precompute C[a] = binomial(n+a-1, a) mod P for a = 0, ..., 29
Process small primes (p <= B):
For each bucket value v in decreasing order
Process large primes (p > B) via segmented sieve:
For each such prime p
Return hi[1]   (which stores S[floor(m/1)] = S[m] = F(m,n))

Complexity Analysis

Time: $O(m)$ , dominated by the large-prime phase. The small-prime phase costs $O(B^2 / \log B) = O(m / \log m)$ . Each large prime $p$ updates at most $\lfloor m/p \rfloor \le B$ buckets; summing over large primes gives $O(\sum_{p > B} m/p) = O(m)$ .
Space: $O(\sqrt{m})$ for the two bucket arrays and the segmented sieve buffer.

Answer

\boxed{345558983}

C++ project_euler/problem_452/solution.cpp

/*
 * Project Euler Problem 452 — Long Products
 *
 * F(m,n) = number of n-tuples (a1,...,an) with ai >= 1 and product <= m.
 * Compute F(10^9, 10^9) mod 1234567891.
 *
 * F(m,n) = sum_{k=1}^{m} d_n(k) where d_n is the n-th Piltz divisor function.
 * d_n is multiplicative: d_n(p^a) = C(n+a-1, a).
 *
 * Algorithm:
 *   - Maintain partial sums S[v] for O(sqrt(m)) bucket values v = floor(m/k).
 *   - Process small primes (p <= sqrt(m)) with full exponent handling.
 *   - Process large primes (p > sqrt(m)) via segmented sieve, a=1 only.
 *
 * Compile: g++ -O2 -o solution solution_final.cpp
 */

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#include <ctime>

using namespace std;
typedef long long ll;

const ll MOD = 1234567891LL;

ll power(ll base, ll exp, ll m) {
    ll r = 1; base %= m;
    while (exp > 0) {
        if (exp & 1) r = r * base % m;
        base = base * base % m;
        exp >>= 1;
    }
    return r;
}

ll modinv(ll a, ll m) { return power(a, m - 2, m); }

ll binom_mod(ll n_val, int a) {
    if (a == 0) return 1;
    ll num = 1, den = 1;
    for (int i = 1; i <= a; i++) {
        num = num * ((n_val + a - i) % MOD) % MOD;
        den = den * ((ll)i % MOD) % MOD;
    }
    return num * modinv(den, MOD) % MOD;
}

ll solve(ll m, ll n) {
    if (m <= 0) return 0;
    if (m == 1) return 1;

    int sqm = (int)sqrt((double)m);
    while ((ll)(sqm + 1) * (sqm + 1) <= m) sqm++;
    int sz = sqm + 1;

    vector<ll> lo(sz, 0), hi(sz, 0);
    for (int i = 1; i < sz; i++) { lo[i] = 1; hi[i] = 1; }

    auto getS = [&](ll v) -> ll {
        return (v <= sqm) ? lo[(int)v] : hi[(int)(m / v)];
    };

    int max_a = (m >= 2) ? (int)(log((double)m) / log(2.0)) + 2 : 1;
    vector<ll> bc(max_a + 1);
    for (int a = 0; a <= max_a; a++) bc[a] = binom_mod(n, a);

    // Sieve small primes
    vector<int> small_primes;
    {
        vector<bool> ip(sqm + 1, true);
        ip[0] = ip[1] = false;
        for (int i = 2; (ll)i * i <= sqm; i++)
            if (ip[i])
                for (int j = i * i; j <= sqm; j += i)
                    ip[j] = false;
        for (int i = 2; i <= sqm; i++)
            if (ip[i]) small_primes.push_back(i);
    }

    printf("  sqm=%d, small primes: %zu\n", sqm, small_primes.size());

    // Step 1: Process small primes
    for (int p : small_primes) {
        // hi buckets (v > sqm), processed in decreasing v order (increasing k)
        for (int k = 1; k < sz; k++) {
            ll v = m / k;
            if (v <= sqm) break;
            if (v < p) break;
            ll pk = p; int a = 1;
            while (pk <= v) {
                hi[k] = (hi[k] + bc[a] * getS(v / pk)) % MOD;
                if (pk > v / p) break;
                pk *= p; a++;
            }
        }
        // lo buckets
        for (int v = sqm; v >= p; v--) {
            ll pk = p; int a = 1;
            while (pk <= v) {
                lo[v] = (lo[v] + bc[a] * getS(v / pk)) % MOD;
                if (pk > v / p) break;
                pk *= p; a++;
            }
        }
    }

    // Step 2: Process large primes via segmented sieve
    // For p > sqm: only a=1. Only hi buckets affected.
    // Process primes in DECREASING order within each segment.
    ll bc1 = bc[1];
    int seg_size = 1 << 19;  // 512K
    vector<bool> seg(seg_size);

    ll total_large = 0;
    for (ll seg_lo = (ll)sqm + 1; seg_lo <= m; seg_lo += seg_size) {
        ll seg_hi_val = min(seg_lo + seg_size - 1, m);
        int seg_len = (int)(seg_hi_val - seg_lo + 1);
        fill(seg.begin(), seg.begin() + seg_len, true);

        for (int p : small_primes) {
            ll start = ((seg_lo + p - 1) / p) * p;
            if (start < (ll)p * p) start = (ll)p * p;  // composites only
            for (ll j = start; j <= seg_hi_val; j += p)
                seg[(int)(j - seg_lo)] = false;
        }

        // Process primes in decreasing order
        for (int i = seg_len - 1; i >= 0; i--) {
            if (!seg[i]) continue;
            ll p = seg_lo + i;
            total_large++;
            for (int k = 1; k < sz && m / k >= p; k++) {
                ll vp = (m / k) / p;
                hi[k] = (hi[k] + bc1 * lo[(int)vp]) % MOD;
            }
        }
    }

    printf("  large primes processed: %lld\n", total_large);
    return (getS(m) % MOD + MOD) % MOD;
}

int main() {
    printf("=== Project Euler #452: Long Products ===\n\n");

    // Verification
    struct { ll m; ll n; ll expected; } tests[] = {
        {10, 10, 571},
        {100, 100, 613658361},
        {1000, 1000, 1229737624},
        {10000, 10000, 1027202788},
        {100000, 100000, 602287876},
    };
    printf("Verification:\n");
    bool all_ok = true;
    for (auto& t : tests) {
        ll r = solve(t.m, t.n);
        bool ok = (r == t.expected);
        if (!ok) all_ok = false;
        printf("  F(%lld, %lld) = %lld  (expected %lld) [%s]\n",
               t.m, t.n, r, t.expected, ok ? "OK" : "FAIL");
    }
    printf(all_ok ? "All tests passed.\n\n" : "SOME TESTS FAILED!\n\n");

    // Main computation
    printf("Computing F(10^9, 10^9) mod %lld ...\n", MOD);
    clock_t t0 = clock();
    ll ans = solve(1000000000LL, 1000000000LL);
    double elapsed = (double)(clock() - t0) / CLOCKS_PER_SEC;
    printf("\n*** ANSWER: F(10^9, 10^9) mod %lld = %lld ***\n", MOD, ans);
    printf("Time: %.2f seconds\n", elapsed);

    return 0;
}

Python project_euler/problem_452/solution.py

#!/usr/bin/env python3
"""
Project Euler Problem 452 — Long Products

F(m, n) = number of n-tuples (a1, ..., an) of positive integers with product <= m.
Compute F(10^9, 10^9) mod 1234567891.

Key insight:
  F(m, n) = sum_{k=1}^{m} d_n(k)

where d_n(k) = number of ordered factorizations of k into exactly n positive
integer factors.  d_n is multiplicative, and for a prime power p^a:
  d_n(p^a) = C(n + a - 1, a)

Since m = 10^9, the maximum exponent a for any prime is 29 (since 2^30 > 10^9).
C(n+a-1, a) is a polynomial of degree a in n, computable mod the prime modulus.

Strategy — recursive prime-power decomposition:
  We iterate over primes p up to m via a sieve.  For each prime power p^a <= m,
  we accumulate the contribution.  We use a "multiplicative prefix sum" approach:
  process primes one by one, maintaining partial sums over integers composed only
  of primes already processed.

  Concretely, let S(x, pi) = sum of d_n(k) over all k <= x whose smallest prime
  factor is > p_pi (the pi-th prime).  Then:
    F(m, n) = S(m, 0)
  with the recursion:
    S(x, pi) = 1 + sum_{p = next prime after p_pi}^{x} sum_{a=1}^{floor(log_p(x))}
               d_n(p^a) * S(x / p^a, index_of(p))

  The "1" accounts for k = 1 (all ones tuple, d_n(1) = 1).

  This is essentially Lucy_Hedgehog's method adapted for a general multiplicative
  function, implemented via memoized recursion with the key being x (since
  all values of x that appear are of the form floor(m / something), of which
  there are O(sqrt(m)) distinct values).

For m = n = 10^9 and MOD = 1234567891, this runs in roughly O(m^{2/3}) time
with appropriate sieve limits and memoization.

However, for the full 10^9 case, pure Python is too slow. This script solves
the smaller test cases and demonstrates the algorithm. The C++ solution handles
the full problem.
"""

import sys
import math
import time

MOD = 1234567891

def solve_small(m, n, mod=None):
    """
    Direct computation of F(m, n) for small m.
    Factorize each k in [1, m], compute d_n(k), sum them up.
    """
    # Sieve smallest prime factor
    spf = list(range(m + 1))
    for i in range(2, int(m**0.5) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, m + 1, i):
                if spf[j] == j:
                    spf[j] = i

    def d_n(k, n_val):
        """Compute d_n(k) = product of C(n+a-1, a) over prime powers p^a || k."""
        if k == 1:
            return 1
        result = 1
        while k > 1:
            p = spf[k]
            a = 0
            while k % p == 0:
                k //= p
                a += 1
            # C(n_val + a - 1, a)
            binom = 1
            for i in range(1, a + 1):
                binom = binom * (n_val + a - i) // i
            if mod:
                binom %= mod
            result = result * binom
            if mod:
                result %= mod
        return result

    total = 0
    for k in range(1, m + 1):
        total += d_n(k, n_val=n)
        if mod:
            total %= mod
    return total if not mod else total % mod

def modinv(a, m):
    """Modular inverse using extended Euclidean algorithm."""
    g, x, _ = extended_gcd(a % m, m)
    if g != 1:
        return None
    return x % m

def extended_gcd(a, b):
    if a == 0:
        return b, 0, 1
    g, x, y = extended_gcd(b % a, a)
    return g, y - (b // a) * x, x

def binom_mod(n_val, a, mod):
    """Compute C(n_val + a - 1, a) mod `mod`. a is small (<=29)."""
    # C(n_val + a - 1, a) = product_{i=1}^{a} (n_val + a - i) / i
    num = 1
    den = 1
    for i in range(1, a + 1):
        num = num * ((n_val + a - i) % mod) % mod
        den = den * i % mod
    return num * modinv(den, mod) % mod

def solve_medium(m, n, mod=MOD):
    """
    Solve F(m, n) mod `mod` using prime-by-prime multiplicative prefix sum.

    We process primes in increasing order.  We maintain an array/dict where
    for each "bucket" value v (which is always of the form floor(m/k)),
    we store the current partial sum of d_n over all integers <= v composed
    only of primes we have already processed.

    Initially (before processing any prime), only k=1 qualifies, so every
    bucket has value 1.

    When we process prime p, for each bucket v (in decreasing order), and
    for each exponent a = 1, 2, ... while p^a <= v, we add:
      binom(n, a) * bucket_value[floor(v / p^a)]
    where binom(n, a) = C(n+a-1, a) is the d_n contribution of p^a.
    """
    if m <= 1:
        return 1 if m == 1 else 0

    # Generate all distinct values of floor(m / k)
    # These are O(2*sqrt(m)) values
    sqrt_m = int(m**0.5)
    vals = []
    for i in range(1, sqrt_m + 1):
        vals.append(m // i)
        vals.append(i)
    vals = sorted(set(vals))

    # Map value -> index
    # Use two arrays: for v <= sqrt_m, index directly; for v > sqrt_m, use m//v
    # which is <= sqrt_m.
    # S[v] = current partial sum for bucket v
    # We store S in a dict for simplicity (Python version)
    S = {}
    for v in vals:
        S[v] = 1  # initially only k=1 contributes

    # Sieve primes up to sqrt(m) (we need primes up to m, but primes > sqrt(m)
    # contribute only at exponent a=1)
    limit = m
    # For a full sieve up to m we need memory; use segmented or just sieve up
    # to sqrt(m) for the "small prime" part, then handle large primes.

    # Actually for this approach we sieve primes up to m.
    # For m up to ~10^6 this is fine in Python.  For 10^9 use C++.
    if m > 5_000_000:
        print(f"WARNING: m={m} is large for Python. Use C++ for full solution.")
        return None

    sieve = bytearray(b'\x01') * (m + 1)
    sieve[0] = sieve[1] = 0
    for i in range(2, int(m**0.5) + 1):
        if sieve[i]:
            for j in range(i * i, m + 1, i):
                sieve[j] = 0
    primes = [i for i in range(2, m + 1) if sieve[i]]

    # Precompute binom values: C(n+a-1, a) mod for a = 1..29
    max_a = int(math.log2(m)) + 1 if m >= 2 else 1
    binom_cache = [0] * (max_a + 1)
    for a in range(1, max_a + 1):
        binom_cache[a] = binom_mod(n, a, mod)

    # Process each prime
    for p in primes:
        # We iterate over bucket values in DECREASING order to avoid
        # counting the same prime twice (similar to knapsack)
        for v in reversed(vals):
            if v < p:
                break
            # For each exponent a = 1, 2, ... while p^a <= v
            pk = p  # p^a
            a = 1
            while pk <= v:
                # floor(v / p^a)
                v_div = v // pk
                S[v] = (S[v] + binom_cache[a] * S[v_div]) % mod
                if pk > v // p:
                    break
                pk *= p
                a += 1

    return S[m]

def verify():
    """Verify against known values."""
    print("Verifying F(10, 10)...")
    result = solve_small(10, 10)
    print(f"  F(10, 10) = {result}  (expected 571)")
    assert result == 571, f"FAILED: got {result}"

    print("Verifying with medium solver...")
    result2 = solve_medium(10, 10, MOD)
    print(f"  F(10, 10) via medium = {result2}  (expected 571)")
    assert result2 == 571, f"FAILED: got {result2}"

    # Test F(100, 100)
    r1 = solve_small(100, 100, MOD)
    r2 = solve_medium(100, 100, MOD)
    print(f"  F(100, 100) = {r1} (small) vs {r2} (medium)")
    assert r1 == r2, f"MISMATCH: {r1} vs {r2}"

    # Test F(1000, 1000)
    r1 = solve_small(1000, 1000, MOD)
    r2 = solve_medium(1000, 1000, MOD)
    print(f"  F(1000, 1000) = {r1} (small) vs {r2} (medium)")
    assert r1 == r2, f"MISMATCH: {r1} vs {r2}"

    print("All verifications passed!")

def create_visualization():
    """Create visualization of the problem structure."""