7
\$\begingroup\$

The problem

I have lately been working Project Euler: 231:

The prime factorisation of binomial coefficients

The binomial coefficient \$ ^{10}C_3 = 120 \$.

\$ 120 = 2^3 × 3 × 5 = 2 × 2 × 2 × 3 × 5 \$, and \$ 2 + 2 + 2 + 3 + 5 = 14 \$.

So the sum of the terms in the prime factorisation of \$^{10}C_3\$ is 14.

Find the sum of the terms in the prime factorisation of \$ ^{20000000}C_{15000000} \$.

Naive approach

Using some built in functions I was able to hack together a brute-force solution fairly quickly:

from primefac import primefac

def slow_factors(n, k):
    total = 0
    nom = range(n - k+1, n+1)
    denom = range(2, k+1)
    for num in nom:
        total += sum(primefac(num))
    for num in denom:
        total -= sum(primefac(num))
    return total

This solution uses that:

\$ ^n C_k = \binom{n}{k} = \frac{n!}{k!(n-k)!} = \frac{\overbrace{n \cdot (n-1) \cdots p}^{k \text{ times}}}{k\cdot (k-1) \cdots 2} \$

For a concrete example of my algorithm:

\$ \binom{10}{3} = \frac{10!}{3!(10-3)!} = \frac{10 \cdot 9 \cdot 8}{3 \cdot 2 \cdot 1} \$

This is respectively nominator and denominator in my code. From here I simply calculate the prime factorization of every number in the denominator. The answer is this minus the sum of the prime factorization of the denominator. I then realized that:

\$ \binom{20000000}{15000000} \$

Is incredibly large! Wolfram Alpha says that it is approximately \$2 \approx 10^{4884377}\$, far bigger than \$9000! \approx 10^{31681}\$. So I was in need of a smarter solution.

Some improvements

The first improvements uses that:

\$ \binom{n}{k} = \binom{n}{n-k} \$

We compute values \$n \cdot (n-1) \cdots\$ until \$k\$. However because of the symmetry we could also compute \$n \cdot (n-1) \cdots\$ until \$n-k\$. We choose based on which gives us the fewest values to compute:

k = min(k, n - k)

It's also inefficient because it calculates the prime factorization of each number separately. My new idea was to make a dict containing a list of sums of the factors under some LIMIT. Then for new numbers I would divide them by primes until I found a match in my list. This is what the factor function does.

Now my code can calculate:

\$ \binom{200000}{150000} \$

In a little over 2 minutes. Since this is 100 times smaller than what was asked for in the problem. I estimated that my code would use well over 3 hours to finish. This is of course still too slow. Are there any other obvious speed improvements I have forgotten? I really think my idea is good, although my code surely could be improved.

Question

I'm not looking for feedback on the function slow_factors this is only included to test the accuracy of the faster function. I'm mainly looking for speed improvements, as well as smarter approaches.

from primefac import primefac, isprime
from collections import Counter
from primesieve import generate_n_primes as primes

LIMIT = 10**4
PRIMS = primes(1000)
sum_factors = Counter()
for i in range(LIMIT):
    sum_factors[i] = sum(primefac(i))


def factors(num):
    primelist = []
    for prime in PRIMS:
        while num % prime == 0:
            primelist.append(prime)
            num //= prime
            if sum_factors[num] > 0:
                val = sum_factors[num]
                return num, val, primelist
    return num, sum(primefac(num)), primelist


def factor_lst(lst):
    total = 0
    for num in lst:
        if sum_factors[num] > 0:
            total += sum_factors[num]
        elif isprime(num) == True:
            sum_factors[num] = num
            total += num
        else:
            # Reiterates over all the newly found nums
            in_factors, val, primelist = factors(num)
            for p in primelist:
                in_factors *= p
                val += p
                sum_factors[in_factors] = val
            total += val
    return total


def binom_factor(n, k):
    total = 0
    k = min(k, n - k)
    nom = xrange(n - k+1, n+1)
    denom = xrange(2, k+1)
    nom_sum = factor_lst(nom)
    denom_sum = factor_lst(denom)
    return nom_sum - denom_sum


def slow_factors(n, k):
    total = 0
    k = min(k, n - k)
    nom = range(n - k+1, n+1)
    denom = range(2, k+1)
    for num in nom:
        total += sum(primefac(num))
    for num in denom:
        total -= sum(primefac(num))
    return total


if __name__ == '__main__':

    print binom_factor(76430, 4321)
    print slow_factors(76430, 4321)
    power = 4
    print binom_factor(20*10**power, 15*10**power)
\$\endgroup\$
4
\$\begingroup\$

The strategy in the post is to iterate over all the numbers that are multiplied together in the numerator and denominator of the binomial coefficient, and factorize each one. But this is very time-consuming, because no-one knows an efficient algorithm for factorization. We would much rather work out the factorization in some other way.

