Problem 320: Factorials Divisible by a Huge Integer

Project Euler Problem 320 Solution

EulerSolve provides an optimized solution for Project Euler Problem 320, Factorials Divisible by a Huge Integer, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $N(i)$ as the smallest integer $n$ such that $$n! \text{ is divisible by } (i!)^{1234567890}.$$ Then $$S(u)=\sum_{i=10}^{u}N(i).$$ We are given $$S(1000)=614538266565663,$$ and must compute $$S(1000000)\bmod 10^{18}.$$ Mathematical Approach 1) Divisibility of factorials is a prime-valuation problem. Write $$i!=\prod_{p\le i}p^{e_{i,p}},\qquad e_{i,p}=v_p(i!).$$ Then $$(i!)^K=\prod_{p\le i}p^{K e_{i,p}},\qquad K=1234567890.$$ Therefore $$(i!)^K\mid n!\iff v_p(n!)\ge K e_{i,p}\quad\text{for every prime }p\le i.$$ So the whole problem reduces to matching required prime exponents inside $n!$. 2) Legendre's formula tells us the exponent of one prime in $n!$. For a fixed prime $p$, $$v_p(n!)=\sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor.$$ This function is monotone increasing in $n$. So if we want the smallest $n$ satisfying $$v_p(n!)\ge t,$$ we can find it by binary search. 3) Per-prime bottlenecks. For each prime $p$, define the current target $$t_p=K\,v_p(i!).$$ Let $n_p$ be the smallest integer with $$v_p(n_p!)\ge t_p.$$ Then the smallest $n$ satisfying all prime conditions simultaneously is simply $$N(i)=\max_{p\le i} n_p.$$ The largest per-prime obstruction determines the answer. 4) Incremental update from $i-1$ to $i$. This is the key optimization....

Detailed mathematical approach

Problem Summary

Define $N(i)$ as the smallest integer $n$ such that

$$n! \text{ is divisible by } (i!)^{1234567890}.$$

Then

$$S(u)=\sum_{i=10}^{u}N(i).$$

We are given

$$S(1000)=614538266565663,$$

and must compute

$$S(1000000)\bmod 10^{18}.$$

Mathematical Approach

1) Divisibility of factorials is a prime-valuation problem.

Write

$$i!=\prod_{p\le i}p^{e_{i,p}},\qquad e_{i,p}=v_p(i!).$$

Then

$$(i!)^K=\prod_{p\le i}p^{K e_{i,p}},\qquad K=1234567890.$$

Therefore

$$(i!)^K\mid n!\iff v_p(n!)\ge K e_{i,p}\quad\text{for every prime }p\le i.$$

So the whole problem reduces to matching required prime exponents inside $n!$.

2) Legendre's formula tells us the exponent of one prime in $n!$.

For a fixed prime $p$,

$$v_p(n!)=\sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor.$$

This function is monotone increasing in $n$. So if we want the smallest $n$ satisfying

$$v_p(n!)\ge t,$$

we can find it by binary search.

3) Per-prime bottlenecks.

For each prime $p$, define the current target

$$t_p=K\,v_p(i!).$$

Let $n_p$ be the smallest integer with

$$v_p(n_p!)\ge t_p.$$

Then the smallest $n$ satisfying all prime conditions simultaneously is simply

$$N(i)=\max_{p\le i} n_p.$$

The largest per-prime obstruction determines the answer.

4) Incremental update from $i-1$ to $i$.

This is the key optimization. Suppose

$$i=\prod p^{\alpha_p}.$$

Then

$$v_p(i!)=v_p((i-1)!)+\alpha_p,$$

so the required target updates only for primes dividing $i$:

$$t_p\leftarrow t_p+K\alpha_p.$$

All other primes keep exactly the same target. Therefore, when we move from $i-1$ to $i$, we only need to recompute the primes appearing in the factorization of $i$.

5) Why the global maximum never decreases.

Every target $t_p$ is nondecreasing as $i$ grows. Hence every minimal witness $n_p$ is also nondecreasing. So

$$N(i)=\max_p n_p$$

can never go down. This is why the code stores a single running value current_max and only raises it when a changed prime produces a larger $n_p$.

6) Lower bound for binary search.

Legendre's formula implies

$$v_p(n!)\le \frac{n}{p-1}.$$

So any $n$ meeting $v_p(n!)\ge t_p$ must satisfy

$$n\ge t_p(p-1).$$

