Category Archives: assembly

О реализации целочисленного аналога функции sqrt

Год назад я написал пост по обращению матриц определённого вида. Он возник не на пустом месте. Исходной посылкой послужило желание сконструировать способ качественного определения быстродействия маленьких быстро выполняющихся функций в условиях, «приближенных к реальности». На тот момент я пытался сравнивать быстродействие различных вариантов реализации функции вычисления целочисленного квадратного корня.

При вычислении sqrt(x) в плавающей точке обычным методом является вычисление 1/sqrt(x) методом Ньютона и затем домножение на x (основная хитрость практической реализации данного метода заключается в способе выбора начального приближения). Эффективно применять этот метод, используя только целочисленные операции, вряд ли возможно, поэтому я решил посмотреть, насколько хорошо (по крайней мере, на конкретном процессоре) в этой задаче работает прямой лобовой метод — двоичный поиск.

Но даже двоичный поиск можно реализовать по-разному. И естественно хочется уметь сравнивать быстродействие разных вариантов в разных условиях.

Неожиданным открытием при первом простом тестировании стала зависимость скорости выполнения функции от её положения в asm-файле (для ассемблерных вариантов). Директива align перед циклами, напротив, не продемонстрировала устойчивого эффекта. Наиболее вероятным «виновником» в данном случае является предсказатель ветвлений и особенности его работы для конкретного CPU. Также может сказываться работа L1I и TLB (здесь, в частности, может удачно или неудачно подействовать директива align).

В случае «простого» тестирования вроде многократного вызова одной функции в цикле есть немалый шанс натолкнуться на проявление таких особенностей во всей их красе. Поэтому я решил поступить несколько хитрее: выполнить псевдослучайное перемешивание разных вызовов. В каждом прогоне — разное количество вызовов разных функций. Выполнив столько же прогонов, сколько у нас есть функций, получим вектор суммарных времён и матрицу количеств вызовов. Решив линейное уравнение, получим оценки времени вызова каждой функции. Можно сделать и больше прогонов и получить решение с помощью МНК (здесь не реализовано). В итоге я всё же не стал использовать аналитически обращаемую матрицу, а задействовал библиотеку uBLAS из состава Boost. См. также о влиянии фрагментации кучи на результаты тестирования.

Использовалась среда Microsoft Visual C++ 2015, компиляция под x86-64 с O2. Причёсывать код однолетней давности я уже поленился — там есть, что выкинуть, но если бы я этим занялся, то… В общем, год уже и так прошёл.

Замечание. В MSVC все исходники .c, .cpp, .asm должны иметь разные имена, так как при компиляции каждого из них порождается одноимённый .obj файл, которые сваливаются в одну общую папку и могут затирать друг друга с печальным результатом во время компоновки.

Итак, исследуемые варианты:

  1. cpp1, asm1 — цикл + ветвление;
  2. cpp2, asm2 — цикл + cmov;
  3. cpp3, asm3 — цикл + ветвление + ilog2 (целочисленный логарифм позволяет выбрать хорошее начальное приближение, для его реализации на x86 используется команда bsr или lzcnt);
  4. cpp4, asm4 — цикл + cmov + ilog2;
  5. asm5 — развёрнутый вариант asm3;
  6. asm6 — развёрнутый вариант asm4 (computed goto);
  7. fpu — использует команды SSE для вычисления isqrt (стандартный метод);
  8. null — ничего не делает, используется для оценки времени вызова.

Итак, результаты тестирования представлены на графике ниже. Числа получены следующим образом: в каждом прогоне использовалось 100’000 псевдослучайных чисел-аргументов, равномерно распределённых от нуля до числа, указанного по оси 0x (1k — 1024, 32k — 32’768, 1M — 220, 1G — 230), выполнялось 1100 повторений, из которых верхние и нижние 5% полученных оценок времён отбрасывались, а от оставшихся 1000 бралось арифметическое среднее. Всего выполнялось по 10 таких прогонов, из результатов которых для каждой функции выбиралось минимальное значение и затем вычиталось соответствующее минимальное значение функции null. В конце все результаты нормировались по функции fpu, так что по оси 0y графиков мы видим, во сколько раз все варианты медленнее вычисления isqrt на fpu.