A binomial coefficient is computed using factorials: $${n \choose k} = {n! \over k!(n-k)!}$$ So let's have a look at some factorizations of small factorials: $$\eqalign{2! &= 2 \\ 3! &= 2·3 \\ 4! &= 2^3·3 \\ 5! &= 2^3·3·5 \\ 6! &= 2^4·3^2·5 }$$ It should be clear that the factorization of \$n!\$ includes every prime \$≤n\$. So instead of iterating over the numbers up to \$n\$ and trying to factorize them, we can iterate over the primes \$p ≤ n\$ and count how many times each one appears in the factorization. We can iterate over the primes efficiently using a sieve. (I'm sure you've implemented several sieves if you've got as far as problem 231, but see this answer for some suggested implementations.)

So how do we compute the number of times each prime appears in the factorization of a factorial? Let's have a look at a concrete example: how many times does 2 appear in the factorization of 10 factorial? Well, every even number contributes a 2: $$ \mathbf{10}·9·\mathbf{8}·7·\mathbf{6}·5·\mathbf{4}·3·\mathbf{2}·1 $$ Every multiple of 4 contributes an extra 2: $$ 10·9·\mathbf{8}·7·6·5·\mathbf{4}·3·2·1 $$ And every multiple of 8 contributes a third 2: $$ 10·9·\mathbf{8}·7·6·5·4·3·2·1 $$ So in the factorization of \$10!\$, the prime 2 appears $$ \eqalign{ \nu_2(10!) &= \bigg\lfloor {10 \over 2} \bigg\rfloor + \bigg\lfloor {10 \over 4} \bigg\rfloor + \bigg\lfloor {10 \over 8} \bigg\rfloor \\ &= 5 + 2 + 1 \\ &= 8 } $$ times. In general, the power of \$p\$ in the factorization of \$n!\$ is $$ \nu_p(n!) = \bigg\lfloor {n \over p} \bigg\rfloor + \bigg\lfloor {n \over p^2} \bigg\rfloor + \bigg\lfloor {n \over p^3} \bigg\rfloor + \dots $$ This gives us an effective way to compute the power of a prime in a factorial:

def power_in_factorial(p, n):
    """Return the exponent of the prime p in the factorization of n!"""
    result = 0
    while True:
        n //= p
        if not n:
            break
        result += n
    return result

Now, the sum of the factors of \$n!\$ is $$ \sum_{p \text{ prime}} p · \nu_p(n!) $$ where \$p\$ ranges over the primes \$≤n\$, and so the sum of the factors of \$n \choose k\$ is $$ \sum_{p \text{ prime}} p \left( \nu_p(n!) - \nu_p(k!) - \nu_p((n - k)!) \right) $$ That is:

def euler231(n, k):
    """Return the sum of the terms in the prime factorisation of n choose k."""
    f = power_in_factorial
    return sum(p * (f(p, n) - f(p, k) - f(p, n - k)) for p in primes(n + 1))

Using sieve8 from here to implement primes, this takes about 5 seconds.

But so long as we're using NumPy, why not vectorize the whole computation?

def sum_factorial_factors(primes, n):
    """Return the sum of terms in the prime factorization of n!"""
    p = primes
    n = np.full_like(p, n)
    result = r = np.zeros_like(p)
    while True:
        n //= p
        if n[-1] == 0:       # some primes not contributing any more?
            l = n.argmin()   # number of primes still contributing
            if l == 0:
                break
            n, r, p = n[:l], r[:l], p[:l]
        r += n
    result *= primes
    return result.sum()

def euler231(n, k):
    """Return the sum of the terms in the prime factorisation of n choose k."""
    f = sum_factorial_factors
    p = primes(n + 1)
    return f(p, n) - f(p, k) - f(p, n - k)

This takes about half a second.

\$\endgroup\$
2
\$\begingroup\$

Project Euler says I solved this one about seven years ago, but that code is nowhere to be found. It wasn't all that hard to figure it out again...

It should be obvious that performing several million factorizations is going to be very slow. Very often in Project Euler the solutions have to do with reformulating the problem in a way that avoids the obvious, slow operation. In your case, since you have contiguous ranges of numbers, it is relatively easy to figure out how many of those are multiples of a certain number:

def multiples_in_range(start, stop, factor):
    """Returns number of multiples of factor in [start, stop).
    """
    return (stop - 1) // factor - (start - 1) // factor

In your case, you want to take into account extra contributions from multiples of the factor squared, cubed, etc... This is relatively straightforward to do also:

def factors_in_product_of_range(start, stop, factor):
    """Returns exponent of factor in the factorization of the product of
    all numbers in [start, stop).
    """
    exponent = 0
    test_factor = factor
    while test_factor < stop:
        exponent += multiples_in_range(start, stop, test_factor)
        test_factor *= factor
    return exponent

All that's left now is to take this function for a ride over all relevant primes:

def project_euler_231(m, n):
    """Solves Project Euler 231.
    """
    solution = 0
    for prime in primes(m+1):
        solution += (factors_in_product_of_range(m-n+1, m+1, prime) -
                     factors_in_product_of_range(1, n+1, prime)) * prime
    return solution

where primes(n) is a function returning all primes smaller than n.

With a very straightforward implementation of the sieve of Erathostenes, my laptop takes about 6 seconds to compute all primes less than 20,000,000, and about 2 more seconds to iterate over them to give the right answer to the problem.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.