The implementation uses this as a cheap lower bound, and it also keeps the previous solution for that prime as another lower bound because $n_p$ never decreases.

7) Example of the incremental logic.

When we move from $i=9$ to $i=10$, only the prime factors of $10$ matter:

$$10=2\cdot 5.$$

So the required exponents update as

$$t_2\leftarrow t_2+K,\qquad t_5\leftarrow t_5+K,$$

while the targets for $3,7,\dots$ stay unchanged. Therefore we recompute only $n_2$ and $n_5$, and then compare them to the existing global maximum.

8) Fast factorization of every $i$.

To make the incremental update cheap, the program precomputes the smallest prime factor (SPF) of every number up to $10^6$. Then each $i$ is factorized in essentially logarithmic time by repeatedly dividing by its SPF.

Algorithm

1) Build an SPF sieve up to $u$.

2) Maintain, for every prime $p$, the required target $t_p$ and its minimal witness $n_p$.

3) For each $i=2,3,\dots,u$, factorize $i$.

4) For each prime $p^{\alpha_p}$ in $i$, update

$$t_p\leftarrow t_p+K\alpha_p,$$

then recompute $n_p$ by binary search on Legendre's formula.

5) Update the running maximum $N(i)$.

6) Add $N(i)$ to the sum whenever $i\ge10$.

Complexity Analysis

The sieve costs roughly linear time in $u$. Each $i$ changes only the primes in its factorization, which is very small on average. Each changed prime needs a binary search, and each midpoint evaluation computes

$$v_p(n!)=\sum_{j\ge1}\left\lfloor\frac{n}{p^j}\right\rfloor$$

in $O(\log_p n)$. This is vastly faster than recomputing all prime constraints from scratch for every $i$.

Checks And Final Result

The source checks

$$S(1000)=614538266565663.$$

For the full problem, it computes

$$S(1000000)\bmod 10^{18}=278157919195482643.$$

Problem 320 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 K = 1234567890ULL;
constexpr u64 OUTPUT_MOD = 1000000000000000000ULL;

u64 valuation_factorial(u64 n, const int p) {
    u64 sum = 0;
    const u64 prime = static_cast<u64>(p);
    while (n > 0) {
        n /= prime;
        sum += n;
    }
    return sum;
}

u64 minimal_n_for_target(const int p, const u64 target, const u64 previous) {
    if (target == 0ULL) {
        return 0ULL;
    }

    const u64 factor = static_cast<u64>(p - 1);
    const u64 lower = std::max(previous, static_cast<u64>(static_cast<u128>(target) * factor));

    u64 hi = lower;
    while (true) {
        const u64 have = valuation_factorial(hi, p);
        if (have >= target) {
            break;
        }
        hi = static_cast<u64>(static_cast<u128>(hi) + static_cast<u128>(target - have) * factor);
    }

    u64 lo = lower;
    while (lo < hi) {
        const u64 mid = lo + (hi - lo) / 2ULL;
        if (valuation_factorial(mid, p) >= target) {
            hi = mid;
        } else {
            lo = mid + 1ULL;
        }
    }
    return lo;
}

