Prob 0021 ‐ Amicable Numbers - maccergit/Euler GitHub Wiki

It turns out this problem is highly sensitive to the speed of the algorithm used to find the divisors (not factors) of a number, and there are several ways to find the divisors of a number. Before getting into the solution for this problem, let's look at different approaches to getting the divisors of a number :

00.01

The simplest way is to do trial division of each integer smaller than the number. We can take a small optimization by realizing that 1 is always a divisor, and that 2 is the next possible divisor - and if the number divides by 2, then the other divisor left after dividing by 2 is the largest divisor for the number - we don't need to look at numbers larger than $\frac{N}{2}$. So, the range of trial divisors is $2\dots\frac{N}{2}$. This is easy to code as a list comprehension, and is fairly fast - the linear growth is reasonable, but very large numbers can take a bit of time, and this time can add up when they are done repeatedly :

00.02

The trial division approach can be further optimized by realizing that once we find a divisor, we can divide the number by that divisor to yield another divisor - and we can do that as we keep moving through trial candidates one by one. This means we are finding the divisors from smallest up, along with those from largest down, at the same time. They also won't meet at $\frac{N}{2}$, but will meet at $\sqrt{N}$, so the range is now $2\dots\sqrt{N}$. This is three orders of magnitude faster - and the scalability curve is a better curve (it's concave down, instead of linear). This means the algorithm will get even more efficient compared to a linear algorithm as the limit increases :

00.03

Another way of finding the divisors is to look at all the combinations of products of the factors of the number. We already have a fast factorization routine (pyprimesieve) - all we need to do is find all the combinations of pairs, triples, etc... of the factors and compute their products. The set of all possible combinations of elements taken 1 at a time, 2 at a time, 3 at a time, and so on is called the "power set". There is also a Python library called "itertools" that are fast implementations of common iteration needs. In this case, itertools documentation includes a recipe for creating a power set by combining the itertools "combinations" and "chain.fom_iterable" functions. Let's see if this is faster :

Even though this is a couple orders of magnitude faster than the initial trial division method, it is still much slower than the optimized trial division approach. This was a surprise, as usually the approaches that use factoring end up being more efficient than the direct methods.

01

"What's a few milliseconds among friends?"

And so we now look at the solutions to be problem. The direct approach is to iterate over the range of candidates, getting the sum of the divisors (note they need to be "proper divisors", so routines that return the number itself need to remove it from the sum), and then getting the sum of the divisors again to see if the result matches the original number.

For our first try, let's use the non-optimal way of finding divisors. It works, and even returns the answer in a couple seconds :

However, the quadratic growth is concerning - especially when we start with a timing measured in seconds.

02

Using the optimized trial division approach to finding divisors, we get much better results :

03

A little research reveals that the sum of the divisors can be computed without getting all of the divisors, but instead can be computed directly from the prime factors and their exponents :

\displaystyle \sum{\text{divisors}} = \prod{\frac{p^{power + 1} - 1}{p - 1}}

Doing this with our basic factorization library, we still get a performance improvement over the best trial division approach :

04

If we combine the previous solution with our fastest factorization library, we get another improvement in performance :