Home

Discoverers

Statistics

News

Articles

Submit new pairs

List of articles

Amicable pairs: searching for them 10 times faster by Sergei Chernykh. November 6th, 2026
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

---

Amicable pairs: searching for them 10 times faster

Sergei Chernykh
November 6th, 2026

Contents

1. Algorithmic improvements
   1.1. n is a linear function of p
   1.2. Sieving instead of trial division
   1.3. The sieve window
   1.4. How deep to sieve
   1.5. Ranges without amicable pairs
   1.6. Checking the survivors
   1.7. Near-degenerate ranges: the band sieve
   1.8. Large primes

2. GPU-specific improvements
3. CPU-specific improvements
4. Source code


Preface
This article continues "Amicable pairs: searching for them efficiently", so please read that one first.
It describes a new algorithm for the exhaustive search of amicable pairs that was developed to push the search limit up to 10²².
The main idea is simple: the larger members of potential amicable pairs are not checked one by one anymore. Instead, they're sieved, just like in the sieve of Eratosthenes, and only very few of them need further checks.

Both the CPU and the GPU versions of the program use the new algorithm. Compared to the previous versions, on typical work units of the 10²² search:

• The GPU version is 11-44 times faster (GeForce RTX 5070)
• The CPU version is 12-31 times faster (Ryzen 9 9950X, 32 threads). Part of it comes from a fix of an old bug, see section 3.
• Some "difficult" work units became more than 100 times faster

The GPU version is still about 5 times faster than the CPU version on the same work units.
Both versions were checked against all 9,359,869 known amicable pairs below 10²² and found all of them.

Terminology
N, m, n, σ(m) - the same as in the previous article
p - the largest prime factor of m
m₀ - m without its largest prime factor: m = m₀ × p. The new algorithm is for the case when p doesn't divide m₀ (99.998% of all known amicable pairs). Numbers divisible by p² or p³ are still searched with the old algorithm.
range - all numbers m = m₀ × p for a fixed m₀, where p goes through all primes greater than the largest prime factor of m₀ (and m ≤ N)
q - a small prime used for sieving, Q - the largest such prime
g = gcd(m₀, σ(m₀))

1. Algorithmic improvements

1.1. n is a linear function of p

As in the previous article, we iterate from m₀ to m = m₀ × p, and σ(m) is known right away. But let's look at the whole range at once:

σ(m) = σ(m₀) × (p + 1)
n = σ(m) - m = σ(m₀) × (p + 1) - m₀ × p = (σ(m₀) - m₀) × p + σ(m₀)

n = a × p + b, where a = σ(m₀) - m₀ and b = σ(m₀) are the same for the whole range.

So n is a linear function of p! The target value is linear in p too: (m, n) is an amicable pair if and only if σ(n) = σ(m) = b × (p + 1).

Now take any prime power q^k. If q doesn't divide a, then q^k divides n if and only if

p ≡ -b × a^-1 (mod q^k)

This means that all p for which q^k divides n form an arithmetic progression with step q^k. This is exactly what a sieve needs!

If q divides a, there are two cases. If q doesn't divide b, then q never divides n. If q divides both a and b, then q divides gcd(a, b) = gcd(σ(m₀) - m₀, σ(m₀)) = gcd(m₀, σ(m₀)) = g.
So n = g × (a/g × p + b/g) for every p in the range: g is a known factor of n, and the sieve works with a/g × p + b/g.

1.2. Sieving instead of trial division

The old algorithm takes every n of the range and does trial divisions (with all the tricks from the previous article) until it can prove that σ(n) ≠ σ(m).
It's fast, but it's still at least a few divisions per number, and a range can have billions of numbers.

The new algorithm processes the whole range at once. It takes an array with one element per p (I call them positions) and, for every prime q ≤ Q and every power q^k, walks through the arithmetic progression of p where q^k divides n, adding a value to every position it hits.
Just like in the sieve of Eratosthenes, it takes only a couple of additions per position on average and no divisions at all, and every position gets the full information about the small prime factors of its n.

