Чисельна перевірка 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)

Читайте також  Кращі емулятори Android на Windows

Результати

Трійок 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 я так і не дочекався, час обчислення становить більше години (можливо я десь помилився з оптимізацією алгоритму можна зробити краще).

Читайте також  Як створити свій сайт з нуля — 10 кращих безкоштовних конструкторів сайтів

Ще цікавіше подивитися результати графічно:

В принципі, цілком очевидно, що залежність кількості можливих трійок від N помітно зростає повільніше самого N, і цілком ймовірно, що результат буде сходитися до якогось конкретного числа для кожного ε (ризикну висловити гіпотезу, що в даному випадку воно не перевищить 256:). До речі, при збільшенні ε число «трійок» помітно скорочується, наприклад при ε=0.4 маємо всього 2 рівності при N<100000 (1 + 4374 = 4375 і 343 + 59049 = 59392). Так що в цілому, схоже що теорема виконується (ну і напевно її вже перевіряли на комп’ютерах потужніший, і можливо, все це вже давно пораховано).

Бажаючі можуть поекспериментувати самостійно, якщо у когось будуть результати для чисел 10000000 і вище, я з задоволенням додам їх до статті. Зрозуміло, було б цікаво «дорахувати» до того моменту, коли безліч «трійок» зовсім перестане рости, але це може зайняти реально довгий час, швидкість розрахунку схоже залежить від N як N*N (а може і N^3), і процес дуже довгий. Але тим не менш, дивовижне поруч, і бажаючі можуть долучитися до пошуку.

Степан Лютий

Обожнюю технології в сучасному світі. Хоча частенько і замислююся над тим, як далеко вони нас заведуть. Не те, щоб я прям і знаюся на ядрах, пікселях, коллайдерах і інших парсеках. Просто приходжу в захват від того, що може в творчому пориві вигадати людський розум.

You may also like...

Залишити відповідь

Ваша e-mail адреса не оприлюднюватиметься. Обов’язкові поля позначені *