Численная проверка abc-гипотезы (да, той самой)

На geektimes habr было уже несколько статей про 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 = %.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;       if (c > pow(valRads[a]*valRads[b]*valRads[c], S) && not_mutual_primes(a,b,c)) {         printf("%d + %d = %d\n", a,b,c);         cnt += 1;       }     }   }   printf("Done, cnt=%d\n", cnt);      time(&end_t);   float diff_t = difftime(end_t, start_t);   printf("N=%lld, T = %.2fs\n", N_MAX, diff_t); } 

Для тех кому лень устанавливать компилятор С++, приведена слегка оптимизированная Python-версия, запустить которую можно в любом онлайн редакторе (я использовал https://repl.it/languages/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.001c
1 + 8 = 9

N=100: 2 «тройки», время выполнения <0.001c
1 + 8 = 9
1 + 80 = 81

N=1000: 8 «троек», время выполнения <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 «тройки», время выполнения 2с

N=100000: 53 «тройки», время выполнения 50c

При N=1000000 имеем всего лишь 102 «тройки», полный список приведен под спойлером.

Увы, программа работает все равно медленно, результатов для N=10000000 я так и не дождался, время вычисления составляет больше часа (возможно я где-то ошибся с оптимизацией алгоритма, и можно сделать лучше).

Еще интереснее посмотреть результаты графически:

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

Желающие могут поэкспериментировать самостоятельно, если у кого будут результаты для чисел 10000000 и выше, я с удовольствием добавлю их к статье.