What value is added? In the end, every position must have the sum

S(p) = Σ log(σ(q^e)/q^e)

over all primes q ≤ Q, where q^e is the exact power of q dividing n.
So the progression of q^k adds log(σ(q^k)/q^k) - log(σ(q^k-1)/q^k-1): when q^e divides n, the values added by q, q², ..., q^e add up to log(σ(q^e)/q^e).
The logarithms are scaled and rounded to integers, so the sieve only adds integers.

Note that the sieve doesn't need to know which p are primes. It simply sieves all p not divisible by 2, 3 and 5 (or by 2, 3, 5 and 7), and only the few p that survive the sieve are tested for primality.
So the search doesn't even need a table of all primes up to 10¹¹ anymore.

1.3. The sieve window

After the sieve, we know S(p) = log(σ(n₁)/n₁) for every position, where n = n₁ × n₂, n₁ - factored portion of n (all prime factors ≤ Q) and n₂ - unfactored portion of n (all prime factors > Q).
If (m, n) is an amicable pair, then σ(n₁) × σ(n₂) = σ(m), which can be written as

σ(n₁)/n₁ = (σ(m)/n) / (σ(n₂)/n₂)

All prime factors of n₂ are greater than Q, so n₂ can have only a few of them, and its abundancy is barely above 1: 1 ≤ σ(n₂)/n₂ ≤ ρ, where ρ is the largest abundancy of a number below n with all prime factors greater than Q.
For example, if Q = 32749 and n < 10²², then n₂ has at most 4 prime factors and ρ < 1.00013.

Taking logarithms, we get:

log(σ(m)/n) - log(ρ) ≤ S(p) ≤ log(σ(m)/n)

This is the sieve window. It's very narrow (log(ρ) < 0.00013 in the example above), and σ(m)/n = b × (p + 1)/(a × p + b) changes very slowly with p, so the window is calculated once for a whole block of positions (using the smallest and the largest σ(m)/n of the block, plus a small margin for the rounding errors).
Only p with S(p) inside the window can give an amicable pair. I call them survivors.
Typically less than 1% of all positions survive. All other numbers are rejected with just one comparison, without ever looking at them individually!

1.4. How deep to sieve

A bigger Q makes the window narrower and gives fewer survivors, but it also makes the sieve slower: every prime q costs a modular inverse per range (to find where its progression starts) plus its hits.
So there's no single best Q. The program has several sieve depths, from Q = 1021 to Q = 32749, and chooses one for every range: short ranges are sieved with a small Q, long ranges with a bigger one.

The best depth also depends on how much abundancy n still needs. Let K be the known part of n: g, multiplied by 2 if n is always even. Then σ(n)/n ≥ σ(K)/K for every p, and the difference

W = log(σ(m)/n) - log(σ(K)/K) (for the largest p of the range)

tells how much abundancy the unknown factors of n must provide. I call it the window height.
If W is smaller than the width of the window, then every p where n has no small prime factors except those of K falls into the window, and the number of survivors explodes. A bigger Q makes the window narrower, so ranges with a small W are sieved deeper.

1.5. Ranges without amicable pairs

σ(m)/n = b × (p + 1)/(a × p + b) grows with p, and K divides n for all p of the range.
So if W < 0, i.e. σ(K)/K > σ(m)/n even for the largest p of the range, then σ(n)/n ≥ σ(K)/K > σ(m)/n for all p. This means σ(n) > σ(m), and the whole range can be skipped without sieving it!
It's the same idea as the condition (σ(g) - g) × (σ(m) - m) < g × m from the previous article ("the more abundant m is, the more deficient n must be"), but applied to a whole range at once.

This can be taken further with the check sum divisibility idea from the previous article.
Take a prime power q^e that exactly divides K. If one more factor q would make the abundancy of n too big for every p of the range, then q^e is the exact power of q in n for every possible amicable p. σ(q^e) must divide σ(n) = σ(m) = b × (p + 1), so

σ(q^e)/gcd(σ(q^e), b) divides p + 1

