WARNING: The text below provides a guidance for a Project Euler problem.

  • If you pride yourself in tackling problems entirely on your own, you are encouraged to exit this post right now.

  • If you are simply looking for a hint for the problem, please visit the »Project Euler’s official forum«. You may find some hints here, but if you don’t want your problem spoiled, scroll cautiously.

  • If you are looking for a readily available copy/paste answer, you will not find it here. All of the code snippets in the post have been stripped of some crucial lines (these are clearly marked). It is my aim that anyone able to replicate the answer based on such snippets will have to understand the solution conceptually. At that point, by virtue of learning something new, I believe he deserves it anyway.

  • Guidance for future problems will not be published before 100 people have solved the problem and there are at least two more recent problems.

Divisors of Binomial Product

The problem number 650 is less demanding conceptually and harder in terms of finding an efficient implementation. Denoting the sum of divisors function, we define a slightly nauseating chain of quantities:

The task is to calculate . Knowing the prime factorization of a number, recall the basic formula for calculating the sum of divisor function:

In order to calculate , the only information we need to know about is its prime factorization. Let factor represent the prime factorization of a number as a dictionary-like object. It is then easy to find the factorization of a product of numbers by simply noting that: , where represents “addition” of dictionaries (merge the keys, sum up the values for each key). Similarly, let represent the -fold “addition” of dictionaries (multiply the values in by ; this corresponds to taking powers) and represent “subtraction” of dictionaries (this corresponds to division). To find the prime factorization of , we may use the following recursive relations:

The following piece of Julia code gives answer in about 20 seconds.

  using Nemo, Primes, Memoize
  R = ResidueRing(ZZ, 10^9+7)
  factor = Primes.factor

  @memoize fact(n) = n == 1 ? factor(1) : merge(+, factor(n), fact(n-1))
  @memoize B(n) = n == 1 ? factor(1) : merge( # REMOVED
  sigma(f) =  prod( (R(p)^(e+1)-R(1)) * R(p-1)^-1 for (p, e) in f )
  S(n) = R(1) + sum(sigma(B(k)) for k=2:n)

  @time S(20000)