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 ofn(numbers less thannwhich divide evenly inton).

If d(a) =band d(b) =a, wherea≠b, thenaandbare an amicable pair and each ofaandbare called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

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

int sumAmicible = 0; int factorsi, factorsj; for (int i = 2; i <= UpperLimit; i++) { factorsi = sumOfFactors(i); if (factorsi > i && factorsi <= UpperLimit) { factorsj = sumOfFactors(factorsi); if (factorsj == i) { sumAmicible += i + factorsi; } } } [/csharp] All we need to figure out now, is how to calculate the sum of factors, and the first try is as I mentioned a simple brute force approach, which in c# can be implemented as [csharp] private int sumOfFactors(int number) { int sqrtOfNumber = (int)Math.Sqrt(number); int sum = 1; //If the number is a perfect square //Count the squareroot once in the sum of factors if (number == sqrtOfNumber * sqrtOfNumber) { sum += sqrtOfNumber; sqrtOfNumber--; } for (int i = 2; i <= sqrtOfNumber; i++) { if (number % i == 0) { sum = sum + i + (number / i); } } return sum; } [/csharp] It is not the fastest solution in the world, but it works well enough for our purpose. The result is [sourcecode] The largest sum of all amicable pairs less than 100000: 852810 Solution took 109 ms [/sourcecode] <h2>Prime Factorisation</h2> 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. There are several sources on the internet dealing with this problem. The best explanation I have found so far come from the <a href="http://mathschallenge.net/faq">FAQ on the mathchallenge.net</a> site, which as far as I understand is the predecessor of Project Euler. It is also treated under <a href="http://en.wikipedia.org/wiki/Divisor_function">divisor function on Wikipedia</a>. 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 <p style="text-align: center;"></p> <p style="text-align: left;">Where the notation <em>d|n</em> means the divisors <em>d</em> of the number <em>n</em>.</p> <h3>Sum of Divisors Function for a Prime</h3> 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 p<sup>a</sup>, where a is an integer then we have <p style="text-align: center;"></p> and if we mulitply by p we have <p style="text-align: center;"></p> we can subtract these two formulas from each other and get <p style="text-align: center;"></p> Which can be rearranged to <p style="text-align: center;"></p> <h3 style="text-align: left;">Multiplicative property of the function</h3> <p style="text-align: left;">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 <a href="http://http://en.wikipedia.org/wiki/Coprime">coprimes</a>. 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.</p> <p style="text-align: left;">Let both m,n be positive integers and coprime. Assuming that the function is multiplicative we have that</p> <p style="text-align: center;"></p> <p style="text-align: left;">The <a href="http://www.mathacademy.com/pr/prime/articles/fta/index.asp">Fundamental Theorem of Arithmetic</a> (<a href="http://planetmath.org/encyclopedia/ProofOfFundamentalTheoremOfArithmetic.html">and a page with a proof</a>) states that each divisors of mn is a product of two unique positive coprime integers d<sub>1</sub> and d<sub>2</sub> dividing m and n. Thus we can write</p> <p style="text-align: center;"></p> <p style="text-align: left;">The second to last rearrangement follows from the multiplicativeness of scalars.</p> <h3 style="text-align: left;">The algorithm</h3> <p style="text-align: left;">Hence for coprimes the function is multiplicative, so we can expand the function from a prime to any number as</p> <p style="text-align: left;"></p> <p style="text-align: center;"></p> <p style="text-align: left;">where the product runs over the set of distinct prime factors of n. <a href="http://mathworld.wolfram.com/DivisorFunction.html">Wolfram</a>, which gives another source for divisor functions refers to <a href="http://www.amazon.com/Ramanujans-Notebooks-Bruce-C-Berndt/dp/0387961100/">Berndt 1985</a> 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.</p> <p style="text-align: left;">If a prime has the exponent 1 then the formula becomes</p> <p style="text-align: center;"></p> <p style="text-align: left;">Something which we will see once we implement the function. In the<a title="Project Euler – Problem 3" href="http://www.mathblog.dk/project-euler-problem-3/"> solution to problem 3</a> we made a prime factorisation for the first time. We have used it in the <a title="Project Euler – Problem 5" href="http://www.mathblog.dk/project-euler-problem-5/">answer to problem 5 </a>as well. Lets modify those solutions to generate the sum of factors instead. The C# implementation looks like</p> private int sumOfFactorsPrime(int number, int[] primelist) { int n = number; int sum = 1; int p = primelist[0]; int j; int i = 0; while (p * p <= n && n > 1 && i < primelist.Length) { p = primelist[i]; i++; if (n % p == 0) { j = p * p; n = n / p; while (n % p == 0) { j = j * p; n = n / p; } sum = sum * (j - 1) / (p - 1); } } //A prime factor larger than the square root remains if (n > 1) { sum *= n + 1; } return sum - number; }

When we want to use this function we need a list of primes, which we can easily generate with the Sieve of Eratosthenes which we constructed in Problem 10. So the main routine is changed to

int sumAmicible = 0;

int factorsi, factorsj;

int[] primelist = ESieve((int)Math.Sqrt(UpperLimit) + 1);

for (int i = 2; i <= UpperLimit; i++) { factorsi = sumOfFactorsPrime(i, primelist); if (factorsi > i && factorsi <= UpperLimit) { factorsj = sumOfFactorsPrime(factorsi, primelist); if (factorsj == i) { sumAmicible += i + factorsi; } } } [/csharp]And the result is[sourcecode] The largest sum of all amicable pairs less than 100000: 852810 Solution took 20 ms [/sourcecode]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

int sumAmicible = 0;

int factors;

int[] primelist = ESieve((int)Math.Sqrt(UpperLimit) + 1);

Dictionary

for (int i = 2; i <= UpperLimit; i++) { factors = sumOfFactorsPrime(i, primelist); if (factors > i) {

dic.Add(i, factors);

} else if (factors < i) { if (dic.ContainsKey(factors)) { if (dic[factors] == i) { sumAmicible = sumAmicible + i + factors; } } } } [/csharp]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.

Hi Kristian, nice code.

Just a small correction: the upper limit in the problem is 10k, not 100k.

Hi Guy

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.

/Kristian

Here is 10x faster (for the 100000 test) version:

Hi Cat_Baxter

Sorry for the long reply time. I had time to check your method now. Basically you use a sieve to generate all the divisor sums. A pretty clever approach. Thanks for sharing that.

/Kristian

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 :).

Hi Aaron

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.

/Kristian

Table 1 :

Comparison of different solutions : Simple, Advanced and AdvancedWithCache.

http://img15.hostingpics.net/pics/498821PE21tableone.jpg

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).

___

Sources:

1) Code of Project Euler – Problem 21

http://www.mathblog.dk/files/euler/Problem21.cs

2) Microsoft Visual C# 2010 Express

3) Microsoft Office Excel (2007)

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.

Thanks for your code! It is really nice. One thing that I didn’t understand is why 6 is not an amicable number since it relies in the rules.

The sum of divisors d(6) = 6, so it is not an amicable number since it requires that sum to be larger than the number. It is however a perfect number.

Sources: https://www.dropbox.com/s/g6pusinrdiausos/21.zip

I know this is kind of late, but seeing as this page is the first hit when googling this problem, I thought I’d mention I saw something which might be a bug.

if (factorsi > i && factorsi <= UpperLimit) {

This discards amicable numbers whose complement is greater than 10,000.

It’s a minor thing, but I was concerned when I saw my answer was different. I got 861468 for the 100,000 upper bound.

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.