All amicable p of the range must be in this arithmetic progression. And for p in this progression, the exact powers of some other small primes in n can be known in advance too, so K gets bigger, W gets smaller, and the same steps can be repeated.
If W becomes negative, the range has no amicable pairs. This proves many ranges empty: in some regions of the search, whole work units consist of such ranges.

The same check works inside the sieve: if a prime q divides n exactly once, then q + 1 = σ(q) must divide σ(m) = b × (p + 1), i.e. (q + 1)/gcd(q + 1, b) must divide p + 1.
The GPU version removes the positions where this fails in the ranges with the narrowest windows, which halves the number of survivors there.

1.6. Checking the survivors

A survivor still needs a full check: p must be a prime, and σ(n) must be equal to σ(m). The trial division from the previous article can do it, but the sieve already knows a lot about n, so the check can be much faster.

Small factors. S(p) is the exact sum of log(σ(q^e)/q^e) over the small prime factors of n. If the factors found so far explain all but δ of S(p), then every small prime factor q that is left satisfies log(1 + 1/q) ≤ δ, i.e. q ≥ 1/δ - 1/2, and the search for the next factor can jump right to that prime. When δ = 0, n₁ is completely factored.

The rest of n. Now n = n₁ × n₂, where n₁ is fully known and n₂ has no prime factors ≤ Q. σ(n₂) must be equal to σ(m)/σ(n₁), so we also know the value

D = σ(m)/σ(n₁) - n₂ - 1

which must be the sum of all divisors of n₂ except 1 and n₂ itself (and if D < 0, the check is over right away). D tells a lot about n₂.
Every prime factor q of n₂ gives a divisor n₂/q, so n₂/q ≤ D, i.e.

all prime factors of n₂ are ≥ n₂/D

and trial division can start at n₂/D instead of Q. For the hardest numbers D is tiny compared to n₂, so trial division skips thousands of primes.

D also limits the number of prime factors of n₂:

• 4 or more prime factors: D ≥ n₂^3/4 (the largest divisor n₂/q₁ alone, where q₁ is the smallest prime factor)
• 3 different prime factors q₁×q₂×q₃: D > q₁×q₂ + q₁×q₃ + q₂×q₃ ≥ 3×n₂^2/3
• q₁²×q₂: D = (q₁ + 1)×(q₁ + q₂) ≥ 1.8899×n₂^2/3
• q³: D = q + q²

So if D < 3×n₂^2/3, then n₂ can only be a prime (D = 0), q² (D = q), q₁×q₂ (D = q₁ + q₂, the quadratic equation from section 2.3 of the previous article), q³ or q₁²×q₂.
Each of these cases can be checked directly, without any trial division at all.

1.7. Near-degenerate ranges: the band sieve

Some ranges are much harder than the others. If W is very close to 0, the known part K of n has almost exactly the abundancy n needs.
Then almost every p for which n has no other small prime factors survives the sieve, no matter how big Q is, and these survivors are hard to check: n₂ has no small factors, and D is tiny.
One example is the ranges with m₀ = 2×5×7×11×23×41×2083 and its multiples: the old algorithm needed more than 15 minutes for the first such range alone.

But the previous section gives a very strong bound here: all prime factors of n₂ are greater than L = n₂/D.
It can also be shown that the smallest one is at most j×L + j + 2, where j ≤ log(n₂)/log(L) is the number of prime factors of n₂.
And in such a range D/n₂ < e^W - 1 for every number, so L > 1/(e^W - 1) for all of them: we know in advance where to look for the smallest prime factor of n₂!

So these ranges get a second sieve, the band sieve. It sieves the primes of a "band" around L (again, every prime gives an arithmetic progression of p) and marks every position with the smallest band prime that divides its n.
If a survivor's smallest prime factor must be inside the band, but no band prime divides its n, the survivor is rejected immediately. This rejects more than 90% of the survivors in these ranges without any trial division, and the others start their trial division right at the band prime.
The band sieve made such work units more than 10 times faster on both CPU and GPU.

1.8. Large primes

