After a few exercises with the focus on other areas, we are in Problem 21 of Project Euler back to focusing on number theory and factorisation. The problem reads

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n). If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

Calculating the proper divisors of one rather small number should not pose much of a problem, but doing it for all numbers below 10.000 and for some of them possibly multiple times means that we have to do a bit of thinking.

I will provide you with an incremental solution, starting from a rather simple brute force approach, with only a few tricks, over to a more advanced calculation of the sum of proper divisors with caching of the results.

Brute Force approach

I have seen solutions where they make two loops from 1 to 10000 checking all combinations, I have chosen to skip directly to a brute force approach which is a bit more intelligent.

The approach runs through all numbers and calculate the sum of factors. If the sum of proper divisors is larger than the number, then we calculate the sum of proper divisors for the just found sum of proper divisors, to check if they are amicable numbers. Since we iterate over all numbers, we only need to check if the sum of proper divisors is larger than the number it self, since otherwise it will already be checked. Those sentences contained a whole lot of “sum of proper divisors”, so let me show you the C# code instead

It is not the fastest solution in the world, but it works well enough for our purpose. The result is

The largest sum of all amicable pairs less than 100000: 852810
Solution took 109 ms

Prime Factorisation

Earlier we have discovered that finding the divisors of a number can be done very efficiently using prime factorisation. This section will provide the method to deduce the sum of divisors based on a prime factorisation of a number.

In the problem d(n) is defined as the sum of proper divisors. However, in this section we will be finding a method for deducing all divisors, so lets call that function t(n) = d(n) + n. So as you see there is a simple relationship between them,

In general this function can be described as

Where the notation d|n means the divisors d of the number n.

Sum of Divisors Function for a Prime

For any prime p, this function would be t(p) = 1 + p, since primes only have the factors one and the number it self.

If we now consider pa, where a is an integer then we have

and if we mulitply by p we have

we can subtract these two formulas from each other and get

Which can be rearranged to

Multiplicative property of the function

So far we have deduced a function which works for a prime factorisation consisting of one prime, or a multiple of that prime. In order to generalise the function to any number, we will need to show that the function is multiplicative for coprimes. Since if it is multiplicative for coprimes, then it will work for any prime factorisation, since any number and a prime are coprimes. None of the aforementioned sources provide a prove of the multiplicative property, but that shouldn’t stop us from deducing it.

Let both m,n be positive integers and coprime. Assuming that the function is multiplicative we have that

The second to last rearrangement follows from the multiplicativeness of scalars.

The algorithm

Hence for coprimes the function is multiplicative, so we can expand the function from a prime to any number as

where the product runs over the set of distinct prime factors of n. Wolfram, which gives another source for divisor functions refers to Berndt 1985 as a source of this function. In other words we can find the sum of factors for an arbitrary number, if we have the the prime factorisation.

If a prime has the exponent 1 then the formula becomes

Something which we will see once we implement the function. In the solution to problem 3 we made a prime factorisation for the first time. We have used it in the answer to problem 5 as well. Lets modify those solutions to generate the sum of factors instead. The C# implementation looks like

The largest sum of all amicable pairs less than 100000: 852810
Solution took 20 ms

A factor five improvement in execution time. Once again a prime factorisation saved us a lot of time.

Caching Results

The last thing we can do in order to save operations is to store the results. If the divisor sum is larger than the number it self, we can store it for later instead of checking that larger number. That also means we have to check the storage if the sum of factors is less than the number.

I have chosen a dictionary, since that is ideal when you have a key and a data object, and we have exactly that. The key is the sum of divisors, and the object is the number it self. The C# implementation looks like

The improvement of running this code is rather small, I was saving 1 millisecond on it, so I wouldn’t deem it worthwhile, but I thought I would mention it anyway.

Closing Remarks

As usual you can find the source code in C#, to play around with, so you don’t have to reconstruct it from the snippets I show you. I am unsure if I have explained the multiplicative property of the sum of divisors well enough, so don’t be shy to add to the explanation or ask questions.

14 Comments For This Post I'd Love to Hear Yours!

Thanks for the compliment, and thanks for noticing. I think I sometime during a test cranked it up to 100.000 to test the timing and then forgot to change it back. I will update the post sometime this week with the correct numbers and timings for the proble.

Thanks for the post, I actually didn’t quite understand the question, so this gave me a better grasp of what it was asking. Reading the answer and understanding it is also nice when it’s a bit beyond my mathematics level :).

Thanks for the comment. Well, I didn’t except that I would help people understand the problem, but I am glad I did. You also seemed to have learned something from blog post, so my mission is accomplished.

Table 1 shows that as the upper limit is growing from 100 000 to 1 000 000, the speed to get the results is also growing from about 7 to about 12 times faster when using Advanced or AdvancedWithCache solutions. All results obtained with Kristian’s algo (1). ___

Hey,thanks for all your solutions. I am really a big fan of yours. I am actually trying on my own to solve the project euler problems. But you give us some 2,3 different methods and clever optimizations. Because of this site,we are able to visualize the problems from different angles. Thanks again. I hope more people visit this site.

Correction on my previous comment — I was incorrectly counting the d(a)==a case. I got 852810, which is the same as yours.

Nonetheless, I’m curious, does your version count a if d(a) > 10,000, or is there a check I’m not seeing. It wouldn’t matter anyway, your solution is correct regardless as there are no amicable pairs straddling the finish line.