22. Горн Гавриила
В математике есть некоторые парадоксы, связанные с бесконечностью.
Так называемый “Горн Гавриила” (Gabriel's Horn) - это фигура, образованная функцией y = 1/x, если создать из нее 3х-мерную фигуру:
Интерес данного объекта в том, что “горн” имеет конечный объем - но при этом бесконечную площадь поверхности. Если бы такая фигура существовала в реальности, ее можно было бы заполнить краской, а вот покрасить целиком не получилось бы.
Действительно, формулы площади (V) и объема (A) “горна” легко найти в Википедии:
Из них видно, что при a = ∞, V становится равным ℼ, а A действительно стремится к бесконечности.
23. Построение графиков функций
Очень многие закономерности гораздо легче исследовать графически. Очень просто и эффективно строить графики с помощью языка Python и библиотеки matplotlib.
Выведем графики функций y = x2 и у = 0.5x2.
import matplotlib.pyplot as plt
plt.title('Graph')
plt.xlabel('X')
plt.ylabel('Y')
x_values = range(0,10)]
y_values1 = map(lambda x: 0.5*x*x, x_values)
y_values2 = map(lambda x: x*x, x_values)
plt.axis([0, 10, 0, 20]) # xmin, xmax, ymin, ymax
plt.plot(x_values, y_values1, 'ro', label='0.5*x*x')
plt.plot(x_values, y_values2, marker='o', linestyle='--', label='x*x')
plt.legend()
plt.show()
Результат работы программы выглядит так:
Как можно видеть, мы имеем массив входных данных x_values = range(0,20), и 2 массива выходных данных y_values1 и y_values2, вычисляемых по соответствующим формулам. Все остальное за нас делает библиотека matplotlib. Команда plt.plot выводит график на экран, причем можно задать разные варианты оформления. Команда выводит подсказку (“легенду”), и является опциональной. Команды plt.title, plt.xlabel и plt.ylabel задают названия графика и осей. Команда plt.show выводит окно с графиком на экран.
Рассмотрим более сложный пример.
Выведем график так называемой “функции распределения простых чисел”, показывающий количество простых чисел меньше либо равных числу x.
Программа вывода графика выглядит так:
import matplotlib.pyplot as plt
import math
def is_prime(n):
if n % 2 == 0 and n > 2:
return False
for i in xrange(3, int(math.sqrt(n)) + 1, 2):
if n % i == 0:
return False
return True
def get_primes(n):
cnt = 0
for i in xrange(1, n):
if is_prime(i):
cnt += 1
return cnt
N = 400
x_values = range(2, N)
y_values = map(lambda n: get_primes(n), x_values)
plt.axis([0, N, 0, N/4])
plt.plot(x_values, y_values, 'ro', markersize=1)
plt.show()
Кстати, согласно Википедии, еще в 18м веке Гауссом и Лежандром было высказано предположение что функция распределения простых чисел выглядит так:
Это легко проверить графически, добавив строчку вывода второго графика:
plt.plot(x_values, map(lambda n: n/math.log(n), x_values), marker='o', linestyle='--')
Окончательный график выглядит так:
С помощью библиотеки numpy строить графики еще проще.
График функции sin(t) в диапазоне 0..2𝜋:
import matplotlib.pyplot as plt
import numpy as np
t = np.arange(0.0, 2.0, 0.01)
s = np.sin(2*np.pi*t)
plt.plot(t, s)
plt.show()
Важно заметить, что t здесь - не одно число, а целый массив - в данном примере используется возможность numpy работать непосредственно с массивами.
И наконец, с помощью matplotlib можно строить даже 3х-мерные графики функции двух переменных, для этого используется модуль mplot3d:
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import math
fig = plt.figure()
ax = fig.gca(projection='3d')
# Axis
ax.set_xlabel('X')
ax.set_xlim(0, 20)
ax.set_ylabel('Y')
ax.set_ylim(0, 20)
ax.set_zlabel('Z')
ax.set_zlim(0, 20)
# Graph
for x in range(0,20):
for y in range(0,20):
z = 10*math.sin(0.1*x + 0.1*y)
ax.scatter(x, y, z, c='r', s=1)
plt.show()
Результат показан на рисунке:
Интересно отметить кросплатформенность кода на Python - приведенная выше программа работает и на Windows и на OSX без каких-либо изменений.
Приложение 1 - Вычисления с помощью видеокарты
Еще 20 лет назад, во времена процессоров 80386, пользователям приходилось покупать математический сопроцессор, позволяющий быстрее выполнять вычисления с плавающей точкой. Сейчас такой сопроцессор покупать уже не надо - благодаря прогрессу в игровой индустрии, даже встроенная видеокарта компьютера имеет весьма неплохую вычислительную мощность. Например, даже бюджетный видеочип Intel Graphics 4600 имеет 20 вычислительных блоков, что превышает количество ядер “основного” процессора. Разумеется, каждое ядро GPU по отдельности слабее CPU, но здесь как раз тот случай, когда количество дает преимущество над качеством. Вычисления с помощью GPU сейчас очень популярны - от майнинга биткоинов до научных расчетов, диапазон ценовых решений также различен, от “бесплатной” встроенной видеокарты до NVIDIA Tesla ценой более 100тыс рублей. Поэтому интересно посмотреть, как же это работает.
Есть две основные библиотеки для GPU-расчетов - NVidia CUDA и OpenCL. Первая обладает большими возможностями, однако работает только с картами NVIDIA. Библиотека OpenCL работает с гораздо большим числом графических карт, поэтому мы рассмотрим именно ее.
Основной принцип GPU-расчетов - параллельность вычислений. Данные, хранящиеся в “глобальной памяти” (global & constant memory) устройства, обрабатываются модулями (каждый модуль называется “ядром”), каждый из которых работает параллельно с другими. Модуль имеет и свою собственную память для промежуточных данных (private memory). Так это выглядит в виде блок-схемы:
Таким образом, если задача может быть разбита на небольшие блоки, параллельно обрабатывающие небольшой фрагмент блока данных, такая задача может эффективно быть решена на GPU.
Рассмотрим пример: необходимо проверить, какие числа в массиве являются простыми. Массив может быть большим, например миллион элементов. Такая задача идеальна для распараллеливания: каждое число может быть проверено независимо от предыдущего.
Для решения такой задачи с помощью OpenCL необходимо выполнить ряд шагов.
Написать код микроядра (kernel):
Этот код будет запускаться непосредственно на графических процессорах видеокарты. Код пишется на языке C. В данном примере мы для упрощения храним код прямо в виде строки в программе.
const char *KernelSource = "\n" \
"__kernel void primes( \n" \
" __global unsigned int* input, \n" \
" __global unsigned int* output) \n" \
"{ \n" \
" unsigned int i = get_global_id(0); \n" \
" //printf(\"Task-%d\\n\", i); \n" \
" output[i] = 0; \n" \
" unsigned int val = input[i]; \n" \
" for(unsigned int p=2; p<=val/2; p++) { \n" \
" if (val % p == 0) \n" \
" return; \n" \
" } \n" \
" output[i] = 1; \n" \
"} \n" \
"\n";
Суть кода проста. Массив input хранит числа, которые нужно проверить, функция get_global_id возвращает индекс задачи, которую выполняет данное ядро. Мы берем число с нужным индексом, проверяем его на простоту, и записываем 0 или 1 в зависимости от результата, в массив output.
Инициализировать подготовку вычислений:
int gpu = 1;
clGetDeviceIDs(NULL, gpu ? CL_DEVICE_TYPE_GPU : CL_DEVICE_TYPE_CPU, 1, &device_id, NULL);
cl_context context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);
cl_command_queue commands = clCreateCommandQueue(context, device_id, 0, &err);
На этом этапе можно выбрать где будут производиться вычисления, на основном процессоре или на GPU. Для отладки удобнее основной процессор, окончательные расчеты быстрее на GPU.
Подготовить данные
#define DATA_SIZE 1024
cl_uint *data = (cl_uint*)malloc(sizeof(cl_uint) * DATA_SIZE);
cl_uint *results = (cl_uint*)malloc(sizeof(cl_uint) * DATA_SIZE);
Загрузить данные и программу из основной памяти в GPU
cl_program program = clCreateProgramWithSource(context, 1, (const char **) & KernelSource, NULL, &err);
clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
cl_kernel kernel = clCreateKernel(program, "primes", &err);
cl_mem output = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(cl_uint) * count, NULL, NULL);
clEnqueueWriteBuffer(commands, input, CL_TRUE, 0, sizeof(cl_uint) * count, data, 0, NULL, NULL);
clSetKernelArg(kernel, 0, sizeof(cl_mem), &output);
clGetKernelWorkGroupInfo(kernel, device_id, CL_KERNEL_WORK_GROUP_SIZE, sizeof(local), &local, NULL);
Запустить вычисления на GPU и дождаться их завершения
global = DATA_SIZE;
clEnqueueNDRangeKernel(commands, kernel, 1, NULL, &global, &local, 0, NULL, NULL);
clFinish(commands);
Загрузить результаты обратно из GPU в основную память
clEnqueueReadBuffer( commands, output, CL_TRUE, 0, sizeof(cl_uint) * count, results, 0, NULL, NULL );
Освободить данные
free(data);
free(results);
clReleaseMemObject(input);
clReleaseMemObject(output);
clReleaseProgram(program);
clReleaseKernel(kernel);
clReleaseCommandQueue(commands);
clReleaseContext(context);
Как можно видеть, процесс довольно-таки громоздкий, но оно того стоит. Для примера, проверка простоты 250000 чисел заняла на процессоре Core i5 около 6 секунд. И всего лишь 0.5 секунд заняло выполнение вышеприведенного кода на встроенной видеокарте. Для дешевого нетбука с процессором Intel Atom этот же код выполнялся 34 секунды на основном процессоре, и 6 секунд на GPU. Т.е. разница весьма прилична.
Разумеется, еще раз стоит повторить, что “игра стоит свеч” лишь в том случае, если задача хорошо распараллеливается на небольшие блоки, в таком случае выигрыш будет заметен.
Владельцы видеокарт NVIDIA (особенно игровых и достаточно мощных) могут также посмотреть в сторону библиотеки NVIDIA CUDA, расчеты с ее помощью должны быть еще быстрее.
Приложение 2 - Быстродействие языка Python
Язык Python очень удобен своей краткостью и лаконичностью, возможностью использования большого количества сторонних библиотек. Однако, один из его минусов, который может быть ключевым для математических расчетов - это быстродействие. Python это интерпретатор, он не создает exe-файл, что разумеется, сказывается на скорости выполнения программы.
Рассмотрим простой пример: рассчитаем сумму квадратов чисел от 1 до 1000000. Также выведем время выполнения программы.
Программа на языке Python выглядит так:
import time
start_time = time.time()
s = 0
for x in xrange(1,1000001):
s += x*x
print("Sum={}, T={}s".format(s, time.time() - start_time))
Результаты работы:
Sum = 333333833333500000, T = 0.47s
Учитывая, что чисел всего миллион, не так уж и быстро. Попробуем ускорить программу, для этого по возможности используем функции встроенных библиотек. Они зачастую написаны на С, и работают быстрее.
import time
start_time = time.time()
l = xrange(1000001)
s = sum(x*x for x in l)
print("Sum = {}, T = {}s".format(s, time.time() - start_time))
Результаты работы:
Sum = 333333833333500000, T = 0.32s
Быстрее, но лишь чуть-чуть.
И наконец, призываем “тяжелую артиллерию”: перепишем программу на языке Cи.
Код выглядит так:
#include <stdio.h>
#include <time.h>
int main()
{
clock_t start = clock();
unsigned long long int sum = 0, i;
for(i=1; i<1000001; i++) {
sum += i*i;
}
clock_t end = clock();
printf("Sum = %llu, T = %fs", sum, (float)(end - start)/CLOCKS_PER_SEC);
return 0;
}
Как можно видеть, он ненамного сложнее python-версии. Перед запуском программы, ее надо скомпилировать, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest". Результат очевиден: T = 0.007 секунд. И еще чуть-чуть: добавляем флаг оптимизации по скорости, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest" -O3. Результат: 0.0035 секунд, разница в быстродействии более 100 раз!