Логарифм на SSE по методу Буля – How to?

Изменено: . Разместил: FadeDemon | Меню
Количество просмотров3`121 / FD1766-6-0

История эта взяла своё начало в том самом Множестве Мандельброта, которому удалось чрезвычайно глубоко затронуть мою программерскую мысль. Дело началось с того, что до меня дошла информация о существовании алгоритма, позволяющего рисовать фрактал Мандельброта так, чтобы карта получаемых значений не имела резких переходов, а была исключительно гладкой, чего на самом деле очень-очень не хватало в моей программе. Алгоритм этот, кстати, даже название своё нехитрое имеет: Smooth coloring algorithm. Фракталы, сгенерированные с его помощью иногда просто поражают своим изяществом и красотой, достигнуть такого-же уровня эффектов классическим алгоритмом если и можно, то очень, и очень непросто. И всё бы кажется с новым нашим алгоритмом замечательно, вот только…


\[ \nu(z) = n - \log_P {\log|z_n| \over \log(N)}, \; z_n \in \mathbb C, P = 2, N = 32 \]
Smooth coloring. Суть

Вот только для реализации наших сглаженных цветов, нам всего-то нужно сделать такую мелочь, как научиться считать логарифмы. Казалось бы, а что тут сложного-то? Но на самом деле если начать вдумываться… Как мы помним, код построения фрактала целиком и полностью написан на SSE, и почти полностью исключает любые обращения к памяти, а SSE, к огромному нашему сожалению, логарифмы считать не умеет — не учили его таким вещам. Логарифмы в нашей системе по-умолчанию умеет считать только FPU, с помощью непонятной на вид инструкции FYL2X, которая делает ни что иное как берет из стека X, ищет его логарифм по основанию 2, и умножает его на Y, который тоже берет из стека. Как можно заметить, сделано это довольно удобно, поскольку таким способом можно легко посчитать, например, натуральный логарифм, подсунув вместо Y константу log(2), но только вот нам этот метод не подходит. И дело даже не в его поразительной медлительности. Да, команда FYL2X работает так долго, что время ее работы сравнимо с ≈10-ю обращениями к памяти (кешируемой), и с ≈25-ю командами SSE (выявлено экспериментально), но это только пол-проблемы. Дело тут в том, что не существует быстрого метода перекинуть содержимое XMM-регистра в FPU стек, кроме как корячиться через память. Да, есть например команда MOVDQ2Q, позволяющая пересылать значения из XMM в MMX, а регистры последнего, как известно, расположены прямо в стеке FPU, но во-первых, FPU стек использует формат чисел long double (C++), или Extended (Delphi), который требует 10 байт на число, а размер MMX регистров всего 8 байт, во вторых не существует я не нашел способа безболезненно вернуться из режима MMX в режим FPU помимо команды EMMS, которая добросовестно помечает все регистры FPU стека пустыми, забывая таким образом любую информацию, которую до этого кто-либо пытался запихать в MMX-регистры. Посему, возможность использования старого доброго x87 чтобы считать наши замечательные логарифмы, пришлось отбросить. Ну, по крайней мере отложить в самый-самый дальний угол головы. Чтобы дать простор творческой мысли. А мысль творческая, разгулявшись, зашла тааак далеко…


Вспоминаем матанализ


В общем, поставленная задача быстро обрела новый облик, выйдя далеко за рамки джедайского ассемблерного программирования. Всё, чему нам нужно было теперь научиться, это научиться считать значение функции натурального логарифма в произвольной точке, используя только бумажку, ручку, и 4 стандартные арифметические операции, которые мы все хорошо знаем. И ладно бы просто считать, дак ещё и делать это быстро, желательно намного быстрее чем делает это FPU. Но вот кое-в чём у нас есть нехилый допуск: в точности. Поскольку карта значений фрактала рано или поздно всё-равно подвергнется растеризации, то-есть преобразованию в дискретные значения цветов, которые могут принимать только целые значения от 0 до 255, то нетрудно догадаться, что требуемая точность, при которой практически полностью будут отсутствовать погрешности, лежит в пределах значения константы 1 / 255 (что приблизительно составляет ≈0.003), и очень неплохим её приближением будет значение 10−3, то-есть, от нашего метода требуется точность, которая всегда обеспечит всего 3 точных цифры после запятой. Можно, конечно, больше, но как видим, совсем не обязательно. А этот факт в некоторых методах, на самом деле, является ключевым, и определяющим эффективность метода… Но не всё сразу.


\[ \log{x \over {x-1}} = \sum_{n=1}^\infty {1 \over {n x^n}} = {1 \over x}+ {1 \over {2x^2}} + {1 \over {3x^3}} + \dots \]
Преобразованный ряд Тейлора логарифма

Исследования мои начались с обыкновенного ряда Тейлора, в который довольно легко разложить логарифм. Но чёрт! Никогда, запомните, никогда в нашей жизни небывает всё так просто… Суть проблемы в том, что ряд расходиться при |x| > 1, сколько его в этих точках не считай. С помощью комплексного анализа можно даже объяснить, почему так происходит, но мы на сложности такие заморачиваться не будем. Мы лучше воспользуемся преобразованной формулой, тактично с3.14зженной из Википедии, которая приведена чуть выше, и демонстрирует разложение сией функции в тот-же самый ряд Тейлора. Но этот ряд уже не глючит, и покорно сходится на всём множестве вещественных чисел. Но Господи! КАК он сходиться… Допустим, значение в точке x = 2 с хотимой точностью 10−3 можно посчитать всего через 5-6 членов. В точке x = 4, количество членов возрастает уже до 20. Точка x = 15 же, требует подсчёта уже ни много ни мало 75 членов ряда, что само по себе совсем не спешит радовать. Но когда мы приближаемся к точке 100… Нет, дорогие мои. Вычисление 500 членов ряда не лезет ну совсем ни в какие ворота, даже очень большие. Эпик, эпик фэйл. Забыть как кошмарный сон.


Но ряд Тейлора, это далеко не единственное средство, которое может помочь нам сосчитать наш логарифм. В одном страшном сне, мне случайно привиделся способ посчитать его используя небезызвестный метод Ньютона, которым обычно решают извращенские трансцендентные уравнения вида F(x) = 0, в которых нельзя явно выразить x. Подобным образом, кстати, некоторые люди любят вычислять квадратный корень, да и кстати, не обязательно квадратный. Конкретно в нашем случае, мне пришла в голову мысль использовать вот такое вот нехитрое уравнение:


\[ x = \log(y) \Rightarrow e^x = y \Rightarrow e^x - y = 0 \]
Выражаем логарифм экспонентой

Ну гениально же, правда? Просто подставляем обе части первого уравнения в экспоненту, и получаем вполне себе симпатичное показательное уравнение, с вполне себе трансцендентным основанием e… Осталось только продифференцировать полученную функцию по x, и напомнить, что из себя представляет метод Ньютона, а кстати вся его суть кроется в такой вот рекуррентной формуле:


\[ x_{n+1} = x_n - {F(x_n) \over F'(x_n)} \]
Метод Ньютона

В эту формулу только мы подставляем x0 равный какой-нибудь белиберде, желательно так-же трансцендентной (исключает вероятность, что метод заглючит), например (π + e) × π × e + 123, тут можно пофантазировать, вычисляем x1, полученное значение опять ставим на место x0, снова что-то вычисляем, опять подставляем это что-то в формулу… И короче делаем вот так, пока не надоест не будет достигнута нужная точность, которую можно отследить, например, по значению модуля дроби в формуле. Таким интересным образом можно решить абсолютно всё. Кроме, правда, проблемы плохих дорог в России. Наше уравнение решается за милую душу. Причем за довольно скромное число итераций, от 3-х до 10, в зависимости от выбора начального значения. И знаете, я бы, наверное, использовал именно этот метод, если бы…


Если бы не экспонента. Да, чёрт возьми! Экспоненту, как и логарифм, так просто на SSE не посчитаешь! Правда тут поправочка, её можно очень даже быстро посчитать рядом Тейлора, который для неё выглядит весьма бесхитростно и непринуждённо:


\[ \mathrm{e}^{x} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots = \sum^{\infin}_{n=0} \frac{x^n}{n!}, \, x\in\mathbb{R} \]
Ряд Тейлора экспоненты

Но не так быстро, как может показаться. Экспоненте нашей, для приемлемой точности, хватало обычно 11 членов. Но проблема в том, что вычислять-то её надо для каждой итерации метода Ньютона! А итераций может быть 10. Сколько всего членов придётся считать? А ведь тут вам не просто циферки складывать, тут вообще-то полноценный полином надо считать… В общем, метод Ньютона, несмотря на свой шарм и очарование, тоже оказался слишком медлительным. С ним можно было ещё поэкспериментировать в плане выбора начального значения, но я не стал. Всё равно было бы слишком медленно. Очень медленно, по сравнению со следующим методом, до которого моя голова случайно додумалась в совсем другом, не менее страшном сне…


Интегрируем численно


Суть безумной мысли была в том, что логарифм, естественным образом можно было выразить через неопределённый интеграл от функции 1 / x, и через вполне себе определённый интеграл — как площадь под кривой графика этой функции начиная с точки x = 0, до точки, в которой нам надо что-то посчитать. Формально это можно записать, например так:


\[ \int\nolimits_1^x {1 \over t} dt = \log(x), \; x > 1 \]
Выражаем логарифм интегрально

Почему именно от единицы, догадаться в принципе не сложно, ведь именно в этой точке функция логарифма обращается в ноль, и начиная именно с неё начинает возвращать положительные значения. А отрицательные, нам, кстати, и не очень-то нужны, по крайней мере, для конкретной поставленной задачи — в ней вообще аргумент всегда будет больше двойки (и меньше некоторого числа), если его не поделить. Вообще, встраивая получившийся код в программу, я как раз таки сделал деление и еще более 100500 преобразований, позволивших избавиться от корня в модуле, завуалировать деление, и вообще свести все к 4-м операциям, но об этом в другой раз. Но логарифмы всё-таки классная штука. Суть тут в том, что наша задача свелась к тому, чтобы уметь считать логарифм всего-то для промежутка от 1 до 1024. 1024 тут, между делом, есть 322, который возникает потому, что именно 32 я решил выбрать в качестве значения N в формуле сглаженных цветов. Об этом подробнее в следующей статье.


Давайте теперь поразмыслим, во что мы превратили нашу задачу. Задача-то вроде как и не сильно-то упростилась, нам нужно научиться считать площадь под кривой функции 1 / x, что если поразмыслить, не так-то и просто. Ну вот какие мы знаем методы интегрирования? Метод трапеций, метод Симпсона, метод Гаусса… И всё? А вот ни фига не всё. Исследуя богатые возможности Wolfram Alpha, мне удалось наткнуться на весьма интересный метод численного интегрирования, одновременно очень точный, и простой как 3 копейки. И звался метод этот, методом Буля по 5-ти точкам или Boole's rule в англоязычной литературе. Давайте взглянем на его сущность:


\[ \begin{aligned} & \int\limits_{x_1}^{x_5} f(x) dx \approx {2h \over 45} \left (7f(x_1)+32f(x_2)+12f(x_3)+32f(x_4)+7f(x_5)\right ) \\ & x_k = x_1 + h(k - 1), \; h = {x_5 - x_1 \over 4 } \end{aligned} \]
Метод Буля, без остаточного члена

Здесь x2 … x4 — точки, стоящие между конечными точками x1 и x5 на равном друг от друга шаге h, который можно вычислить, как 1/4 разности x5 и x1, что в дальнейшем сыграло нам очень неплохую службу. Вся эта формула, кстати, описывает только один интервал интегрирования, а как вы понимаете, интервалов таких можно использовать сколько угодно. Но нам сколько угодно интервалов не пригодилось. Конечно, если в лоб считать площадь под кривой 1/x от 1 до, скажем, 1000, то там потребуется очень, очень много промежутков, чтобы получить приемлемую точность, минимум 15, но идею считать в лоб я быстро отбросил. Сейчас поясним, почему.


Идём мы значит в Wolfram Alpha, и набираем там следующие тексты:


Код: DOS
integrate 1/x from 1 to 2 using boole's rule
integrate 1/x from 2 to 4 using boole's rule
integrate 1/x from 4 to 8 using boole's rule
...
integrate 1/x from 512 to 1024 using boole's rule

И что же мы видим? А мы видим сразу несколько замечательных вещей. Во первых, метод Буля вне конкуренции по точности. Во вторых, на каждом из этих промежутков он даёт точность порядка 10−4, что для нас очень и очень даже неплохо. Ну а в третьих, делает наш метод всё это всего за один интервал. Не за 10, не за 5, а за 1 интервал. Наглядно всё это выглядит примерно так:


Интегральный метод Буля
Рис. №900. Интегральный метод Буля

Исчерпывающее сравнение методов интегрирования, чуть дальше:


Сравнение методов численного интегрирования
Рис. №901. Сравнение методов численного интегрирования

Наглядно видим, что метод Буля ну просто вне конкуренции. Чтобы добить вас окончательно, приведём еще формулу, в которую превращается метод Буля в нашем случае интегрирования функции f(x) = 1 / x. Here it is:


\[ \int\limits_{x_1}^{x_5} {1 \over t} dt \approx {2h \over 45} \left ({7 \over x_1} + {32 \over x_2} + {12 \over x_3} + {32 \over x_4} + {7 \over x_5}\right ) \approx \log(x_5) - \log(x_1) \]
Выражаем логарифм методом Буля

Из этой формулы, на самом деле, становится ясно очень многое, кстати правая её часть, есть ни что иное, как следствие из всем известной теоремы Ньютона-Лейбница о, скажем так, связи значения определённого интеграла и первообразной. Вообще говоря, справедливость этого тождества даже наглядно заметна на первом скриншоте вольфрама — точный результат интегрирования там всегда есть log(2), что есть ни что иное, как например log(8) − log(4), или допустим log(1024) − log(512), я надеюсь, все помнят, как находить разности логарифмов. Но это только начало. Мы теперь подобрались к самой вкусноте…


Бинарный поиск – панацея от всех болезней


Давайте напряжём мозг, и подумаем, как мы формулу выше можем применить. Мы заметили, что нам, вообще говоря, очень выгодно интегрировать промежутки функции, лежащие между соседними степенями двойки. Здесь мы сделаем важное, и довольно очевидное утверждение, о том, что если мы попытаемся интегрировать промежуток, который является подмножеством промежутка между двумя соседними степенями двойки, то для него, наш метод Буля выдаст результат с не меньшей точностью, а, что вполне разумно, даже большей. Но что нам всё это даёт? А это даёт гениально простую идею, что мы ведь можем посчитать точные значения логарифма в точках, являющимися степенями двойки, запомнить их, а дальше действовать по такому вот нехитрому сценарию:

1. Найти для параметра x такие значения k, что выполнялось бы неравенство 2k ≤ x < 2k + 1;
2. Вычислить определённый интеграл методом Буля на промежутке 2k … x;
3. К полученному результату прибавить log(2k), или, с позволения сказать, k ⋅ log(2) (хотя тут запись – не принципиальна);
4. Всё! Натуральный логарифм вычислен. Можем уже себя хвалить, и удивляться простоте действа.

Продолжая думать головой, мы сразу можем сказать, что сложность поиска по степеням двойки уже составит O(log2(x)), поскольку степеней этих, как догадались вы уже, самих по себе не более, чем логарифм по основанию 2 аргумента x. Но это в случае брутального, линейного поиска… А что в ближайшем заголовке-то написано? Да да! Вообще говоря, log2(1024), где 1024 — максимально возможный аргумент в нашей задаче, уже не превышает равен 10, а если учесть ещё и логарифмическую сложность бинарного поиска, то асимптотика нашей задачи принимает вот такой, хоть и весьма наукообразный, но такой сладкий и хорошенький вид:


\[ O(\log_2(\log_2(x_{\max})) + C) \]
Сложность метода

А значит это, что к примеру для нашей задачи, количество абстрактных операций (обращений к памяти, например), необходимых для нахождения нужного промежутка 2k ≤ x < 2k + 1 никогда не превысит логарифма по основанию 2 числа 10, то-есть никогда не будет больше 4. А метод Буля вообще всегда, и при любой погоде будет работать за константу времени, да ещё какую константу! О такой константе можно только мечтать, прямо как о идеальной женщине…


На этом месте, мы покидаем образ математиков-мыслителей, и снова становимся программистами-джедаями. Мы вспоминаем, что в нашем арсенале есть такой мощное оружие, как SSE вплоть до версии 4, и кэшируемая память. Мы готовы сорвать покровы. Поехали.

• Мы помним, требуемая точность невысока. Следовательно использовать будем 4-х байтовые float, что даёт нам 4 операнда на 1 действие. А в формуле Буля всего 5 членов.

• Мы видим, что константу 2/45 в формуле Буля можно внести внутрь, перемножить с другими константами, и получить вот такой вектор констант
\[ \left [ {2 \cdot 32 \over 45}, {2 \cdot 12 \over 45}, {2 \cdot 32 \over 45}, {2 \cdot 7 \over 45} \right ] \]
, который мы, понятное дело, можем заранее посчитать, и на всю жизнь запомнить в выровненной памяти.

• Мы замечаем, что в векторе выше нету первой семёрки. Оказывается, это неспроста. Дело всё в том, что исходя из нехитрых соображений, x1 в первом члене, всегда бует равен 2k, из чего следует, что для каждой из степеней двоек, мы так-же можем посчитать по одной персональной константе, которая будет равна
\[ {2 \cdot 7 \over 45 \cdot 2^k} \equiv {7 \over 45 \cdot 2^{k - 1} } \]
, причём как именно считать, не принципиально, ибо сделано это будет один раз в жизни программы. В расчётах, остаётся только эту константу множить на h, и прибавлять к результату. Красота.

• Мы ощущаем, что слушающие переживают по тому поводу, что несмотря на то, что мы уже знаем значения x1 и x5, все остальные иксы остаются для нас пока загадкой. Развеем эти страхи, заметив, что значения x2, x3, x4 и x5 можно вычислить с помощью вектора констант приблизительно вот так
\[ \overline{X}_{2 \mathellipsis 5} = \left (x - 2^k\right ) \cdot \left [{1 \over 4}, {1 \over 2}, {3 \over 4}, 1\right ] + \left [2^k, 2^k, 2^k, 2^k \right ] \]
. Я думаю уже стало понятно, почему все векторы содержат именно 4 компонента (если понятно не всем, то напомню, что в XMM регистр помещается именно 4 float'а), как и тот простой факт, что 4-вектор с одинаковыми значениями из одного числа можно получить великой и могучей командой SHUFPS, с нулевым последним параметром. Остаётся напомнить, что вектор с дробями надо запомнить. В выровненной памяти. На всю жизнь.

• Мы помним, что
\[ h = {x_5 - x_1 \over 4} = {x - 2^k \over 4} = x_2 - 2^k \]
. А вы, посмотрю я, уже начали переживать… Заметим, что всё до ужаса упрощается, если сразу после умножения на вектор, мы будем сохранять первую его компоненту как h. Не прибавляя к нему 2k. Так, собственно, и сделаем.

• Мы просто возьмём, и возложим вычитание 2k из x на команду SUBSS, умножение двух векторов на одну команду MULPS, сложение оных на ADDPS, делить на ноль x1 … x4 заставим команду DIVPS — поверьте, она не сможет нам в этом отказать, ну а чтобы всё это умножить на вектор из h, и феерично между собой сложить (если точнее, вычислить скалярное произведение двух векторов), мы привлечём тяжелую артиллерию из возможностей SSE4 — одну единственную, брутальную и устрашающую команду DPPS, которую именно для этой цели и придумали. После её работы, останется только прибавить к результату константу первого члена, умноженную на h, добавить к ней log(2k), выпекать 15 минут на медленном огне, и всё. Мы только что посчитали натуральный логарифм.

• Мы ещё отметим, что для бинарного поиска будем использовать команду UCOMISS для прямого сравнения нашего x с уже посчитанными значениями 2k, хранящимися в недалёкой памяти. Сам бинарный поиск реализуем без рекурсии, с помощью двух границ, с выбором средней точки. Продолжать будем, пока границы не сойдутся.

Теперь мы знаем почти всю теорию процесса. Тут надо сразу сказать, что метод с бинарным поиском, хоть и показал не самые плохие результаты, но на самом деле, оказался далеко не идеальным. Но об этом чуть позже. Сейчас мы перейдём от слов к делам, и покажем первую часть кода, который подготавливает для нашего метода все многочисленные константы, которые могут ему потребоваться. Here it is:


Код: Delphi
type
  TSingleArray = array[0..0] of Single;
  PSingleArray = ^TSingleArray;

var
  perz, perz_o:PSingleArray;
  w:Single;

begin
    aligned_get_mem(Pointer(perz_o), Pointer(perz), 256 * 4, 16);

    w := 2.0 / 45.0;

    for i:=0 to 15 do begin
      perz[i] := 1 shl i;                 // 2^x
      perz[i + 16] := ln(perz[i]);        // ln(2^x)
      perz[i + 32] := w * (7.0 / perz[i]);// первый член: 2/45 * 7*f(x1)
    end;

    i := 64;

    perz[i] := 0.25;
    perz[i + 1] := 0.5;
    perz[i + 2] := 0.75;
    perz[i + 3] := 1.0;

    perz[i + 4] := 32.0 * w;
    perz[i + 5] := 12.0 * w;
    perz[i + 6] := 32.0 * w;
    perz[i + 7] := 7.0 * w;
end;

Если вы знакомы с предыдущей статьёй из этого цикла, то вы уже знаете, как работает, и зачем применяется aligned_get_mem, но а если не знакомы, то поясним, функция это занимается отведением блока памяти так, чтобы адрес на него указывающий оказался выровненным. Оригинальный адрес блока после вызова GetMem также сохраняется, с целью блок потом корректно освободить. Просто SSE штука такая, что без выровненной памяти она смотрит на нас, как на несмышлёных школьников… Ну а в остальном тут все просто. Память мы выделяем под массив, и в разные его части для всех x записываем сначала 2x, затем log(2x), и заодно вычисляем коэффициент первого члена формулы. В нашем примере x пробегает значения от 0 до 15, чего нам хватит аж для чисел порядка 215, то-есть до 65535. В конце мы ещё вычисляем 2 вектора констант, первый для вычисления иксов, второй — коэффициенты. Адреса на начала этих блоков, как можно отметить, тоже получаются выровненные.


Теперь, самое время взглянуть на более интересную часть кода, которая, собственно, занимается бинарным поиском, и делает все вычисления. В данном куске кода, в качестве указателя на таблицу предподсчитанных констант используется регистр eax, переменные vin и vout объявлены как Single, и служат для алгоритма входными и выходными данными, соответственно. То-есть, из первой переменной алгоритм читает x, а во вторую записывает посчитанное значение логарифма. Собственно, как всё это выглядит:


Код: Delphi
//pointer to data table
mov eax, perz

// ecx всегда в старших битах ноль
xor ecx, ecx

//input X value
movss xmm0, vin

//binary search closest left 2^x
//DL : Left bound
//DH : Right bound
//dh = 12, dl = 0
mov dx, 3072

// собстно, бинпоиск
@xloop:

  //M = floor((R + L) / 2)
  mov cl, dh
  add cl, dl
  //делим на 2
  sar cl, 1

  //сравниваем с очереджной степенью двойки
  ucomiss xmm0, dword ptr [eax + ecx * 4]

  //обновляем границы
  jc @xless

	mov dl, cl
	inc dl

	//while L < R
	cmp dl, dh
	// JNZ значительно быстрее JL - факт
	jnz @xloop
	jmp @xexit

  @xless:

	mov dh, cl

	//while L < R
	cmp dl, dh
	// JNZ значительно быстрее JL - факт
	jnz @xloop

@xexit:

// store found index into ECX
mov cl, dl
// calculate the effective address
lea ecx, [eax + ecx * 4]

// XMM7 = 2^x (closest)
movss xmm7, dword ptr [ecx - 4]

// XMM0 = X - 2^x
subss xmm0, xmm7

// заполняем регистр первым элментом
shufps xmm0, xmm0, 0
shufps xmm7, xmm7, 0

// XMM0 = X2,X3,X4,X5
mulps xmm0, dqword ptr [eax + $100]

// XMM2 = H (in Boole's rule), H = (X - 2^n) / 4
movss xmm2, xmm0

// XMM5 = Коэффициенты при X2..X5
movaps xmm5, dqword ptr [eax + $110]

// XMM0 += 2^x
addps xmm0, xmm7

shufps xmm2, xmm2, 0

// XMM5 = (2/45 * Coeffs) / X2,X3,X4,X5
divps xmm5, xmm0

// XMM5[0] = (F(X2) + F(X3) + F(X4) + F(X5)) * H
dpps xmm5, xmm2, $F1
// для более древних компилеров
// db $66, $0F, $3A, $40, $EA, $F1

  (*********************************
  более медленная, как ни странно, альтернатива DPPS
  безнадежно устарела, не успев появиться на свет
  movaps xmm7, xmm5
  shufps xmm5, xmm5, 27
  addps xmm5, xmm7
  movaps xmm7, xmm5
  shufps xmm5, xmm5, 177
  addss xmm5, xmm7
  *********************************)

// X1 = CoeffX1 * H
mulss xmm2, dword ptr [ecx + $80 - 4]

// XMM5[0] = C*F(X1) + TEMP(X2..X5) + LN(2^x)
addss xmm5, dword ptr [ecx + $40 - 4]
addss xmm5, xmm2

movss vout, xmm5

Собственно, всё что мы тут видим, мы детально описали в теории выше. В коде, хоть и используются некоторые ухищрения (о которых мы расскажем чуть позже), всё равно присутствует довольно много не оптимальных мест, которые к концу статьи мы все (или почти все) искореним. Ведь наша задача — получить максимально быстродействующий вариант, не так-ли? :) Ах да, быстродействие… Вот что касается последнего, то у этой версии, оно не слишком радует.


Для определения быстродействия кусков подопытного кода, мне не пришло ничего в голову умнее, чем прогонять этот код в цикле из 100 млн. итераций, и банально замерять время с помощью QueryPerformanceCounter(), хотя понятное дело, что метод этот весьма неточный, шаткий и плавучий, но я вас уверяю, приблизительную оценку получить он всё-таки позволяет, а нам, в общем-то, большего и не требуется. Алгоритм подвергся отдельным испытаниям, в общей сложности, на 4-х различных тестах, чтобы в совокупности хоть как-то имитировать боевые условия использования. Тесты были следующие: в первом, 100 млн. раз вычислялся логарифм константы 431 + π, во втором на вход подавались совершенно случайные значения, в третьем значение некоторым образом зависело от номера итерации, и в четвертом вычислялось что-то вроде сходящегося ряда — каждое входящее значение зависело от предыдущего результата. Псевдокод алгоритмов будут чуть ниже, пока заметим, что этим тестам сначала подвергся и обычный логарифм из FPU, то-есть команда FYL2X, с целью сравнения результатов, а именно, за эталон логарифма на базе FPU взялся вот такой нехитрый код:


Код:
fldln2
fld vin
fyl2x
fstp vout

Всё, как видим, честно, здесь присутствуют точно такие-же обращения к памяти (переменным vin и vout), как и в нашем коде. Разве что, выглядит он весьма менее объёмистым, но поверьте, это совершенно пока ни о чём не говорит :) Что-либо сказать нам смогут тесты, а вот, кстати, и результаты оных:


Тестируем подсчёт логарифма бинарным поиском
Рис. №902. Тестируем подсчёт логарифма бинарным поиском

Честно говоря, результаты эти меня слегка огорчили. Нет, ну правда, я ожидал большего… Получается, выигрыш мы получили только на константе, и на повторяющихся аргументах, а на случайных значениях наша функция лажает, и иногда работает даже медленнее FPU. А во фрактале, значения в основном, как раз-таки довольно-таки случайные, хотя и стоящие довольно близко друг к другу… В общем быстро я понял, что что-то тут надо менять. И вот, в один прекрасный момент, когда меня посетило Его Величество Вдохновение™…


Так стоп, а что мы вообще можем тут предпринять и оптимизировать? Мне очень не нравилась громоздкость задачи бинарного поиска, хотя она и была довольно-таки неплохо оптимизирована. Замена коротких регистров на 32-битные, кстати (и еще некоторые мелочи), дало прирост примерно в 50 миллисекунд на всех тестах, но этого было мало. Нужно было пересматривать концепцию… И как не странно, отказываться от громоздкого двоичного поиска.


На самом деле, давайте слегка пересмотрим задачу поиска ближайшей степени двойки. Как например её можно перефразировать? А просто ведь. Если рассматривать целые числа в двоичной системе счисления, то ближайшая меньшая степень двойки, будет ни чем иным, как индексом первого значащего бита справа. И всё. Никаких двоичных поисков не надо. Но вы скажете, что у нас не целые числа. А мы возразим, что если переконвертировать вещественное число в целое с отсечением дробной части (командой CVTTSS2SI, например), то ближайшая его степень двойки от этого измениться не может. Тогда что же нам остаётся сделать? Верно, просто взять, и найти ближайший ненулевой бит справа, и запомнить его индекс. Что я и сделал. Но поделка моя работала на маленьких числах еще медленнее, чем двоичный поиск, хотя и давала небольшой выигрыш на больших числах. Удивляться тут в общем-то нечего, ибо сложность при таком подходе переставала быть логарифмической, и для, допустим, числа π, и максимуме равном 215, мы имели, ни много ни мало, около 14 сравнений. Пусть чуть более быстрых чем при бинарном поиске, но, как видим, не намного. На этом месте мне показалось, что это тупик, но вдруг, в одну прекрасную бессонную ночь, читая мануал от Intel, я вдруг обнаружил…


BSR — Bit Scan Reverse


Чёрт! Я вдруг обнаружил, что совершенно незнаком со многими, нечасто используемыми командами архитектуры x86… Среди которых я случайно заметил интересную такую команду BSR, которая, судя по описанию… Да. Решала нашу задачу поиска первого значащего бита одной командой. Одной! Без всяких циклов, и прочих свистелок. Я конечно, слегка офигел, но придя в себя, по-быстрому прикрутил команду к одной из версий своего кода. И знаете — помогло.


Считаем логарифм на SSE с помощью BSR
Рис. №903. Считаем логарифм на SSE с помощью BSR

При самой первой проверке, скорость работы подскочила примерно в 1.6 раз, в сравнении с бинарным поиском, но, как оказалось, я исчерпал не все возможности низкоуровневой оптимизации. Взглянем сразу на то, что в итоге получилось. В коде ниже, в качестве указателя на таблицу с константами, как можно заметить, используется регистр EBX:


Код: Delphi
//input X value
movss xmm0, vin
cvttss2si edx, xmm0

//ищем эту вашу степень двойки
bsr edx, edx

// xmm3 = 2^x (closest)
movss xmm3, dword ptr [ebx + edx * 4]

// XMM0 = X - 2^x
subss xmm0, xmm3

// заполняем регистр первым элментом
shufps xmm0, xmm0, 0
shufps xmm3, xmm3, 0

// XMM0 = X2,X3,X4,X5
mulps xmm0, dqword ptr [ebx + $100]

// XMM2 = H (in Boole's rule), H = (X - 2^n) / 4
movss xmm2, xmm0

// XMM5 = Коэффициенты при X2..X5
movaps xmm5, dqword ptr [ebx + $110]

// XMM0 += 2^x
addps xmm0, xmm3

shufps xmm2, xmm2, 0

// XMM5 = (2/45 * Coeffs) / X2,X3,X4,X5
divps xmm5, xmm0

// XMM5[0] = (F(X2) + F(X3) + F(X4) + F(X5)) * H
dpps xmm5, xmm2, $F1

// X1 = CoeffX1 * H
mulss xmm2, dword ptr [ebx + edx * 4 + $80]

// XMM5[0] = C*F(X1) + TEMP(X2..X5) + LN(2^x)
addss xmm5, dword ptr [ebx + edx * 4 + $40]
addss xmm5, xmm2

movss vout, xmm5

И всё. По сравнению с кодом предыдущим, тут используется всего 2 регистра общего назначения, и раза в 2 меньше команд. Тот случай, когда решающую роль играет знание железа, а не алгоритмики… Кстати, между словом, инициализируются необходимые для этого кода константы точно так-же, как и для предыдущего.


Наверно, можно было пойти в оптимизации и дальше, но меня устроил и такой результат. Прирост скорости на случайных числах, как видно из картинки, составил почти 90%, и это ведь еще при том, что тесты включают в себя время вычисления собственно самого рандома, которое, на самом деле, не такое уж и маленькое. В общем нам такой скорости точно хватит.


Напоследок, расскажем немного об основных методах низкоуровневой оптимизации ассемблерного кода, по крайней мере тех, которые мне помогли сэкономить в общей сложности почти 200 миллисекунд.

• Для начала, нужно хорошенько продумать алгоритм, и по возможности выкинуть лишние команды. Например, если надо сложить значение регистра со значением переменной в памяти, то совсем не обязательно его сначала загружать в другой регистр. Можно, и нужно, это делать одной командой сложения, в качестве второго операнда указывая сразу адрес в памяти. Хотя вообще, это должно быть очевидно.

• Вычисления адресов в командах наподобие addss xmm5, dword ptr [ebx + edx * 4 + $40], оказывается, делаются значительно быстрее, чем предварительное их вычисление с сохранением в регистр, предположим командой LEA. На этом удалось сэкономить около 20-30 миллисекунд. Именно поэтому из финального кода LEA я убрал.

• Короткие регистры — зло. Нет, ну правда, заменив DL и DH на полноценные EAX и EDX, я получил прирост почти в 50 миллисекунд. Касается, похоже, абсолютно всех команд. Видимо это связано с особенностями устройства железа процессора — ну удобнее ему с 32 битами работать…

• Команды условных переходов, оказывается, довольно сильно различаются по времени исполнения. Например, JC работает значительно быстрее JL (чуть-ли не в 2 раза), а JZ работает чуть-чуть быстрее JC. Учет этой особенности позволяет неплохо экономить наносекунды.

• Можно еще экономит на прыжках по другому, а именно, неплохо помогает дублирование кода :). В первом коде с бинарным поиском это как можно заметить — сэкономило 10-20 миллисекунд.

• Наш враг — команды, зависимые друг от друга, наш друг — суперскалярность! Если вы хотите, чтобы ваш код работал быстро, старайтесь команды группировать небольшими блоками так, чтобы в этих блоках команды друг от друга не зависели, то-есть, например, не использовали одни и те-же регистры на запись, одни и те-же адреса памяти и т.д, тогда, есть приличная вероятность того, что ваш суперскалярный процессор их будет выполнять практически параллельно! Что, конечно-же положительно скажется на скорости. Так положительно, что до позеленения переставляя местами команды SSE в коде выше, мне удалось выжать из него прирост почти в 75-100 миллисекунд. Эпик вин, в общем! обратите, кстати, внимания, как они расположены: почти каждая команда не зависит от другой, идут как-бы попарно, выполняя при этом две задачи практически параллельно.

• Сто́ит ещё напомнить, что такие часто встречающиеся операции с целыми числами, как умножение и деление на 2 (в том числе, и вычисление степеней двойки), можно, и нужно выполнять с помощью битовых сдвигов, а именно, командами SAR и SAL (по скорости работы не отличаются от SHR/SHL). Но думаю эти базовые вещи известны всем.

• Ну и последний совет, на примере команды BSR: читайте мануалы от Intel (если мы говорим про x86-64). Это действительно бывает очень полезно. Только оттуда я узнал, как-же на самом деле работает CPUID, и что, оказывается, существует команда такая — VBROADCASTSS, с помощью которой, существуй её аналог под SSE, можно было бы из кода убрать все SHUFPS X, X, 0 (но не судьба — это уже частички расширения AVX).


На этой мажорной ноте, наше маленькое исследование, наверное, можно считать оконченным, результат мы, в общем-то, получили неплохой… Хотя ещё как сказать. Будем, всё-таки считать, что он неплохой :). Напоследок еще надо сказать, что метод наш можно с таким-же успехом использовать для вычисления отрицательных логарифмов, то-есть когда аругмент меньше единицы, и больше нуля. В этом случае, можно единичку поделить на аргумент, и для полученного числа использовать ту-же тактику поиска ближайшей степени двойки, только на этот раз верхней. Ну, и в конце поделить 1 на найденную степень. Ну а дальше вычесть и проинтегрировать. Знак только не забыть сменить. Всё и здесь довольно просто…


Полученный в результате код, кстати, уже встроен в мою программу-генератор фракталов, и отлично там работает. В следующей статье, надеюсь, основательно займёмся параллельными вычислениями, прогрессированием, и синхронизацией. А сейчас… До скорых встреч, друзья!

Ключевые слова математика, Программизм, Множество Мандельброта, FYL2X, логарифм, SSE, метод Ньютона, интеграл, Boole's rule, метод Буля, DPPS, бинарный поиск, BSR, оптимизация
QR-код
Интересное
А теперь, внимание, вопрос: а можно ли обменять значение двух переменных так, чтобы третью явно или косвенно не использовать? Оказывается, можно.

Читать »»
Случайные фото