Как реализовать atoi с помощью SIMD?

Я хотел бы попытаться написать реализацию atoi с помощью инструкций SIMD, которые будут включены в RapidJSON (библиотека чтения / записи C ++ JSON). В настоящее время он имеет некоторые оптимизации SSE2 и SSE4.2 в других местах.

Если это коэффициент усиления скорости, несколько результатов atoi могут быть выполнены параллельно. Строки изначально поступают из буфера данных JSON, поэтому функция multi-atoi должна выполнять любые требуемые swizzling.

Алгоритм, который я придумал, таков:

Я могу инициализировать вектор длины N следующим образом: [10 ^ N..10 ^ 1]
Я конвертирую каждый символ в буфер в целое число и помещаю его в другой вектор.
Я беру каждое число в векторе значимых цифр и умножаю его на число совпадений в векторе чисел и суммируем результаты.

Я нацелен на архитектуры x86 и x86-64.

Я знаю, что AVX2 поддерживает три операнда Fused Multiply-Add, поэтому я смогу выполнить Sum = Number * Significant Digit + Sum.
Вот где я до сих пор.
Правильно ли мой алгоритм? Есть ли способ лучше?
Существует ли эталонная реализация для atoi с использованием набора инструкций SIMD?

Теперь алгоритм и его реализация завершены. Он завершен и (умеренно) протестирован ( Обновлен для менее постоянного использования памяти и переносится плюс-char ).

Свойства этого кода следующие:

Работает для int и uint , от MIN_INT=-2147483648 до MAX_INT=2147483647 и от MIN_UINT=0 до MAX_UINT=4294967295
Ведущий символ '-' указывает отрицательное число (как разумное), символ '+' игнорируется
Ведущие нули (с символом или без знака) игнорируются
Переполнение игнорируется – большие числа просто обходятся
Строки нулевой длины приводят к значению 0 = -0
Недействительные символы распознаются, и преобразование завершается с первым недопустимым символом
Не менее 16 байт после последнего начального нуля должны быть доступны и возможные последствия для безопасности чтения после EOS оставлены вызывающему абоненту
Требуется только SSE4.2

Об этой реализации:

Этот пример кода может быть запущен с помощью GNU Assembler ( as ) с использованием .intel_syntax noprefix в начале
Объем данных для констант – 64 байта (4 * 128 бит XMM), равный одной строке кэша.
Кодовый след – 46 инструкций с 51 микро-операцией и задержкой в 64 цикла
Один цикл для удаления начальных нhive, в противном случае нет переходов, кроме обработки ошибок, поэтому …
Сложность времени – 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**) PONMLKJI | HGFEDCBA -> HGFE | DCBA (XMM0) PONM | LKJI (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 ! - instruction not supported, was not accounted in Analysis комментарии ! - instruction not supported, was not accounted in Analysis . Тем не менее, хотя мне нравится сокращение OpCodes и сокращение от (2uops, 2c latency) до (1 uop, 3c latency - «хуже латентности, но все еще один m-op на 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, но я не знаю, существует ли такая инструкция для целых чисел.
Если на шаге 4 не было найдено значения больше 9 , перейдите к шагу 2.
Верните аккумулятор.

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

Имейте в виду, что этот алгоритм читает за конец строки и, следовательно, не является заменой для atoi .

Самый быстрый способ сделать горизонтальную векторную сумму float на x86

Почему оставить «mov esp, ebp» в сборке x86?

Можно ли сообщить предсказателю ветви, насколько вероятно, что он должен следовать за веткой?

Как написать привет мир в ассемблере под Windows?

Какие целые операции с дополнением 2 можно использовать без обнуления высоких бит в входах, если требуется только низкая часть результата?

Почему этот код SSE в 6 раз медленнее без VZEROUPPER на Skylake?

Что означает «DS: » означает в сборке?

Почему x86 маленький endian?

В чем смысл «не временных» обращений к памяти в x86

Почему вы не можете установить указатель инструкции напрямую?

Наблюдение за устаревшей инструкцией по x86 с самомодифицируемым кодом