Процессор, на котором производился тест — AMD K10. Итак, мы видим, что «бодаться» с FPU на «больших» процессорах (пусть даже старых) бессмысленно. Пиррова победа в режиме 1k у asm4/cpp4. Применение ilog2 в комбинации с cmov даёт весьма неплохой эффект, который может оказаться ещё выше на простых процессорах без хорошего предсказателя переходов. А вот разворачивание цикла здесь не помогает (впрочем, опять же на «простом» процессоре ситуация может быть другой).

GLU-1: SSE2 реализация Vec2

Введение

Старый вариант GLUPoly переименован в GLUPoly-1.0.7z. Ниже описаны отличия GLUPoly-1.1.7z.

Написана реализация Vec2 поверх SSE2, использующая интринсики (помещена в Vec2sse2.h). Старый Vec2.h переименован в Vec2std.h. Новый Vec2.h включает Vec2sse2.h, если компилируется MSVC под x64 или с параметром /arch:SSE2 (или /arch:AVX), иначе включает Vec2std.h.

Некоторые изменения были внесены в Bbox2.h (ускорение построения Bbox2 по набору Vec2).

Проблема 1

Для эффективного использования SSE следует выравнивать данные по 16-байтной границе, поэтому значения типа __m128d (пара 64-битных чисел с плавающей точкой, упакованная в 128-битное значение) выравниваются по 16-байт. Это свойство автоматически передаётся всем структурам/классам, включающим в себя __m128d.

Однако стандартные средства выделения памяти не гарантируют выравнивания нестандартных размеров (т.е., например, new int[100] будет выровнен по 4-байтной границе для 4-байтного int, а вот new __m128d[100] по 16-байтной границе уже может быть и не выровнен).

Простое решение этой проблемы достигается перегрузкой operator new и operator delete (попарно всех нужных форм, избыточные варианты delete требуются из-за того, что в случае бросания исключения из конструктора автоматически вызываются только «симметричные» delete, поэтому в случае отсутствия такого delete память не будет освобождена). Эти перегруженные операторы могут выделять память системно-зависимыми функциями, гарантирующими нужное выравнивание.

Недостаток такого решения — невозможность использования с контейнерами, использующими rebind::other аллокатора (например, std::list).

Принятое решение основано на двух дополнительных классах, размещённых в новом файле AlignedAllocator.h.

AlignedNewDeleteImpl<Align> реализует набор operator new/operator delete, выравнивающих по Align байт (Align == 16 для __m128d и, соответственно, SSE2-версии Vec2). Класс, которому нужно обеспечить выравнивание своих объектов, может наследовать от AlignedNewDeleteImpl. С другой стороны, AlignedNewDeleteImpl является обёрткой системно-зависимых функций выделения памяти (только отсюда они вызываются напрямую).

Замечание. Vec2 для простоты использует public-наследование. С точки зрения принципов использования наследования следует использовать private-наследование, однако оказалось, что using AlignedNewDeleteImpl::operator new и т.д. в public-секции не работает. Пока мне не ясно, с чем это связано. Конечно, можно это обойти, написав явную переадресацию вызовов.

AlignedAllocator<T, Align> реализует STL-совместимый аллокатор, для выделения объектов T средствами AlignedNewDeleteImpl. Его можно использовать непосредственно (Align по умолчанию устанавливается равным __alignof T) в качестве параметра шаблона STL-контейнера. «Правильный» аллокатор решает проблему rebind::other.

Кроме того, для удобства можно определить std::allocator для своего класса, чтобы каждый раз не приходилось указывать аллокатор (т.к. стандартный аллокатор всё равно не подходит, и его использование может оказаться латентной ошибкой), что и сделано для Vec2 путём наследования std::allocator<Vec2> от AlignedAllocator<Vec2>.

