# Numerical test of abc hypothesis (yes, the one)

Hi, Habr.

On Geektimes Habr there were already several articles about the abc-hypothesis (for example, in 2013 and in 2018 ). The story itself about the theorem, which for many years cannot be proved for many years, and then cannot be checked for as many years, certainly deserves at least a feature film. But in the shadow of this wonderful story, the theorem itself is considered too superficially, although it is no less interesting. Even by the fact that the abc hypothesis is one of the few unsolved problems of modern science, the task of which can be understood even by a fifth-grader. If this hypothesis is indeed true, then the proof of other important theorems easily follows from it, for example, the proof of Fermat's theorem .

Without claiming Mochizuki’s laurels, I also decided to try and decided to check with the computer how much the equality promised in the hypothesis was fulfilled. Actually, why not - modern processors are not only to play games - why not use a computer in its main (compute - compute) purpose ...

Who cares what happened, I ask for cat.

## Formulation of the problem

Start over. What is the theorem about? As Wikipedia says (the wording in the English version is a little more understandable), for mutually simple (without common divisors) numbers a, b and c, such that a + b = c, for any ε> 0 there is a limited number of triples a + b = c, such that: The function rad is called the radical , and denotes the product of prime factors of a number. For example, rad (16) = rad (2 * 2 * 2 * 2) = 2, rad (17) = 17 (17 is a prime number), rad (18) = rad (2 * 3 * 3) = 2 * 3 = 6, rad (1,000,000) = rad (2 ^ 6 ⋅ 5 ^ 6) = 2 * 5 = 10.

Actually, the essence of the theorem is that the number of such triples is quite small. For example, if we take at random ε = 0.2 and equality 100 + 27 = 127: rad (100) = rad (2 * 2 * 5 * 5) = 10, rad (27) = rad (3 * 3 * 3) = 3, rad (127) = 127, rad (a * b * c) = rad (a) * rad (b) * rad (s) = 3810, 3810 ^ 1.2 is clearly greater than 127, the inequality is not satisfied. But there are exceptions, for example, for equality 49 + 576 = 625, the condition of the theorem is fulfilled (those who wish can check it themselves).

The next key moment for us is of these equalities, according to the theorem, a limited number. Those. This means that you can just try them all on your computer. In the end, this gives us the Nobel Prize is quite an interesting programming problem.

So let's get started.

## Source

The first version was written in Python, and although this language is too slow for such calculations, writing code on it is easy and simple, which is convenient for prototyping.

Getting the radical : decompose the number into prime factors, then remove the repetitions, converting the array into a set. Then we just get the product of all the elements.

``````defprime_factors(n):
factors = []
# Print the number of two's that divide nwhile n % 2 == 0:
factors.append(int(2))
n = n / 2# n must be odd at this point so a skip of 2 ( i = i + 2) can be usedfor i in range(3, int(math.sqrt(n)) + 1, 2):
# while i divides n , print i ad divide nwhile n % i == 0:
factors.append(int(i))
n = n / i
# Condition if n is a prime number greater than 2if n > 2:
factors.append(int(n))
return set(factors)
result = 1for num in prime_factors(n):
result *= num
return result
``````

Mutually simple numbers : decompose numbers into factors, and just check the intersection of sets.

``````defnot_mutual_primes(a,b,c):
fa, fb, fc = prime_factors(a), prime_factors(b), prime_factors(c)
return len(fa.intersection(fb)) == 0and len(fa.intersection(fc)) == 0and len(fb.intersection(fc)) == 0``````

Check : we use already created functions, everything is simple.

``````defcheck(a,b,c):
S = 1.2# Eps=0.2if c > (rad(a)*rad(b)*rad(c))**S and not_mutual_primes(a, b, c):
print("{} + {} = {} - PASSED".format(a, b, c))
else:
print("{} + {} = {} - FAILED".format(a, b, c))
check(10, 17, 27)
check(49, 576, 625)
``````

Those interested can experiment on their own by copying the above code into any online Python editor. Of course, the code works as expectedly slowly, and going through all the triples to at least a million would be too long. Below is a optimized version under the spoiler, it is recommended to use it.

