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

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

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

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

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?