Проблема 2

Компилятор может «отказаться» передавать объекты с выравниванием поверх стандартных встроенных типов в функции как копии. Например, void foo (const __m128d&); компилируется, а void foo (__m128d); — нет.

К сожалению, в std::vector нашлась такая функция, а именно — resize (size_type, T), принимающая заполняющее значение как копию. Можно по-разному пытаться обойти эту проблему. В данном случае я выбрал решение «в лоб»: написать свою реализацию динамического массива Vec2: Contour (Contour.h и Contour.cpp). Данная реализация также использует SSE2 (для копирования /заполнения диапазонов), если был включен Vec2sse2.h.

UPDATE 08/12. В стандартной библиотеке C++11 эта функция изменена и теперь принимает ссылку. Таким образом, можно использовать std::vector из современных дистрибутивов STL (например, в составе Visual Studio 2012).

Тестирование

Solution теперь содержит два проекта: исходный GLUPoly и вспомогательный Tests (содержит только один файл test.cpp). Цель второго проекта — аккумуляция unit-тестов, проверяющих код из GLUPoly. На данный момент тестами покрыты Vec2, Bbox2 и Contour. Тесты построены на библиотеке Boost.Test.

Своевременное написание тестов позволяет быстро найти и исправить множество ошибок, способствует повышению производительности труда и качества кода.

Ещё

Добавлены FusedFunctor.h и FusedFunctorTuple.h, но о них я намерен написать в следующий раз (планируется внести некоторые изменения).

ReBits 32 x64

В ReBits 32 и ReBits 32 x4 использовался компилятор под IA-32 (часто обозначаемую как «x86»). Однако, почти все x86-совместимые ЦП, выпускаемые в настоящее время, поддерживают 64-битные расширения (часто называемые «x86-64» или, сокращённо, «x64»), предложенные AMD в 2000 (первый процессор вышел в 2003). ОС, работающие в 64-битном режиме всё шире распространяются на домашних ПК. Поэтому, на мой взгляд, именно x64, а не IA-32 следует рассматривать как основную целевую платформу для новых приложений, особенно, если они должны будут обрабатывать большие объёмы данных.

Система команд x64 включает в себя обязательную поддержку SSE2, причём количество регистров xmm удвоено (до 16). Количество РОНов также удвоено, что, в частности, позволяет ускорить вызов функций (использовать больше регистров для передачи параметров). Посмотрим, как будут работать полученные ранее решения в 64-битном режиме.

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

Для условной компиляции кода, имеющего разные версии под IA-32 и x64, MSVC предлагает макрос _M_X64 (определён для целевой платформы x64).

Рассматривались следующие варианты: N0, H0, H1 и цикл H5, как написанные на C и не требующие изменений перед перекомпиляцией под x64. Были переписаны H5 (как лучший доступный вариант для обработки одного 32-битного блока), циклы с использованием SSE2 и SSSE3. MSVC предлагает следующие заголовочные файлы:

  • intrin.h — набор всех интринсиков (включает все ниже перечисленные), новый вариант H5 использует _rotl/_rotr (вместо rol/ror: циклический сдвиг влево/вправо) и _byteswap_ulong (вместо bswap: обращение порядка байт);
  • mmintrin.h — команды MMX;
  • xmmintrin.h — команды SSE, включает mmintrin.h;
  • emmintrin.h — команды SSE2, включает xmmintrin.h;
  • pmmintrin.h (Prescott) — команды SSE3, включает emmintrin.h;
  • tmmintrin.h (Tejas) — команды SSSE3, включает pmmintrin.h;
  • smmintrin.h (?) — команды SSE4.1, включает tmmintrin.h;
  • nmmintrin.h (Nehalem) — команды SSE4.2, включает smmintrin.h;
  • wmmintrin.h (Westmere) — команды AES NI и CLMUL, включает nmmintrin.h;
  • immintrin.h (?) — команды AVX, включает wmmintrin.h;
  • ammintrin.h — команды, предложенные AMD (бывший SSE5).

