Когда-то давно я решил «потрогать» Fortran. Единственную задачу которую я придумал — генерация фракталов (заодно и OpenMP в Fortran'е можно было бы попробовать). В процессе написания я часто сталкивался с проблемами, решение которых приходилось додумывать самому (например в интернете не так много примеров использования чисел двойной точности или бинарной записи в файл). Но рано или поздно все проблемы решились, и я хочу написать этот текст, который возможно кому-нибудь поможет.
Писать я буду на диалекте Fortran 90, но с GNU расширениями (те же числа двойной точности).
Немного теории о фракталах
Фрактал (fractus – дробленый, сломанный, разбитый) — в узком смысле, сложная геометрическая фигура, обладающая свойством самоподобия, то есть составленная из нескольких частей, каждая из которых подобна всей фигуре целиком.
Существует большая группа фракталов — алгебраические фракталы. Это название вытекает из принципа построения таких фракталов. Для их построения используют рекурсивные функции. Например алгебраическими фракталами являются: Множество Жюлиа, Множество Мандельброта, Бассейн(фрактал) Ньютона.
Построение алгебраических фракталов
Один из методов построения представляет собой итерационный расчет 
Поведение функции для разных точек комплексной плоскости может иметь разное поведение:
- Стремится к бесконечности
- Стремится к 0
- Принимает несколько фиксированных значений и не выходит за их пределы
- Хаотичное поведение. Отсутствие каких либо тенденций
Один из простейших (как для понимания, так и для реализации) алгоритмов построения алгебраических фракталов известен как «escape time algorithm». Если кратко, то производится итеративный расчет числа для каждой точки комплексной плоскости.
Было доказано, что если 
Я рассмотрю два простых фрактала – Множество Жюлиа и Множество Мандельброта. Расчитываются они по одной и той-же функции, а отличаются лишь константой, где у Жюлиа это постоянная константа, а у Мандельброба константа зависит от точки комплексной плоскости.
Построение Множества Жюлиа
Начало и конец программы простой и тривиальный:
program frac
end
Дальше нам нужно ввести несколько констант:
INTEGER, PARAMETER :: iterations = 2048 !Количество итераций. Чем больше, тем глубже будет идти просчет
REAL(8), PARAMETER :: mag_ratio = 1.0_8 !Приближение
REAL(8), PARAMETER :: x0 = 0.0_8 !Смещение комплексной плоскости
REAL(8), PARAMETER :: y0 = 0.0_8 !
INTEGER, PARAMETER :: resox = 1024 !Разрешение получившегося изображения
INTEGER, PARAMETER :: resoy = resox
REAL(8), PARAMETER :: xshift = (-1.5_8 / mag_ratio) + x0 !Смещение центра комплексной оси в
REAL(8), PARAMETER :: yshift = (-1.5_8 / mag_ratio) + y0 !центр изображения
REAL(8), PARAMETER :: CXmin = -1.5_8 / mag_ratio !Следующие константы растягивают
REAL(8), PARAMETER :: CXmax = 1.5_8 / mag_ratio !комплексную плоскость до размеров resox x resoy
REAL(8), PARAMETER :: CXshag = (DABS(CXmin) + DABS(CXmax)) / resox !Связываем значение одного пикселя и точки на комплексной плоскости
REAL(8), PARAMETER :: CYmin = -1.5_8 / mag_ratio
REAL(8), PARAMETER :: CYmax = 1.5_8 / mag_ratio
REAL(8), PARAMETER :: CYshag = (DABS(CYmin) + DABS(CYmax)) / resoy
COMPLEX(8), PARAMETER :: c = DCMPLX(0.285_8, 0.01_8) !Наша основная константа
И вот тут нужно кое-что объяснить. REAL и COMPLEX это число с плавающей точкой одинарной точности и комплексное число, состоящее из двух REAL. REAL(8) это число двойной точности, а COMPLEX(8), соотвественно, комплексное число, состоящее из двух REAL(8). Так-же к стандартным функциям добавляется литера D (как в случае с ABS), что указывает на использование чисел двойной точности.
Дальше нам нужно ввести несколько переменных:
CHARACTER, DIMENSION(:), ALLOCATABLE :: matrix !Массив точек
COMPLEX(8) :: z
INTEGER :: x, y, iter
REAL(8) :: iter2 !Понадобится нам для сглаживания
ALLOCATE(matrix(resox * resoy * 3))
Использование ALLOCATABLE обязательно! Т.к. в ином случае:
CHARACTER, DIMENSION(0:resox * resoy * 3) :: matrix !Массив точек
Память у нас выделится на стеке, а это не приемлемо при использовании в нескольких потоках. Поэтому мы выделяем память на куче.
Так-же я не использую двумерные массивы, т.к. при выделении на куче массив будет занимать больше места, да и записать его в файл будет сложнее.
Дальше нам нужно указать количество потоков, порождаемых OpenMP:
call omp_set_num_threads(16)
А тут начинаются непосредственно вычисления. В Fortran'е директивы для OpenMP указываются через коментарии (в C/C++ для это есть специльный макрос #pragma).
!$OMP PARALLEL DO
do x = 0, resox-1
do y = 0, resoy-1
iter = 0 !Обнуляем количество итераций
z = DCMPLX(x * CXshag + xshift, y * CYshag + yshift) !Задаем Z начальное значение
do while ((iter .lt. iterations) .and. (CDABS(z) .le. 4.0_8)) !Обычно за bailout берут 2, но я взял 4, т.к. результат лучше сглаживается
z = z**2 + c
iter = iter + 1
end do
iter2 = REAL(iter) - DLOG(DLOG(CDABS(z))) / DLOG(2.0_8) !Формула для сглаживания iter-ln(log2(|Z|))
matrix((x + y * resox) * 3 + 0) = CHAR(NINT(DMOD(iter2 * 7.0_8, 256.0_8))) !Один из способов расскрашивания
matrix((x + y * resox) * 3 + 1) = CHAR(NINT(DMOD(iter2 * 14.0_8, 256.0_8)))
matrix((x + y * resox) * 3 + 2) = CHAR(NINT(DMOD(iter2 * 2.0_8, 256.0_8)))
end do
end do
!$OMP END PARALLEL DO
Самая трудоемкая часть закончена. Дальше нам остается вывести информацию в удобном для восприятия виде — изображение. Я буду использовать формат PNM, т.к. это самый простой формат изображений.
PNM состоит из нескольких секций:
P6
#Комментарий
1024 1024
255
Первая строчка это указание формата записи информации о пикселях:
- P1/P4 — черно-белое изображение
- P2/P5 — серое изображение
- P3/P6 — цветное изображение
Первая P это вариант, где цвет пикселя записывается ASCII символом, а вторая P дает возможность записывать цвет пикселя в бинарном виде (что значительно экономит место на диске).
Следующей строкой идет комментарий, дальше разрешение изображения и количество цветов на пиксель. После количества цветов идет непосредственно информация о пикселях.
Начинаем записывать изображение в файл:
open(unit=8, file = trim("test.pnm"))
write(8, '(a)') "P6" !(a) это текстовая строка
write(8, '(a)') ""
write(8, '(I0, a, I0)') resox, ' ', resoy
write(8, '(I0)') 255
write(8, *) matrix
close(8)
DEALLOCATE(matrix)
DEALLOCATE мы выполняем для приличия, ибо при выходе из программы, ОС все равно вернет занятую память системе.
Для сборки программы можно использовать компилятор от GNU — gfortran:
gfortran -std=gnu frac.f90 -o frac.run -fopenmp
Должно получится следующее изображение:
Как сгенерировать Множество Мандельброта
Это сделать несложно, достаточно заменить c и изменить несколько констант. В данном случаем убираем константу c и добавляем переменную c:
COMPLEX(8) :: z, c
z = DCMPLX(x * CXshag + xshift, y * CYshag + yshift) !Задаем Z начальное значение
c = z
А так-же меняем старые константы:
INTEGER, PARAMETER :: MAIN_RES = 1024
INTEGER, PARAMETER :: resox = (MAIN_RES / 3) * 3
INTEGER, PARAMETER :: resoy = (resox * 2) / 3
REAL(8), PARAMETER :: xshift = (-2.0_8 / mag_ratio) + x0
REAL(8), PARAMETER :: yshift = (-1.0_8 / mag_ratio) + y0
REAL(8), PARAMETER :: CXmin = -2.0_8 / mag_ratio
REAL(8), PARAMETER :: CXmax = 1.0_8 / mag_ratio
REAL(8), PARAMETER :: CXshag = (DABS(CXmin) + DABS(CXmax)) / resox
REAL(8), PARAMETER :: CYmin = -1.0_8 / mag_ratio
REAL(8), PARAMETER :: CYmax = 1.0_8 / mag_ratio
REAL(8), PARAMETER :: CYshag = (DABS(CYmin) + DABS(CYmax)) / resoy
Должно получится следующее изображение:
Используемая литература
en.wikipedia.org/wiki/Fractal
en.wikipedia.org/wiki/Julia_set
en.wikipedia.org/wiki/Mandelbrot_set
en.wikipedia.org/wiki/OpenMP
openmp.org/wp/openmp-specifications/
gcc.gnu.org/onlinedocs/gfortran/
Автор: BratSinot
