Algorithm to express n! as product of powers of prime

algorithm, c

Solution

The most efficient way to find the prime factorisation of `n!` that I know of is counting how often each prime appears as a factor in `n!`. Obviously, no prime larger than `n` appears, so let `p <= n` be a prime.

Among the numbers `1, ..., n`, there are

q1 = floor(n/p)

multiples of `p`. Among these, there are

q2 = floor(q1/p) = floor(n/p²)

multiples of `p²`. Among these, there are

q3 = floor(q2/p) = floor(n/p³)

multiples of `p³`. Etc. So the exponent of `p` in the prime factorisation of `n!` is

q1 + q2 + q3 + ...

(An `a = p^k*b`, with `b` not divisible by `p` contributes `k` to the exponent, and appears in the `k` lists corresponding to the counts `q1, ..., qk`.) We can succinctly write a function for that:

unsigned long long factorial_exponent(unsigned long long n, unsigned long long p) {
    unsigned long long exponent = 0;
    do {
        n /= p;
        exponent += n;
    }while(n);
    return exponent;
}

That uses `floor(log n/log p) + 1` divisions, so, if the primes not exceeding `n` are known, that contributes approximately

log n * ∑ (1/log p + 1) ≈ 2n/log n
       p≤n

divisions and additions to the total work. (Note: since most of the primes `≤ n` are `> √n`, and for primes `p > √n` obviously `q2 = 0`, it is faster to calculate their exponent directly: `n/p`, that reduces the number of divisions needed by about half.)

Finding the primes not exceeding `n` is best done with a sieve, if you already have a good implementation, the Sieve of Atkin does it in O(n) or O(n/log log n) operations (depends on the implementation, but for feasible ranges, `log log n` can be considered a constant), otherwise use the Sieve of Eratosthenes, it's simple to implement and finds the primes not exceeding `n` in - again depending on the implementation - O(n*log log n) or O(n) operations.

The total work for that algorithm is dominated by finding the primes (but for feasible `n`, the contribution of the determination of the exponents is still not negligible).

On the other hand, the work needed for finding the prime factorisation of each `k ≤ n` of course depends on the algorithm used for that. Using trial division for that would result in a total work of about `c*n^1.5/log n` - I haven't done anything to determine the constant `c`, and depending on details, you may have a factor of `log n` in the numerator or denominator, but it's basically `n^1.5`. A better method of finding the factorisations would be first finding the smallest prime factor [or any prime factor] with a modification of the Sieve of Eratosthenes, again in O(n*log log n) operations, and then use that to find the factorisation. You can store the factorisations and then, when processing `k` with the known prime divisor `p`, look up the factorisation of `k/p`, or generate the factorisation on the fly by recursively looking up the known prime factor `q` of `k/p`, then of `k/(p*q)` etc. until the factorisation is complete - that is much simpler to handle if the known prime factor is always the smallest.

On average, the prime factorisation of `k` contains `≈ log log k` terms, so that method would lead to O(n*log log n) overall complexity. But the constant factors in this method are considerably larger than in the first, so even if the prime-finding gives the same complexity, the first is faster.

Problem

Original source