Problem 926: Total Roundness

View on Project Euler

Project Euler Problem 926 Solution

EulerSolve provides an optimized solution for Project Euler Problem 926, Total Roundness, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer \(m=\prod p_i^{a_i}\), its roundness is the largest integer \(r\ge 1\) such that \(m\) is a perfect \(r\)th power. For every nontrivial \(m\), that is exactly $$\rho(m)=\gcd(a_1,a_2,\dots).$$ The total roundness of \(N\) is defined by summing \(\rho(d)\) over every divisor \(d\mid N\) with \(d>1\): $$T(N)=\sum_{\substack{d\mid N\\ d>1}}\rho(d).$$ Problem 926 asks for \(T(10^7!)\bmod 10^9+7\). The divisor \(1\) is excluded for an important reason: it is a perfect \(k\)th power for every \(k\), so any layer-by-layer counting formula would otherwise contain an infinite extra contribution. Mathematical Approach The implementations never enumerate divisors of \(10^7!\), and they never compute \(\gcd\) values divisor by divisor. Instead, they describe every divisor through the prime exponents of the factorial, convert roundness into a sum over perfect-power layers, and then maintain the layer counts with a sweep over the breakpoints of floor quotients....

Detailed mathematical approach

Problem Summary

For a positive integer \(m=\prod p_i^{a_i}\), its roundness is the largest integer \(r\ge 1\) such that \(m\) is a perfect \(r\)th power. For every nontrivial \(m\), that is exactly

$$\rho(m)=\gcd(a_1,a_2,\dots).$$

The total roundness of \(N\) is defined by summing \(\rho(d)\) over every divisor \(d\mid N\) with \(d>1\):

$$T(N)=\sum_{\substack{d\mid N\\ d>1}}\rho(d).$$

Problem 926 asks for \(T(10^7!)\bmod 10^9+7\). The divisor \(1\) is excluded for an important reason: it is a perfect \(k\)th power for every \(k\), so any layer-by-layer counting formula would otherwise contain an infinite extra contribution.

Mathematical Approach

The implementations never enumerate divisors of \(10^7!\), and they never compute \(\gcd\) values divisor by divisor. Instead, they describe every divisor through the prime exponents of the factorial, convert roundness into a sum over perfect-power layers, and then maintain the layer counts with a sweep over the breakpoints of floor quotients.

Prime-Exponent Description of the Divisors of \(n!\)

Write

$$n!=\prod_{p\le n} p^{e_p},$$

where each exponent is given by Legendre's formula

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

Every divisor of \(n!\) is then uniquely determined by choosing exponents

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

For \(d>1\), only the positive exponents matter, so

$$\rho(d)=\gcd\{a_p:\ a_p>0\}.$$

This is the real problem-specific object behind the task: the answer depends only on the multiset of factorial exponents \(\{e_p\}\), not on the divisors themselves as standalone integers.

Counting Perfect-Power Layers Instead of Individual \(\gcd\) Values

A divisor with roundness \(g\) is a perfect \(k\)th power for exactly the integers \(k\) with \(1\le k\le g\). Equivalently,

$$\rho(d)=\sum_{k\ge 1}\mathbf{1}\!\left(d\text{ is a perfect }k\text{th power}\right).$$

Summing that identity over all nontrivial divisors gives

$$T(n!)=\sum_{k\ge 1}\#\left\{d\mid n!:\ d>1,\ d\text{ is a perfect }k\text{th power}\right\}.$$

For a fixed \(k\), the divisor \(d=\prod p^{a_p}\) is a perfect \(k\)th power exactly when each exponent is divisible by \(k\), so we can write \(a_p=k b_p\) with

$$0\le b_p\le \left\lfloor\frac{e_p}{k}\right\rfloor.$$

That makes the number of perfect \(k\)th-power divisors

$$D_k-1,\qquad D_k=\prod_{p\le n}\left(\left\lfloor\frac{e_p}{k}\right\rfloor+1\right).$$

The factor \(D_k\) counts all choices of the scaled exponents \(b_p\), including the all-zero choice that produces the divisor \(1\); the final \(-1\) removes that forbidden divisor.

