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

May 10th, 2016

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

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

- 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:

If p from (1) is integer and prime, then we've found a breeder (u×p, a). We know that p > p

(u×p

or

(u×p

which limits the set of values of 'a' to check against. The larger p

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

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

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:

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:

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.

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

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

It's easy to see that D=1 for a=2

For example, for a=2

D=2

if D=2 then p = 2

if D=4 then p = 2

...

if D=2

so we can take

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:

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

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:

where b = S(a)/g

d divides a/g

for a=2

Applying (8) and (9) to a=2

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:

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:

Exhaustive search is still manageable, and I was able to find most (3, 1) amicable pairs with a < 10

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

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

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

for some very small

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

21 Euler 1750

11498355=3

12024045=3

u=2581=29×89

a=4455=3

(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

41 TeRiele 1982

3489655839073826955942266055=3

3650550738768150894849977145=3

Another example:

u=397×479×141403=26889618689

a=3

(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

61 Chernykh 2016

3623171640394727379955351072502990492618130497038578416351325756458993451782755208025=3

3639906844276458222356992416948731788196043339981019887096597468959958779504615339175=3

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

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

December 20th, 2015

2. Generating new breeders from known breeders

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

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.

and

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

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.

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

- 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^{6}...10^{7}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.

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

- Step 1: a ≤ 10
^{13}, 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
^{13}. It took a couple of days to finish this step. - Step 3 took ~15 minutes.

Source code for the described algorithm: BreederSearch_2015.zip

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):

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

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

If t = u + S(u) is a prime number then for n = 1,2,3,...

A=au×t

B=a×t

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

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.

November 11th, 2015

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

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

and then to find over 1 million new amicable pairs above 10

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

σ(m) - sum of all divisors of m, including 1 and m.

- 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

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

For example, if m=p

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
_{1}, ..., p_{n}in a loop to get factor powers (k_{1}, ..., k_{n}). Integer division is slow.

Divisions can be eliminated using a more tricky version of this method (see the link above).

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
^{2}, ..., p^{n}

- Iterate from m to m × p where p doesn't divide m.

Then m_{1}= m × p and σ(m_{1}) = σ(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_{1}= (m × p^{k}) × p - one multiplication

σ(m_{1}) = σ(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_{1}).

As a result, we get m and σ(m) already calculated. Perfect!

Knowing that, we can try to find ways to prove that σ(n) ≠ σ(m) way before completing factorization.

Indeed, we know that n = n

Hence σ(n) = σ(n

Lower bound for σ(n): σ(n) = σ(n

Upper bound for σ(n): σ(n) = σ(n

- n
_{2}= p_{1}^{k1}×...×p_{n}^{kn}where p_{0}≤ p_{1}< p_{2}< ... < p_{n}, p_{0}- smallest prime factor we haven't tried yet. - n
_{2}≤ t for some number t, t≥n/n_{1}

Herman te Riele in "Computation of All the Amicable Pairs Below 10

He uses analytical estimation of max(σ(n

This table is of course has an upper limit for p

My implementation calculates it for all p

What about cases where n

I derived these formulas:

max(σ(n

max(σ(n

n

If σ(n

This technique prunes numbers very effectively. It discards most candidates after finding just a few smallest prime factors.

Indeed, n

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

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

x

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

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

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

Let's try to transform it to the form F

σ(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

This formula has another equivalent form (pointed to me by David Einstein):

But it involves floating-point divisions and thus can be unreliable when dealing with huge numbers.

The general meaning of both formulas is

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

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

Download the source code: Amicable.zip

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