Многие слышали о великом и ужасном быстром преобразовании Фурье (БПФ / FFT — fast fourier transform) — но как его можно применять для решения практических задач за исключением JPEG/MPEG сжатия и разложения звука по частотам (эквалайзеры и проч.) — зачастую остается неясным вопросом.Недавно я наткнулся на интересную реализацию игры «Жизнь» Конвея, использующую быстрое преобразование Фурье — и надеюсь, оно поможет вам понять применимость этого алгоритма в весьма неожиданных местах.
Правила
Вспомним правила классической «Жизни» — на поле с квадратными клетками, живая клетка погибает если у неё больше 3 или меньше 2 соседей, и если у пустой клетки ровно 3 соседей — она рождается. Соответственно, для эффективной реализации алгоритма нужно быстро считать количество соседей вокруг клетки.Алгоритмов для этого существует целая куча (в том числе и моя JS реализация). Но есть у задачи и математическое решение, которое может давать хорошую скорость для плотно заполненных полей, и быстро уходит в отрыв с ростом сложности правил и площади/объема суммирования — например в Smoothlife-подобных (1,2), и 3D вариантах.
Реализация на БПФ
Идея алгоритма следующая:- Формируем матрицу суммирования (filter), где 1 стоят в ячейках, сумму которых нам нужно получить (8 единиц, остальные нули). Выполняем над матрицей прямое преобразование Фурье (это нужно сделать только 1 раз).
- Выполняем прямое преобразование Фурье над матрицей с содержимым игрового поля.
- Перемножаем каждый элемент результата на соответствующий элемент матрицы «суммирования» из пункта 1.
- Выполняем обратное преобразование Фурье — и получаем матрицу с нужной нам суммой количества соседей для каждой клетки.
Реализация на C++
//Author: Mark VandeWettering http://brainwagon.org/2012/10/11/crazy-programming-experiment-of-the-evening/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <complex.h> #include <unistd.h> #include <fftw3.h> #define SIZE (512) #define SHIFT (18) fftw_complex *filter ; fftw_complex *state ; fftw_complex *tmp ; fftw_complex *sum ; int main(int argc, char *argv[]) { fftw_plan fwd, rev, flt ; fftw_complex *ip, *jp ; int x, y, g ; srand48(getpid()) ; filter = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ; state = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ; tmp = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ; sum = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ; flt = fftw_plan_dft_2d(SIZE, SIZE, filter, filter, FFTW_FORWARD, FFTW_ESTIMATE) ; fwd = fftw_plan_dft_2d(SIZE, SIZE, state, tmp, FFTW_FORWARD, FFTW_ESTIMATE) ; rev = fftw_plan_dft_2d(SIZE, SIZE, tmp, sum, FFTW_BACKWARD, FFTW_ESTIMATE) ; /* initialize the state */ for (y=0, ip=state; y<SIZE; y++) { for (x=0; x<SIZE; x++, ip++) { *ip = (fftw_complex) (lrand48() % 2) ; } } /* initialize and compute the filter */ for (y=0, ip=filter; y<SIZE; y++, ip++) { for (x=0; x<SIZE; x++) { *ip = 0. ; } } #define IDX(x, y) (((x + SIZE) % SIZE) + ((y+SIZE) % SIZE) * SIZE) filter[IDX(-1, -1)] = 1. ; filter[IDX( 0, -1)] = 1. ; filter[IDX( 1, -1)] = 1. ; filter[IDX(-1, 0)] = 1. ; filter[IDX( 1, 0)] = 1. ; filter[IDX(-1, 1)] = 1. ; filter[IDX( 0, 1)] = 1. ; filter[IDX( 1, 1)] = 1. ; fftw_execute(flt) ; for (g = 0; g < 1000; g++) { fprintf(stderr, "generation %03d\r", g) ; fftw_execute(fwd) ; /* convolve */ for (y=0, ip=tmp, jp=filter; y<SIZE; y++) { for (x=0; x<SIZE; x++, ip++, jp++) { *ip *= *jp ; } } /* go back to the sums */ fftw_execute(rev) ; for (y=0, ip=state, jp=sum; y<SIZE; y++) { for (x=0; x<SIZE; x++, ip++, jp++) { int s = (int) round(creal(*ip)) ; int t = ((int) round(creal(*jp))) >> SHIFT ; if (s) *ip = (t == 2 || t == 3) ; else *ip = (t == 3) ; } } /* that's it! dump the frame! */ char fname[80] ; sprintf(fname, "frame.%04d.pgm", g) ; FILE *fp = fopen(fname, "wb") ; fprintf(fp, "P5\n%d %d\n%d\n", SIZE, SIZE, 255) ; for (y=0, ip=state; y<SIZE; y++) { for (x=0; x<SIZE; x++, ip++) { int s = ((int) creal(*ip)) ; fputc(255*s, fp) ; } } fclose(fp) ; } fprintf(stderr, "\n") ; return 0 ; }
Для сборки нужна библиотека FFTW. Ключи для сборки в gcc:
gcc life.cpp -lfftw3 -lm -lstdc++в Visual Studio нужны изменения в работе с комплексными числами.
Изображения во время выполнения алгоритма:
Исходная матрица фильтра (начало координат — левый верхний угол, поле «закольцовано»):
Действительная часть фильтра после прямого БПФ:
Исходное поле — 1 глайдер:
БПФ исходного поля, действительная и мнимая части:
После умножения на матрицу фильтра:
После обратного БПФ — получаем суммы:
Результат вполне ожидаемый:
«Закольцовывание» поля получается автоматически из-за БПФ.
Плюшки БПФ
- Вы можете суммировать любое количество элементов с любыми коэффициентами — а время работы остается фиксированным, N2logN. Т.е. если для классической жизни — обычные алгоритмы на заполненных полях все еще достаточно быстрые, то с увеличением площади/объема суммирования — они становятся все медленнее, а скорость работы БПФ остается фиксированной.
- БПФ — уже написан, отлажен и оптимизирован идеально — с использованием всех возможностей процессора: sse, sse2, avx, altivec и neon.
- Вы легко можете использовать все процессоры и видеокарты взяв готовые многопроцессорные и CUDA/OpenCL реализации БПФ. Опять же, об оптимизации БПФ вам заботится не нужно. FFTW кстати поддерживает многопроцессорное выполнение (--enable-openmp во время сборки), чем больше поле — тем меньше потери на синхронизацию.
- Тот же подход применим и к 3D пространству.
Сравниваем скорость работы с «наивной» реализацией
FFTW собран в режиме одинарной точности (--enable-float, дает прирост скорости примерно в 1.5 раза). Процессор — Core2Duo E8600. gcc 4.7.2 -O3Исходный текст наивной реализации
for (int i=8; i<1000; i++) { delta[i][0] = (lrand48() % 21)-10; delta[i][1] = (lrand48() % 21)-10; } #define NUMDELTA 18 for(int frame=0;frame<10;frame++) { double start = ticks(); for (y=0; y<SIZE; y++) for (x=0; x<SIZE; x++) { int sum=0; for(int i=0;i<NUMDELTA;i++) sum+=array[readpage][(SIZE+y+delta[i][0]) % SIZE][(SIZE+x+delta[i][1]) % SIZE]; if(array[readpage][y][x]) array[1-readpage][y][x] = (sum>3 || sum<2)?0:1;else array[1-readpage][y][x] = (sum==3)?1:0; } readpage = 1-readpage; printf("NaiveDelta: %10.10f\n", ticks()-start); }
| Конфигурация теста | При каком количестве суммируемых клеток скорость FFT и наивной реализации равна |
| 1 ядро, 512x512 | 29 |
| 2 ядра, 512x512 | 18 |
| 1 ядро, 4096x4096 | 88 |
| 2 ядра, 4096x4096 | 65 |
Update:
Как оказалось, FFTW при указании флага FFTW_ESTIMATE использует далеко не оптимальные планы вычисления FFT. При указании FFTW_MEASURE — скорость сильно выросла, и ситуация выглядит уже радостнее (при суммировании менее 18 клеток — наивная реализация резко становится в 3 раза быстрее, потому тут все упирается в 18):
| Конфигурация теста | При каком количестве суммируемых клеток скорость FFT и наивной реализации равна |
| 1 ядро, 512x512 | 18 |
| 2 ядра, 512x512 | 18 |
| 1 ядро, 4096x4096 | 28 |
| 2 ядра, 4096x4096 | 19 |
Чтобы ни говорили — математика в программировании иногда может пригодиться в самых неожиданных местах.
И да пребудет с вами сила FFT!