Since the largest factorial exponent is \(e_2=v_2(n!)\), the sum stops at

$$E=\max_{p\le n} e_p=v_2(n!).$$

So the whole task becomes

$$\boxed{T(n!)=\sum_{k=1}^{E}\left(\prod_{p\le n}\left(\left\lfloor\frac{e_p}{k}\right\rfloor+1\right)-1\right)\pmod{10^9+7}.}$$

At \(k=1\), this product is simply the divisor count \(\tau(n!)\), so the first layer already contributes \(\tau(n!)-1\).

Compressing Equal Factorial Exponents

Rebuilding the product over all primes for every \(k\) would still be too slow, because \(n=10^7\) has hundreds of thousands of primes. The key compression used by all three implementations is to group primes by the value of their factorial exponent. Define

$$c_e=\#\left\{p\le n:\ v_p(n!)=e\right\}.$$

Then the layer count becomes

$$D_k=\prod_{e\ge 1}\left(\left\lfloor\frac{e}{k}\right\rfloor+1\right)^{c_e}.$$

Now the problem is no longer indexed by primes, but by distinct exponent values \(e\). Many primes share the same \(e\), so this representation is dramatically smaller than the raw prime list while remaining exactly equivalent.

Where the Product Changes, and Why a Sweep Works

For one fixed exponent value \(e\), only the factor

$$f_e(k)=\left\lfloor\frac{e}{k}\right\rfloor+1$$

matters. This function is piecewise constant in \(k\). If at some left endpoint \(l\) we have

$$q=\left\lfloor\frac{e}{l}\right\rfloor,$$

then the same quotient persists up to

$$r=\left\lfloor\frac{e}{q}\right\rfloor.$$

So one quotient value controls the entire interval \(l\le k\le r\). The implementations exploit exactly those interval boundaries. They start from

$$D_1=\prod_{e\ge 1}(e+1)^{c_e},$$

and whenever a group with multiplicity \(c_e\) changes from an old factor \(t_{\text{old}}\) to a new factor \(t_{\text{new}}\), the current product is updated by

$$\left(\frac{t_{\text{new}}}{t_{\text{old}}}\right)^{c_e} \pmod{10^9+7}.$$

Because the modulus is prime, those divisions are implemented with modular inverses. When \(k=e+1\), the factor for that group becomes \(1\) permanently, so no further updates from that group are needed.

A concrete one-group example makes the event logic visible. For \(e=8\), the values of \(\left\lfloor 8/k\right\rfloor+1\) are

$$9\text{ at }k=1,\quad 5\text{ at }k=2,\quad 3\text{ at }k=3,4,\quad 2\text{ at }k=5,6,7,8,\quad 1\text{ afterwards}.$$

So this single exponent group generates updates exactly at the breakpoints \(k=2,3,5,9\). The global algorithm does the same for every distinct exponent value, sorts all update events, and sweeps upward through \(k\).

Worked Example: \(10!\)

For \(n=10\), Legendre's formula gives

$$10!=2^8\cdot 3^4\cdot 5^2\cdot 7^1.$$

So the factorial exponents are \(8,4,2,1\), and the perfect-power layer counts are

$$\begin{aligned} D_1&=(8+1)(4+1)(2+1)(1+1)=270,\\ D_2&=(4+1)(2+1)(1+1)(0+1)=30,\\ D_3&=(2+1)(1+1)(0+1)(0+1)=6,\\ D_4&=(2+1)(1+1)(0+1)(0+1)=6,\\ D_5&=D_6=D_7=D_8=2. \end{aligned}$$

Therefore

$$T(10!)=(270-1)+(30-1)+(6-1)+(6-1)+4\cdot(2-1)=312.$$

This example shows the exact meaning of the layer formula. The first layer counts all nontrivial divisors, the second layer counts all nontrivial square divisors, the third layer counts all nontrivial cube divisors, and so on. Summing those layers reproduces the sum of roundness values without ever enumerating the divisors themselves.

How the Code Works

Computing the Exponent Groups

