Problem 646: Bounded Divisors

Project Euler Problem 646 Solution

EulerSolve provides an optimized solution for Project Euler Problem 646, Bounded Divisors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write the factorial as $$n!=\prod_{p\le n}p^{a_p},\qquad a_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$$ The quantity to compute is the signed divisor sum $$S(n;L,H)=\sum_{\substack{d\mid n!\\L\le d\le H}}\lambda(d)\,d,\qquad \lambda(d)=(-1)^{\Omega(d)},$$ reduced modulo $10^9+7$. The interval $[L,H]$ is enormous in the target case, so the solution cannot enumerate divisors naively; it must exploit the multiplicative structure of $n!$ and the fact that only divisors up to $H$ can matter. Mathematical Approach The main idea is to represent each divisor by its prime exponents, split those primes into two independent groups, and turn the interval condition into two prefix queries. Step 1: Prime Exponents of the Factorial For every prime $p\le n$, Legendre's formula gives the exponent $a_p$ of $p$ in $n!$. Therefore every divisor of $n!$ has a unique form $$d=\prod_{p\le n}p^{e_p},\qquad 0\le e_p\le a_p.$$ So the problem is equivalent to iterating over all allowable exponent vectors $(e_p)$, but doing so intelligently enough that the upper bound $H$ removes most impossible branches before they are fully built....

Detailed mathematical approach

Problem Summary

Write the factorial as

$$n!=\prod_{p\le n}p^{a_p},\qquad a_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$$

The quantity to compute is the signed divisor sum

$$S(n;L,H)=\sum_{\substack{d\mid n!\\L\le d\le H}}\lambda(d)\,d,\qquad \lambda(d)=(-1)^{\Omega(d)},$$

reduced modulo $10^9+7$. The interval $[L,H]$ is enormous in the target case, so the solution cannot enumerate divisors naively; it must exploit the multiplicative structure of $n!$ and the fact that only divisors up to $H$ can matter.

Mathematical Approach

The main idea is to represent each divisor by its prime exponents, split those primes into two independent groups, and turn the interval condition into two prefix queries.

Step 1: Prime Exponents of the Factorial

For every prime $p\le n$, Legendre's formula gives the exponent $a_p$ of $p$ in $n!$. Therefore every divisor of $n!$ has a unique form

$$d=\prod_{p\le n}p^{e_p},\qquad 0\le e_p\le a_p.$$

So the problem is equivalent to iterating over all allowable exponent vectors $(e_p)$, but doing so intelligently enough that the upper bound $H$ removes most impossible branches before they are fully built.

Step 2: The Signed Weight is Completely Multiplicative on Prime Powers

If $d=\prod p^{e_p}$, then

$$\lambda(d)\,d=(-1)^{\Omega(d)}\prod p^{e_p}=\prod (-p)^{e_p}.$$

Without the interval restriction, the total signed sum over all divisors would factor as

$$\sum_{d\mid n!}\lambda(d)\,d=\prod_{p\le n}\left(\sum_{e=0}^{a_p}(-p)^e\right).$$

The interval $L\le d\le H$ destroys that simple product formula, because different primes now interact through the size of the full divisor. That is exactly why the implementation switches to a meet-in-the-middle strategy.

Step 3: Split the Prime Factors into Two Independent Groups

Partition the prime powers of $n!$ into two disjoint groups. For the first group define partial states

$$u=\prod_{p\in A}p^{e_p},\qquad \alpha=\prod_{p\in A}(-p)^{e_p},$$

and for the second group define

$$v=\prod_{p\in B}p^{e_p},\qquad \beta=\prod_{p\in B}(-p)^{e_p}.$$

Because the groups use disjoint primes, every divisor corresponds to exactly one pair of partial states, with

$$d=u\,v,\qquad \lambda(d)\,d=\alpha\beta.$$

This reduces one very large search over all prime exponents to a merge of two precomputed lists.

Step 4: Enumerate Only Partial Products that Stay at Most $H$

