# Finding pairs of Amicable Numbers, up to 10000

I have this program to find pairs of amicable numbers in C:

#include<stdio.h>

int main()
{
printf("|---------- PROGRAM FOR AMICABLE NUMBERS.----------|");
int num1,num2,sum=0;
for(num1=1; num1<=10000; num1++)
{
for(num2=1; num2<=10000; num2++)
{
if ((num1==sum_of_divisors(num2)) && (num2==sum_of_divisors(num1)) && num1!=num2)
{
printf("\n%d\t\t%d", num1,num2);
}
}
}
return 0;
}

int sum_of_divisors(int n)
{
int sum=0,i;
for(i=1; i<n; i++)
{
if(n%i==0)
{
sum=sum+i;
}
}
return sum;
}


I'd like to improve the performance - it's taking around 10 minutes to find even the first 4 pairs.

• You'd get a slight improvement in readability if you used C99 features such as declaring the loop index variables inside the init-statements of the for statements. – Snowbody Jan 9 '18 at 14:04
• Since nobody mentioned it yet: format your source code consistently (indentation, spacing around operators, etc.). Also, don’t waste vertical space: put opening braces on the preceding line. This particular advice might be somewhat controversial but I maintain it objectively improves readability. – Konrad Rudolph Jan 9 '18 at 17:16
• This question is about Project Euler, Problem 21. – Olivier Grégoire Jan 11 '18 at 10:20

    printf("|---------- PROGRAM FOR AMICABLE NUMBERS.----------|");


This kind of banner message makes it harder to use the output of your program in a pipeline. I'd suggest removing it (the user has chosen to run it; trust them to know what they're doing!).

    int num1,num2,sum=0;


What's sum for? It doesn't seem to be used.

    for(num1=1; num1<=10000; num1++)


Where does the magic value 10000 come from? It should be a named constant (or, better, specifiable as a command-line argument).

        for(num2=1; num2<=10000; num2++)


