Conference rusure::math

An 'abundant' number is one the sum of whose divisors is greater than itself,
in case you didn't already know that...
An example is 24, whose divisors are 1,2,3,4,6,8,12, which sum to 36.
(To be precise, I guess you have to specify aliquot divisors, which excludes 
the number itself as a divisor.) - Lynn

I have such a working program -- the problem is that it'll be a few months
before it completes.  Determining which numbers are abundant is easy; the
hard part is determining the number of ways a number < 10^6 can be expressed
as the sum of two of them.  The straight-forward approach is "order N squared",
where N about 10^12.

Any suggestions for fast algorithms ot do this?  Thanks.

Although it doesn't matter, I've assumed that X+Y and Y+X are two different
ways of expressing the sum.

The number 997920 can be expressed as the sum of two abundant numbers
in 234159 different ways.  This is the maximum over all integers < 10^6.
This is a little surprising, as there are only 247543 abundant numbers < 10^6
(and only 247023 abundant numbers < 997920).

Determining that 997920 yields the maximum of 234159 required checking 999999
through 937884 (there are only 234158 abundant numbers less than 937884).

This was a good way to burn 21 hours of CPU time.

One trick (to make the program run faster) involved recognizing the scarcity
of odd abundant numbers (only about 2000 out of the 247543).  I'll add a
longer note later.

					- Gilbert

The program now runs in under 4 CPU minutes.

The approach is as follows.

First, let abund[i] be a boolean that is true iff i is an abundant number.
These are computed as a simple function of the factorization of each i.

Set modcnt[i] to the number of abundant numbers <= 10^6 that equal i, mod 840.
Let skip[i] be an 840 element bitvector, initially all false.
Initialize max_number_of_ways to 0, and best_n to 0.

Now, varying n from 10^6 downto 1,

    Adjust modcnt[] to be the number of abundant numbers < n.

    If skip[n mod 840] is zero, n cannot produce a new max_number_of_ways,
	so go try the next lower n.

    Cheaply compute an upper bound on the number of ways n can be expressed
	as the sum of two abundant numbers, by summing (over i = 0 to 839)
	    max (modcnt[i], modcnt[(n - i) mod 840]).
	If this sum is less than max_number_of_ways, no number < n that
	    is equal to n (mod 840) will have a larger sum, so set
	    skip[n mod 840] to true, and go try the next lower n.
	    Additionally, if the skip[] are now all true, we are done.

    Now compute the precise number of ways that n can be expressed as
	the sum of two abundant numbers, by counting the number of times
	    that abund[i] and abund[n-i] are both true, over i = 0 to n.

    If this sum is > max_number_of_ways, we have a new maximum.
	Set max_number_of_ways to the sum, and set best_n to n.

When we are done, best_n and max_number_of_ways give the solution.


The reason this works well is because abundant numbers tend to have many
small factors, so the values in modcnt[] are poorly distributed (e.g.,
the values in modcnt[60m] are abnormally high, and values for modcnt[60m�1]
are abnormally low), and the cheap computation eliminates many candidates
for being the new best.

The reason for using 840 rather than 60 (which also has many small factors)
is that it made the program run a bit faster.  If, instead of 10^6, the
problem was for 10^7, a value even larger than 840 would be warranted, since
it becomes more worthwhile to avoid the full computation.