Как реализовать atoi с помощью SIMD?
Я хотел бы попробовать написать реализацию atoi с использованием инструкций SIMD, которая будет включена в RapidJSON (библиотека для чтения / записи C++ JSON). В настоящее время он имеет некоторые оптимизации SSE2 и SSE4.2 в других местах.
Если это увеличение скорости, несколько atoi результаты могут быть сделаны параллельно. Первоначально строки поступают из буфера данных JSON, поэтому многопользовательская функция должна выполнять любые необходимые действия.
Алгоритм, который я придумал, следующий:
- Я могу инициализировать вектор длины N следующим образом: [10^N..10^1]
- Я конвертирую каждый символ в буфере в целое число и помещаю их в другой вектор.
- Я беру каждое число в векторе значащих цифр и умножаю его на соответствующее число в векторе чисел и суммирую результаты.
Я ориентируюсь на архитектуры x86 и x86-64.
Я знаю, что AVX2 поддерживает три операнда Fused Multiply-Add, так что я смогу выполнить Sum = Number * Значительная цифра + Sum.
Вот где я так далеко.
Мой алгоритм правильный? Есть ли способ лучше?
Существует ли эталонная реализация для atoi с использованием каких-либо SIMD-инструкций?
2 ответа
Решение
Алгоритм и его реализация завершены. Он завершен и (умеренно) протестирован (обновлен для менее постоянного использования памяти и допускает использование плюсовых символов).
Свойства этого кода следующие:
- Работает на
intа такжеuint, отMIN_INT=-2147483648вMAX_INT=2147483647и изMIN_UINT=0вMAX_UINT=4294967295 - Ведущий
'-'char обозначает отрицательное число (как разумное), ведущее'+'символ игнорируется - Ведущие нули (со знаком или без знака) игнорируются
- Переполнение игнорируется - большие числа просто переносятся
- Строки нулевой длины приводят к значению
0 = -0 - Неверные символы распознаются, и преобразование заканчивается на первом неверном символе
- По крайней мере 16 байтов после последнего ведущего нуля должны быть доступны, и возможные последствия для безопасности чтения после EOS должны быть оставлены вызывающей стороне
- Нужен только SSE4.2
Об этой реализации:
- Этот пример кода можно запустить с помощью GNU Assembler (
as) с помощью.intel_syntax noprefixв начале - Объем данных для констант составляет 64 байта (4*128 бит XMM), что соответствует одной строке кэша.
- Код занимает 46 инструкций с 51 микрооперацией и задержкой 64 цикла
- Один цикл для удаления начальных нулей, в противном случае никаких переходов, кроме обработки ошибок, так что...
- Временная сложность O(1)
Подход алгоритма:
- Pointer to number string is expected in ESI
- Check if first char is '-', then indicate if negative number in EDX (**A**)
- Check for leading zeros and EOS (**B**)
- Check string for valid digits and get strlen() of valid chars (**C**)
- Reverse string so that power of
10^0 is always at BYTE 15
10^1 is always at BYTE 14
10^2 is always at BYTE 13
10^3 is always at BYTE 12
10^4 is always at BYTE 11
...
and mask out all following chars (**D**)
- Subtract saturated '0' from each of the 16 possible chars (**1**)
- Take 16 consecutive byte-values and and split them to WORDs
in two XMM-registers (**2**)
P O N M L K J I | H G F E D C B A ->
H G F E | D C B A (XMM0)
P O N M | L K J I (XMM1)
- Multiply each WORD by its place-value modulo 10000 (1,10,100,1000)
(factors smaller then MAX_WORD, 4 factors per QWORD/halfXMM)
(**2**) so we can horizontally combine twice before another multiply.
The PMADDWD instruction can do this and the next step:
- Horizontally add adjacent WORDs to DWORDs (**3**)
H*1000+G*100 F*10+E*1 | D*1000+C*100 B*10+A*1 (XMM0)
P*1000+O*100 N*10+M*1 | L*1000+K*100 J*10+I*1 (XMM1)
- Horizontally add adjacent DWORDs from XMM0 and XMM1 to XMM0 (**4**)
xmmDst[31-0] = xmm0[63-32] + xmm0[31-0]
xmmDst[63-32] = xmm0[127-96] + xmm0[95-64]
xmmDst[95-64] = xmm1[63-32] + xmm1[31-0]
xmmDst[127-96] = xmm1[127-96] + xmm1[95-64]
- Values in XMM0 are multiplied with the factors (**5**)
P*1000+O*100+N*10+M*1 (DWORD factor 1000000000000 = too big for DWORD, but possibly useful for QWORD number strings)
L*1000+K*100+J*10+I*1 (DWORD factor 100000000)
H*1000+G*100+F*10+E*1 (DWORD factor 10000)
D*1000+C*100+B*10+A*1 (DWORD factor 1)
- The last step is adding these four DWORDs together with 2*PHADDD emulated by 2*(PSHUFD+PADDD)
- xmm0[31-0] = xmm0[63-32] + xmm0[31-0] (**6**)
xmm0[63-32] = xmm0[127-96] + xmm0[95-64]
(the upper QWORD contains the same and is ignored)
- xmm0[31-0] = xmm0[63-32] + xmm0[31-0] (**7**)
- If the number is negative (indicated in EDX by 000...0=pos or 111...1=neg), negate it(**8**)
И пример реализации в GNU Assembler с синтаксисом Intel:
.intel_syntax noprefix
.data
.align 64
ddqDigitRange: .byte '0','9',0,0,0,0,0,0,0,0,0,0,0,0,0,0
ddqShuffleMask:.byte 15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0
ddqFactor1: .word 1,10,100,1000, 1,10,100,1000
ddqFactor2: .long 1,10000,100000000,0
.text
_start:
mov esi, lpInputNumberString
/* (**A**) indicate negative number in EDX */
mov eax, -1
xor ecx, ecx
xor edx, edx
mov bl, byte ptr [esi]
cmp bl, '-'
cmove edx, eax
cmp bl, '+'
cmove ecx, eax
sub esi, edx
sub esi, ecx
/* (**B**)remove leading zeros */
xor eax,eax /* return value ZERO */
remove_leading_zeros:
inc esi
cmp byte ptr [esi-1], '0' /* skip leading zeros */
je remove_leading_zeros
cmp byte ptr [esi-1], 0 /* catch empty string/number */
je FINISH
dec esi
/* check for valid digit-chars and invert from front to back */
pxor xmm2, xmm2
movdqa xmm0, xmmword ptr [ddqDigitRange]
movdqu xmm1, xmmword ptr [esi]
pcmpistri xmm0, xmm1, 0b00010100 /* (**C**) iim8=Unsigned bytes, Ranges, Negative Polarity(-), returns strlen() in ECX */
jo FINISH /* if first char is invalid return 0 - prevent processing empty string - 0 is still in EAX */
mov al , '0' /* value to subtract from chars */
sub ecx, 16 /* len-16=negative to zero for shuffle mask */
movd xmm0, ecx
pshufb xmm0, xmm2 /* broadcast CL to all 16 BYTEs */
paddb xmm0, xmmword ptr [ddqShuffleMask] /* Generate permute mask for PSHUFB - all bytes < 0 have highest bit set means place gets zeroed */
pshufb xmm1, xmm0 /* (**D**) permute - now from highest to lowest BYTE are factors 10^0, 10^1, 10^2, ... */
movd xmm0, eax /* AL='0' from above */
pshufb xmm0, xmm2 /* broadcast AL to XMM0 */
psubusb xmm1, xmm0 /* (**1**) */
movdqa xmm0, xmm1
punpcklbw xmm0, xmm2 /* (**2**) */
punpckhbw xmm1, xmm2
pmaddwd xmm0, xmmword ptr [ddqFactor1] /* (**3**) */
pmaddwd xmm1, xmmword ptr [ddqFactor1]
phaddd xmm0, xmm1 /* (**4**) */
pmulld xmm0, xmmword ptr [ddqFactor2] /* (**5**) */
pshufd xmm1, xmm0, 0b11101110 /* (**6**) */
paddd xmm0, xmm1
pshufd xmm1, xmm0, 0b01010101 /* (**7**) */
paddd xmm0, xmm1
movd eax, xmm0
/* negate if negative number */
add eax, edx /* (**8**) */
xor eax, edx
FINISH:
/* EAX is return (u)int value */
Результат анализа пропускной способности Intel-IACA для 32-разрядной версии Haswell:
Throughput Analysis Report
--------------------------
Block Throughput: 16.10 Cycles Throughput Bottleneck: InterIteration
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 9.5 0.0 | 10.0 | 4.5 4.5 | 4.5 4.5 | 0.0 | 11.1 | 11.4 | 0.0 |
---------------------------------------------------------------------------------------
N - port number or number of cycles resource conflict caused delay, DV - Divider pipe (on port 0)
D - Data fetch pipe (on ports 2 and 3), CP - on a critical path
F - Macro Fusion with the previous instruction occurred
* - instruction micro-ops not bound to a port
^ - Micro Fusion happened
# - ESP Tracking sync uop was issued
@ - SSE instruction followed an AVX256 instruction, dozens of cycles penalty is expected
! - instruction not supported, was not accounted in Analysis
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 0* | | | | | | | | | | xor eax, eax
| 0* | | | | | | | | | | xor ecx, ecx
| 0* | | | | | | | | | | xor edx, edx
| 1 | | 0.1 | | | | | 0.9 | | | dec eax
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | CP | mov bl, byte ptr [esi]
| 1 | | | | | | | 1.0 | | CP | cmp bl, 0x2d
| 2 | 0.1 | 0.2 | | | | | 1.8 | | CP | cmovz edx, eax
| 1 | 0.1 | 0.5 | | | | | 0.4 | | CP | cmp bl, 0x2b
| 2 | 0.5 | 0.2 | | | | | 1.2 | | CP | cmovz ecx, eax
| 1 | 0.2 | 0.5 | | | | | 0.2 | | CP | sub esi, edx
| 1 | 0.2 | 0.5 | | | | | 0.3 | | CP | sub esi, ecx
| 0* | | | | | | | | | | xor eax, eax
| 1 | 0.3 | 0.1 | | | | | 0.6 | | CP | inc esi
| 2^ | 0.3 | | 0.5 0.5 | 0.5 0.5 | | | 0.6 | | | cmp byte ptr [esi-0x1], 0x30
| 0F | | | | | | | | | | jz 0xfffffffb
| 2^ | 0.6 | | 0.5 0.5 | 0.5 0.5 | | | 0.4 | | | cmp byte ptr [esi-0x1], 0x0
| 0F | | | | | | | | | | jz 0x8b
| 1 | 0.1 | 0.9 | | | | | | | CP | dec esi
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | | movdqa xmm0, xmmword ptr [0x80492f0]
| 1 | | | 0.5 0.5 | 0.5 0.5 | | | | | CP | movdqu xmm1, xmmword ptr [esi]
| 0* | | | | | | | | | | pxor xmm2, xmm2
| 3 | 2.0 | 1.0 | | | | | | | CP | pcmpistri xmm0, xmm1, 0x14
| 1 | | | | | | | 1.0 | | | jo 0x6e
| 1 | | 0.4 | | | | 0.1 | 0.5 | | | mov al, 0x30
| 1 | 0.1 | 0.5 | | | | 0.1 | 0.3 | | CP | sub ecx, 0x10
| 1 | | | | | | 1.0 | | | CP | movd xmm0, ecx
| 1 | | | | | | 1.0 | | | CP | pshufb xmm0, xmm2
| 2^ | | 1.0 | 0.5 0.5 | 0.5 0.5 | | | | | CP | paddb xmm0, xmmword ptr [0x80492c0]
| 1 | | | | | | 1.0 | | | CP | pshufb xmm1, xmm0
| 1 | | | | | | 1.0 | | | | movd xmm0, eax
| 1 | | | | | | 1.0 | | | | pshufb xmm0, xmm2
| 1 | | 1.0 | | | | | | | CP | psubusb xmm1, xmm0
| 0* | | | | | | | | | CP | movdqa xmm0, xmm1
| 1 | | | | | | 1.0 | | | CP | punpcklbw xmm0, xmm2
| 1 | | | | | | 1.0 | | | | punpckhbw xmm1, xmm2
| 2^ | 1.0 | | 0.5 0.5 | 0.5 0.5 | | | | | CP | pmaddwd xmm0, xmmword ptr [0x80492d0]
| 2^ | 1.0 | | 0.5 0.5 | 0.5 0.5 | | | | | | pmaddwd xmm1, xmmword ptr [0x80492d0]
| 3 | | 1.0 | | | | 2.0 | | | CP | phaddd xmm0, xmm1
| 3^ | 2.0 | | 0.5 0.5 | 0.5 0.5 | | | | | CP | pmulld xmm0, xmmword ptr [0x80492e0]
| 1 | | | | | | 1.0 | | | CP | pshufd xmm1, xmm0, 0xee
| 1 | | 1.0 | | | | | | | CP | paddd xmm0, xmm1
| 1 | | | | | | 1.0 | | | CP | pshufd xmm1, xmm0, 0x55
| 1 | | 1.0 | | | | | | | CP | paddd xmm0, xmm1
| 1 | 1.0 | | | | | | | | CP | movd eax, xmm0
| 1 | | | | | | | 1.0 | | CP | add eax, edx
| 1 | | | | | | | 1.0 | | CP | xor eax, edx
Total Num Of Uops: 51
Результат анализа задержки Intel-IACA для 32-битного Haswell:
Latency Analysis Report
---------------------------
Latency: 64 Cycles
N - port number or number of cycles resource conflict caused delay, DV - Divider pipe (on port 0)
D - Data fetch pipe (on ports 2 and 3), CP - on a critical path
F - Macro Fusion with the previous instruction occurred
* - instruction micro-ops not bound to a port
^ - Micro Fusion happened
# - ESP Tracking sync uop was issued
@ - Intel(R) AVX to Intel(R) SSE code switch, dozens of cycles penalty is expected
! - instruction not supported, was not accounted in Analysis
The Resource delay is counted since all the sources of the instructions are ready
and until the needed resource becomes available
| Inst | Resource Delay In Cycles | |
| Num | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | FE | |
-------------------------------------------------------------------------
| 0 | | | | | | | | | | | xor eax, eax
| 1 | | | | | | | | | | | xor ecx, ecx
| 2 | | | | | | | | | | | xor edx, edx
| 3 | | | | | | | | | | | dec eax
| 4 | | | | | | | | | 1 | CP | mov bl, byte ptr [esi]
| 5 | | | | | | | | | | CP | cmp bl, 0x2d
| 6 | | | | | | | | | | CP | cmovz edx, eax
| 7 | | | | | | | | | | CP | cmp bl, 0x2b
| 8 | | | | | | | 1 | | | CP | cmovz ecx, eax
| 9 | | | | | | | | | | CP | sub esi, edx
| 10 | | | | | | | | | | CP | sub esi, ecx
| 11 | | | | | | | | | 3 | | xor eax, eax
| 12 | | | | | | | | | | CP | inc esi
| 13 | | | | | | | | | | | cmp byte ptr [esi-0x1], 0x30
| 14 | | | | | | | | | | | jz 0xfffffffb
| 15 | | | | | | | | | | | cmp byte ptr [esi-0x1], 0x0
| 16 | | | | | | | | | | | jz 0x8b
| 17 | | | | | | | | | | CP | dec esi
| 18 | | | | | | | | | 4 | | movdqa xmm0, xmmword ptr [0x80492f0]
| 19 | | | | | | | | | | CP | movdqu xmm1, xmmword ptr [esi]
| 20 | | | | | | | | | 5 | | pxor xmm2, xmm2
| 21 | | | | | | | | | | CP | pcmpistri xmm0, xmm1, 0x14
| 22 | | | | | | | | | | | jo 0x6e
| 23 | | | | | | | | | 6 | | mov al, 0x30
| 24 | | | | | | | | | | CP | sub ecx, 0x10
| 25 | | | | | | | | | | CP | movd xmm0, ecx
| 26 | | | | | | | | | | CP | pshufb xmm0, xmm2
| 27 | | | | | | | | | 7 | CP | paddb xmm0, xmmword ptr [0x80492c0]
| 28 | | | | | | | | | | CP | pshufb xmm1, xmm0
| 29 | | | | | | 1 | | | | | movd xmm0, eax
| 30 | | | | | | 1 | | | | | pshufb xmm0, xmm2
| 31 | | | | | | | | | | CP | psubusb xmm1, xmm0
| 32 | | | | | | | | | | CP | movdqa xmm0, xmm1
| 33 | | | | | | | | | | CP | punpcklbw xmm0, xmm2
| 34 | | | | | | | | | | | punpckhbw xmm1, xmm2
| 35 | | | | | | | | | 9 | CP | pmaddwd xmm0, xmmword ptr [0x80492d0]
| 36 | | | | | | | | | 9 | | pmaddwd xmm1, xmmword ptr [0x80492d0]
| 37 | | | | | | | | | | CP | phaddd xmm0, xmm1
| 38 | | | | | | | | | 10 | CP | pmulld xmm0, xmmword ptr [0x80492e0]
| 39 | | | | | | | | | | CP | pshufd xmm1, xmm0, 0xee
| 40 | | | | | | | | | | CP | paddd xmm0, xmm1
| 41 | | | | | | | | | | CP | pshufd xmm1, xmm0, 0x55
| 42 | | | | | | | | | | CP | paddd xmm0, xmm1
| 43 | | | | | | | | | | CP | movd eax, xmm0
| 44 | | | | | | | | | | CP | add eax, edx
| 45 | | | | | | | | | | CP | xor eax, edx
Resource Conflict on Critical Paths:
-----------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
-----------------------------------------------------------------
| Cycles | 0 0 | 0 | 0 0 | 0 0 | 0 | 0 | 1 | 0 |
-----------------------------------------------------------------
List Of Delays On Critical Paths
-------------------------------
6 --> 8 1 Cycles Delay On Port6
Альтернативная обработка, предложенная в комментариях Питером Кордесом, заменяет последние два add+xor инструкции от imul, Эта концентрация OpCodes, вероятно, будет выше. К сожалению, IACA не поддерживает эту инструкцию и выдает ! - instruction not supported, was not accounted in Analysis комментарий. Тем не менее, хотя мне нравится сокращение OpCodes и сокращение с (2 моп, 2 с задержкой) до (1 моп, 3 с задержкой - "худшая задержка, но все еще одна операция на AMD"), я предпочитаю оставить это разработчику, который способ выбрать. Я не проверял, достаточен ли следующий код для анализа любого числа. Это только упомянуто для полноты, и модификации кода в других частях могут быть необходимы (особенно обработка положительных чисел).
Альтернативой может быть замена двух последних строк на:
...
/* negate if negative number */
imul eax, edx
FINISH:
/* EAX is return (u)int value */
Я хотел бы подойти к этой проблеме следующим образом:
- Инициализируйте аккумулятор до 0.
- Загрузите следующие четыре символа строки в регистр SSE.
- Вычтите значение
'0'от каждого персонажа. - Найти первое значение в векторе, значение которого без знака больше
9, - Если значение найдено, сместите компоненты вектора вправо, чтобы значение, найденное на предыдущем шаге, было просто сдвинуто.
- Загрузите вектор, содержащий степени десяти (
1000,100,10,1) и умножить с этим. - Вычислить сумму всех записей в векторе.
- Умножьте аккумулятор на соответствующее значение (в зависимости от количества смен на шаге 5) и добавьте вектор. Вы можете использовать инструкцию FMA для этого, но я не знаю, существует ли такая инструкция для целых чисел.
- Если значение не превышает
9был найден в шаге четыре, перейдите к шагу 2. - Верните аккумулятор.
Вы можете упростить алгоритм, обнулив все записи, начиная с неправильного в шаге 5, вместо сдвига, а затем разделив на подходящую степень десяти в конце.
Пожалуйста, имейте в виду, что этот алгоритм читает после конца строки и поэтому не является заменой atoi,