Do we really need to consider every pair in both directions? We could simply iterate while num2 < num1 (which also saves ourselves a comparison later). But see below why this loop isn't needed at all.

                if ((num1==sum_of_divisors(num2))


There's no prototype visible for sum_of_divisors() - enable the relevant compiler warnings and provide the prototype (e.g. by moving the definition ahead of main()).

        for(num2=1; num2<=10000; num2++)
{
if ((num1==sum_of_divisors(num2))
&& (num2==sum_of_divisors(num1))
&& num1!=num2)


In the inner loop, sum_of_divisors(num1) is a constant, so we could save it in the outer loop, and make that the first comparison (so that short-circuit && will then evaluate sum_of_divisors(num2) just once). Then observe that we have a loop which tests whether its index is equal to a particular value - that means that we replace the inner loop with a single test:

    for (num1 = 1;  num1 <= 10000;  num1++) {
num2 = sum_of_divisors(num1);
if (num1 < num2 && sum_of_divisors(num2) == num1) {
/* we found a pair */
}
}


The num2 < num1 test there saves us re-finding pairs we've already seen. num2 > num1 would also work, of course - I chose this version so that we print each pair lowest-first, as that seems to be the convention.

                    printf("\n%d\t\t%d", num1,num2);


It's better to consistently end your output with newline, rather than beginning with newline. This works better in programs that might emit errors or other diagnostics, and plays nicely with line-buffered output such as interactive terminals (where the newline causes a flush).

    for(i=1; i<n; i++) if(n%i==0) sum += i;


This is a really inefficient way to accumulate factors. You can reduce it to the range (1, √n) with sum += i + n/i (adjusted a little to avoid double-counting √n when n is square); it may be better instead to generate all the prime factors and use those to generate the composite factors.

After making the above improvements, I found the runtime was below my measurement threshold. Increasing the limit to 1 million gives a runtime of about 2 seconds - that might be acceptable; if not, the problem parallelizes well (other than outputting the results).

# Modified version

#include <stdio.h>

#if USE_LONG_TYPE
typedef unsigned long Number;
#define FMT "%lu"
#else
typedef unsigned int Number;
#define FMT "%u"
#endif

static const Number MAX_VALUE = 1000000;

Number sum_of_factors(Number n)
{
Number sum = 1;
Number i = 2;
for (;  i < n/i;  ++i) {
if (n/i*i == n) {
sum += i + n/i;
}
}
if (i*i == n)
sum += i; /* add square root only once */
return sum;
}

int main()
{
//#pragma omp parallel for
for (Number num1 = 1;  num1 <= MAX_VALUE;  ++num1) {
Number num2 = sum_of_factors(num1);
if (num2 > num1 && num1 == sum_of_factors(num2)) {
printf(FMT "\t\t" FMT "\n", num1, num2);
}
}
}

• Is it possible that a pair (n, n) is skipped because of the test num2 > num1 ? (I would have written num2 >= num1 or handled equality as a special case) – SylvainD Jan 9 '18 at 22:15
• i*i < n is a problem on overflow, E.g. a long loop with 32-bit unsigned and n == 4294967291, a large prime. Safe alternative i < n/i – chux - Reinstate Monica Jan 10 '18 at 3:48
• The 2nd i*i == n is tricker to avoid OF. I think if (n%i == 0) should work. – chux - Reinstate Monica Jan 10 '18 at 4:18
• @Josay: "Is it possible that a pair (n, n) is skipped...?" Such an (n, n) is not a "pair", it's a single number! Such a number n (whose sum-of-divisors is equal to the number n itself) is called "perfect." Perfect numbers are (by definition) not amicable numbers, so it's correct to skip them. – Quuxplusone Jan 10 '18 at 8:03
• @chux, I've updated my version to avoid the overflow in the loop constraint. As a bonus, it increases the speed, because we can re-use the n/i in the divisibility test (surprisingly, I got an almost 2× speedup, despite n/i and n%i being results of the same division - I thought that modern x86 processors produced quotient and remainder together? Maybe GCC doesn't spot the connection?). – Toby Speight Jan 10 '18 at 10:06
    for(num1=1; num1<=10000; num1++)
{
for(num2=1; num2<=10000; num2++)
{
if (... && (num2==sum_of_divisors(num1)) && ...)


There's an obvious optimisation here which will give you roughly a 10000 times speedup...

int sum_of_divisors(int n)
{
int sum=0,i;
for(i=1; i<n; i++)
{
if(n%i==0)
{
sum=sum+i;


This is not the fastest way to calculate it.

Optimisation 1: only loop up to i*i <= n and use the fact that if n / i = j then n / j = i.

Optimisation 2: what is the formula for the sum of divisors given the prime factorisation? How can you find the prime factorisation efficiently?

Optimisation 3: you don't actually want the sum of divisors of one number, but of all of them up to N. How can you adapt the sieve of Eratosthenes to calculate the sum of divisors of all numbers up to N efficiently?

• Do you mean i <= n/2? Because i*i <= n is giving wrong results. Suppose n=10 then loop will stop at i=4 but 5 is a factor of 10. – Rahul Jan 9 '18 at 16:02
• That sentence contains two halves linked by and. If you implement the first half and not the second half, things will go wrong. – Peter Taylor Jan 9 '18 at 16:07
• @CrisLuengo, if you think that was hinting at merely pulling sum_of_divisors(num1) out of the loop then you have seriously missed the point, which is that the inner loop does not need to be a loop at all. – Peter Taylor Jan 9 '18 at 21:24
• @PeterTaylor, yes I missed that. It's obvious once you see it. :) – Cris Luengo Jan 9 '18 at 21:48
• I saw your "Optimization 3" hint the other day and it caused me to scratch my head. I finally had the AH-HA moment. Well done! – brian_o Jan 11 '18 at 16:18

Try the following (MUCH faster):

int num, num2, num3;
for (num1 = 1; num1<=10000; num1++)
{
num2 = sum_of_divisors(num1);
if (num1 < num2)    /* at least one of any amicable pair is larger */
{
num3 = sum_of_divisors(num2);
if (num3 == num1)    /* FOUND AN AMICABLE PAIR! */
... (record the pair) ...
}
}

• That's super fast. Even calculated till 8th pair in under 5 seconds. – Rahul Jan 9 '18 at 18:27
• Upvoting because you were the first one to write the code for this solution. It was only hinted by Peter Taylor, and Toby Speight only wrote it in a later edit. – Cœur Jan 10 '18 at 2:12
• Welcome to Code Review! You have presented an alternative solution, but haven't reviewed the code. Please edit it to explain your reasoning (how your solution works and how it improves upon the original) so that everyone can learn from your thought process. – Toby Speight Jan 10 '18 at 10:08
for(num1=1; num1<=10000; num1++)
{
for(num2=1; num2<=10000; num2++)
{
if ((num1==sum_of_divisors(num2)) && (num2==sum_of_divisors(num1)) && num1!=num2)


To add a minor note that the other 2 answers didn't mention:

You're calling sum_of_divisors() about 20000x as many times as you should be.

You're calling it many times with the same input value and it's a fairly expensive function.

If you calculate sum_of_divisors() for the numbers 1 to 10000 once, right at the beginning of your program and store the output in an array and look at that whenever you need the value then you won't end up repeating the same work ~10,000x as many times as you need to and you can instead just do it once.

• Actually everything there is in my answer, just in the form of hints. – Peter Taylor Jan 9 '18 at 14:51
• Yes I got that when Peter mentioned 10000 times speedup. Thanks for explaining though. Will be helpful for future viewers. – Rahul Jan 9 '18 at 15:44
• @PeterTaylor your hint made me think of some alternative approaches, I wasn't sure whether you were referring to precomputation or to a few of the other possible optimizations like limiting direction of comparison and/or moving one of the calls into a higher loop or cutting out pairs that can never be amicable. – Murphy Jan 9 '18 at 16:17
• This observation is also in my answer; admittedly it's easily missed: "you can reduce the second loop to a single test". – Toby Speight Jan 9 '18 at 18:03

One thing to consider is the sum_of_divisors function.

It is a pure function, meaning that its return value depends only on the input parameter i.

In addition, it is called 100,000,000 times, but with only values 1 to 1,000,000.

It is ripe for memoization. This basically means that, when the function is called, the return value is cached with the input parameter being the key. Subsequent calls with the same input parameter just return the cached value.

Thus,

int sum_of_divisors(int n)
{
int sum=0,i;
for(i=1; i<n; i++)
{
if(n%i==0)
{
sum=sum+i;
}
}
return sum;
}


becomes

int sums;
int sum_of_divisors(int n)
{
if(sums[n]) {
return sums[n];
}
int sum=0,i;
for(i=1; i<n; i++)
{
if(n%i==0)
{
sum=sum+i;
}
}
sums[n] = sum;
return sum;
}

• Good idea. Make the array static inside the function though. – Peter - Reinstate Monica Jan 10 '18 at 8:34
• Further clarification: if (sums[n]) { for unknown ns works, because int[]array elements get initialized with 0 and 0 is equivalent to false in this case. – Minix Jan 10 '18 at 13:07
• Memoization is a good idea, yet this sum_of_divisors() is O(n) and should be O(sqrt(n)). That O() change has a much bigger impact. – chux - Reinstate Monica Jan 10 '18 at 18:59

I think this is one of the best solutions to your problem. It works in around 1 second and finds all pairs between 1 and 20 millions.

First, we need to calculate the smallest divisor of each number with this method: Sieve of Eratosthenes.

Then we use that sigma(n) - sum of all divisors of n is Multiplicative function. For each number i we know its smallest divisor - p. Find greatest k that i is divisible by p^k. Now, sigma(i) is sigma(p^k) * sigma(i / p^k)

Now, we can calculate sum of proper divisors of each number. sum_proper[i] = sigma[i] - i.

Finally, let's iterate over i and check whether sigma[sigma[i]] = i and print the answer.

# The code

#include <stdio.h>

const int N = 20000000;
const int CNT_PRIMES = N; // upper bound on number of primes within range from 1 to N

int main()
{
int smallest_divisor[N];
int sigma[N];
int sum_proper[N];
int primes[CNT_PRIMES];
int curcnt = 0;
for (int i = 2; i < N; ++i)
smallest_divisor[i] = i;
for (int i = 2; i < N; ++i)
{
if (smallest_divisor[i] == i)
primes[curcnt++] = i;
for (int j = 0; j < curcnt && primes[j] <= smallest_divisor[i] && primes[j] * i < N; ++j)
smallest_divisor[primes[j] * i] = primes[j];
}
sigma = 1;
for (int i = 2; i < N; ++i)
{
int cur = i;
int p = smallest_divisor[i];
int k = 0;
while (cur % p == 0)
cur /= p, ++k;
int pk = i / cur;
sigma[i] = (pk * p - 1) / (p - 1); // sigma(p^k) = 1 + p + ... + p^k = (p^(k+1) - 1) / (p - 1)
sigma[i] *= sigma[cur];
sum_proper[i] = sigma[i] - i;
}
for (int i = 2; i < N; ++i)
{
int j = sum_proper[i];
if (j < N && j > i && sum_proper[j] == i)
printf("%d\t\t%d\n", i, j);
}
}


Ideone Demo

• I gave it a quick check and it's throwing too many exceptions. Mostly variably modified smallest_divisor at file scope. Same for other variables declared outside main. – Rahul Jan 9 '18 at 21:24
• @Rahul what is your compilation arguments? I compiled it with gcc a.cpp -o a -O2 – Alexander Morozov Jan 9 '18 at 21:28
• What insight does this provide about the code in the question? – greybeard Jan 9 '18 at 21:49
• @Rahul You're compiling C code with a compiler for a different language. Of course it does not work as expected. – pipe Jan 10 '18 at 10:30
• @pipe: Dev C++ works both for C and C++. It has MinGW compiler. – Rahul Jan 10 '18 at 16:55

This is actually Project Euler's problem 21.

If you want pure performance, just get rid of the divisions altogether.

I don't know C well enough to actually review the code, but I wrote an efficient solution for Java a few years back, and wanted to share it. The algorithm is rather simple and requires only some changes in the way of thinking about what is a divisor. You'll see, I'm basically only using the addition (except at one time, for... doubling).

// Note: this is Java, not C
public static int projectEuler21(final int limit) {

int[] sumOfDivisors = new int[limit];
int sum = 0;

// For each number below the limit
for (int number = 1; number < limit; number++) {

// Add it as divisor to each of its multiples
for (int multiple = number * 2; multiple < limit; multiple += number) {
sumOfDivisors[multiple] += number;
}

// sum both number and sumOfDivisors[number] if they're amicable, but only if number is the larger of the pair
if (sumOfDivisors[number] < number && number == sumOfDivisors[sumOfDivisors[number]]) {
sum += number + sumOfDivisors[number];
}
}
return sum;
}


The code above is performing extremely well as it runs in the order of 2 or 3 milliseconds. It runs between 5 and 6 seconds if you change the limit to 20,000,000.