When p is very large, m₀ = m/p is small. The search handles these numbers separately: for p > 10¹¹ (with N = 10²²), it takes a precomputed list of all numbers m₀ that can be a part of an amicable pair m₀ × p (about 2.4 billion numbers) and checks them against all primes p from an interval [P₁, P₂].

The sieve works here too: every m₀ from the list is a range of p from P₁ to P₂.
Since all of them share the same interval of p, the primes of [P₁, P₂] are sieved only once, into a bitmap, and the survivors don't need a primality test at all.
Only the GPU version uses the sieve for large primes so far.

2. GPU-specific improvements

At first sight, a sieve is a bad fit for a GPU: it's all about writing to scattered places in memory, and GPUs prefer the opposite.
The first straight-forward version of the GPU sieve was 3.6 times slower than the old GPU code! Here's what made it fast:

• The sieve runs entirely in the fast on-chip memory. A work-group of 1024 threads sieves 8×1024 positions at a time in its local memory, and only the survivors ever leave the chip.
• Every thread gets a fixed set of arithmetic progressions to sieve. The program plans this at start-up, so that all threads of a warp do similar amounts of work and no warp has to wait for the others: the busiest warp now does only about 10% more work than an average one (it was 70% at first).
• The GPU doesn't need a table of primes at all, so the start-up time dropped from about 20 seconds to less than 1 second.
• Survivors are checked with the old GPU code, extended with the shortcuts from section 1.6. The numbers are sorted by the amount of work they need before the checks, because all threads of a warp wait for the slowest one.
• p is tested for primality only after the first trial divisions of n, and all these tests run in parallel in a separate kernel. For p < 2³², one strong probable prime test to base 2 is enough, together with a table of the 2256 base-2 strong pseudoprimes below 2³² that are not divisible by 3, 5 or 7.
• The CPU prepares the next batch of ranges while the GPU is busy, and the GPU kernels are chained, so the CPU waits for the GPU only once per batch.
• After compiling its kernels, the program sieves around a few dozen known amicable pairs and stops with an error if it doesn't find all of them. A buggy driver or GPU can't silently miss amicable pairs.

On GeForce RTX 5070, the new GPU version is 11-44 times faster than the old one on typical work units of the 10²² search, and up to 266 times faster on work units with near-degenerate ranges.

3. CPU-specific improvements

The CPU version uses the same algorithm (only for ranges with at least 512 primes: in shorter ranges the setup of the sieve costs more than it saves) and a few CPU-specific tricks:

• A bug fix in the old code: the 128-bit version of the first check of n didn't check for D < 0, so up to 0.1% of the numbers went through trial division by all 1.36 million primes of its table. With N = 10²², this bug made the old version up to 17 times slower (1830 instead of 107 nanoseconds per number in one test), so part of the speed-up comes from this fix.
• The table of primes up to 10¹¹ is calculated by all threads in parallel: the start-up time dropped from 11.5 to 0.9 seconds.
• Long ranges are split into chunks that all threads share. Before, one thread had to check a whole range alone (up to 4 billion primes p), and at the end of a work unit all other threads could wait for it for a long time.
• The sieve depth adapts to the real work: if checking the survivors takes more than 40% of the time of a chunk, the next chunk of the range is sieved with a bigger Q, and if it takes less than 10%, with a smaller one.
• Long ranges skip all p divisible by 7 too (48 of every 210 numbers are sieved instead of 8 of every 30), short ranges use the simpler layout that is faster to set up.
• The small factors of the survivors are found with the help of the sieve sums (section 1.6), not by trial division by all small primes.
• The hottest loops use AVX2 or AVX-512 when the CPU has them (chosen at run time). For example, the trial division of n₂ tests 8 primes per AVX-512 instruction, using double-precision numbers where all products are exact. On ARM and RISC-V CPUs the program uses portable code.

On Ryzen 9 9950X (all 32 threads), the new CPU version is 12-31 times faster than the old one on typical work units of the 10²² search, and more than 95 times faster on work units with near-degenerate ranges.

4. Source code

The source code of the CPU version
The source code of the GPU version

---

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