The C++, Python, and Java implementations first generate the primes up to \(n\) with a sieve. For each prime \(p\), they evaluate \(v_p(n!)\) by repeated division, which is exactly Legendre's formula in iterative form. While doing that, they build the multiplicities \(c_e\) and record the maximum exponent \(E\).

Turning Quotient Drops into Multiplicative Events

Once the groups \((e,c_e)\) are known, the implementations compute the initial product \(D_1\). Then each distinct exponent value is processed independently. Using the interval rule \(r=\lfloor e/q\rfloor\), the code locates every \(k\) where \(\lfloor e/k\rfloor\) changes and emits one event containing that sweep position and the multiplicative ratio that updates the contribution of the whole group. After all groups have emitted their events, those events are sorted by \(k\).

Sweeping \(k\) While Maintaining the Invariant \( \text{current}=D_k \)

The final sweep runs from \(k=1\) to \(k=E\). Before adding the contribution for a given \(k\), the implementation applies every event scheduled at that position. The maintained invariant is simple and exact: after processing all events at \(k\), the running product equals \(D_k\) modulo \(10^9+7\). The answer then adds \(D_k-1\), again subtracting the trivial divisor \(1\). That is the whole algorithm.

Complexity Analysis

Let \(E=v_2(n!)\), and let \(B\) be the total number of update events produced by all distinct exponent groups. The prime sieve costs \(O(n\log\log n)\) time and \(O(n)\) memory. Computing all factorial exponents costs \(O\!\left(\sum_{p\le n}\log_p n\right)\) arithmetic steps.

After that, event generation costs \(O(B)\), sorting costs \(O(B\log B)\), and the sweep costs \(O(E+B)\). In practice, \(B\) is far smaller than the naive \(E\cdot\pi(n)\) recomputation cost because each exponent group changes only at quotient breakpoints rather than at every value of \(k\). Memory usage is \(O(n+B)\), dominated by the sieve storage and the event list.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=926
  2. Legendre's formula: Wikipedia - Legendre's formula
  3. Perfect power: Wikipedia - Perfect power
  4. Divisor function: Wikipedia - Divisor function
  5. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
  6. Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse

Problem 926 source code

C++

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

namespace {

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

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

u64 add_mod(u64 a, u64 b) {
    a += b;
    if (a >= kMod) {
        a -= kMod;
    }
    return a;
}

u64 mul_mod(u64 a, u64 b) {
    return static_cast<u64>((static_cast<u128>(a) * b) % kMod);
}

u64 pow_mod(u64 base, u64 exp) {
    u64 result = 1;
    while (exp > 0) {
        if (exp & 1ULL) {
            result = mul_mod(result, base);
        }
        base = mul_mod(base, base);
        exp >>= 1ULL;
    }
    return result;
}

u64 inv_mod(u64 x) {
    return pow_mod(x, kMod - 2);
}

struct Event {
    int k;
    u64 mul;