Удвоенное количество регистров позволяет хранить инвертированные маски (и сократить таким образом количество инструкций), а также обрабатывать не одну четвёрку, а, например, две четвёрки за итерацию цикла. При этом вторая загрузка происходит после первого блока вычислений. Теоретически такое «смешивание» команд должно минимизировать простой XMM блоков ЦП. Итак, переписанные на интринсиках варианты: для SSE2, для SSSE3. Сравним их результаты (за 100% приняты соответственно H1 и цикл на H5, обозначенный H5x4).

Можно заметить, что SSE2 улучшил свои относительные показатели. N0 не показан из-за его низкой скорости (чтобы не портить масштаб), кроме того, для N0 переход на x64 ничего не дал (на Deneb даже стало медленнее на 11%). Ниже представлена диаграмма, показывающая отношение времени, затраченного x64 версией, ко времени IA-32 версии. Для большинства рассмотренных вариантов переход на x64 оказался полезен.

ReBits 32 x4

В ReBits 32 рассматривалась задача обращения одного 32-битного целого в рамках набора команд P6. Обычно, код, который стоит оптимизировать, используется внутри циклов с большим числом повторений, и обрабатывает большое количество данных (если вычисления достаточно просты – как в нашем случае). SIMD расширения архитектуры x86 позволяют обрабатывать 2 (MMX), 4 (SSE) или 8 (AVX) 32-битных целых «за раз» одной командой. MMX не сулит большого преимущества относительно SSE (под SSE подразумеваются все семейства команд вплоть до SSE4.2), а AVX не поддерживается процессорами, имеющимися у меня в наличии. Поэтому далее рассматривается SSE (фактически два варианта: для SSE2 и SSSE3).

Ни трюк с lea, ни циклический сдвиг недоступны в рамках SSE, поэтому отталкиваться будем от варианта H1. Сравнивать будем с развёрнутым циклом, вызывающим H5. Тест заключался в 200-кратном прогоне массива в 32Мб.

Итак, вариант для SSE2.

  1. Перемещение данных по 16 байт происходит эффективнее, если адрес выровнен по 16 (остаток от деления на 16 равен 0), поэтому кусочки массива, выдающиеся за границы выравнивания, обрабатываются отдельными циклами, вызывающими H5.
  2. К регистрам SSE (xmm0 – xmm7) нельзя применять константы, заданные непосредственно в коде, поэтому маски стоит разместить в регистрах перед входом в основной цикл, делается это путём помещения 32 битов в РОН, из которого они перемещаются в младшую четверть xmm-регистра, а затем копируются в оставшиеся четверти (broadcast).
  3. Инвертированные маски хранить необязательно, так как есть команда pandn (and-not).
  4. Последний обмен (16-битных половинок) осуществляется комбинацией команд pshuflw/pshufhw.
  5. Команда prefetchnta в начале цикла «намекает» процессору, что полезно начать подгружать данные по указанному адресу, которые мы будем читать и писать. Это может увеличить скорость обработки больших массивов (сотни элементов и больше).

SSSE3 позволяет сократить и ускорить код благодаря наличию команды pshufb, позволяющей произвольно переставлять байты внутри 128-битного xmm-регистра. Уровень поддержки процессором тех или иных SIMD расширений можно проверить с помощью инструкции cpuid непосредственно или используя соответствующий интринсик, как в этом примере. Со стороны пользователя проверка выполняется, например, так

if (doesCpuSupport().SSE3) { /* do SSE3 code */ }

См. также, например, AMD cpuid reference. Результат тестирования представлен на диаграмме (Deneb не поддерживает SSSE3).

ReBits 32

Введение

