Во-вторых, доказать эту гипотезу чрезвычайно сложно. На момент написания книги, только один японский математик Синъити Мотидзуки попытался доказать ее, но это доказательство занимает 500 страниц, и его корректность пока что никто не смог подтвердить (у некоторых ученых были сомнения в его корректности, да и разобраться в таком объеме непросто).
Также считается, что на основе этой гипотезы могут быть доказаны другие известные теоремы, так что ее доказательство было бы важным для математики в целом.
Разумеется, мы здесь не претендуем на доказательство гипотезы, но можем проверить некоторые значения на языке Python.
Для начала, возьмем функцию нахождения простых множителей числа N.
def prime_factors(n):
factors = []
while n % 2 == 0:
factors.append(int(2))
n = n / 2
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
if n > 2:
factors.append(int(n))
return factors
Как можно видеть, функция последовательно делит числа, и если остаток от деления равен 0, результат заносится в массив. Отдельно проверяются два частных случая - если число четное (делится на 2) и если число простое (ни на что не делится кроме себя). В результате функция возвращает массив чисел, напримем для 8 это будет [2,2,2] (2*2*2 = 8).
Следующий шаг - напишем функцию вычисления радикала. Для этого мы воспользуемся множеством (set), которое хранит неповторяющиеся элементы массива prime_factors.
def rad(n):
result = 1
for num in set(prime_factors(n)):
result *= num
return result
Третья функция - проверка на то, что 3 числа являются взаимно простыми, и у них нет общих делителей.
def not_mutual_primes(a,b,c):
fa, fb, fc = set(prime_factors(a)), set(prime_factors(b)), set(prime_factors(c))
return len(fa.intersection(fb)) == 0 and len(fa.intersection(fc)) == 0 and
len(fb.intersection(fc)) == 0
Здесь мы пользуемся функцией нахождения пересечения множеств языка Python.
Написанных 3х функций достаточно, чтобы написать функцию вывода троек чисел A,B,C:
def calculate(N):
S = 1.2
cnt = 0
for a in range(1, N):
for b in range(1, N):
if b < a: continue
c = a+b
if not_mutual_primes(a, b, c):
if c > (rad(a*b*c))**S:
print("{} + {} = {}".format(a, b, c))
cnt += 1
print("N: {}, CNT: {}".format(N, cnt))
return cnt
Как написано в условии задачи, количество таких троек должно быть ограниченно. Мы можем вывести график количества троек от максимального значения N, тогда будет примерно видно, как быстро он растет.
import matplotlib.pyplot as plt
x_values = range(1,2000,25)
y_values = list(map(calculate, x_values))
plt.plot(x_values, y_values, 'ro', label='count(N), S=1.2')
plt.legend()
plt.show()
Как можно видеть, количество троек реально мало - в диапазоне до 2000 их всего 10:
1 + 8 = 9
1 + 80 = 81
1 + 242 = 243
1 + 288 = 289
1 + 512 = 513
3 + 125 = 128
5 + 1024 = 1029
13 + 243 = 256
49 + 576 = 625
81 + 1250 = 1331
На графике это хорошо видно, более того, видно что рост графика замедляется, и число троек скорее всего, действительно ограничено для любого диапазона чисел.
Приложение 1 - Вычисления с помощью видеокарты
Еще 20 лет назад, во времена процессоров 80386, пользователям приходилось покупать математический сопроцессор, позволяющий быстрее выполнять вычисления с плавающей точкой. Сейчас такой сопроцессор покупать уже не надо - благодаря прогрессу в игровой индустрии, даже встроенная видеокарта компьютера имеет весьма неплохую вычислительную мощность. Например, даже бюджетный видеочип Intel Graphics 4600 имеет 20 вычислительных блоков, что превышает количество ядер “основного” процессора. Разумеется, каждое ядро GPU по отдельности слабее CPU, но здесь как раз тот случай, когда количество дает преимущество над качеством. Вычисления с помощью GPU сейчас очень популярны - от майнинга биткоинов до научных расчетов, диапазон ценовых решений также различен, от “бесплатной” встроенной видеокарты до NVIDIA Tesla ценой более 100тыс рублей. Поэтому интересно посмотреть, как же это работает.
Есть две основные библиотеки для GPU-расчетов - NVidia CUDA и OpenCL. Первая обладает большими возможностями, однако работает только с картами NVIDIA. Библиотека OpenCL работает с гораздо большим числом графических карт, поэтому мы рассмотрим именно ее.
Основной принцип GPU-расчетов - параллельность вычислений. Данные, хранящиеся в “глобальной памяти” (global & constant memory) устройства, обрабатываются модулями (каждый модуль называется “ядром”, или “kernel”), каждый из которых работает параллельно с другими. Модуль имеет и свою собственную память для промежуточных данных (private memory). Так это выглядит в виде блок-схемы:
Таким образом, если задача может быть разбита на небольшие блоки, параллельно обрабатывающие небольшой фрагмент блока данных, такая задача может эффективно быть решена на GPU.
Для запуска OpenCL-программы из Python, необходимо поставить библиотеку pyopencl, скачать ее дистрибутив можно со страницы https://wiki.tiker.net/PyOpenCL/Installation/Windows. Для установки библиотеки под Windows необходимо скачать файл и ввести команду pip install pyopencl-2018.1.1+cl12-cp27-cp27m-win_amd64.whl для 64-разрядной версии Python, или pip install pyopencl‑2018.1.1+cl12‑cp27‑cp27m‑win32.whl для 32-разрядной.
Для тестирования pyopencl можно запустить программу, выводящую информацию о системе:
import pyopencl as cl
for plat in cl.get_platforms():
print 'Platform: {}'.format(plat.name)
print 'Version: ' + plat.version
devices = plat.get_devices(cl.device_type.ALL)
print 'Devices:'
for dev in devices:
print('\t{} ({})'.format(dev.name, dev.vendor))
flags = [('Version', dev.version),
('Type', cl.device_type.to_string(dev.type)),
('Memory (global), MB', str(dev.global_mem_size/(1024*1024))),
('Memory (local), KB', str(dev.local_mem_size/1024)),
('Max work item dims', str(dev.max_work_item_dimensions)),
('Max work group size', str(dev.max_work_group_size)),
('Max compute units', str(dev.max_compute_units)),
('Driver version', dev.driver_version),
('Device available', str(bool(dev.available))),
('Compiler available', str(bool(dev.compiler_available)))]
for name, flag in flags:
print '\t\t{0:<25}{1:<10}'.format(name + ':', flag)
Если pyopencl и драйвер видеокарты установлен корректно, мы увидим на экране примерно такой вывод (программа запускалась на Macbook Pro):
Platform: Apple
Version: OpenCL 1.2 (Sep 12 2017 16:28:17)
Devices:
Iris Pro (Intel)
Version: OpenCL 1.2
Type: GPU
Memory (global), MB: 1536
Memory (local), KB: 64
Max work item dims: 3
Max work group size: 512
Max compute units: 40
Driver version: 1.2 (Oct 4 2017 01:28:36)
Device available: True
Compiler available: True
Программа удобна своей кросс-платформенностью, один и тот же код может работать и на OSX и на Windows без каких-либо изменений.
Рассмотрим пример: сформировать массив простых чисел от 1 до N. Для решения такой задачи можно использовать алгоритм “решето Эратосфена”, но мы в учебных целях будем решать задачу напрямую, проверяя каждое число отдельно.
Код программы:
import pyopencl as cl
import pyopencl.array as cl_array
import numpy, time
type = cl.device_type.GPU # cl.device_type.CPU
platform = cl.get_platforms()[0]
devices = platform.get_devices(device_type=type)
ctx = cl.Context(devices)
queue = cl.CommandQueue(ctx)
prg = cl.Program(ctx, """
__kernel void primes(__global int *results)
{
int gid = get_global_id(0);
if (gid < 2) return;
int num = gid + 1;
for(int p=2; p<=num/2 + 1; p++) {
if ((num % p) == 0) {
return;
}
}
int val = atomic_add(&results[0], 1);
results[val+1] = num;
}
""").build()
N = 500000
results = numpy.zeros(N, dtype=numpy.int32)
dest_dev = cl_array.to_device(queue, results)
size = (results.shape[0]-1,)
prg.primes(queue, size, None, dest_dev.data)
res = dest_dev.get()
count = res[0]
primes = numpy.sort(res[1:count+1])
print "Found prime numbers:", count
# for p in primes:
# print p
Разберем текст программы подробнее.
1. Мы создаем контекст устройства , передав ему в качестве параметра тип устройства device_type.GPU. Переменная типа CommandQueue хранит очередь команд, которые будут посланы на устройство.
2. Класс cl.Program получает в качестве параметра программу ядра (kernel) и компилирует ее с помощью вызова функции build. Функция primes написана на языке Си, и будет выполняться параллельно на всех ядрах видеокарты.
3. С помощью функции numpy.zeros мы создаем заполненный нулями массив размера N, затем копируем его на видеокарту с помощью вызова dest_dev = cl_array.to_device. Важно понимать, что существуют две разных копии массива - один в памяти компьютера, второй на видеокарте, с которым и будут выполняться вычисления.
4. Вызовом функции prg.primes выполняется параллельный расчет на видеокарте. При этом в качестве параметра функция получает массив, общий для всех экземпляров функции. Библиотека OpenCL сама распараллеливает вызовы функций, общее количество вызовов равно числу элементов массива, а для получения конкретного идентификатора мы используем функцию get_global_id. В нашем случае идентификатор соответствует числу, которое мы хотим проверить.
5. Если простое число найдено, мы кладем его в глобальный массив results. Но т.к. несколько ядер работают с одним массивом одновременно, мы должны быть уверены что данные будут сохранены корректно, для этого мы используем функцию atomic_add. Первую ячейку массива мы таким образом используем для хранения количества найденных простых чисел.
6. Когда выполнение программы завершено, мы загружаем данные обратно с помощью функции dest_dev.get(). Т.к. ядра запускались параллельно, то простые числа лежат в массиве вперемешку, чтобы их отсортировать, мы вызываем функцию numpy.sort.
Как можно видеть, процесс довольно-таки громоздкий, но оно того стоит. Однопоточная программа, написанная на Си с таким же кодом, выполнялась 14 секунд. Вышеприведенная программа, запускаемая на видеокарте, выполнила расчет за 1.2 секунды.
Разумеется, еще раз стоит повторить, что “игра стоит свеч” лишь в том случае, если задача хорошо распараллеливается на небольшие блоки, в таком случае выигрыш будет заметен.
Рассмотрим теперь тот же код, но написанный с помощью библиотеки pycuda. Разумеется, выполнить его смогут только владельцы видеокарт NVIDIA.
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule