Amicable Numbers (Part II)

In the first part we discussed how to compute the sum of proper divisors function s(n) from the sum of divisors function, \sigma(n). We saw that if we factorize the number n then compute the list divisors from prime factor, it was much faster (on average) than testing each possible numbers from 1 to n (or \sqrt{n}, since we get one large divisor for every small divisors).

fleuron2

So now, let’s all put that together to get an efficient program that enumerates amicable pairs. First, we obtain a list of primes (not missing any) where the largest prime is O(\sqrt{m}), the largest number we will be testing—we need it to compute (up to) s(m). If we were to search for all amicable numbers, the list would need to be unbounded, maybe we would need to produce it as we go along; but if we look for amicable pairs up to, say, about 2^{32}, we only need primes up to 2^{16}=65536 (and there’sn’t that many of them). For the sake of efficiency, we also need a way of not testing a number twice.

This is readily done by memoizing the tests, and by pruning. Let’s say we want to find out if k is a member of an amicable pair. We first check if k is in the “future reject” list (we’ll explain how a number ends up in it in a moment). If it’s not, we compute s(k). If s(k)=k, then we have a perfect number (a self-amicable number), but we do not care much for those, and we reject k. If s(k)>k, we can test if s(s(k))=k. If they are equal, we have an amicable pair. We then store s(k) in a “future reject” list. The list maintains numbers already tested, successfully or not.

The ideal storage for a reject list is a std::set. It offers O(\lg n) insertion/deletion/testing. A simple implementation might look something like:

int main()
 {
  std::list<bigint> known_primes;
  std::set<bigint> seen;

  // load primes
  get_primes(known_primes,"32bits-primes.txt");

  // bigint is a class for large numbers
  // (really uint64_t here, but could be
  // an arbitrary-precision integer class)
  //
  for (bigint n=1;
       n<std::numeric_limits<unsigned>::max()
       ;n++)

   if (!seen.count(n)) // not on reject list?
    {
     bigint m=s(n,known_primes);

     // check if not tested yet,
     // or not already seen
     //
     if (m>n)
      {
       bigint nn=s(m,known_primes);
       if (n==nn)
        // found a pair!
        std::cout
         << "(" << n << "," << m << ")"
         << std::endl;

       seen.insert(m); // seen!
      }
     else
      ; // perfect, or m already tested!


     // once in a while
     if ((n % flush_interval)==0)
      {
       // prevents the seen list
       // growing to much: we remove
       // numbers we pass by
       std::set<bigint>::const_iterator u=seen.upper_bound(n);
       size_t before=seen.size();
      }
    }
  }

This program takes about 1h to enumerate all 904 amicable pairs with the smallest of the friends smaller than 2^{32}.

*
* *

One first interesting thing we find out is how large the excess candidates are compared to the actual larger friend in an amicable pair. That is, how large s(k) gets (whether or not (k,s(k)) is an amicable pair). The following plot shows us that the largest s(k) is about four time as big as k.

excess

However, amicable pairs aren’t very different, on average. If we plot pairs, we get the following plot:

amicables

Plotting the distribution of the ratio of friends (the ratio of the largest to the smallest) we find this plot:

ratios

The plot gives us a likely criterion to further prune search. When s(k)>k, we have to test s(s(k))=k to decide whether or not k and s(k) form an amicable pair. But the distribution of ratio hints that we could use ratio to prune the search. If s(k):k is much grater than 2:1 (the largest ratio is 1.36:1, and the average 1.08:1), then we are reasonably sure they do not form an amicable pair. Of course, I did not know that until I got enough data, and I’m still to show that this ratio test holds for all amicable pairs. It’s valid for the first few thousands, but what happens later on?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: