Home

Discoverers

Statistics

News

Articles

Submit new pairs

List of articles

A collection of methods for finding type (i, 1) amicable pairs by Sergei Chernykh. May 10th, 2016
A million new amicable pairs every day: some new methods for finding amicable pair breeders by Sergei Chernykh. December 20th, 2015
Amicable pairs: searching for them efficiently by Sergei Chernykh. November 11th, 2015

---

A collection of methods for finding type (i, 1) amicable pairs/breeders

Sergei Chernykh
May 10th, 2016

Contents

1. Exhaustive methods
   1.1. Exhaustive search
   1.2. Exhaustive search with unknown largest prime factor

2. Constructive methods
   2.1. The general formula to construct type (i, 1) amicable pairs/breeders
   2.2. Type (2, 1) amicable pairs/breeders
   2.3. Type (3, 1) amicable pairs/breeders
   2.4. Type (4, 1) amicable pairs/breeders
   2.5. Type (5+, 1) amicable pairs/breeders
   2.6. Other constructive methods

3. Conclusion

1. Exhaustive methods

1.1. Exhaustive search

- Take two sets of numbers (a and u)
- Sort one of sets by ratio: (S(a)-a)/a or (u-1)/S(u)
- Iterate through the second set and find all pairs which satisfy (S(a)-a)/a = (u-1)/S(u)

For details refer to the previous article.

1.2. Exhaustive search with unknown largest prime factor

It allows to find many more breeders in the same amount of time as the regular exhaustive search.

- Take two sets of numbers (a and u)
- Sort first set by ratio: (S(a)-a)/a
- Iterate through the second set and find all pairs which satisfy (S(a)-a)/a = (u×p-1)/(S(u)×(p+1)) for some prime p (p is the largest prime factor in u×p)

For each (u, a) pair: (S(a)-a)/a=(u×p-1)/(S(u)×(p+1))

After simplyfying:

p=(a+D)/(a×u-D) where D=(S(a)-a)×S(u) (1)

If p from (1) is integer and prime, then we've found a breeder (u×p, a). We know that p > p₀ (largest factor in u) which puts restrictions on the value of x=(u×p-1)/(S(u)×(p+1))=(S(a)-a)/a

(u×p₀-1)/(S(u)×(p₀+1)) < x < u/S(u)
or
(u×p₀-1)/(S(u)×(p₀+1)) < (S(a)-a)/a < u/S(u)

which limits the set of values of 'a' to check against. The larger p₀ is, the smaller this set will be, so it makes sense to explicitly check first few primes p_i > p₀ until this set gets small enough. A rule of thumb: if 'k' is how many values (p_i > p₀, i = 1, 2, ..., k) you check explicitly, then you end up with

(u×p_k-1)/(S(u)×(p_k+1)) < (S(a)-a)/a < u/S(u) (2)

You should then choose k such that the number of values of 'a' satisfying (2) is less than k and stop as soon as it gets less than k. This reduces asymptotic complexity to check each value of 'u' from O(N) to O(sqrt(N)).

This method proved to be effective for generating lots of small breeders (10-30 digits for u). Small breeders have the highest chances to generate other breeders/amicable pairs using te Riele's rule or any other constructive rule.

When used with parameteres u = p₁×p₂×p₃×p₄×p₅ ≤ 10¹⁵, a ≤ 10⁶, it found (among many others) the following breeder:

u=71×73×1061×3019×4211×65521001603
a=2×5×13×59

which in turn generated the first type (8, 1) amicable pair using te Riele's rule:

81 Chernykh 2016
14438796350376953678174322791271636288778421212339464839728409858893107308551670=2×5×13×59×71×73×1061×3019×4211×65521001603×4841941814273029742117557×84877110421196263915758493
14865176913743789894174101450879657763077918214583046351774081022511427255902730=2×5×13×59×1938093469849255527271721180036461246815895464743552327480323470992363397119

2. Constructive methods

2.1. The general formula to construct type (i, 1) amicable pairs/breeders

If we take any pair of numbers (u, a) then (u×p×q, a) is a breeder if:

S(a)/a=(u×pq+S(u)×(pq+p+q+1)-1)/S(u)/(pq+p+q+1)

for some prime numbers p and q. After simplifying we get this BDE:

(E×p-D)(E×q-D)=a×(E+D×u) (3)

where D=(S(a)-a)×S(u), E=u×a-D

If (u, a) is already a breeder then (3) becomes this after further simplifying:

(p-(u-1))(q-(u-1))=u×(u-1)+1 (4)

This is te Riele's rule to generate new breeders from existing breeders. Other rules for generating new breeders from existing breeders were given in the previous article.

2.2. Type (2, 1) amicable pairs/breeders

These have the form (p×q, a), so u = 1 in this case. Substituting (1, a) into (3) and simplifying:

(D×p-E)(D×q-E)=a² (5)

where D = a×2 - S(a), E = S(a)-a

This BDE is very convenient: it doesn't require factorization of the right side so we can just iterate all posible factorizations of 'a' in given range. I solved this BDE for all a ≤ 1.2 × 10¹³ and for many a ≤ 10¹⁸ with small prime factors. The most suitable choice of 'a' for BDE (5) is such that D(a)=1 or D(a)/g(a)=1 where g(a)=GCD(a, S(a)). It turns (5) into:

(p - E/g)(q - E/g) = a²/g² (6)

Let's find as many such values of 'a' as possible.

It's easy to see that D=1 for a=2ⁿ. Only a few amicable pairs of the form (2ⁿ×p×q, 2ⁿ×r) are known however because 2ⁿ grows very quickly and demolishes all chances for p, q, r to be prime at the same time. But it's possible to find many more relatively small values of 'a' such that D(a)/g(a)=1.

For example, for a=2ⁿ×p:

D=2ⁿ⁺¹×p - (2ⁿ⁺¹-1)×(p+1) = 2ⁿ⁺¹×p - 2ⁿ⁺¹×p - 2ⁿ⁺¹ + p + 1 = p+1 - 2ⁿ⁺¹

if D=2 then p = 2ⁿ⁺¹ + 1, S(p)=2ⁿ⁺¹ + 2, g(a) = GCD(a, S(a)) = 2, D/g=1
if D=4 then p = 2ⁿ⁺¹ + 3, S(p)=2ⁿ⁺¹ + 4, g(a) = GCD(a, S(a)) = 4, D/g=1
...
if D=2ⁿ⁺¹ then p=2ⁿ⁺²-1, S(p)=2ⁿ⁺², g(a) = GCD(2ⁿ⁺¹×p, (2ⁿ⁺²-1)×2ⁿ⁺²) = 2ⁿ⁺¹, D/g=1

so we can take

a = 2ⁿ × (2ⁿ⁺¹ + 2^k - 1) (7)

where k = 1, 2, ..., n+1 and the number in brackets is prime.

Then comes the case where we can take any 'a' with D(a)/g(a) = A ≥ 1 and find p and q such that D(a×p×q)/g(a×p×q) = k:

The first case is when p doesn't divide g(a×p):

a×2×pq - S(a)×(pq+p+q+1) = k×g×d, d divides a/g

simplifying it, we get:

(A×p-b)(A×q-b) = b×(b+A) + A×k×d (8)
where A = (a×2-S(a))/g
b = S(a)/g
d divides a/g
k=A generates new values with the same D(a)/g(a)
k=1 generates new values with D(a)/g(a)=1
"A×k×d" generates new values with D(a)/g(a)=k, 1 ≤ k ≤ A

for a=2ⁿ we have n+1 equations (0 ≤ k ≤ n):

(p-(2ⁿ⁺¹-1))(q-(2ⁿ⁺¹-1))=(2ⁿ⁺¹-1)×2ⁿ⁺¹+2^k (8.1)

The second case is when p divides g(a×p). Let's consider this case only for D(a)/g(a) = A = 1:

(a×2-S(a))×pq-S(a)×(p+q+1)=g×d×p

simplifying it, we get:

(p-b)(q-b-d))=b×(b+d+1) (9)
where b = S(a)/g
d divides a/g

for a=2ⁿ we have n+1 equations (0 ≤ k ≤ n):

(p-(2ⁿ⁺¹-1))(q-((2ⁿ⁺¹-1)+2^k))=(2ⁿ⁺¹-1)×(2ⁿ⁺¹+2^k) (9.1)

Applying (8) and (9) to a=2ⁿ, a=2ⁿ×p found before produced a lot of new values a=2ⁿ×p₁×p₂×...×p_k. (8) and (9) were applied to these new values again and again as long as it was manageable to factorize right sides of these equations. I also ran a big exhaustive search for 'a' such that D(a)/g(a) ≤ 10⁴, a ≤ 10¹⁸ and substituted found values into (8) and (9) with k = 1. As a result, I found almost 170,000 values of 'a' with D(a)/g(a) = 1 and solved (6) for all of them. It produced hundreds of new type (2, 1) amicable pairs. The record for the largest type (2, 1) amicable pair was also increased from 329 digits (Garcia, 1998) to 545 digits: 29 new pairs were larger than the previous record.

2.3. Type (3, 1) amicable pairs/breeders

These have the form (p×q×r, a), so u = p in this case. Substituting it into (3) we get

D = (S(a)-a)×(p+1)
E = p×a-D = p×a-(S(a)-a)×(p+1) = p×a-S(a)×(p+1)+a×(p+1) = p×a-S(a)×p-S(a)+a×p+a = p×(a×2-S(a))-(S(a)-a)

E must be > 0 which sets a lower bound for p:

p > (S(a)-a)/(a×2-S(a)) (10)

A rule of thumb: for each 'a' try just the first few primes larger than the right side of (10) to get the highest chances to find a solution of (3)

There is also a higher bound for p:

p < 3×(S(a)-a)/(a×2-S(a)) (11)

Exhaustive search is still manageable, and I was able to find most (3, 1) amicable pairs with a < 10¹². Values of 'a' with D(a)/g(a) = 1 are also suitable here and generated hundreds of new amicable pairs as well.

2.4. Type (4, 1) amicable pairs/breeders

u = p×q×r×s, so an exhaustive search is not feasible here. But trying first few possible values for p and q and solving (3) for r and s is very fruitful. (10) still holds for p, and we have this lower bound for q:

q > (S(a)-a)×(p+1)/(a×p - (S(a)-a)×(p+1)) (12)

I tried p = smallest prime satisfying (10) and q = first few primes (up to a hundred) satisfying (12) for given p, and vice versa: p = first few primes (up to a hundred) satisfying (10) and q = smallest prime satisfying (12) for given p. In each case, I considered equation (3) for r and s, and if E was relatively small after dividing out common factors (≤ 10⁴ in general or ≤ 10⁶ in the case of smallest possible p and q), it was solved.

This approach resulted in finding over 1,000 new type (4, 1) amicable pairs and over 25,000 breeders of the form (p×q×r×s, a). More than 10⁹ equations (3) were solved in the process.

2.5. Type (5+, 1) amicable pairs/breeders

It's possible to use the same approach as in 2.4. for type (5, 1), but it doesn't produce as many new pairs. The smallest (5+, 1) amicable pairs were found by some form of an exhaustive search (1.1 or 1.2), the other pairs were found by te Riele's rule or other constructive methods.

2.6. Other constructive methods

A rule to generate breeders of the form (u×q, a×p) from (u, a) was described in the previous article. A Thabit rule for amicable pairs/breeders was also described there.

There is also a small modification of te Riele's rule that can generate some more pairs/breeders: if (u, a) is a breeder, then (3) becomes (4), i.e. E becomes equal to 1. But if (u, a) is "almost" a breeder, then E becomes very small and (3) is still likely to have a solution. By "almost" a breeder I mean pairs (u, a) which satisfy

S(a)/a = (u+S(u)-1)/S(u) + ε (13)

for some very small ε. Breeder (u, a) produces "almost" breeders of these types:

(u, a×p) where p > S(u)×S(a)/a and is prime (14)
(u×p, a) where p > u and is prime (15)

For example, take one of Euler's amicable pairs:

21 Euler 1750
11498355=3⁴×5×11×29×89
12024045=3⁴×5×11×2699

u=2581=29×89
a=4455=3⁴×5×11

(14) gives lower bound for p: p > 2700×8712/4455 = 5280. If we take p = 5281 for (14) then solve (3), we'll generate a breeder (29×89×13674611×4202577551, 3⁴×5×11×5281) which is also a type (4, 1) amicable pair found by te Riele (of course) in 1982:

41 TeRiele 1982
3489655839073826955942266055=3⁴×5×11×5281×29×89×13674611×4202577551
3650550738768150894849977145=3⁴×5×11×5281×155165267043476524799

Another example:

u=397×479×141403=26889618689
a=3²×5²×13×29=84825

(15) gives lower bound for p: p > u = 26889618689. If we take p = 26889618697 for (15) then solve (3), we'll generate a breeder (397×479×141403×26889618697×80339065939437225367×735307463850521298183292182550219351727, 3²×5²×13×29) which is also a type (6, 1) amicable pair:

61 Chernykh 2016
3623171640394727379955351072502990492618130497038578416351325756458993451782755208025=3²×5²×13×29×397×479×141403×26889618697×80339065939437225367×735307463850521298183292182550219351727
3639906844276458222356992416948731788196043339981019887096597468959958779504615339175=3²×5²×13×29×42910779183925236927285498578823834815161135749849924987876185899911096722718719

3. Conclusion

The methods described in this article were successfully used to generate thousands new amicable pairs of type (i, 1) and over 70,000 new breeders.
Borho's rule was applied to most of these new amicable pairs/breeders which resulted in bringing the total number of known amicable pairs to 1,000,000,000.

The most productive breeder was this one:

u=163×1042759×773144273090265929×3071110954662140597
a=3×5³×7²×101×109×193×269×563

It generated 17,010,927 new amicable pairs of type (5, 2).

---

A million new amicable pairs every day: some new methods for finding amicable pair breeders

Sergei Chernykh
December 20th, 2015

Contents

1. An almost exhaustive search for breeders
2. Generating new breeders from known breeders
3. Thabit rule for generating new huge breeders from known breeders

Preface
For a definition of what breeder is and how it can be used to generate amicable pairs, refer to
"Amicable pairs, a survey" by M. Garcia, J.M. Pedersen, H.J.J. te Riele - section 5, "Searches of amicable pairs of a special form", definition 5.4.

1. An almost exhaustive search for breeders

If (au, a) is a breeder for Borho's rule then by definition there is a natural number x such that:

au + ax = S(au) (1)
and
S(au) = S(a)(x + 1) (2)

Let's assume that gcd(a, u) = 1. It cuts some breeders, but simplifies a lot of things:

(2) becomes S(a)×S(u) = S(a)×(x+1), so x = S(u) - 1
(1) then becomes a×u + a×(S(u) - 1) = S(a)×S(u)

divide both sides by a: S(a)/a×S(u) = u+S(u)-1
divide both sides by S(u) and we get

S(a)/a = (u+S(u)-1)/S(u) (3)

Searching for breeders is now just finding all pairs of numbers (a, u) for which (3) holds.

• Observation 1: "a" is the denominator on the left side of (3) whereas S(u) is the denominator on the right side of (3).
  Divisor sums tend to have only small prime factors in their factorizations, so it's reasonable to use "increasing largest prime factor order" (see the previous article) for values of "a".
  
  
• Observation 2: (3) can be rewritten as (S(a)-a)/a = (u - 1)/S(u) which gives (S(a)-a)/a < 1.
  The other known limitation is S(a)/a > 3/2, or (S(a)-a)/a > 1/2.
  
  
• Observation 3: since both a and u will be used to generate amicable pairs, skip "a" which are divisible by 6.
  Also consider only "u" coprime to 6, since extremely few breeders have "u" divisible by 2 or 3.

I used the following algorithm to find (almost) all such pairs in the specified range.

1. Generate as many values of (S(a)-a)/a as can fit in the available physical memory.
   Use "increasing largest prime factor order" and consider only values of "a" for which 0.5 < (S(a)-a)/a < 1 holds.
   We know that they're all between 0.5 and 1, so their binary representation will be "0.1..."
   We can store 64 bits which come after "0.1" for each (S(a)-a)/a value, and round the last bit to nearest,
   so maximal absolute error will not exceed 2^-66. Such high precision gives very few false positives when searching,
   and they all can be filtered by checking the exact values in (3), rewritten as S(a)×S(u) = (u+S(u)-1)×a to avoid rounding errors because of division.
   
   
2. Iterate over all values of u, starting from 1. Use modified Sieve of Eratosthenes to get all values of S(u) in range and avoid doing factorizations for each u in range.
   Typical range size would be around 10⁶...10⁷ for optimal performance. For each (u, S(u)) pair calculate (u-1)/S(u) and if it's between 0.5 and 1,
   search for this value in the set of (S(a)-a)/a values generated in step 1. I stored (S(a)-a)/a values in ascending order which allowed binary search,
   plus I used some tweaks like cache-friendly binary search. Hash map would be faster, but it wastes memory: we need to store as many values as possible in all available memory, every byte counts here.
   For each "u" value which gave a match, save (u, S(u)) pair to a file on disc.
   
   
3. Generate the same (a, S(a)) values as in step 1 again and match them against the set of (u, S(u)) pairs found in step 2.
   This time, use the exact formula S(a)×S(u) = (u+S(u)-1)×a for matching.

I ran the algorithm described above with the following parameters:

• Step 1: a ≤ 10¹³, 62 GB of memory: 8,321,499,136 (S(a)-a)/a values used for searching. It took ~30 minutes including sorting the array.
• Step 2: u ≤ 10¹³. It took a couple of days to finish this step.
• Step 3 took ~15 minutes.

As a result of this search and a number of smaller searches, I found 3884 breeders, most of them were new.

Source code for the described algorithm: BreederSearch_2015.zip

2. Generating new breeders from known breeders

The rule for generating new breeders described below is new. It wasn't published anywhere before.

Given a pair (a, u) which corresponds to a breeder (au, a), find prime numbers p, q such that (a×p, u×q) give a new breeder (ap×uq, ap).

Substitute (a×p, u×q) into (3):

S(a)×(p+1)/(a×p) = (u×q+S(u)×(q+1)-1)/S(u)/(q+1)

After some boring calculations and getting rid of "a" by substituting (S(a)-a)/a = (u - 1)/S(u) and S(a)/a×S(u)=(u+S(u)-1) we get this Bilinear Diophantine Equation (BDE):

(p-(u+S(u)-1))×(q-u) = (u+S(u)-1)×(u+1) (4)

Example:
a=8, u=253, S(u)=288
(4) becomes (p-540)×(q-253) = 137160
which has a solution p=541, q=137413
We've successfully generated (8×541, 253×137413) from (8, 253).

Applying this rule to the list of 3884 breeders found before, I got 480 breeders in the second generation and 27 breeders in the third generation.
It directly follows from the left side of (4) that q > u which means u×q > u², so each child breeder is at least twice as long as its parent breeder.
Indeed, first generation breeders were up to 13 digits long, second generation - up to 38 digits, third generation - up to 75 digits.
I also applied this rule to many of the previously known breeders extracted from known amicable pairs.
My list grew to 18035 breeders as a result, 563 of them are 50 digits or bigger.

Second and third generation breeders produced an awful lot of new amicable pairs - more than 36 millions in just 3 weeks!

Source code is not provided since it's just a straight-forward BDE solver for equation (4).

3. Thabit rule for generating new huge breeders from known breeders

Borho and Hoffman in "Breeding Amicable Numbers in Abundance" [1986] give a Thabit rule for generating amicable pairs from breeders of the form (au, a):

If t = u + S(u) is a prime number then for n = 1,2,3,...
A=au×tⁿ×[tⁿ(u+1)-1]
B=a×tⁿ×[tⁿ(u+1)(t-u)-1]
are an amicable pair if both numbers in square brackets are prime.

It's easy to see that (A, B) is an amicable pair of type (i, 1) which means that the same rule can be modified to find new breeders:

If t = u + S(u) is a prime number then for n = 1,2,3,...
A=a×tⁿ
U=u×[tⁿ(u+1)-1]
are a breeder if the number in square brackets is prime.

This rule generates very huge breeders due to quickly growing tⁿ.
Applying it to breeders in my list produced 1378 new 50+ digit breeders, the largest having 4852 decimal digits for "u".
This breeder can potentially produce 15000+ digit amicable pairs. I tried some smaller breeders and found 4 new 1000+ digit amicable pairs.

---

Amicable pairs: searching for them efficiently

Sergei Chernykh
November 11th, 2015

Contents

0. Straight-forward algorithm
1. Avoiding factorization of m
2. Trimming down factorization of n
2.1. Check sum bounds
2.2. Check sum divisibility
2.3. The finish line of trial division
3. Reducing search space
4. Conclusion and source code


Preface
This article summarizes all techniques for an exhaustive search of amicable pairs currently known to me.
The C++ source code implementing these techniques is provided.
This source code implements the leading-edge algorithm for running an exhaustive search of amicable pairs.
It was successfully used to push the exhaustive search limit from 10¹⁴ up to 10¹⁷ (finding 282,197 new amicable pairs)
and then to find over 1 million new amicable pairs above 10¹⁷.

Terminology
N - the limit for the search. The algorithm must find all amicable pairs with smaller member m ≤ N
m - smaller number in a pair
n - larger number in a pair, or an arbitrary natural number when used in enumerations (i.e. p₁, ..., p_n)
σ(m) - sum of all divisors of m, including 1 and m.

0. Straight-forward algorithm

The following basic algorithm can be derived from the very definition of amicable pairs:

• Iterate over all numbers between 1 and N
• For each number m, calculate n=σ(m)-m
• If n > m, then calculate σ(n)
• If σ(n) = σ(m), then (m, n) - amicable pair

This is very inefficient because it requires 1 and in many cases 2 full factorizations per iteration.

1. Avoiding factorization of m

Straight-forward algorithm iterates from m to m + 1 on each step.

One way to eliminate factorization of m is to store factor table of sufficient size.
It was first used by David and Paul Moews: A SEARCH FOR ALIQUOT CYCLES BELOW 10¹⁰ (1991)
For example, if m=p₁^k₁×...×p_n^k_n then we know that p₁ divides m+p₁, p₂ divides m+p₂ and so on.
We only need to store factors for all numbers between m and m + sqrt(N) at each step.
If we store this data from the beginning of search, we'll effectively eliminate factorization of m at each step.

This approach has several downsides:

• The factor table will grow up to sqrt(N) elements in the end of the search. It will not fit in CPU cache.
• Each parallel thread must have its own factor table, making memory requirements and CPU cache contention even worse.
• You'll still need to divide m by p₁, ..., p_n in a loop to get factor powers (k₁, ..., k_n). Integer division is slow.
  Divisions can be eliminated using a more tricky version of this method (see the link above).

A better way to eliminate factorization of m is to change iteration order!
Instead of iterating from m to m + 1, we'll iterate from m to m × p where p is a prime.
Thus we'll know the factorization for each m right away. We don't even need to store the actual factorization, what we need to store is σ(m).
This idea was first used by David Einstein in the beginning of 2000's.

One more benefit from this iteration order is that amicable numbers (more precisely, smaller members of amicable pairs) tend to have only small prime factors in their factorization.
This iteration order finds more than a half of all amicable pairs below N in only a small fraction of the total search time.

A simple recursive algorithm implements this iteration order:

• On the first level of recursion iterate over all primes from 2 to N/20 (there are no amicable numbers of the form k × p where k < 20 and p is prime)
• On the second level of recursion iterate over all primes from 2 to p (but not including p), where p is a prime from the first level of recursion
• and so on
• at each recursion level iterate over powers of current prime: p, p², ..., pⁿ

There are two cases here:

• Iterate from m to m × p where p doesn't divide m.
  Then m₁ = m × p and σ(m₁) = σ(m) × (p + 1).
  Two multiplications, one addition.
  
  
• Iterate from m × p^k to m × p^k+1, k > 0.
  We already have σ(m), p^k, m × p^k, σ(m) × p^k, σ(m × p^k) at our disposal here.
  The last 4 values can be saved when calculating m × p^k and σ(m × p^k).
  
  m₁ = (m × p^k) × p - one multiplication
  
  σ(m₁) = σ(m) × σ(p^k+1)
  = σ(m) × (1 + p + ... + p^k+1)
  = σ(m) × (1 + p + ... + p^k) + σ(m) × p^k+1
  = σ(m × p^k) + σ(m) × p^k+1
  = σ(m × p^k) + (σ(m) × p^k) × p
  
  Both σ(m × p^k) and σ(m) × p^k can be taken from a previous step, so it takes one multiplication and one addition to calculate σ(m₁).
  

Both cases require only two multiplications, one addition, no divisions at all and only a fixed and small amount of memory to iterate from m to m₁.
As a result, we get m and σ(m) already calculated. Perfect!

2. Trimming down factorization of n

One can notice that factorization of n is not what's needed. What's actually needed is to check that σ(n) = σ(m).
Knowing that, we can try to find ways to prove that σ(n) ≠ σ(m) way before completing factorization.

2.1. Check sum bounds

First and obvious thing is to check that min(σ(n)) ≤ σ(m) ≤ max(σ(n)) after each trial division step.
Indeed, we know that n = n₁ × n₂, where n₁ - factored portion of n and n₂ - unfactored portion of n.
Hence σ(n) = σ(n₁) × σ(n₂). We know σ(n₁) and we know that n₂=p₁^k₁×...×p_n^k_n where p₀ ≤ p₁ < p₂ < ... < p_n, p₀ - smallest prime factor we haven't tried yet.

Lower bound for σ(n): σ(n) = σ(n₁) × σ(n₂) ≥ σ(n₁) × (n₂ + 1) (σ(n₂) = n₂ + 1 when n₂ is prime)
Upper bound for σ(n): σ(n) = σ(n₁) × σ(n₂) ≤ σ(n₁) × max(σ(n₂)) where max(σ(n₂)) is a maximal sum of divisors for all numbers satisfying these conditions:

• n₂ = p₁^k₁×...×p_n^k_n where p₀ ≤ p₁ < p₂ < ... < p_n, p₀ - smallest prime factor we haven't tried yet.
• n₂ ≤ t for some number t, t≥n/n₁

We can pre-calculate a look-up table which gives max(σ(n₂)) for each pair (p₀, t), where both p₀ and t are taken from fixed sets of numbers.
Herman te Riele in "Computation of All the Amicable Pairs Below 10¹⁰" (1986) gives more details about this table.
He uses analytical estimation of max(σ(n₂)) though, but it's possible (and it's actually implemented in the code) to calculate the exact values of max(σ(n₂)) for each given (p₀, t).

This table is of course has an upper limit for p₀, because we don't have infinite amounts of memory.
My implementation calculates it for all p₀ < 2¹⁶, t is a primorial number < 2⁶⁴, so it can be used for all numbers < 2⁶⁴ and containing 4 or more prime factors.

What about cases where n₂ has 3 or 2 prime factors? The table can't be used here, but it's possible to get formulas for max(σ(n₂)) for such simple cases.

I derived these formulas:

