Project Euler #549: Divisibility of factorials

Question

This is the problem:

Calculate

$$\sum_{i=2}^{10^8} s(i)$$

where \$s(n)\$ is the smallest \$m\$ such that \$n\$ divides \$m!\$.

Quite mathematical, I've found a better way than brute force by using largest prime factor and greatest common divisor. It works fine when n = 100 but still very slow when n gets larger and far too slow to solve the problem. So of course there are better ways.

Here is the code in Python:

import math

def maxPrimeFactor(n):
    p = 2
    while (p <= n/p):
        if (n%p):
            p += 1
        else:
            n/=p
    return n

def gcd(a,b):
    c = 1
    while (b):
        c = b
        b = a % b
        a = c
    return a

def sn(n):
    solution = prod = 1
    p = maxPrimeFactor(n)
    pFac = math.factorial(p)
    if((pFac%n) == 0 and pFac>=n):
        solution = p
    else:
        rest = n / gcd(pFac, n)
        solution = p+1
        prod = p + 1
        while(prod < rest):
            solution += 1
            prod *= solution
        while (prod % rest):
            solution += 1
            prod *= solution
    return solution

sum = 0
for i in range(2,100000001):
    sum += sn(i)
print sum

Answer 1

Caveat: at the moment I have neither the time nor the knowledge for solving Problem Euler #549 properly - all I can offer is more effective ways of applying brute force.

The biggest obstacle to solving this problem by brute force using factorials is the sheer size of the numbers involved. Stirling's approximation says that 100000000! has more than 2 * 10^9 bits; with a good MP library like GMP such numbers might be just barely tractable but it's going to be painfully slow.

Moreover, these factorials would need to be factored, which under the circumstances would boil down to trial division by up to pi(10^8) = 5,761,455 primes - and this would give a whole new meaning to the phrase 'painfully slow'.

A much better way is to avoid factorials entirely, working (almost) exlusively in the realm of factors. The most important insight - on which all solutions presented here are based - is this:

s(p^k * m) = max(s(p^k), s(m)) if gcd(p, m) == 1

In other words (and formulated without recursion): given the factorisation of a number n, the value s(n) is the maximum s(p^k) over all distinct prime factors p.

This gives a viable recipe for computing s(n), since factorisation of small integers is trivial and the computation of s(p^k) is tractable. Here's a rendering of the latter that's not very mathematical but easy to verify and code (using C# for clarity, since I don't speak Python well enough):

// computes s(p^e) for p prime
static int s_for_prime_power (int p, int e)
{
    int k = 0;

    while (e > p)
    {
        k += p;
        e -= p + 1;

        for (int t = k; (t /= p) % p == 0; )
            --e;
    }

    return (k + Math.Max(0 , e)) * p;
}

With the help of that function the computation of s() is easily formulated:

static int s (int n)
{
    return s(factorisation(n));
}

static int s (List<int> factors, int lo = 0)
{
    int p = factors[lo], hi = lo;

    while (++hi < factors.Count && factors[hi] == p)
        ;

    int result = s_for_prime_power(p, hi - lo);

    if (hi < factors.Count)
        result = Math.Max(result, s(factors, hi));

    return result;
}

This function implements the recursion shown earlier: it determines the power of the least prime factor (by counting the number of occurrences), computes the s value for that, and then returns the maximum of that value and the s value for the rest.

I'm glossing over some unimportant things - like factorisation and so on - since Python has everything that's needed; the C# file is Euler0549_simple.linq on PasteBin. This is a script for use with the (free) LINQPad with which I do almost all of my algorithmic prototyping and Eulering/SPOJing etc. But the file is also valid .cs that can be compiled and run unchanged with Mono, for example.

This code performs suprisingly well (considering), which is nice because that makes it easier to test more advanced solutions:

# benching 'Euler0549_simple' ...

S(       10):                 39        0,0 ms
S(      100):               2012        0,0 ms
S(     1000):             136817        0,1 ms
S(    10000):           10125843        1,3 ms
S(   100000):          793183093       16,1 ms
S(  1000000):        64938007616      216,5 ms
S( 10000000):      5494373412573    3.538,4 ms
S(100000000):    #45a57e1327a02a   66.756,3 ms

In the last row some hex digits from the MD5 of the decimal result are printed, in order to avoid showing the actual result. The increase in time is higher than linear because the factorisations get longer as numbers get bigger. One major point of inefficiency here is the fact that numbers are enumerated as numbers and then factored, instead of enumerating factorisations directly.

There is a type an algorithm that effectively enumerates an initial segment of natural numbers in the form of least factor decompositions p * m where p is not greater than the least prime factor of m; see Linear time Euler's Totient Function calculation for example.

This algorithm is ideally suited for computing tables of functions where

f(p * m) = g(f(p), f(m))

which is - almost - the case here. 'Almost' because we need a decomposition of the form p^k * m with gcd(p, m) == 1; however, for a first draft this can be arranged by the simple expedient of dividing p out of m if necessary.

static int[] s;

static long E549_via_lpf_decomposition (int limit)
{
    int small_prime_limit = (int)Math.Sqrt(limit);
    var small_primes = new List<int>();
    int m = 2;

    s = new int[limit + 1];

    for (int half = limit / 2; m <= half; ++m)
    {
        if (s[m] == 0)  // prime
        {
            s[m] = m;

            if (m <= small_prime_limit)
                small_primes.Add(m);
        }

        int s_m = s[m], threshold = limit / m;

        foreach (int p in small_primes)
        {
            if (p > threshold)  // i.e. if (p * m > limit)
                break;

            if (m % p != 0)  // not a factor of m yet
            {
                // p is a new least factor -> s(p) < s(m) -> s(p * m) = s(m)
                s[p * m] = s_m;             
            }
            else // p is already a (least) factor of m
            {
                int e = 2, q = m;

                while ((q /= p) % p == 0)
                    ++e;

                s[p * m] = Math.Max(s_for_prime_power(p, e), s[q]);
                break;
            }
        }
    }

    for ( ; m <= limit; ++m)
        if (s[m] == 0)
            s[m] = m;

    long sum = 0;
    for (int i = 2; i <= limit; ++i)
        sum += s[i];

    return sum;
}

This can be made more efficient in various ways but I think in this form it nicely shows the logic of the algorithm. The timings are much more linear in n, since the algorithm itself is linear. Only s_for_prime_power() and my fumbling with factoring out powers of the lpf add a tiny logarithmic influence on top:

# benching 'Euler0549_lpf_decomp' ...

S(       10):                 39        0,0 ms
S(      100):               2012        0,0 ms
S(     1000):             136817        0,0 ms
S(    10000):           10125843        0,3 ms
S(   100000):          793183093        1,9 ms
S(  1000000):        64938007616       19,5 ms
S( 10000000):      5494373412573      196,1 ms
S(100000000):    #45a57e1327a02a    1.990,3 ms

Although this is based on raw firepower instead of mathematics, it is almost decent already... Another option would be a modified Sieve of Eratosthenes (using the s[] vector like a sieve, instead of of a bit vector). Or, the devious application of advanced mathematics, but that is well beyond my capabilities ATM. ;-)