The final version was rewritten in C ++ using multithreading and some optimization (working on C with intersection of sets would be too hardcore, although probably faster). The source code is under a spoiler, it can be compiled in the free g ++ compiler, the code works under Windows, OSX and even on Raspberry Pi.

C ++ Code
``````// To compile: g++ abc.cpp -O3 -fopenmp -oabc#include<string.h>#include<math.h>#include<stdbool.h>#include<stdint.h>#include<stdio.h>#include<vector>#include<set>#include<map>#include<algorithm>#include<time.h>typedefunsignedlongint valType;
typedefstd::vector<valType> valList;
typedefstd::set<valType> valSet;
typedef valList::iterator valListIterator;
std::vector<valList> valFactors;
valList factors(valType n){
valList results;
valType z = 2;
while (z * z <= n) {
if (n % z == 0) {
results.push_back(z);
n /= z;
} else {
z++;
}
}
if (n > 1) {
results.push_back(n);
}
return results;
}
valList unique_factors(valType n){
valList results = factors(n);
valSet vs(results.begin(), results.end());
valList unique(vs.begin(), vs.end());
std::sort(unique.begin(), unique.end());
return unique;
}
valList f = valFactors[n];
double result = 1;
for (valListIterator it=f.begin(); it<f.end(); it++) {
result *= *it;
}
return result;
}
boolnot_mutual_primes(valType a, valType b, valType c){
valList res1 = valFactors[a], res2 = valFactors[b], res3; // = valFactors[c];
valList c12, c13, c23;
set_intersection(res1.begin(),res1.end(), res2.begin(),res2.end(), back_inserter(c12));
if (c12.size() != 0) returnfalse;
res3 = valFactors[c];
set_intersection(res1.begin(),res1.end(), res3.begin(),res3.end(), back_inserter(c13));
if (c13.size() != 0) returnfalse;
set_intersection(res2.begin(),res2.end(), res3.begin(),res3.end(), back_inserter(c23));
return c23.size() == 0;
}
intmain(){
time_tstart_t, end_t;
time(&start_t);
int cnt=0;
double S = 1.2;
valType N_MAX = 10000000;
printf("Getting prime factors...\n");
valFactors.resize(2*N_MAX+2);
for(valType val=1; val<=2*N_MAX+1; val++) {
valFactors[val] = unique_factors(val);
}
time(&end_t);
printf("Done, T = %.2fs\n", difftime(end_t, start_t));
printf("Calculating...\n");
#pragma omp parallel for reduction(+:cnt)for(int a=1; a<=N_MAX; a++) {
for(int b=a; b<=N_MAX; b++) {
int c = a+b;
printf("%d + %d = %d\n", a,b,c);
cnt += 1;
}
}
}
printf("Done, cnt=%d\n", cnt);
time(&end_t);
floatdiff_t = difftime(end_t, start_t);
printf("N=%lld, T = %.2fs\n", N_MAX, diff_t);
}
``````