max(σ(n₂)) ≤ n₂ + p₀ + n₂/p₀ + 1 when n₂ has at most 2 factors, both ≥ p₀
max(σ(n₂)) ≤ n₂ + (n₂/p₀ + p₀² + p₀) × 2 when n₂ has at most 3 factors, all of them ≥ p₀

2.2. Check sum divisibility

The next observation is that if σ(n) = σ(n₁) × σ(n₂) and σ(n) = σ(m) then σ(n₂) = σ(m) / σ(n₁).
n₂ is an integer, σ(n₂) is of course also an integer which means that σ(m) / σ(n₁) must be an integer.
If σ(n₁) doesn't divide σ(m) then we can immediately conclude that σ(n) ≠ σ(m) and stop factorization.
This technique prunes numbers very effectively. It discards most candidates after finding just a few smallest prime factors.

2.3. The finish line of trial division

If in spite of all previous filters we got to a case where σ(n) = σ(n₁) × σ(n₂) where n₂ can have at most two factors, then we can just calculate these factors instead of doing trial divisions!
Indeed, n₂ can be either a prime, or a perfect square, or a product of two distinct prime numbers.

First two cases are simple to check, the last case gives two equations (implying that we want to check if "σ(n) = σ(m)" is true)

p×q = n₂
(p+1)×(q+1)=σ(m)/σ(n₁)

Simplifying it, we'll get this:

p×q = n₂
p+q = σ(m)/σ(n₁)-n₂-1

which is a Vieta's formula for this quadratic equation:

x² - (σ(m)/σ(n₁)-n₂-1)×x + n₂ = 0

Find roots, check that they're distinct positive integers satisfying equations for p and q, and you're done!

3. Reducing search space

The algorithm can be even faster if it doesn't need to iterate over all numbers between 1 and N, right? Some numbers can never be a member an amicable pair.

For example, all numbers divisible by 6 can't be a part of an even-even amicable pair.
They can still be a part of an even-odd amicable pair, but such pairs must have a very special form which can be easily and quickly searched for separately.

Another consideration: m must be abundant, i.e. σ(m) > m × 2

The greatest common divisor of m and n must be deficient because n must be deficient: σ(g) < g×2 where g = gcd(m,σ(m))
In fact, we can use even stronger restrictions: if (m,n) is an amicable pair, m < n, then σ(m)=σ(n)=m+n

Dividing it by m and n respectively, we get:
σ(m)/m-1 = n/m = x
σ(n)/n-1 = m/n = 1/x

1/x = σ(n)/n-1 > σ(g)/g-1 because n = g × n₁, n₁ > 1
σ(g)/g - 1 < 1/x

Let's try to transform it to the form F₁(g, m) < F₂(g, m) so we can use it for filtering values of m:

σ(g)/g - 1 < 1/x
σ(g)/g - 1 < 1/(σ(m)/m-1)
σ(g) - g < g/(σ(m)/m-1)
σ(g) - g < (g × m)/(σ(m)-m)
(σ(g) - g)×(σ(m) - m) < g × m

(σ(g) - g) × (σ(m) - m) < g × m

This formula has another equivalent form (pointed to me by David Einstein): m/σ(m) + g/σ(g) > 1
But it involves floating-point divisions and thus can be unreliable when dealing with huge numbers.

The general meaning of both formulas is "the more abundant m is, the more deficient n (and g) must be".

So, the final set of requirements for m is:

• m is not divisible by 6
• m is abundant
• (σ(g) - g) × (σ(m) - m) < g × m where g = gcd(m, σ(m))

Very few numbers satisfy all three requirements combined, so it significantly speeds up the search.

4. Conclusion and source code

My implementation of the algorithm described above finds all amicable pairs < 10¹³ in just ~24 minutes on Core i7-4770K (8 parallel threads).
The asymptotic complexity of the algorithm seems to be O(N×log(log(N))).
I have no strict proof of it, but run times for N=10¹⁰...10¹⁷ align almost perfectly with C×N×log(log(N)).

Download the source code: Amicable.zip
The latest source code is also available on GitHub: https://github.com/SChernykh/Amicable