Finding pairs of Amicable Numbers, up to 10000

Question

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.

Answer 1

    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);
        }
    }
}

Answer 2

    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?

Answer 3

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

Answer 4

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.

Answer 5

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;
}

Answer 6

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

Answer 7

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.