Предположим, что перед нами поставлена простая задача: дано 32-битное целое, поменять порядок бит в нём на обратный (т.е. обменять биты 0-31, 1-30, 2-29 и т.д.). Я не мог решить, как лучше назвать эту операцию: «reflect bits» или «reverse bits», поэтому сократил оба варианта до «rebits».

Говорят, что обычно существует по крайней мере три алгоритма решения той или иной задачи: «наивный», «умный» и «гениальный». Гениального придумать не получилось, поэтому рассмотрим наивный (N*) и умный (H*) алгоритмы обращения порядка бит в 32-битном целом.

Наивный алгоритм

Итак, напишем цикл, перебирающий биты в исходном числе и размещающий их в отражённую позицию в результирующем числе. На стандартном Си это можно сделать так: N0 (код будет работать не только для 32 бит). Какова его эффективность? Для замера скорости работы я прогонял все варианты в цикле (228 вызовов) и вычитал время прогона функции-заглушки, возвращающей свой аргумент (способ непрофессиональный, но простой). Использовался 32-битный компилятор Microsoft Visual C++ 2010SP1 (был задан /Ox и ряд других оптимизационных параметров). Итак, приведённый код требует в среднем более 100 тактов ЦП на обработку одного 32-битного целого. Не впечатляет.

Вставки на ассемблере

Компилятор Microsoft не является самым агрессивным и эффективно оптимизирующим, тем не менее, признан довольно хорошим. Соревноваться с ним в умении производить быстро работающий код на ассемблере обычно не имеет смысла. Тем не менее, существуют ситуации, когда приходится спускаться с «небес» языка высокого уровня на «землю» машинных инструкций. Компиляторы C/C++ обычно предоставляют три способа сделать это:

  1. Написать реализации «ускоряемых» функций полностью на ассемблере, ассемблировать их в объектные файлы и компоновать с кодом на С, код для разных целевых платформ размещается в отдельных файлах.
  2. Использовать ключевое слово asm или его аналог (__asm в Visual C++) для размещения машинных инструкций прямо в теле функций на С, пролог/эпилог вызова функции компилятор сделает сам.
  3. Использовать встроенные функции (intrinsics; у каждого компилятора свои), являющиеся «оболочками» (оформленными в виде функций С) инструкций ЦП, не доступных в языке высокого уровня непосредственно.

Принято считать, что (1) предпочтительнее, чем другие варианты, потому что в этом случае код на C/C++ не загрязняется нестандартным платформозависимым кодом. В свою очередь, (3) считается предпочтительнее (2), так как при использовании интринсиков код обычно получается короче, понятнее и допускает дополнительную оптимизацию компилятором. Пока что будем использовать способ (2), потому что он даёт максимальный контроль над получаемым кодом.

На всякий случай, приведу прописную истину: в реальной работе необходимо, в первую очередь, выбирать оптимальный алгоритм решения задачи и только потом, если потребуется, проводить низкоуровневую оптимизацию. Далее мы увидим, что самый медленный вариант «умного» алгоритма работает быстрее самого быстрого варианта «наивного» алгоритма.

Итак, перепишем N0 на ассемблере: N1. Простой код, никаких изысков, проверим оптимизатор? Удивительно, но этот ассемблерный код действительно быстрее, хотя и незначительно (примерно на 5%).

В наборе инструкций x86 есть команды, которые могут пригодиться в процессе перебора бит. Например, bsf/bsr (bit search forward/reverse), которые определяют номер самого младшего/старшего установленного бита. Может быть, они позволят ускорить цикл? Проверим: N2. Увы, этот код медленнее N0 (на 2-5%). Команды манипуляции битами «сложные» и требуют несколько тактов на выполнение.

Как ещё можно улучшить исходный цикл? Попробуем использовать «двигающуюся маску»: в регистр ecx положим 231, а в edi – 20, будем тестировать бит первым регистром и, в зависимости от результата (используя команду cmov), добавлять 0 или edi к результату. На каждой итерации будем сдвигать ecx на бит вправо, а edi на бит влево. Итак: N3, особого прироста скорости опять не получилось (на 8% быстрее N0 на AMD K10), зато он привёл меня к мысли о том, как устанавливать не один, а два бита за итерацию цикла, используя циклические сдвиги (команды rol/ror). На примере байта:

abcd efgh (исходные биты)
hgfe dcba (результат)
bcde fgha (цикл. сдвиг влево на 1)
defg habc (цикл. сдвиг влево на 2)
fgha bcde (цикл. сдвиг влево на 2)
habc defg (цикл. сдвиг влево на 2)

Метод легко обобщается на 32 бита: N4. N4 вдвое быстрее N0.

Иногда оптимизация по размеру кода может быть важнее оптимизации по скорости. Команды циклического сдвига позволяют организовать из регистра нечто вроде стека. Возможно, это наиболее короткий код на x86 ассемблере, выполняющий поставленную задачу: N5. Увы, он и самый медленный: вдвое медленнее, чем N0. Именно по этой причине я не использовал команду loop. Если заменить её на пару инструкций sub/jnz [N6], то результат будет быстрее N0 на 5-30% (в зависимости от процессора).

Умный алгоритм

Обменивать по биту за итерацию довольно медленно. Используя поразрядные операции над регистрами, и благодаря тому, что размер регистра в битах есть степень двойки, можно обратить его содержимое за O(log N) операций: обменять соседние биты, затем пары бит, четвёрки и т.д. Итак, вариант на Си для 32-битного числа [H1] уже в 7-8 раз быстрее N0. Алгоритм можно представить и в виде, не привязанном к числу 32: H0, хотя такая реализация не отличается скоростью. H1, переписанный «в лоб» на ассемблере (H2; приводить не буду), работает на 0-6% быстрее H1 (конечно, здесь уже мало манёвра для оптимизатора…). Нетрудно заметить, что последний обмен пар байт можно сделать одной командой циклического сдвига: H3. Результат почти на треть быстрее H1.

Начиная с Intel 486, процессоры x86 поддерживают команду bswap, обращающую порядок байт в регистре. Команда была введена для ускорения обмена данными между little-endian процессорами x86 и big-endian процессорами других архитектур. В нашем случае ею можно заменить и циклический сдвиг и обмен соседних байт, что даёт код, вдвое более быстрый, чем H1. Однако, это ещё не всё. Комбинации вида mov edx, eax + shl edx, 1 можно заменить одной быстрой командой lea edx, [eax+eax], что позволяет выжать ещё 10-20%. Окончательный вариант: H4.

Можно ли ещё улучшить производительность? Посмотрим, что делает оптимизатор. Выражение вида

(a & mask) | (b & ~mask) он заменяет на

((a ^ b) & mask) ^ b,

видимо, с целью избежать двух длинных команд с 32-битной константой. Это интересно, но вряд ли полезно. Обратим внимание, что мы сдвигаем одно и то же число на равное количество бит в разные стороны. Эти два линейных сдвига можно заменить на один циклический (соответственно на 2, 4 и 8 бит) с последующим компенсационным сдвигом в обратную сторону на 7 бит: H5. Более того, сдвигаемые регистры можно чередовать (практика показала, что это полезно). H5 чуть-чуть быстрее H4 (примерно на 3%) и короче на одну инструкцию.

Итак, подведём итоги. Ниже приведены результаты замеров рассмотренных вариантов для двух процессоров (AMD Phenom II 955 [Deneb] 3.2ГГц и Intel Celeron U3400 [Arrandale] 1.07ГГц). На верхней диаграмме за 100% принято время N0, на нижней – H1.

Ради интереса можно сравнить эффективность выполнения этого кода на использовавшихся ЦП AMD vs Intel (усреднённая оценка (IPC на Arrandale) / (IPC на Deneb) * 100%). Зелёным выделены варианты эффективнее выполняемые на ЦП AMD, голубым – на ЦП Intel. Вариант «nothing» отражает эффективность выполнения вызова функции по указателю внутри тесного цикла.