    bool operator<(const Event& other) const {
        return k < other.k;
    }
};

std::vector<int> primes_up_to(int n) {
    std::vector<bool> composite(n + 1, false);
    std::vector<int> primes;
    primes.reserve(n / 10);

    for (int i = 2; i <= n; ++i) {
        if (!composite[i]) {
            primes.push_back(i);
            if (static_cast<i64>(i) * i <= n) {
                for (int j = i * i; j <= n; j += i) {
                    composite[j] = true;
                }
            }
        }
    }
    return primes;
}

int vp_factorial(int n, int p) {
    int e = 0;
    while (n > 0) {
        n /= p;
        e += n;
    }
    return e;
}

u64 r_from_exponents_direct(const std::vector<int>& exponents) {
    int max_e = 0;
    for (int e : exponents) {
        max_e = std::max(max_e, e);
    }

    u64 total = 0;
    for (int k = 1; k <= max_e; ++k) {
        u64 d = 1;
        for (int e : exponents) {
            d = mul_mod(d, static_cast<u64>(e / k + 1));
        }
        total = add_mod(total, d - 1);
    }
    return total;
}

std::vector<int> factor_exponents_u64(u64 n) {
    std::vector<int> exponents;
    for (u64 p = 2; p * p <= n; ++p) {
        if (n % p != 0) {
            continue;
        }
        int c = 0;
        while (n % p == 0) {
            n /= p;
            ++c;
        }
        exponents.push_back(c);
    }
    if (n > 1) {
        exponents.push_back(1);
    }
    return exponents;
}

u64 r_of_number_direct(u64 n) {
    if (n <= 1) {
        return 0;
    }
    return r_from_exponents_direct(factor_exponents_u64(n));
}

u64 r_of_factorial_optimized(int n) {
    if (n < 2) {
        return 0;
    }

    const std::vector<int> primes = primes_up_to(n);

    std::unordered_map<int, int> count_by_exp;
    count_by_exp.reserve(primes.size() * 2);

    int max_e = 0;
    for (int p : primes) {
        const int e = vp_factorial(n, p);
        ++count_by_exp[e];
        max_e = std::max(max_e, e);
    }

    std::vector<std::pair<int, int>> groups;
    groups.reserve(count_by_exp.size());
    for (const auto& kv : count_by_exp) {
        groups.push_back(kv);
    }
    std::sort(groups.begin(), groups.end());

    u64 d1 = 1;
    for (const auto& [e, c] : groups) {
        d1 = mul_mod(d1, pow_mod(static_cast<u64>(e + 1), static_cast<u64>(c)));
    }

    std::vector<Event> events;
    events.reserve(groups.size() * 300);

    for (const auto& [e, c] : groups) {
        u64 prev_term = static_cast<u64>(e + 1);

        int l = 2;
        while (l <= e) {
            const int q = e / l;
            const int r = e / q;
            const u64 term = static_cast<u64>(q + 1);
            if (term != prev_term) {
                const u64 ratio = mul_mod(term, inv_mod(prev_term));
                events.push_back({l, pow_mod(ratio, static_cast<u64>(c))});
                prev_term = term;
            }
            l = r + 1;
        }

        if (e + 1 <= max_e && prev_term != 1) {
            const u64 ratio = inv_mod(prev_term);
            events.push_back({e + 1, pow_mod(ratio, static_cast<u64>(c))});
        }
    }

    std::sort(events.begin(), events.end());

    u64 curr_d = d1;
    u64 answer = curr_d - 1;

    std::size_t idx = 0;
    for (int k = 2; k <= max_e; ++k) {
        while (idx < events.size() && events[idx].k == k) {
            curr_d = mul_mod(curr_d, events[idx].mul);
            ++idx;
        }
        answer = add_mod(answer, curr_d - 1);
    }

    return answer;
}

std::vector<int> factorial_exponents_direct(int n) {
    const std::vector<int> primes = primes_up_to(n);
    std::vector<int> exponents;
    exponents.reserve(primes.size());
    for (int p : primes) {
        exponents.push_back(vp_factorial(n, p));
    }
    return exponents;
}

void run_validations() {
    assert(r_of_number_direct(20) == 6);
    assert(r_of_factorial_optimized(10) == 312);

    for (int n : {2, 3, 5, 8, 10, 20, 50, 100}) {
        const u64 optimized = r_of_factorial_optimized(n);
        const u64 direct = r_from_exponents_direct(factorial_exponents_direct(n));
        assert(optimized == direct);
    }
}

}  // namespace

int main() {
    run_validations();
    constexpr int kN = 10'000'000;
    std::cout << r_of_factorial_optimized(kN) << '\n';
    return 0;
}

Python

import math
kMod = 1000000007

def primes_up_to(n):
    composite = [False] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if not composite[i]:
            primes.append(i)
            if i * i <= n:
                for j in range(i * i, n + 1, i):
                    composite[j] = True
    return primes

def vp_factorial(n, p):
    e = 0
    while n > 0:
        n //= p
        e += n
    return e