For those who are too lazy to install the C ++ compiler, a slightly optimized Python version is shown, which can be run in any online editor (I used https://repl.it/languages/python ).

Python version
``````from __future__ import print_function
import math
import time
import multiprocessing
prime_factors_list = []
defprime_factors(n):
factors = []
# Print the number of two's that divide nwhile n % 2 == 0:
factors.append(int(2))
n = n / 2# n must be odd at this point so a skip of 2 ( i = i + 2) can be usedfor i in range(3, int(math.sqrt(n)) + 1, 2):
# while i divides n , print i ad divide nwhile n % i == 0:
factors.append(int(i))
n = n / i
# Condition if n is a prime number greater than 2if n > 2:
factors.append(int(n))
return factors
result = 1for num in prime_factors_list[n]:
result *= num
return result
defnot_mutual_primes(a,b,c):
fa, fb, fc = prime_factors_list[a], prime_factors_list[b], prime_factors_list[c]
return len(fa.intersection(fb)) == 0and len(fa.intersection(fc)) == 0and len(fb.intersection(fc)) == 0defcalculate(N):
S = 1.2
cnt = 0for a in range(1, N):
for b in range(a, N):
c = a+b
print("{} + {} = {}".format(a, b, c))
cnt += 1
print("N: {}, CNT: {}".format(N, cnt))
return cnt
if __name__ == '__main__':
t1 = time.time()
NMAX = 100000
prime_factors_list = *(2*NMAX+2)
for p in range(1, 2*NMAX+2):
prime_factors_list[p] = set(prime_factors(p))
calculate(NMAX)
print("Done", time.time() - t1)
``````

## results

There are very few a, b, c tracks.

Some results are shown below:
N = 10 : 1 “troika”, execution time <0.001c
1 + 8 = 9

N = 100 : 2 “triples”, execution time <0.001c
1 + 8 = 9
1 + 80 = 81

N = 1000 : 8 "triples", execution time <0.01c
1 + 8 = 9
1 + 80 = 81
1 + 242 = 243
1 + 288 = 289
1 + 512 = 513
3 + 125 = 128
13 + 243 = 256
49 + 576 = 625

N = 10000 : 23 "triples", runtime 2s

Triplets A, B, C up to 10,000
```1 + 8 = 9 1 + 80 = 81 1 + 242 = 243 1 + 288 = 289 1 + 512 = 513 1 + 2400 = 2401 1 + 4374 = 4375 1 + 5831 = 5832 1 + 6560 = 6561 1 + 6655 = 6656 1 + 6859 = 6860 3 + 125 = 128 5 + 1024 = 1029 10 + 2187 = 2197 11 + 3125 = 3136 13 + 243 = 256 49 + 576 = 625 1331 + 9604 = 10935 81 + 1250 = 1331 125 + 2187 = 2312 243 + 1805 = 2048 289 + 6272 = 6561 625 + 2048 = 2673```

N = 100000 : 53 "triples", run time 50c

Three A, B, C to 100,000
```1 + 8 = 9 1 + 80 = 81 1 + 242 = 243 1 + 288 = 289 1 + 512 = 513 1 + 2400 = 2401 1 + 4374 = 4375 1 + 5831 = 5832 1 + 6560 = 6561 1 + 6655 = 6656 1 + 6859 = 6860 1 + 12167 = 12168 1 + 14336 = 14337 1 + 57121 = 57122 1 + 59048 = 59049 1 + 71874 = 71875 3 + 125 = 128 3 + 65533 = 65536 5 + 1024 = 1029 7 + 32761 = 32768 9 + 15616 = 15625 9 + 64000 = 64009 10 + 2187 = 2197 11 + 3125 = 3136 13 + 243 = 256 28 + 50625 = 50653 31 + 19652 = 19683 37 + 32768 = 32805 49 + 576 = 625 49 + 16335 = 16384 73 + 15552 = 15625 81 + 1250 = 1331 121 + 12167 = 12288 125 + 2187 = 2312 125 + 50176 = 50301 128 + 59049 = 59177 169 + 58880 = 59049 243 + 1805 = 2048 243 + 21632 = 21875 289 + 6272 = 6561 343 + 59049 = 59392 423 + 16384 = 16807 507 + 32768 = 33275 625 + 2048 = 2673 1331 + 9604 = 10935 1625 + 16807 = 18432 28561 + 89088 = 117649 28561 + 98415 = 126976 3584 + 14641 = 18225 6561 + 22000 = 28561 7168 + 78125 = 85293 8192 + 75843 = 84035 36864 + 41261 = 78125```

With N = 1,000,000, we have only 102 “triples”, the full list is given under the spoiler.

Three A, B, C up to 1,000,000
```1 + 8 = 9 1 + 80 = 81 1 + 242 = 243 1 + 288 = 289 1 + 512 = 513 1 + 2400 = 2401 1 + 4374 = 4375 1 + 5831 = 5832 1 + 6560 = 6561 1 + 6655 = 6656 1 + 6859 = 6860 1 + 12167 = 12168 1 + 14336 = 14337 1 + 57121 = 57122 1 + 59048 = 59049 1 + 71874 = 71875 1 + 137780 = 137781 1 + 156249 = 156250 1 + 229375 = 229376 1 + 263168 = 263169 1 + 499999 = 500000 1 + 512000 = 512001 1 + 688127 = 688128 3 + 125 = 128 3 + 65533 = 65536 5 + 1024 = 1029 5 + 177147 = 177152 7 + 32761 = 32768 9 + 15616 = 15625 9 + 64000 = 64009 10 + 2187 = 2197 11 + 3125 = 3136 13 + 243 = 256 13 + 421875 = 421888 17 + 140608 = 140625 25 + 294912 = 294937 28 + 50625 = 50653 31 + 19652 = 19683 37 + 32768 = 32805 43 + 492032 = 492075 47 + 250000 = 250047 49 + 576 = 625 49 + 16335 = 16384 49 + 531392 = 531441 64 + 190269 = 190333 73 + 15552 = 15625 81 + 1250 = 1331 81 + 123823 = 123904 81 + 134375 = 134456 95 + 279841 = 279936 121 + 12167 = 12288 121 + 255879 = 256000 125 + 2187 = 2312 125 + 50176 = 50301 128 + 59049 = 59177 128 + 109375 = 109503 128 + 483025 = 483153 169 + 58880 = 59049 243 + 1805 = 2048 243 + 21632 = 21875 289 + 6272 = 6561 338 + 390625 = 390963 343 + 59049 = 59392 423 + 16384 = 16807 507 + 32768 = 33275 625 + 2048 = 2673 864 + 923521 = 924385 1025 + 262144 = 263169 1331 + 9604 = 10935 1375 + 279841 = 281216 1625 + 16807 = 18432 2197 + 583443 = 585640 2197 + 700928 = 703125 3481 + 262144 = 265625 3584 + 14641 = 18225 5103 + 130321 = 135424 6125 + 334611 = 340736 6561 + 22000 = 28561 7153 + 524288 = 531441 7168 + 78125 = 85293 8192 + 75843 = 84035 8192 + 634933 = 643125 9583 + 524288 = 533871 10816 + 520625 = 531441 12005 + 161051 = 173056 12672 + 117649 = 130321 15625 + 701784 = 717409 18225 + 112847 = 131072 19683 + 228125 = 247808 24389 + 393216 = 417605 28561 + 89088 = 117649 28561 + 98415 = 126976 28561 + 702464 = 731025 32768 + 859375 = 892143 296875 + 371293 = 668168 36864 + 41261 = 78125 38307 + 371293 = 409600 303264 + 390625 = 693889 62192 + 823543 = 885735 71875 + 190269 = 262144 131072 + 221875 = 352947 132651 + 588245 = 720896```

Alas, the program still works slowly, I didn’t wait for the results for N = 10,000,000, the calculation time is more than an hour (maybe I was mistaken somewhere with the optimization of the algorithm, and it can be done better).

It is even more interesting to look at the results graphically: In principle, it is quite obvious that the dependence of the number of possible triples on N grows noticeably slower than N itself, and it is quite likely that the result will converge to some specific number for each ε. By the way, as ε increases, the number of “triples” decreases noticeably, for example, when ε = 0.4, we have only 2 equalities for N <100000 (1 + 4374 = 4375 and 343 + 59049 = 59392). So in general, it seems that the theorem is indeed satisfied (well, and probably it has already been tested on computers more powerful, and it is possible that all this has long been calculated).

Those interested can experiment on their own, if anyone has results for numbers 10,000,000 or more, I’m happy to add them to the article. Of course, it would be interesting to “calculate” until the moment when the set of “triples” stops growing completely, but it can take a really long time, the speed of calculation seems to depend on N as N * N (and maybe N ^ 3), and the process very long. Nevertheless, an amazing number, and those who wish may well join the search.

Edit: as suggested in the comments, the table with the results already exists in Wikipedia - in the range N up to 10 ^ 18 the number of "triples" is still growing, so the "end" of the set has not yet been found. The more interesting - the intrigue is still preserved.