During recursion, if the current partial product already exceeds $H$, or if multiplying by one more copy of the current prime would exceed $H$, then no continuation can ever contribute to a divisor in $[L,H]$ or to a prefix sum up to $H$. Every extra prime power only makes the divisor larger.

So each side generates only feasible pairs $(u,\alpha)$ or $(v,\beta)$ with numeric value at most $H$. The sign is carried through the factor $-p$, while the divisor magnitude itself is stored exactly for later ordering and comparison.

Step 5: Turn the Interval into Prefix Queries

Define the bounded prefix function

$$P(X)=\sum_{\substack{d\mid n!\\d\le X}}\lambda(d)\,d.$$

After splitting, this becomes

$$P(X)=\sum_{(u,\alpha)}\alpha\left(\sum_{\substack{(v,\beta)\\v\le X/u}}\beta\right).$$

If the second list is sorted by $v$ and its signed weights are accumulated into prefix sums, then the inner quantity is available immediately once we know how far the threshold $X/u$ reaches. Sorting the first list as well allows a monotone sweep: as $u$ grows, the allowed bound $X/u$ shrinks, so one pointer moves only in one direction.

Step 6: Recover the Required Interval Sum

The target interval sum is just the difference of two prefixes:

$$S(n;L,H)=P(H)-P(L-1)\pmod{10^9+7}.$$

That identity is the final reduction. The algorithm therefore computes two threshold queries, one at $H$ and one at $L-1$, and subtracts them modulo $10^9+7$.

Worked Example: $n=5$, $L=6$, $H=30$

Here

$$5!=2^3\cdot 3\cdot 5.$$

Split the prime factors into $A=\{2^3,3\}$ and $B=\{5\}$. The first side produces the sorted partial states

$$\{(1,1),(2,-2),(3,-3),(4,4),(6,6),(8,-8),(12,-12),(24,24)\}.$$

The second side produces

$$\{(1,1),(5,-5)\},$$

whose prefix sums of weights are $1$ and $-4$.

For $P(30)$, the bound $30/u$ is at least $5$ when $u=1,2,3,4,6$, so the inner sum is $-4$. For $u=8,12,24$, only the state $v=1$ is allowed, so the inner sum is $1$. Thus

$$P(30)=-4+8+12-16-24-8-12+24=-20.$$

For $P(5)$, only $u=1,2,3,4$ contribute, giving

$$P(5)=-4-2-3+4=-5.$$

Therefore

$$S(5;6,30)=P(30)-P(5)=-15\equiv 10^9+7-15 \pmod{10^9+7}.$$

This small example is exactly the same merge logic used in the full computation.

How the Code Works

The C++, Python, and Java implementations first generate all primes up to $n$ and compute their exponents in $n!$ by repeated division. They then split the resulting prime-factor list after the first five primes, which keeps both recursive enumerations manageable for the target instance.

Each side recursively enumerates every partial divisor whose exact value does not exceed $H$. Alongside the exact divisor value, the implementation carries the signed multiplicative weight modulo $10^9+7$. The divisor values themselves use exact wide integers: 256-bit arithmetic in C++ and arbitrary-precision integers in Python and Java, so comparisons against $10^{60}$ remain exact.

After sorting both partial-state lists by divisor value, the implementation builds prefix sums of the signed weights on the second list. A linear sweep over the first list then evaluates $P(H)$ and $P(L-1)$. If $L\le 1$, the lower prefix is zero; otherwise the algorithm evaluates the second threshold normally and subtracts the two results modulo $10^9+7$.

Complexity Analysis

Generating primes up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory. Computing the exponents $a_p$ by repeated division contributes $O(\pi(n)\log n)$, which is small here. If the two recursive enumerations produce $A$ and $B$ feasible partial states after pruning, then enumeration costs $O(A+B)$, sorting costs $O(A\log A+B\log B)$, and each threshold sweep costs $O(A+B)$ because the pointer over the second list moves monotonically. The memory usage is $O(A+B)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=646
Legendre's formula: Wikipedia — Legendre's formula
Liouville function: Wikipedia — Liouville function
Prefix sum: Wikipedia — Prefix sum
Divisor: Wikipedia — Divisor

Problem 646 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using u64 = std::uint64_t;
using boost::multiprecision::uint256_t;

constexpr u64 kMod = 1'000'000'007ULL;

struct Factor {
    int p;
    int a;
};

struct Elem {
    uint256_t v;
    u64 w;
};

uint256_t pow10_u256(int e) {
    uint256_t x = 1;
    for (int i = 0; i < e; ++i) x *= 10;
    return x;
}

std::vector<int> primes_upto(int n) {
    std::vector<char> is_prime(static_cast<std::size_t>(n + 1), 1);
    if (n >= 0) is_prime[0] = 0;
    if (n >= 1) is_prime[1] = 0;
    for (int p = 2; (long long)p * p <= n; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) continue;
        for (int m = p * p; m <= n; m += p) is_prime[static_cast<std::size_t>(m)] = 0;
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

std::vector<Factor> factorial_factors(int n) {
    std::vector<Factor> out;
    const auto primes = primes_upto(n);
    for (int p : primes) {
        int a = 0;
        for (int x = n; x; x /= p) a += x / p;
        out.push_back(Factor{p, a});
    }
    return out;
}

u64 mod_mul(u64 a, u64 b) { return (a * b) % kMod; }

void gen_divs(const std::vector<Factor>& fac, const int idx, const uint256_t& cap, uint256_t v, u64 w,
              std::vector<Elem>& out) {
    if (idx == (int)fac.size()) {
        out.push_back(Elem{v, w});
        return;
    }
    const int p = fac[(std::size_t)idx].p;
    const int a = fac[(std::size_t)idx].a;
    const u64 negp = (kMod - static_cast<u64>(p) % kMod) % kMod;
    uint256_t vv = v;
    u64 ww = w;
    for (int e = 0; e <= a; ++e) {
        gen_divs(fac, idx + 1, cap, vv, ww, out);
        if (e == a) break;
        if (vv > cap / static_cast<unsigned>(p)) break;
        vv *= static_cast<unsigned>(p);
        ww = mod_mul(ww, negp);
    }
}

u64 prefix_sum(const uint256_t& X, const std::vector<Elem>& A, const std::vector<Elem>& B,
               const std::vector<u64>& prefB) {
    std::size_t j = B.size();
    u64 ans = 0;
    for (const auto& a : A) {
        if (a.v > X) break;
        const uint256_t lim = X / a.v;
        while (j > 0 && B[j - 1].v > lim) --j;
        ans += mod_mul(a.w, prefB[j]);
        ans %= kMod;
    }
    return ans;
}

u64 S_factorial_bounded(const int n, const uint256_t& L, const uint256_t& H) {
    const std::vector<Factor> fac = factorial_factors(n);
    const int split = std::min<int>(5, static_cast<int>(fac.size()));
    const std::vector<Factor> fa(fac.begin(), fac.begin() + split);
    const std::vector<Factor> fb(fac.begin() + split, fac.end());

    u64 cntA = 1;
    for (const auto& f : fa) cntA *= static_cast<u64>(f.a + 1);
    u64 cntB = 1;
    for (const auto& f : fb) cntB *= static_cast<u64>(f.a + 1);

    std::vector<Elem> A;
    std::vector<Elem> B;
    A.reserve(static_cast<std::size_t>(cntA));
    B.reserve(static_cast<std::size_t>(cntB));

    gen_divs(fa, 0, H, 1, 1, A);
    gen_divs(fb, 0, H, 1, 1, B);

    std::sort(A.begin(), A.end(), [](const Elem& x, const Elem& y) { return x.v < y.v; });
    std::sort(B.begin(), B.end(), [](const Elem& x, const Elem& y) { return x.v < y.v; });

    std::vector<u64> prefB(B.size() + 1, 0);
    for (std::size_t i = 0; i < B.size(); ++i) prefB[i + 1] = (prefB[i] + B[i].w) % kMod;

    const u64 sumH = prefix_sum(H, A, B, prefB);
    const uint256_t Lm1 = L - 1;
    const u64 sumL = (L == 0 || L == 1) ? 0 : prefix_sum(Lm1, A, B, prefB);
    return (sumH + kMod - sumL) % kMod;
}

}  // namespace

