In the first part we discussed how to compute the sum of proper divisors function from the sum of divisors function, . We saw that if we factorize the number then compute the list divisors from prime factor, it was much faster (on average) than testing each possible numbers from to (or , 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 , the largest number we will be testing—we need it to compute (up to) . 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 , we only need primes up to (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 is a member of an amicable pair. We first check if 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 . If , then we have a perfect number (a self-amicable number), but we do not care much for those, and we reject . If , we can test if . If they are equal, we have an amicable pair. We then store 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 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 .

*

* *

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 gets (whether or not is an amicable pair). The following plot shows us that the largest is about four time as big as .

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 , we have to test to decide whether or not and form an amicable pair. But the distribution of ratio hints that we could use ratio to prune the search. If 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?