# 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. Commented Jan 9, 2018 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. Commented Jan 9, 2018 at 17:16
• This question is about Project Euler, Problem 21. Commented Jan 11, 2018 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) Commented Jan 9, 2018 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 Commented Jan 10, 2018 at 3:48
• The 2nd i*i == n is tricker to avoid OF. I think if (n%i == 0) should work. Commented Jan 10, 2018 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. Commented Jan 10, 2018 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?). Commented Jan 10, 2018 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. Commented Jan 9, 2018 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. Commented Jan 9, 2018 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. Commented Jan 9, 2018 at 21:24
• @PeterTaylor, yes I missed that. It's obvious once you see it. :) Commented Jan 9, 2018 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! Commented Jan 11, 2018 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. Commented Jan 9, 2018 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
Commented Jan 10, 2018 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. Commented Jan 10, 2018 at 10:08
• @TobySpeight Bearing in mind it isn't completely agreed upon by all, there is a justification, albeit short, but it does somewhat "answer the question" (though it doesn't give concrete measurements of how "MUCH faster" it is) Commented Mar 8, 2022 at 22:48
• It would also benefit from a brief summary of why it's faster, too. Commented Mar 9, 2022 at 8:29
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. Commented Jan 9, 2018 at 14:51
• Yes I got that when Peter mentioned 10000 times speedup. Thanks for explaining though. Will be helpful for future viewers. Commented Jan 9, 2018 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. Commented Jan 9, 2018 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". Commented Jan 9, 2018 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[1000000];
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. Commented Jan 10, 2018 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. Commented Jan 10, 2018 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. Commented Jan 10, 2018 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] = 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. Commented Jan 9, 2018 at 21:24
• What insight does this provide about the code in the question? Commented Jan 9, 2018 at 21:49
• A valid observation about the code to review is that it uses brute force, bound to exceed the time limit of any self-respecting challenge. While the approach above is excellent, its presentation is not - if interested in peer views, present it as a question on code review. Commented Jan 10, 2018 at 4:49
• @Rahul You're compiling C code with a compiler for a different language. Of course it does not work as expected.
– pipe
Commented Jan 10, 2018 at 10:30
• @pipe: Dev C++ works both for C and C++. It has MinGW compiler. Commented Jan 10, 2018 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.

• The Project Euler specification is slightly different, although it does require finding amicable pairs. The question code simply prints the pairs, but PE is asking for the sum of all the amicable numbers in a range. Always check the specification before starting to code! Commented Aug 25, 2023 at 8:15
• @TobySpeight The request was about speed improvement, it beats any of the other answers in terms of performance, even your algorithm, but whatever... I guess there's a reason why Stack Overflow and friends are on the decline. Commented Aug 26, 2023 at 17:41
• Sure, just pointing out that you misread the question code, according to your first sentence. The other suggestions are great. Commented Aug 26, 2023 at 19:26