Генерация фракталов с помощью SSE/SSE2
ОГЛАВЛЕНИЕ
В статье описывается способ генерации фракталов с помощью SSE/SSE2
Введение
Вероятно, вы раньше слышали о фракталах. Это красивые изображения, подобные показанному выше. Любой фрактал можно описать с помощью итерационных формул. Фрактал генерируется путем вычисления этих формул и нахождения цвета каждого пикселя. Это тяжелая вычислительная задача, и рисование фрактала требует быстрого процессора и тщательно оптимизированной программы.
Фракталы и их свойства изучаются в нелинейной динамике. Помимо ученых, некоторым людям нравится изучать фракталы, потому что изображения изящные и причудливые.
В этой статье опущены математические основы. Необходимые объяснения есть тут:
• Наборы Джулии и Мандельброта, автор – Дэвид Е. Джойс.
• Что такое набор Мандельброта и другие статьи о фракталах, автор –Фабио Цесари.
• Ваш фрактал, сэр. - Спасибо, Джеймс, автор – q123456789.
• Наборы Джулии и Мандельброта. Теория, автор – Мартин Фингстл.
• Теория фракталов, автор – Cygnus Software.
Последние три статьи также дают некоторый код C для генерации фракталов. Если вы впервые знакомитесь с фракталами, прочтите статьи Фабио Цесари. Его объяснения просты и понятны, с таблицами и картинками, стоящими тысячи слов.
Цели и обзор нескольких генераторов фракталов
Надо найти или разработать учебную программу для изучения фракталов и их свойств. Программа будет применяться для демонстрации ученикам, что такое наборы Мандельброта и Джулии, и будет давать им некоторые упражнения. Ожидается, что ученики будут работать с программой всего один урок, поэтому она должна быть простой и интуитивно понятной.
Главные особенности следующие:
• генерация наборов Джулии и Мандельброта с помощью быстрого алгоритма под управлением SSE;
• просмотр орбиты функции (т.е., последовательности итераций) в виде графика и в виде таблицы;
• демонстрация осей и координат любой точки;
• легкое масштабирование (щелчок левой кнопкой мыши – для увеличения, щелчок правой кнопкой мыши – для уменьшения или нечто подобное);
• генерация набора Джулии с параметром C, соответствующим точке набора Мандельброта.
Последний пункт значит, что можно выбрать значение C для набора Джулии из изображения набора Мандельброта. Вы просто держите кнопку мыши, перемещаете курсор над набором Мандельброта, и набор Джулии меняется, в то время как вы делаете это:
Это называется режимом PnP (изображение в изображении) для краткости.
И несколько слов о тех целях, которые не будут реализованы:
• яркие цвета и палитры (если это сделать, то ученики сосредоточатся на схемах окраски, а не на математике);
• любые типы фракталов, отличные от Джулии и Мандельброта;
• поддержка для 3DNow (расширения SIMD, имеющиеся в старых процессорах AMD; современный Athlon XP использует SSE);
• перенос на другие операционные системы и процессоры (x86 и Win9x/XP было достаточно);
• поддержка нескольких процессоров.
Первой попыткой был Chaos Pro. Это сложное бесплатное приложение, способное генерировать наборы Джулии и Мандельброта и много других фракталов. Можно добавлять собственные формулы фракталов с помощью встроенного языка программирования. Несмотря на то, что программа может создавать превосходные трехмерные изображения и анимации, использовать разные палитры и т.д., в ней нет некоторых важных возможностей. Во-первых, нельзя просмотреть орбиту и оси. Во-вторых, набор Мандельброта генерируется за 1.6 секунды на ноутбуке с процессором Pentium M 1.5 ГГц, а рисование набора Джулии с C = -0.12+0.74i занимает 1.4 секунды. Chaos Pro поддерживает режим PnP, но его неудобно использовать из-за медленной отрисовки. Фракталы нельзя полностью отрисовать в реальном времени (нельзя заставить изображение меняться мгновенно при перемещении курсора). В-третьих (главная причина), в Chaos Pro часто возникает ошибка "нарушение доступа" или другие глюки. Это нормально для бесплатного приложения, написанного любителем, но нужно нечто более надежное.
Также есть FracInt, старая добрая программа DOS для рисования фракталов. Она была разработана несколько лет назад, и сейчас доступна 20.0 версия этого приложения. Программа весьма развитая и мощная, с множеством интересных особенностей, таких как возможность просматривать несколько дюжин типов фракталов, воспроизводить музыку фрактала, использовать режим PnP и создавать стереоскопические изображения. Ее основной недостаток – устаревший пользовательский интерфейс. Масштабирование занимает много времени: вы перетаскиваете мышь, чтобы установить нужный размер прямоугольника, затем перемещаете мышь в точку и, наконец, щелкаете дважды (левой кнопкой – для увеличения, правой кнопкой – для уменьшения). Это было неприемлемо, потому что обычный ученик потратит пол-урока на понимание, как использовать интерфейс масштабирования. Доступны видеорежимы до 1600x1200 (1024x768 VESA должен работать с большинством видеоплат), но при запуске FracInt в Windows XP и нажатии Alt+Tab видеодрайвер дал сбой. Вывод на звуковую плату не работал, и пришлось использовать встроенный динамик. Стоит ли говорить, что это старое приложение не поддерживает MMX и SSE.
После нескольких часов навигации по интернету было найдено «Быстрое плавающее фрактальное веселье», или FFFF. Это очень быстрый генератор набора Мандельброта, поддерживающий SSE, SSE2, 3DNow и программы вершин графического процессора. Если бы он был известен раньше, не была бы начата разработка собственного механизма отрисовки фракталов. Однако у FFFF немного функций помимо масштабирования и переключения между разными механизмами расчета, поэтому если вы решите использовать его, то вам придется добавить генератор набора Джулии и разработать более продуманный графический пользовательский интерфейс.
Алгоритм
Ниже приведен алгоритм для генерации набора Мандельброта в псевдокоде:
комплексные P, Z;
для каждой точки P на комплексной плоскости {
Z = 0;
for i = 0 to ITER {
if abs(Z) > 4.0
установить цвет точки на palette[i];
прервать внутренний цикл и перейти к следующей точке;
Z = Z * Z + P;
}
установить цвет точки на черный;
}
P и Z – комплексные переменные. Можно прочитать общие объяснения о комплексных числах в вышеназванных статьях.
Для каждого пикселя P изображения цикл повторяется ITER раз. Выполнение цикла занимает O(X*Y*ITER) времени, где ITER – количество итераций, а размеры изображения равны X x Y пикселей.
abs(Z) – функция, вычисляющая модуль комплексного числа. Представьте круг с радиусом 2.0 и центром в начале координат. Если точка Z находится вне этого круга, abs(Z) будет больше 4.0. В этом случае внутренний цикл должен прерваться, потому что известно, что каждая точка набора Мандельброта лежит внутри круга. И если Z вышло за пределы круга, то точка P не принадлежит набору Мандельброта. Отметим эту точку цветом, показывающим, как быстро она покидает круг при повторении вычисления Z = Z * Z + P. Количество полных итераций i показывает это. Берется i-й цвет из некоторой палитры и продолжается со следующей точки. Палитра – массив цветов. Например, точки, выходящие за пределы круга после первой итерации, отмечаются темно-зеленым цветом; точки, покидающие круг после второй итерации, отмечаются светло-зеленым, и т.д.
Иначе говоря, если цикл выполнился ITER раз, и точка Z вообще не вышла за пределы круга, то точка принадлежит набору Мандельброта. Она окрашивается черным цветом.
Внутренний цикл многократно вычисляет Z<SUB>N+1</SUB> = Z<SUB>N</SUB>2 + P, итерационную формулу для набора Мандельброта. Z инициализируется нулем, но легко понять, что можно установить Z = P, потому что Z<SUB>1</SUB> = Z<SUB>0</SUB>2 + P = 02 + P = P.
Использование палитры – простейший и быстрейший способ окрашивания в наборе Мандельброта, но полученное изображение будет иметь только ITER+1 цветов. Есть методы отрисовки 24-битных изображений. Например, Chaos Pro использует следующую формулу:
color_index = (i - log(log(abs(Z))) * 1.442695 + 2) / 252.0;
Коэффициенты 252.0 и 2.0 уменьшают значение color_index, чтобы оно вмещалось в диапазон 0.1. Затем color_index переводится в 24-битный цвет. Увы, вычисление логарифма – медленная операция, увеличивающая общее время выполнения на 30-50%. Поэтому от 24-битного цвета отказались ради скорости и выбрали ограниченную палитру цветов.
Чтобы перевести алгоритм в код C, надо помнить формулы для умножения комплексных чисел и нахождения их модулей.
void PaintMandelbrotFPU(int w, int h, long* b, double left, double top,
double right, double bottom) {
int x, y, i; double x2, y2, px, py, zx, zy, dx, dy;
dx = (right - left) / w; // Найти горизонтальное и вертикальное увеличение для каждого пикселя
dy = (bottom - top) / h;
py = top;
for(y = 0; y < h; ++y) {
px = left;
for(x = 0; x < w; ++x) {
zx = px, zy = py;
for(i = 0; i < ITER; ++i) {
x2 = zx * zx, y2 = zy * zy;
if(x2 + y2 > 4.0)
break;
zy = zx * zy * 2 + py; // Вычислить z = z * z + p
zx = x2 - y2 + px;
}
*b++ = a[i];
px += dx; // Перейти к следующему пикселю
}
py += dy; // Перейти к следующей строке пикселей
}
}
Комплексная переменная p превратилась в px (действительная часть) и py (мнимая часть). То же самое происходит для zx и zy. Квадраты zx и zy используются дважды: при вычислении модуля и при нахождении z * z + p. Поэтому для них было сделано две переменных x2 и y2. Переменные zx и zy инициализируются значениями px и py вместо нулей, потому что пропуск умножения 0*0 экономит время.
Обратите внимание на целые переменные x и y, считающие координаты текущего пикселя, и на пару (px, py), являющуюся координатами на комплексной плоскости. px и py увеличиваются независимо от x и y. Это кардинально повышает производительность, потому что не надо переводить экранные координаты (x, y) в координаты комплексной плоскости (px, py) при каждом повторении цикла. Плата за такое увеличение производительности – потеря точности. Числа с плавающей точкой имеют ограниченную точность, и ошибки накапливаются в ходе многократных сложений: px += dx, py += dy. Это вызывает неприятный эффект: при увеличении масштаба изображения оно будет дергаться. Альтернативным решением могло быть вычисление px = (right - left) * x / w, py = (bottom - top) * y / h на каждой итерации, но оно непрактично из-за большой задержки на деление.
Наконец, a – палитра цветов, а b – битовая карта, генерируемая функцией. Битовая карта b позже переносится на экран с помощью функций DirectDraw или GDI(интерфейс графического устройства). Детали этого процесса неважны. Просто помните, что размер b равен w * h пикселей, и что надо заполнить его цветами из палитры.
Когда внутренний цикл прерывается с помощью break, значение i меньше ITER. А когда он не прерывается, i == ITER после цикла. Благодаря данному свойству был создан массив a[ITER+1], и последнему элементу этого массива был присвоен черный цвет. После цикла a[i] ссылается на черный a[ITER]. Более строгий способ сделать это был бы:
// ...
for(i = 0; i < ITER; ++i) {
x2 = zx * zx, y2 = zy * zy;
if(x2 + y2 > 4.0) {
*b++ = a[i];
px += dx;
goto NEXTPIXEL;
}
zy = zx * zy * 2 + py; // Вычислить z = z * z + p
zx = x2 - y2 + px;
}
*b++ = RGB(0,0,0);
px += dx; // Перейти к следующему пикселю
NEXTPIXEL:
// ...
Алгоритм для генерации набора Джулии схож с алгоритмом для набора Мандельброта, поэтому не показан тут. Главное отличие – введение комплексного параметра (PX, PY) в Джулии.
Метод MOV-ADD
Профилирование показало, что функция PaintMandelbrotFPU занимает около 95%-99% полного времени выполнения. Рисование занимает только 5% с GDI или 1% при использовании DirectDraw. Оптимизировать данную функцию гораздо важнее, чем любую иную часть программы.
Один из лучших способов оптимизации - использование SSE (потоковые SIMD-расширения). Эти команды позволяют вычислять значения 4 пикселей одновременно. Следовательно, хотя бы в теории, скорость возрастет в 4 раза (реальные цифры - около 3.2-3.3). Дополнительную информацию о SSE можно найти в Руководстве разработчика ПО под архитектуру Intel, том 1, глава 10, и в Руководстве программиста под архитектуру AMD64, том 1, глава 4. Дальнейшие объяснения предполагают, что вы знаете основы программирования на ассемблере под SSE и x86.
Расположить команды SSE для этой функции можно разными способами, и некоторые из них быстрее остальных. Ниже приведена очень быстрая последовательность для процессора Pentium IV, называемая методом MOV-ADD:
; xmm0 = zx xmm1 = zy xmm2 = tmp xmm3 = tmp
; xmm4 = px xmm5 = py xmm6 = результат xmm7 = 4.0
; eax = tmp ecx = счётчик цикла
MOVAPS xmm0,xmm4
XORPS xmm6,xmm6
MOVAPS xmm1,xmm5
mov ecx, ITER
ILOOP:
; xmm0 = zx xmm1 = zy
MOVAPS xmm2, xmm0
MULPS xmm0, xmm0
MOVAPS xmm3, xmm1
ADDPS xmm1, xmm1
; xmm0 = zx^2 xmm1 = 2 * zy xmm2 = zx xmm3 = zy
MULPS xmm1, xmm2
MOVAPS xmm2, xmm0
MULPS xmm3, xmm3
; xmm0 = zx^2 xmm1 = 2*zy*zx xmm2 = zx^2 xmm3 = zy^2
ADDPS xmm1, xmm5
SUBPS xmm0, xmm3
ADDPS xmm2, xmm3
; xmm0 = zx^2 - zy^2 xmm1=2*zy*zx+py xmm2 = zx^2 + zy^2 xmm3 = zy^2
CMPLEPS xmm2, xmm7
ADDPS xmm0, xmm4
MOVMSKPS eax, xmm2
test eax, eax
jz EXIT
ANDPS xmm2, xmm7 ; xmm6 += (xmm2 < 4.0) ? 4.0 : 0.0;
ADDPS xmm6, xmm2
sub ecx, 1
jnz ILOOP
EXIT:
Здесь показан только внутренний цикл. К началу цикла xmm0 и xmm1 содержат px и py, начальные значения для zx и zy. Регистры xmm4 и xmm5 хранят px и py (две постоянные цикла). Регистр xmm7 содержит константу 4.0 во всех двойных словах.
Программа начинает работу с сохранения xmm0 и xmm1 в буферные регистры xmm2 и xmm3. Затем zx возводится в квадрат, и zy удваивается путем сложения xmm1 с самим собой. Остальные шаги расчета показаны в комментариях выше.
Регистр xmm6 хранит число полных итераций. Заметьте, что число в xmm6 должно считаться независимо для каждого из 4 пикселей. Команда CMPLEPS помещает битовую маску в регистр xmm2, и последовательность MovMskPS-test-jz проверяет, вышли ли все 4 пикселя за пределы круга радиусом 2.0. Если это происходит, то цикл прерывается. В иных случаях программа должна увеличить число полных итераций для пикселей, все еще пребывающих внутри круга. Затем производится операция "логическое и"над маской и 4.0, чтобы получить 4.0 в двойных словах, соответствующих таким пикселям, и 0.0 в других двойных словах. Это значение прибавляется к xmm6. Короче говоря, когда цикл заканчивается, xmm6 будет содержать 4.0 * число полных итераций для каждого из 4 пикселей. Значение позже будет использоваться как сдвиг в палитре.
Теперь вычисляется время выполнения для Pentium IV и находятся узкие места в этом коде. Руководство от Anger Fog описывает, как провести такой анализ; если вы впервые знакомитесь с оптимизацией для Pentium IV, прочитайте руководство. Также смотрите официальные документы от Intel: Справочник по оптимизации для архитектуры Intel и Микроархитектура процессора Pentium 4.
Во-первых, обращение к памяти не является узким местом, потому что вычисление занимает гораздо больше времени, чем запись в память. Для обычного разрешения 1024x768 программа записывает 1024*768*4 байта = 3 Мб. Типичная производительность современной памяти DDR составляет 1-2 Гб/с. Грубая оценка показывает, что обращение к памяти займет 1.5-3 миллисекунды, тогда как реальное время выполнения этого кода - 70-80 мс. Анализ кода дает такой же результат: происходит всего 2 обращения к памяти для 1 пикселя (смотрите код в следующем разделе), но внутренний цикл выполняется ITER раз (обычно около 64 раз), и он включает в себя команды SSE с долгими задержками.
Поэтому разные приемы оптимизации производительности кэша бесполезны для этого кода. То же самое верно для расшифровки команд: циклы столь малы, что они полностью вмещаются в трассировочный кэш и расшифровываются лишь один раз. Реальное узкое место – долгие задержки и ограниченная производительность команд SSE.
Для простоты предположим, что конвейер вначале пустой, и весь код находится в трассировочном кэше. Имена регистров xmm1, xmm2, и т.д. являются сокращениями для x1, x2 чтобы они вошли на рисунок.
Первые три команды (MOVAPS x2, x0, MULPS x0, x0, и MOVAPS x3, x1) отправляются из трассировочного кэша в очередь операций. Кластер выполнения FP/SSE имеет два порта: один - для операций add, mul, div, названный блоком fp в руководстве Agner Fog, а второй порт - для пересылок между регистрами и для сохранений (блок выполнения mov). Эти два блока могут начать выполнять команды MOVAPS и MULPS одновременно.
Во втором тактовом цикле следующие три команды уходят в очередь. Из них MOVAPS x3, x1 и ADDPS x1, x1 могут выполняться параллельно. Взаимная производительность для команды MOVAPS равна 1, следовательно, в этом цикле может запуститься вторая команда.
Задержка команды MOVAPS составляет 6 тактовых циклов, и MULPS x1, x2 зависит от результатов MOVAPS x2, x0 (зависимость показана стрелками на схеме). Следовательно, MULPS может запуститься в шестом цикле (к этому времени очередь должна быть полна операций, ожидающих выполнения).
Следующая команда MOVAPS x2, x0 зависит от MULPS x0, x0. Команда MULPS также имеет дополнительную задержку 1, а значит, ей надо ждать 1 лишний тактовый цикл, если следующая операция уходит в другой блок выполнения. MULPS и MOVAPS родом из разных блоков выполнения fp и mov, и между ними пришлось вставить задержку в один цикл.
Теперь команда MULPS x3, x3. Подблок умножения (умножитель FP в понятиях Intel) может производить одно умножение каждые два цикла, поэтому второй команде MULPS придется ждать один цикл. Она запустится в восьмом цикле.
Следующие пять команд уходят в подблок суммирования FP (сумматор FP). Все они имеют задержку в 4 тактовых цикла. ADDPS x1, x5 может запуститься на 12-м цикле, когда будут готовы результаты умножения MUL x1, x2. То же самое касается SUBPS x0, x3; она запустится в 14-м цикле. Команда ADDPS x2, x3 может запуститься раньше, чем на 16-м цикле, но производительность сумматора FP ограничивает число одновременно выполняющихся команд до 2. Последовательность SUBPS x0, x3; ADDPS x2, x3; ADDPS x0, x4; CMPPS x2, x7 идеально соответствует этому требованию сумматора FP. Две соседних команды взяты из разных цепей зависимости, следовательно, они полностью независимы и могут выполняться параллельно.
Задача – завершить вычисление x0 = x0 ^ 2 - x1 ^ 2 + x4, x1 = 2 * x0 * x1 + x5 как можно быстрее. На графике команды, вычисляющие новые значения для x0 и x1, отмечены разными оттенками серого, а неважные команды белые. Надо включить в график серые команды раньше остальных, потому что следующая итерация цикла не может начаться без новых значений для x0 и x1. Команды ANDPS и MOVMSKPS из предыдущей итерации могут выполняться параллельно со следующей итерацией, при условии, что переход будет предсказан правильно. Стало быть, их можно включить в график позже серых команд.
Проведенные измерения времени подтвердили эти соображения. Сейчас есть долгая задержка между MOVAPS x2, x0 и ADDPS x2, x3 (смотрите схему). Но если поместить ADDPS x2, x3 раньше SUBPS x0, x3, то полное время выполнения увеличится. Значит, задержка не важна, при условии что x0 и x1 вычисляются рано.
Быстродействие кода зависит от задержек команд и производительности, потому что есть три независимых цепи, способных выполняться одновременно. Критический путь отмечен темно-серым цветом. MOVAPS добавляет лишнюю задержку (может быть, если бы архитектура x86 имела трехадресные команды, как RISC(вычисления с сокращённым набором команд), выполнялся бы быстрее, как знать). Подблоки не используются полностью. На первой стадии цикла (циклы 0-12) большая нагрузка ложится на блок MOV и на умножитель FP; позже, в циклах 14-22, перегружается сумматор FP, а умножитель FP вообще не используется. Можно взять команды из двух итераций и перемешать их, чтобы получить более равномерную нагрузку, но для этого надо больше регистров SSE. В AMD64 (или Intel EM64T) имеется 8 дополнительных регистров, следовательно, код можно дополнительно оптимизировать в 64-битном режиме.
После рисования схемы команды были упорядочены по их времени запуска. Например, команды MOVAPS x2, x0 и MULPS x0, x0 запускаются раньше всех, поэтому поставлены в коде первыми. Затем запускаются MOVAPS x3, x1 и ADDPS x1, x1, поэтому они являются следующими командами в исходном коде. Упорядочивание команд сокращает число операций, ждущих в очереди. Теперь процессору не надо переупорядочивать команды при выполнении: они уже упорядочены согласно их задержкам, производительности и зависимостям.
Но как насчет предсказания ветвлений в этом коде? Может ли оно быть узким местом? По сути, потери на непредсказанные ветвления отнимают много времени выполнения для этой программы. Но число полных итераций является результатом, который надо получить. Нельзя точно предсказать, сколько раз выполнится цикл, потому что само это значение является результатом программы. Если точно известно, сколько раз цикл повторяется для каждого пикселя, зачем тогда запускать эту программу?
Были случайно выбраны два увеличенных изображения наборов, затем было вычислено распределение числа итераций для них и для одного большого изображения:
Голубые полосы обозначают большое изображение набора Мандельброта; красные и желтые полосы показывают два увеличенных изображения. ITER в программе равно 64, поэтому диаграмма показывает частоты точек от 0 до 5 итераций, от 6 до 11 итераций, и так далее до 60 итераций и более. Минимальное и максимальное числа итераций составляют крупнейшие части изображения. Для голубой полосы минимальное значение равно 0-5 итераций (59% всех пикселей); для красной полосы - 12-17 (43%); и для желтой полосы - 30-35 (8%) и 36-41 (15%). Минимальное число итераций соответствует огромной фоновой области, закрашенной одним цветом, а максимальное число итераций означает черные точки, принадлежащие набору. Таким образом, 26% пикселей принадлежат набору Мандельброта на большом изображении, а 21% и 62% точек принадлежат набору на двух увеличенных изображениях.
Из этого анализа ясно, что если цикл прерывается, то велика вероятность, что он будет прерван рано, на минимальном числе итераций. Но если цикл не прерывается рано, обычно он вообще не прерывается. Эта статистика интересна, но ей не нашлось применения в коде. Может быть, читатель сможет сделать это.
Ниже сделано несколько выводов о методе MOV-ADD:
• Общая стратегия – как можно раньше вычислить значения, требуемые в следующей итерации (в данном случае, x0 и x1). Эти вычисления составляют критический путь программы, и полное время выполнения зависит от их времени завершения.
• Планирование одновременного выполнения независимых команд улучшает параллелизм на уровне команд (ILP). Например, первые четыре команды выполняют сложение, умножение и две пересылки параллельно.
• Команды упорядочены по их времени запуска, что уменьшает длину очереди операций.
Метод FFFF
FFFF v3.2.2, вышеназванный генератор набора Мандельброта, использует другой метод. FFFF имеет открытый исходный код, поэтому публикация его кода здесь полностью правомерна:
; xmm0 = zx xmm1 = zy xmm2 = tmp xmm3 = tmp xmm4 = px
; xmm5 = py xmm6 = результат xmm7 = 4.0
; eax = tmp ecx = счётчик цикла
MOVAPS xmm0,xmm4
XORPS xmm6,xmm6
MOVAPS xmm1,xmm5
mov ecx, ITER
ILOOP:
; xmm0=zx xmm1=zy
MOVAPS xmm2, xmm0
MULPS xmm2, xmm1
MULPS xmm0, xmm0
MULPS xmm1, xmm1
; xmm0=zx^2 xmm1=zy^2 xmm2=zx*zy
addps xmm2, xmm2
movaps xmm3, xmm0
addps xmm3, xmm1
; xmm0=zx^2 xmm1=zy^2 xmm2=2*zx*zy xmm3=zx^2+zy^2
cmpltps xmm3, xmm7
movmskps eax, xmm3
test eax, eax
jz EXIT
subps xmm0, xmm1
movaps xmm1, xmm2
; xmm0=zx^2-zy^2 xmm1=2*zx*zy xmm3=zx^2+zy^2
ANDPS xmm3, xmm7 ; xmm6 += (xmm3<радиус) ? 4.0 : 0.0
ADDPS xmm6, xmm3
ADDPS xmm0, xmm4
ADDPS xmm1, xmm5
dec ecx
jnz ILOOP
EXIT:
FFFF пытается сравнить x0^2 + x1^2 с 4.0, и только если сравнение не удается, находит новые значения для x0 и x1. Если сравнение возвращает истину, цикл прерывается, и нет нужды вычислять x0 и x1. Хотя такая стратегия может казаться хорошей, эксперименты показывают, что это не так. Посмотрите на длинную цепь: MOVAPS xmm3, xmm0; ADDPS xmm3, xmm1; CMPLTPS xmm3, xmm7; MOVMSKPS eax, xmm3. При выполнении этой цепи есть много возможностей использовать свободные блоки выполнения и параллельно выполнить вторую цепь. Но авторы FFFF не вставили никаких команд между командами этой цепи. Планировщик вынужден сам переупорядочивать команды. Иногда он терпит неудачу, и блоки выполнения остаются неиспользованными. Например, в коде между MULPS xmm1, xmm1 и SUBPS xmm0, xmm1 столь большое расстояние, что вторая команда может не попасть в очереди вовремя. В таком случае между командами добавляется дополнительная задержка. То же самое касается ADDPS xmm2, xmm2 и MOVAPS xmm1, xmm2: маловероятно, что они выполнятся без задержки.
Схема показывает идеальный случай вообще без задержек. Даже в таком случае на вычисление x0 требуется 20 тактовых циклов, а на вычисление x1 требуется 27 тактовых циклов. Для метода MOV-ADD цифры - 16 и 22 тактовых циклов, соответственно. Критический путь (отмеченный темно-серым цветом) включает в себя две команды MOVAPS, имеющие долгие задержки, тогда как метод MOV-ADD использует только одну MOVAPS на своем критическом пути.
Есть дополнительные задержки перед MOVAPS x1, x2 и MOVAPS x3, x0, потому что предыдущая команда выполнилась в другом блоке. ADDPS x3, x1 откладывается на один цикл из-за ограниченной производительности сумматора FP.
Во всем остальном MOV-ADD схож с FFFF. Оба метода используют SSE для вычисления значений 4 пикселей вместе, оба имеют одинаковое ограничение точности и одинаковый эффект дергания при масштабировании. Был использован схожий алгоритм, но команды в нем были расположены иначе для повышения производительности. Другие различия между FFFF и приведенной программой будут показаны позже.
В коде FFFF выше имена регистров изменены на имена, используемые в MOV-ADD (xmm0 для zx, xmm1 для zy, и так далее). Это изменение не влияет на производительность, но облегчает сравнение кода методов.
Метод Vodnaya
Также был разработан метод Vodnaya, оптимизированный для Pentium M. Название лишено смысла, поэтому если вам трудно произносить русские слова, можете называть его Vod-метод или V-метод. Он превосходит по производительности методы MOV-ADD и FFFF на Pentium M и Pentium III. О микроархитектуре Pentium M известно мало, поэтому не удалось проанализировать задержки и производительность. После многих проб и ошибок была найдена последовательность команд, казавшаяся самой быстрой на этих процессорах. Метод также является лидером по числу команд: он на одну команду меньше, чем MOV-ADD или FFFF.
MOVAPS xmm0,xmm4
XORPS xmm6,xmm6
MOVAPS xmm1,xmm5
mov ecx, ITER
ILOOP:
; xmm0 = zx xmm1 = zy
MOVAPS xmm2, xmm1
MULPS xmm2, xmm2
MULPS xmm1, xmm0
; xmm0 = zx xmm1 = zy*zx xmm2 = zy^2
MOVAPS xmm3, xmm2 ; save zy^2
MULPS xmm0, xmm0
ADDPS xmm2, xmm0
; xmm0 = zx^2 xmm1 = zy*zx xmm2 = zx^2+zy^2 xmm3 = zy^2
CMPLEPS xmm2, xmm7
ADDPS xmm1, xmm1
SUBPS xmm0, xmm3
; xmm0 = zx^2-zy^2 xmm1 = zy*zx*2 xmm2 = zx^2+zy^2
ADDPS xmm1, xmm5
ADDPS xmm0, xmm4
MOVMSKPS eax, xmm2
test eax, eax
jz EXIT
ANDPS xmm2, xmm7
ADDPS xmm6, xmm2
sub ecx, 1
jnz ILOOP
EXIT:
Сумма zx^2+zy^2 вычисляется в буферном регистре xmm2. zy^2 копируется сюда во избежание его перезаписи, когда программа вычтет zx^2-zy^2. Новые значения для zx и zy вычисляются на месте в регистрах xmm0 и xmm1. Весь метод довольно изящен. Увы, он не самый быстрый для процессоров Pentium-IV и Athlon XP.
Великая битва: MOV-ADD против FFFF против Vodnaya
Сравнение быстродействия трех вышеназванных методов:
В наборе тестов размер изображения был 1024x768 пикселей; время измерялось с помощью команды RDTSC (смотрите исходный код в файле SSEroutines_test.asm). Тесты были повторены 5 раз, и был взят минимальный результат для каждого метода. Времена на диаграмме выражены в миллионах тактовых циклов. Более низкие времена лучше. После названий процессоров показаны их сигнатуры CPUID. Например, Celeron F13 – процессор семейства 15, модели 1, степпинга 3, то есть сокращенный Pentium IV. Его кэши меньше, чем у Pentium IV, но ядро такое же, оттого результаты этого Celeron должны почти равняться результатам Pentium IV (увы, не удалось достать Pentium IV для тестов). "Rabbit(кролик)" означает набор Джулии с C = -0.12+0.74i.
Видимо, MOV-ADD лучший метод для процессоров Pentium IV и Athlon XP. Vodnaya показывает хорошие результаты на процессорах Intel семейства 6 (в частности, Pentium III и Pentium M). FFFF однозначно проигрывает на всех типах процессоров. Повышение быстродействия MOV-ADD по сравнению с FFFF составляет 12% на Pentium IV, и разница между Vodnaya и FFFF составляет 7-11% для Pentium M и Pentium III.
Имеющаяся готовая программа использует Vodnaya на процессорах Intel семейства 6 и метод MOV-ADD на остальных процессорах. Метод выбирается при выполнении с помощью сигнатуры процессора, извлеченной с помощью команды CPUID.
CvtTss2si
После окончания расчета программа должна перевести результаты в значения цветов. Одна строка кода (*b++ = a[i];) может произвести перевод в C, но на языке ассемблера придется сделать больше работы.
Программа устанавливает xmm6 в ноль и прибавляет к нему 4.0 на каждой итерации (смотрите код любого метода выше). Если внутренний цикл не прерывается, он повторяется ITER раз. При завершении цикла число полных итераций умножается на 4 в регистре xmm6 при условии, что точка не принадлежит набору Мандельброта, и там есть значение ITER * 4 в противном случае. Можно провернуть такой же трюк с последним элементом массива: добавить один элемент к массиву и приравнять его значение к черному. С помощью такого приема можно использовать один код окраски для обработки точек, входящих в набор Мандельброта, и точек, не входящих в него. FFFF v3.2.2 по-разному обрабатывает эти типы точек, используя условную пересылку (это более медленный путь, и команду CMOV не поддерживают некоторые старые процессоры AMD).
Заметьте, что регистр xmm6 хранит число полных итераций, умноженное на 4. FFFF делит его, чтобы получить индекс в палитре. Для деления он использует команду shr, медленно работающую на процессорах Pentium IV с номером модели меньше 3. Но один элемент палитры занимает 4 байта в памяти. Таким образом, FFFF делит число на 4, а затем умножает его на 4, чтобы получить сдвиг в памяти. Какая ужасная потеря времени! Это делается гораздо проще:
CVTTSS2SI eax, xmm6
; xmm6 содержит сдвиг в палитре, т.е. 4*число итераций
mov eax, [esi + eax] ; палитра начинается в [esi]
mov [edi], eax ; текущий пиксель битовой карты находится в [edi]
Команда CvtTss2si преобразует, с усечением, значение с плавающей точкой с одинарной точностью в целое значение в регистре eax. Количество итераций, разумеется, целое число. Оно преобразуется без округления, поэтому режим округления тут не важен. Следующая команда извлекает цвет из палитры, а последняя команда записывает цвет в битовую карту. Битовая карта может находиться в видеопамяти или в ОЗУ (пока не беспокойтесь об этом).
Есть другие команды для преобразования чисел с плавающей точкой (а именно, CvtPs2Pi и CvtPs2Dq), но ни одна из них не записывает результат в регистры общего назначения. Приходится сохранять число в памяти и сразу читать его из того же самого места, что вызывает дополнительные задержки в Pentium IV. Тестирование кода с CvtPs2Dq показало, что он не быстрее кода, использующего CvtTss2si.
Мощь макроса FASM
Написание нескольких вариантов одного и того же кода для набора Мандельброта и для набора Джулии, для чёрно-белого изображения и для цветного изображения, для SSE и для SSE2, требует много рутинной работы. К счастью, большинство ассемблеров имеет специальные инструменты, способные сделать эту работу за вас. Инструменты - макрокоманды и условная компоновка. Макросы позволяют повторять последовательность команд с разными параметрами, а условная компоновка позволяет вставлять команды в зависимости от этих параметров.
Был выбран FASM 1.62 из нескольких ассемблеров, потому что он маленький, быстрый и имеет богатый синтаксис макросов. Руководство по FASM подробно описывает синтаксис. Ниже приведен маленький кусок из метода Vodnaya для иллюстрации использования макросов в программе:
macro JuliaMandelPaint color, type {
; ...
ADDPS xmm1, xmm1
SUBPS xmm0, xmm3
if type = Julia
ADDPS xmm1, dqword[cy1]
ADDPS xmm0, dqword[cx1]
else if type = Mandel
ADDPS xmm1, xmm5
ADDPS xmm0, xmm4
end if
; ...
}
Если аргумент type равен Julia, FASM сгенерирует код, добавляющий cx1 к xmm0 и cy1 к xmm1, иначе он будет использовать значения из регистров xmm4-xmm5. Оба варианта кода будут сгенерированы: первый с меткой JULIA и второй с MANDELBROT. Вы лишь проверяете переменную fractaltype и решаете, какой вариант использовать.
Можно сказать, что есть более легкий способ сделать это:
mov eax, [fractaltype]
test eax, eax
jz MANDELBROT
ADDPS xmm1, dqword[cy1]
ADDPS xmm0, dqword[cx1]
jmp NEXT
MANDELBROT:
ADDPS xmm1, xmm5
ADDPS xmm0, xmm4
NEXT:
Почему бы не использовать этот код? Потому что проверка находится в глубоко вложенных циклах, и сравнение повторится несколько миллионов раз, тогда как его надо делать всего один раз. Нечего и говорить о задержках, возникающих при не предсказании условного перехода. Поэтому выгодно вынести сравнение за пределы цикла. Макросы помогают кратко записать это.
Следующая интересная задача – адаптация методов под SSE2. Команды SSE2 параллельно работают над двумя значениями двойной точности. В иных аспектах они не отличаются от SSE. Даже задержки и производительности для SSE2 те же, что и для SSE.
Надо провести массовый поиск с заменой, чтобы изменить ADDPS на ADDPD, MULPS на MULPD, MOVAPS на MOVAPD, и т.д. Затем надо поправить увеличения указателей и получить две функции: одну – для SSE и одну – для SSE2. Такой подход применяется в FFFF и нескольких других программах. У него есть один серьезный недостаток. При обновлении программы приходится менять один и тот же код в обеих функциях. Это нудно и неприятно. Иногда исправляют одну функцию, но забывают изменить другую.
Идея такова – создать макрос по имени MOVAP, заменяемый на MOVAPS для SSE и на MOVAPD для SSE2. Теперь вместо двух функций пишется только одна, и в ней используются макросы MOVAP, MULP и ADDP. Во время компиляции реальные команды SSE или SSE2 используются вместо этих макросов:
macro MOVAP a, b {
if ssewidth = 4
MOVAPS a, b
else
MOVAPD a, b
end if
}
Не спешите писать этот код! Есть более элегантное решение:
macro defineinstr instr, [args] {
common
if ssewidth = 4
instr#S args
else
instr#D args
end if
}
Макрос defineinstr принимает переменное число аргументов. Например, он может использоваться для SHUFPS с 3 аргументами и для ADDPS с 2 аргументами. Макрос выдает команду с суффиксом S, если ssewidth = 4, или ту же самую команду с суффиксом D в противном случае. Оператор # соединяет имена; common означает, что команда повторится всего один раз.
Можно еще сократить этот код:
macro defineinstr instr {
macro instr [args] \{
\common
if ssewidth = 4
instr#S args
else
instr#D args
end if
\}
Тут есть два вложенных макроса. Внешний макрос defineinstr определяет внутренний макрос instr. Заметьте, что instr также является аргументом defineinstr, поэтому при написании defineinstr MOVAP будет определен макрос по имени MOVAP. Это как раз то, что надо. Фигурные скобки и директиву common надо экранировать обратным слешем из-за каких-то странных проблем парсинга внутри FASM. Теперь определение новой команды SSE/SSE2 очень краткое и простое: пишется defineinstr и имя команды. Число аргументов не имеет значения.
Теперь вы убедились, что макросы FASM очень мощные и изящные. Если нет, попробуйте сделать макрос препроцессора C с переменным числом аргументов, и сравните его с эквивалентом FASM!
Копирование битовых карт с помощью GDI и DirectDraw
Функция SetPixel кажется самым медленным из возможных способов вывода пикселей на экран:
case WM_PAINT: {
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
for(x = 0; x < w; ++x)
for(y = 0; y < h; ++y)
SetPixel(hdc, x,y, b[i]);
EndPaint(hWnd, &ps);
}
Используйте эту функцию только для установки 1 или 2 пикселей; она слишком медленная для больших изображений. Ниже приведена более быстрая упрощенная версия программы:
int w, h; HBITMAP hBMP = 0, holdBMP; HDC hmem;
BITMAPINFOHEADER bi = {sizeof(BITMAPINFOHEADER), 0, 0, 1, 32, BI_RGB};
void ReDraw() {
PaintJuliaSSE(b, ceil(w, 4), h);
SetDIBits(hmem, hBMP, 0, h, b,
(BITMAPINFO*)&bi, DIB_RGB_COLORS);
// Копировать битовую карту в hmem
}
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wParam,LPARAM lParam) {
switch(msg) {
case WM_SIZE: {
HDC hdc = GetDC(hWnd);
if(hBMP) { DeleteObject(SelectObject(hmem, holdBMP)); DeleteDC(hmem); }
w = LOWORD(lParam), bi.biWidth = ceil(w, 4),
h = HIWORD(lParam), bi.biHeight = h;
bi.biSizeImage = bi.biWidth * h * 4;
b = (long*)realloc(b, bi.biSizeImage);
hBMP = CreateCompatibleBitmap(hdc, ceil(w, 4), h);
hmem = CreateCompatibleDC(hdc);
holdBMP = (HBITMAP)SelectObject(hmem, hBMP);
ReleaseDC(hWnd, hdc);
ReDraw();
return 0; }
case WM_PAINT: {
PAINTSTRUCT ps; HDC hdc = BeginPaint(hWnd, &ps);
BitBlt(hdc, 0, 0, w, h, hmem, 0, 0, SRCCOPY); // Копировать hmem на экран
DrawPolyline(hdc);
EndPaint(hWnd, &ps);
return 0; }
case WM_DESTROY:
DeleteObject(SelectObject(hmem, holdBMP)); DeleteDC(hmem);
PostQuitMessage(0); return 0;
default:
return DefWindowProc(hWnd,msg,wParam,lParam);
}
}
Функция ReDraw вызывается при каждом полном обновлении изображения, т.е. при его увеличении или уменьшении. Она использует PaintJuliaSSE для генерации битовой карты фрактала в b. Затем вызывается функция SetDIBits. Она копирует весь массив b в контекст устройства памяти hmem.
В ответ на WM_SIZE программа устанавливает новый размер битовой карты и заново создает ее. Также память для битовой карты перераспределяется, и фрактал полностью перерисовывается. Обработчик для WM_PAINT просто копирует контекст устройства памяти на экран и добавляет ломаную линию, показывающую орбиту функции.
Изображение фрактала хранится в контексте устройства памяти, существующем на протяжении всего времени работы программы. Имея это изображение, можно быстро нарисовать движущуюся ломаную путем копирования сохраненного изображения и рисования нового положения вышеназванной ломаной. Это вызовет небольшое мерцание ломаной, что допустимо. Добавление второго контекста устройства памяти и второй битовой карты может полностью убрать мерцание, но тогда программа будет тормозить и потреблять много памяти.
Команды SSE работают над 4 значениями параллельно, и все вышеприведенные методы SSE требуют, чтобы число пикселей в линии было кратно 4. Поэтому ширина (w) округляется в большую сторону функцией ceil(w). Битовая карта становится немного больше чем надо, и в обработчике WM_PAINT из нее копируется только прямоугольник w * h.
Кроме GDI, программа также поддерживает DirectDraw. Код здесь не показан, потому что он очень длинный и скучный (обработка ошибок занимает много места). По сути, он сходен с кодом GDI. Создается скрытая поверхность; фрактал отрисовывается на нее. В WM_PAINT скрытая поверхность копируется на главную поверхность (на экран) и рисуется ломаная. Современные видеоплаты хранят обе поверхности в видеопамяти, следовательно, копирование с помощью DirectDraw может быть быстрее в 50 раз, чем копирование через GDI (данный результат был получен на ноутбуке с простой встроенной видеоплатой Intel 915 GM). DirectDraw повышает общее быстроту действия на 5-6% по сравнению с GDI.
DirectDraw поддерживается для Windows 98 и выше. Для поклонников Windows 95 было написано несколько строк кода, поэтому программа изящно перейдет на GDI, если DirectDraw отсутствует.
Заключение
Созданная программа является очень быстрым генератором набора Мандельброта и Джулии, использующим SSE/SSE2 и DirectDraw. Она не называется "быстрейшим генератором фракталов для Win32", как сделал автор FFFF, но программа очень быстрая и удобная. Некоторые пункты все еще требуют доработки, например:
• есть более быстрые алгоритмы, "угадывающие" значения соседних пикселей;
• можно дополнительно оптимизировать методы MOV-ADD и Vodnaya, чтобы получить максимальную производительность на 64-битных процессорах;
• может быть, существует какой-то хитрый подход к оптимизации предсказания ветвлений, но его не нашли;
• программа не может сохранить текущее положение и параметр C, чтобы позволить вернуться к ним впредь (следующая версия будет делать это);
• изменение палитры не поддерживается, но такой цели и не было.
Итак, код хороший, но не идеальный.
Можно использовать представленные в данной статье методы во множестве разных графических приложений. Анализ задержек и производительности помогает найти узкие места. Затем можно убрать операции с длинными задержками с критического пути, сделать несколько независимых цепей в коде и увеличить параллелизм на уровне команд путем переупорядочивания команд. Наконец, макросы ассемблера освобождают от такой скучной работы, как переписывание того же самого кода для SSE2.