Чисельна перевірка abc-гіпотези (так, тієї самої)
На Geektimes було вже кілька статей про abc-гіпотезу (наприклад у 2013 та 2018 роках). Сама історія про теорему, яку спочатку багато років не можуть довести, а потім стільки ж років не можуть перевірити, безумовно заслуговує як мінімум, художнього фільму. Але в тіні цієї чудової історії, сама теорема розглядається занадто поверхово, хоча вона не менш цікава. Вже хоча б тим, що abc-гіпотеза — одна з небагатьох невирішених проблем сучасної науки, постановку завдання якої зможе зрозуміти навіть п’ятикласник. Якщо ж ця гіпотеза справді вірна, то з неї легко слід доказ інших важливих теорем, наприклад доказ теореми Ферма.
Не претендуючи на лаври Мотідзукі, я теж вирішив спробувати вирішив перевірити за допомогою комп’ютера, наскільки виконуються обіцяні в гіпотезі рівності. Власне, чому б ні — сучасні процесори адже не тільки для того щоб грати в ігри — чому б не використовувати комп’ютер за своїм основним (compute — обчислювати) призначення…
Постановка задачі
Почнемо з початку. Про що власне, теорема? Як свідчить Вікіпедія (формулювання в англійській версії трохи більш зрозуміла), для взаємно простих (не мають спільних дільників) чисел a, b і з, таких, що a+b=c, для будь-якого ε>0 існує обмежене число трійок a+b=c, таких що:
Функція rad називається радикалом, і позначає добуток простих множників числа. Наприклад, rad(16) = rad(2*2*2*2) = 2, rad(17) = 17 (17 просте число), rad(18) = rad(2*3*3) = 2*3 = 6, rad(1000000) = rad(2^6 ⋅ 5^6) = 2*5 = 10.
Власне, суть теореми у тому, що кількість таких трійок досить мало. Наприклад, якщо взяти навмання ε=0.2 і рівність 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(с) = 3810, 3810^1.2 явно більше 127, нерівність не виконується. Але бувають і винятки, наприклад для рівності 49 + 576 = 625 умова теореми виконується (бажаючі можуть перевірити самостійно).
Наступний ключовий для нас момент — цих рівностей, згідно з теоремою, обмежене число. Тобто це означає, що їх можна просто спробувати перебрати на комп’ютері. У підсумку, це дає нам Нобелівську премію цілком цікаву задачу з програмування.
Отже, приступимо.
Вихідний код
Перша версія була написана на Python, і хоча ця мова занадто повільний для подібних розрахунків, писати код на ньому легко і просто, що зручно для прототипування.
Отримання радикала: розкладаємо число на прості множники, потім прибираємо повтори, перетворюючи масив безліч. Потім просто отримуємо добуток всіх елементів.
def prime_factors(n):
factors = []
# Print the number of two's that divide n
while 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 used
for i in range(3, int(math.sqrt(n)) + 1, 2):
# while i divides n , print i ad divide n
while n % i == 0:
factors.append(int(i))
n = n / i
# Condition if n is a prime number greater than 2
if n > 2:
factors.append(int(n))
return set(factors)
def rad(n):
result = 1
for num in prime_factors(n):
result *= num
return result
Взаємно прості числа: розкладаємо числа на множники, і просто перевіряємо перетин множин.
def not_mutual_primes(a,b,c):
fa, fb, fc = prime_factors(a), prime_factors(b), prime_factors(c)
return len(fa.intersection(fb)) == 0 and len(fa.intersection(fc)) == 0 and len(fb.intersection(fc)) == 0
Перевірка: використовуємо вже створені функції, тут все просто.
def check(a,b,c):
S = 1.2 # Eps=0.2
if 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)
Бажаючі можуть поекспериментувати самостійно, скопіювавши вищенаведений код в будь онлайн-редактор мови Python. Зрозуміло, код працює очікувано повільно, і перебір всіх трійок хоча б до мільйона був би надто довгим. Нижче під спойлером є оптимізована версія, рекомендується використовувати її.
Остаточна версія була переписана на с++ З використанням багатопоточності і деякої оптимізації (працювати на Сі з перетином множин було б занадто хардкорно, хоча ймовірно і швидше). Вихідний код під спойлером, його можна скомпілювати в безкоштовному компіляторі g++, код працює під Windows, OSX і навіть на Raspberry Pi.
Код на С++
// 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>
typedef unsigned long int valType;
typedef std::vector<valType> valList;
typedef std::set<valType> valSet;
typedef valList::iterator valListIterator;
std::vector<valList> valFactors;
std::vector<double> valRads;
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;
}
double rad(valType n) {
valList f = valFactors[n];
double result = 1;
for (valListIterator it=f.begin(); it<f.end(); it++) {
result *= *it;
}
return result;
}
bool not_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) return false;
res3 = valFactors[c];
set_intersection(res1.begin(),res1.end(), res3.begin(),res3.end(), back_inserter(c13));
if (c13.size() != 0) return false;
set_intersection(res2.begin(),res2.end(), res3.begin(),res3.end(), back_inserter(c23));
return c23.size() == 0;
}
int main()
{
time_t start_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);
valRads.resize(2*N_MAX+2);
for(valType val=1; val<=2*N_MAX+1; val++) {
valFactors[val] = unique_factors(val);
valRads[val] = rad(val);
}
time(&end_t);
printf("Done, T = %.2fsn", 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;
if (c > pow(valRads[a]*valRads[b]*valRads[c], S) && not_mutual_primes(a,b,c)) {
printf("%d + %d = %dn", a,b,c);
cnt += 1;
}
}
}
printf("Done, cnt=%dn", cnt);
time(&end_t);
float diff_t = difftime(end_t, start_t);
printf("N=%lld, T = %.2fsn", N_MAX, diff_t);
}
Для тих кому ліньки встановлювати компілятор С++, наведена злегка оптимізована Python-версія, запустити яку можна в будь-якому онлайн редакторі (я використовував https://repl.it/languages/python).
Python-версія
from __future__ import print_function
import math
import time
import multiprocessing
prime_factors_list = []
rad_list = []
def prime_factors(n):
factors = []
# Print the number of two's that divide n
while 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 used
for i in range(3, int(math.sqrt(n)) + 1, 2):
# while i divides n , print i ad divide n
while n % i == 0:
factors.append(int(i))
n = n / i
# Condition if n is a prime number greater than 2
if n > 2:
factors.append(int(n))
return factors
def rad(n):
result = 1
for num in prime_factors_list[n]:
result *= num
return result
def not_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)) == 0 and len(fa.intersection(fc)) == 0 and len(fb.intersection(fc)) == 0
def calculate(N):
S = 1.2
cnt = 0
for a in range(1, N):
for b in range(a, N):
c = a+b
if c > (rad_list[a]*rad_list[b]*rad_list[c])**S and not_mutual_primes(a, b, c):
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 = [0]*(2*NMAX+2)
rad_list = [0]*(2*NMAX+2)
for p in range(1, 2*NMAX+2):
prime_factors_list[p] = set(prime_factors(p))
rad_list[p] = rad(p)
calculate(NMAX)
print("Done", time.time() - t1)
Результати
Трійок a,b,c дійсно дуже мало.
Деякі результати наведені нижче:
N=10: 1 «трійка», час виконання <0.001 c
1 + 8 = 9
N=100: 2 «трійки», час виконання <0.001 c
1 + 8 = 9
1 + 80 = 81
N=1000: 8 «трійок», час виконання <0.01 c
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 «трійки», час виконання 2с
Трійки A,B,C до 100001 + 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 «трійки», час виконання 50c
Трійки A,B,C до 1000001 + 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
При N=1000000 маємо всього лише 102 «трійки», повний список наведено під спойлером.
Трійки A,B,C до 10000001 + 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
На жаль, програма працює все одно повільно, результатів для N=10000000 я так і не дочекався, час обчислення становить більше години (можливо я десь помилився з оптимізацією алгоритму можна зробити краще).
Ще цікавіше подивитися результати графічно:
В принципі, цілком очевидно, що залежність кількості можливих трійок від N помітно зростає повільніше самого N, і цілком ймовірно, що результат буде сходитися до якогось конкретного числа для кожного ε (ризикну висловити гіпотезу, що в даному випадку воно не перевищить 256:). До речі, при збільшенні ε число «трійок» помітно скорочується, наприклад при ε=0.4 маємо всього 2 рівності при N<100000 (1 + 4374 = 4375 і 343 + 59049 = 59392). Так що в цілому, схоже що теорема виконується (ну і напевно її вже перевіряли на комп’ютерах потужніший, і можливо, все це вже давно пораховано).
Бажаючі можуть поекспериментувати самостійно, якщо у когось будуть результати для чисел 10000000 і вище, я з задоволенням додам їх до статті. Зрозуміло, було б цікаво «дорахувати» до того моменту, коли безліч «трійок» зовсім перестане рости, але це може зайняти реально довгий час, швидкість розрахунку схоже залежить від N як N*N (а може і N^3), і процес дуже довгий. Але тим не менш, дивовижне поруч, і бажаючі можуть долучитися до пошуку.