def r_of_factorial_optimized(n):
    if n < 2: return 0
    primes = primes_up_to(n)
    count_by_exp = {}
    max_e = 0
    for p in primes:
        e = vp_factorial(n, p)
        count_by_exp[e] = count_by_exp.get(e, 0) + 1
        max_e = max(max_e, e)
        
    groups = sorted(count_by_exp.items())
    
    d1 = 1
    for e, c in groups:
        d1 = (d1 * pow(e + 1, c, kMod)) % kMod
        
    events = []
    for e, c in groups:
        prev_term = e + 1
        l = 2
        while l <= e:
            q = e // l
            r = e // q
            term = q + 1
            if term != prev_term:
                ratio = (term * pow(prev_term, -1, kMod)) % kMod
                events.append((l, pow(ratio, c, kMod)))
                prev_term = term
            l = r + 1
        if e + 1 <= max_e and prev_term != 1:
            ratio = pow(prev_term, -1, kMod)
            events.append((e + 1, pow(ratio, c, kMod)))
            
    events.sort()
    
    curr_d = d1
    answer = curr_d - 1
    
    idx = 0
    for k in range(2, max_e + 1):
        while idx < len(events) and events[idx][0] == k:
            curr_d = (curr_d * events[idx][1]) % kMod
            idx += 1
        answer = (answer + curr_d - 1) % kMod
        
    return answer

def solve():
    return str(r_of_factorial_optimized(10000000))

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

Java

import java.util.*;

public class Euler926 {

    static final long kMod = 1000000007L;

    static long powMod(long base, long exp) {
        long res = 1;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * base) % kMod;
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return res;
    }

    static long invMod(long x) {
        return powMod(x, kMod - 2);
    }

    static List<Integer> primesUpTo(int n) {
        boolean[] isComp = new boolean[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (!isComp[i]) {
                primes.add(i);
                if ((long) i * i <= n) {
                    for (int j = i * i; j <= n; j += i)
                        isComp[j] = true;
                }
            }
        }
        return primes;
    }

    static int vpFactorial(int n, int p) {
        int e = 0;
        while (n > 0) {
            n /= p;
            e += n;
        }
        return e;
    }

    static class Event implements Comparable<Event> {
        int k;
        long mul;

        Event(int k, long mul) {
            this.k = k;
            this.mul = mul;
        }

        public int compareTo(Event o) {
            return Integer.compare(this.k, o.k);
        }
    }

    static class Pair {
        int e, c;

        Pair(int e, int c) {
            this.e = e;
            this.c = c;
        }
    }

    static long rOfFactorialOptimized(int n) {
        if (n < 2)
            return 0;

        List<Integer> primes = primesUpTo(n);
        Map<Integer, Integer> countByExp = new HashMap<>();
        int maxE = 0;
        for (int p : primes) {
            int e = vpFactorial(n, p);
            countByExp.put(e, countByExp.getOrDefault(e, 0) + 1);
            maxE = Math.max(maxE, e);
        }

        List<Pair> groups = new ArrayList<>();
        for (Map.Entry<Integer, Integer> entry : countByExp.entrySet()) {
            groups.add(new Pair(entry.getKey(), entry.getValue()));
        }
        Collections.sort(groups, (a, b) -> Integer.compare(a.e, b.e));

        long d1 = 1;
        for (Pair p : groups) {
            d1 = (d1 * powMod(p.e + 1, p.c)) % kMod;
        }

        List<Event> events = new ArrayList<>();
        for (Pair p : groups) {
            int e = p.e;
            int c = p.c;
            long prevTerm = e + 1;

            int l = 2;
            while (l <= e) {
                int q = e / l;
                int r = e / q;
                long term = q + 1;
                if (term != prevTerm) {
                    long ratio = (term * invMod(prevTerm)) % kMod;
                    events.add(new Event(l, powMod(ratio, c)));
                    prevTerm = term;
                }
                l = r + 1;
            }

            if (e + 1 <= maxE && prevTerm != 1) {
                long ratio = invMod(prevTerm);
                events.add(new Event(e + 1, powMod(ratio, c)));
            }
        }

        Collections.sort(events);

        long currD = d1;
        long answer = currD - 1;

        int idx = 0;
        for (int k = 2; k <= maxE; ++k) {
            while (idx < events.size() && events.get(idx).k == k) {
                currD = (currD * events.get(idx).mul) % kMod;
                idx++;
            }
            answer = (answer + currD - 1 + kMod) % kMod;
        }

        return answer;
    }

    public static String solve() {
        return Long.toString(rOfFactorialOptimized(10000000));
    }

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