int main() {
    assert(S_factorial_bounded(10, 100, 1000) == 1457);
    assert(S_factorial_bounded(15, 1000, 100000) == (kMod - 107'974ULL));
    assert(S_factorial_bounded(30, 100'000'000, 1'000'000'000'000ULL) == (9766732243224ULL % kMod));

    const uint256_t L = pow10_u256(20);
    const uint256_t H = pow10_u256(60);
    std::cout << S_factorial_bounded(70, L, H) << "\n";
    return 0;
}

Python

import sys

sys.setrecursionlimit(2000)
MOD = 1000000007

class Factor:
    def __init__(self, p, a):
        self.p = p
        self.a = a

class Elem:
    def __init__(self, v, w):
        self.v = v
        self.w = w

def primes_upto(n):
    is_prime = bytearray(n + 1)
    for i in range(len(is_prime)): is_prime[i] = 1
    if n >= 0: is_prime[0] = 0
    if n >= 1: is_prime[1] = 0
    for p in range(2, int(n**0.5) + 1):
        if is_prime[p]:
            for m in range(p * p, n + 1, p):
                is_prime[m] = 0
    return [i for i in range(2, n + 1) if is_prime[i]]

def factorial_factors(n):
    out = []
    primes = primes_upto(n)
    for p in primes:
        a = 0
        x = n
        while x:
            x //= p
            a += x
        out.append(Factor(p, a))
    return out

def mod_mul(a, b):
    return (a * b) % MOD

def gen_divs(fac, idx, cap, v, w, out):
    if idx == len(fac):
        out.append(Elem(v, w))
        return
    p = fac[idx].p
    a = fac[idx].a
    negp = (MOD - (p % MOD)) % MOD
    vv = v
    ww = w
    for e in range(a + 1):
        gen_divs(fac, idx + 1, cap, vv, ww, out)
        if e == a: break
        if vv > cap // p: break
        vv *= p
        ww = mod_mul(ww, negp)

def prefix_sum(X, A, B, prefB):
    j = len(B)
    ans = 0
    for a in A:
        if a.v > X: break
        lim = X // a.v
        while j > 0 and B[j - 1].v > lim:
            j -= 1
        ans = (ans + mod_mul(a.w, prefB[j])) % MOD
    return ans

def S_factorial_bounded(n, L, H):
    fac = factorial_factors(n)
    split = min(5, len(fac))
    fa = fac[:split]
    fb = fac[split:]

    A = []
    B = []
    gen_divs(fa, 0, H, 1, 1, A)
    gen_divs(fb, 0, H, 1, 1, B)

    A.sort(key=lambda x: x.v)
    B.sort(key=lambda x: x.v)

    prefB = [0] * (len(B) + 1)
    for i in range(len(B)):
        prefB[i + 1] = (prefB[i] + B[i].w) % MOD

    sumH = prefix_sum(H, A, B, prefB)
    Lm1 = L - 1
    sumL = 0 if L <= 1 else prefix_sum(Lm1, A, B, prefB)
    return (sumH + MOD - sumL) % MOD

def solve():
    L = 10**20
    H = 10**60
    ans = S_factorial_bounded(70, L, H)
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;

public class Euler646 {
    static final long kMod = 1000000007L;

    static class Factor {
        int p, a;

        Factor(int p, int a) {
            this.p = p;
            this.a = a;
        }
    }

    static class Elem {
        BigInteger v;
        long w;

        Elem(BigInteger v, long w) {
            this.v = v;
            this.w = w;
        }
    }

    static ArrayList<Integer> primesUpto(int n) {
        byte[] is_prime = new byte[n + 1];
        for (int i = 0; i <= n; i++)
            is_prime[i] = 1;
        if (n >= 0)
            is_prime[0] = 0;
        if (n >= 1)
            is_prime[1] = 0;
        for (int p = 2; (long) p * p <= n; ++p) {
            if (is_prime[p] == 1) {
                for (int m = p * p; m <= n; m += p)
                    is_prime[m] = 0;
            }
        }
        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; ++i) {
            if (is_prime[i] == 1)
                primes.add(i);
        }
        return primes;
    }

    static ArrayList<Factor> factorialFactors(int n) {
        ArrayList<Factor> out = new ArrayList<>();
        ArrayList<Integer> primes = primesUpto(n);
        for (int p : primes) {
            int a = 0;
            for (int x = n; x > 0; x /= p)
                a += x / p;
            out.add(new Factor(p, a));
        }
        return out;
    }

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

    static void genDivs(ArrayList<Factor> fac, int idx, BigInteger cap, BigInteger v, long w, ArrayList<Elem> out) {
        if (idx == fac.size()) {
            out.add(new Elem(v, w));
            return;
        }
        int p = fac.get(idx).p;
        int a = fac.get(idx).a;
        long negp = (kMod - (long) p % kMod) % kMod;
        BigInteger vv = v;
        long ww = w;
        BigInteger bigP = BigInteger.valueOf(p);
        for (int e = 0; e <= a; ++e) {
            genDivs(fac, idx + 1, cap, vv, ww, out);
            if (e == a)
                break;
            if (vv.compareTo(cap.divide(bigP)) > 0)
                break;
            vv = vv.multiply(bigP);
            ww = modMul(ww, negp);
        }
    }

    static long prefixSum(BigInteger X, ArrayList<Elem> A, ArrayList<Elem> B, long[] prefB) {
        int j = B.size();
        long ans = 0;
        for (Elem a : A) {
            if (a.v.compareTo(X) > 0)
                break;
            BigInteger lim = X.divide(a.v);
            while (j > 0 && B.get(j - 1).v.compareTo(lim) > 0)
                --j;
            ans += modMul(a.w, prefB[j]);
            ans %= kMod;
        }
        return ans;
    }

    static long sFactorialBounded(int n, BigInteger L, BigInteger H) {
        ArrayList<Factor> fac = factorialFactors(n);
        int split = Math.min(5, fac.size());
        ArrayList<Factor> fa = new ArrayList<>(fac.subList(0, split));
        ArrayList<Factor> fb = new ArrayList<>(fac.subList(split, fac.size()));

        ArrayList<Elem> A = new ArrayList<>();
        ArrayList<Elem> B = new ArrayList<>();

        genDivs(fa, 0, H, BigInteger.ONE, 1, A);
        genDivs(fb, 0, H, BigInteger.ONE, 1, B);

        Comparator<Elem> cmp = new Comparator<Elem>() {
            public int compare(Elem x, Elem y) {
                return x.v.compareTo(y.v);
            }
        };

        Collections.sort(A, cmp);
        Collections.sort(B, cmp);

        long[] prefB = new long[B.size() + 1];
        for (int i = 0; i < B.size(); ++i) {
            prefB[i + 1] = (prefB[i] + B.get(i).w) % kMod;
        }

        long sumH = prefixSum(H, A, B, prefB);
        BigInteger Lm1 = L.subtract(BigInteger.ONE);
        long sumL = (L.compareTo(BigInteger.ONE) <= 0) ? 0 : prefixSum(Lm1, A, B, prefB);
        return (sumH + kMod - sumL) % kMod;
    }

    public static String solve() {
        BigInteger L = BigInteger.TEN.pow(20);
        BigInteger H = BigInteger.TEN.pow(60);
        long ans = sFactorialBounded(70, L, H);
        return Long.toString(ans);
    }

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