std::vector<int> build_spf(const int limit) {
    std::vector<int> spf(static_cast<std::size_t>(limit + 1));
    std::iota(spf.begin(), spf.end(), 0);

    for (int i = 2; static_cast<long long>(i) * i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return spf;
}

u128 solve_sum(const int upper) {
    const std::vector<int> spf = build_spf(upper);

    std::vector<u64> required(static_cast<std::size_t>(upper + 1), 0ULL);
    std::vector<u64> needed_n(static_cast<std::size_t>(upper + 1), 0ULL);

    u64 current_max = 0;
    u128 sum = 0;

    for (int i = 2; i <= upper; ++i) {
        int x = i;
        while (x > 1) {
            const int p = spf[static_cast<std::size_t>(x)];
            int exp = 0;
            while (x % p == 0) {
                x /= p;
                ++exp;
            }

            required[static_cast<std::size_t>(p)] += K * static_cast<u64>(exp);
            const u64 updated = minimal_n_for_target(p, required[static_cast<std::size_t>(p)],
                                                     needed_n[static_cast<std::size_t>(p)]);
            needed_n[static_cast<std::size_t>(p)] = updated;
            if (updated > current_max) {
                current_max = updated;
            }
        }

        if (i >= 10) {
            sum += static_cast<u128>(current_max);
        }
    }

    return sum;
}

bool run_checkpoints() {
    const u128 sample = solve_sum(1000);
    if (sample != static_cast<u128>(614538266565663ULL)) {
        std::cerr << "Checkpoint failed: S(1000) mismatch\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const u128 total = solve_sum(1000000);
    const u64 answer = static_cast<u64>(total % static_cast<u128>(OUTPUT_MOD));
    std::cout << answer << '\n';
    return 0;
}

Python

def solve():
    K = 1234567890
    OUTPUT_MOD = 10**18
    upper = 1000000

    def valuation_factorial(n, p):
        s = 0
        while n > 0:
            n //= p
            s += n
        return s

    def minimal_n_for_target(p, target, previous):
        if target == 0:
            return 0
        factor = p - 1
        lower = max(previous, target * factor)
        hi = lower
        while True:
            have = valuation_factorial(hi, p)
            if have >= target:
                break
            hi += (target - have) * factor
        lo = lower
        while lo < hi:
            mid = lo + (hi - lo) // 2
            if valuation_factorial(mid, p) >= target:
                hi = mid
            else:
                lo = mid + 1
        return lo

    # Build SPF
    spf = list(range(upper + 1))
    for i in range(2, upper + 1):
        if i * i > upper:
            break
        if spf[i] == i:
            for j in range(i*i, upper + 1, i):
                if spf[j] == j:
                    spf[j] = i

    required = [0] * (upper + 1)
    needed_n = [0] * (upper + 1)
    current_max = 0
    total = 0

    for i in range(2, upper + 1):
        x = i
        while x > 1:
            p = spf[x]
            exp = 0
            while x % p == 0:
                x //= p
                exp += 1
            required[p] += K * exp
            updated = minimal_n_for_target(p, required[p], needed_n[p])
            needed_n[p] = updated
            if updated > current_max:
                current_max = updated
        if i >= 10:
            total += current_max

    return str(total % OUTPUT_MOD)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;

public class Euler320 {
    static final long K = 1234567890L;
    static final long OUTPUT_MOD = 1000000000000000000L;

    static long valuationFactorial(long n, int p) {
        long sum = 0;
        long prime = p;
        while (n > 0) {
            n /= prime;
            sum += n;
        }
        return sum;
    }

    static long minimalNForTarget(int p, long target, long previous) {
        if (target == 0)
            return 0;

        long factor = p - 1;
        long lower = Math.max(previous, target * factor);

        long hi = lower;
        while (true) {
            long have = valuationFactorial(hi, p);
            if (have >= target) {
                break;
            }
            hi += (target - have) * factor;
        }

        long lo = lower;
        while (lo < hi) {
            long mid = lo + (hi - lo) / 2;
            if (valuationFactorial(mid, p) >= target) {
                hi = mid;
            } else {
                lo = mid + 1;
            }
        }
        return lo;
    }

    static int[] buildSpf(int limit) {
        int[] spf = new int[limit + 1];
        for (int i = 0; i <= limit; i++)
            spf[i] = i;
        for (int i = 2; i * i <= limit; ++i) {
            if (spf[i] == i) {
                for (int j = i * i; j <= limit; j += i) {
                    if (spf[j] == j)
                        spf[j] = i;
                }
            }
        }
        return spf;
    }

    public static String solve() {
        int upper = 1000000;
        int[] spf = buildSpf(upper);

        long[] required = new long[upper + 1];
        long[] neededN = new long[upper + 1];

        long currentMax = 0;
        BigInteger sum = BigInteger.ZERO;

        for (int i = 2; i <= upper; ++i) {
            int x = i;
            while (x > 1) {
                int p = spf[x];
                int exp = 0;
                while (x % p == 0) {
                    x /= p;
                    ++exp;
                }
                required[p] += K * exp;
                long updated = minimalNForTarget(p, required[p], neededN[p]);
                neededN[p] = updated;
                if (updated > currentMax) {
                    currentMax = updated;
                }
            }
            if (i >= 10) {
                sum = sum.add(BigInteger.valueOf(currentMax));
            }
        }

        BigInteger outMod = BigInteger.valueOf(1000000000000000000L);
        return String.valueOf(sum.mod(outMod).longValue());
    }

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

Problem 320: Factorials Divisible by a Huge Integer

Project Euler Problem 320 Solution

Detailed mathematical approach

Problem Summary

Mathematical Approach

Algorithm

Complexity Analysis

Checks And Final Result

Further Reading

Problem 320 source code

C++

Python

Java