Под машинным способом я главным образом подразумеваю использование битовых операций
Если Вы под этим завуалировали "как получить быстрое решение?", то мой вариант: слепите 4 бита сравнений координат с нулем и проверьте в switch'е
все 16 комбинаций.
Можете, пожалуйста, расписать подробно, так как я не уловил суть подхода?
int Quadrant( int x, int y ){
unsigned a;
//
a = ( x > 0 ) ? 1 : 0;
a |= ( x < 0 ) ? 2 : 0;
a |= ( y > 0 ) ? 4 : 0;
a |= ( y < 0 ) ? 8 : 0;
//
switch( a ){
case 5: return 1;
case 6: return 2;
case 9: return 4; // это не ошибка в расчетах квадранта,
case 10: return 3; // это Декарт виноват :)
default: return -1;
}
}
int Quadrant( int x, int y, int cx, int cy ){
return Quadrant( x - cx, y - cy );
}
Суть в том, чтобы не рассчитывать признаки повторно.
--
Здесь 4 + 4
операции сравнения, 3
- сложения, т.е. для получения результата требуется от 8 до 11
операций. Еще можно попробовать заменить операции 4 + 3
на 1 - SIMD-сравнение
, получится где-то от 2 до 5
операций, но скорее всего, это будет медленнее решения со SO
, там, насколько я понял, всего 4
операции.
--
Если вместо switch
'а задействовать бинарный поиск, то можно сэкономить на одной операции.
--
Можно, кстати, вообще убрать switch
и использовать "a"
в качестве результата, но, конечно, номера будут непопорядку, как и на SO
.
--
Зато быстро: грубо говоря - за 1 SIMD-сравнение
или за 4 + 3
обычных (сравнения + сложения) ))
--
Немного доработал алгоритм, чтобы получить правильный порядок четвертей. Это потребовало еще 1 дополнительную операцию:
int Quadrant( int x, int y ){
unsigned a;
//
a = ( x > 0 ) ? 1 : 0;
a |= ( x < 0 ) ? 2 : 0;
a |= ( y > 0 ) ? 4 : 0;
//a |= ( y < 0 ) ? 8 : 0;
if( y < 0 ){
a ^= 3; // инверсия признаков в I и II полуплоскостях
a |= 8;
}
//
switch( a ){
case 5: return 1;
case 6: return 2;
case 9: return 3; // здесь уже
case 10: return 4; // все верно, по Декарту! :)
default: return -1;
}
}
Если задействовать SIMD
, то придется сначала выполнить SIMD-сравнение
с 0
и за'xor'ить первые два юнита, если 4 юнит не ноль
. Т.е. выходит три операции вместо одной из предыдущего варианта.
Note: при задействовании PCMPGTD (SSE2), такую скорость можно получить только в потоковом режиме, т.е. на наборе точек > 2 штук, т.к., насколько я знаю, нет в этой спецификации такой SIMD-операции, которая может применить к каждому юниту произвольную операцию сравнения (но не факт :), могли уже ввести, на более современном железе).
Тут, к примеру, для первого варианта алгоритма придется в первом SIMD-векторе вычислить { x1 > 0, x2 > 0, y1 > 0, y2 > 0 }, а на втором - { 0 > x1, 0 > x2, 0 > y1, 0 > y2 }. Комбинация этих двух векторов и дает результат для двух точек.
Выходит, 2 операции на 2 точки, а не 1 на 1. Для второго алгоритма все аналогично: 6 операций на 2 точки.
--
Немного доработал алгоритм, чтобы получить еще и правильный номер четвертей:
int Quadrant( int x, int y ){
int a, b;
//
a = ( x > 0 ) ? 1 : 0;
a |= ( x < 0 ) ? 2 : 0;
b = ( y > 0 ) ? 1 : 0;
if( y < 0 ){
a ^= 3;
b |= 2;
}
return a + b;
}
Здесь 4 сравнения 3 сложения 1 - xor, т.е. 8 простых операций, без проверки на равенство 0: если a
или b
равны 0, то точка лежит на базисе, поэтому результат будет неверным.
Может показаться, что SSE
здесь так же может ускорить расчет: 4 сравнения и 2 сложения
заменяются на 1 SIMD-сравнение
, оставляя 1 сложение
и 1 xor
, но увы - только если в SSE
будет введена команда 4-битного сложения, где бит соответствует одному юниту вектора, а такого в спецификации я не наблюдал.
Подвожу свои итоги:
- Самая быстрая реализация первого алгоритма дает результат за одну операцию, что условно в 3 раза быстрее, чем на
SO
. Размер результата 128 бит. Если надо меньше, то можно его сжать до 32 бит, если задействовать еще 1-2 операции упаковки
(сейчас точное число операций упаковки не скажу, надо смотреть документацию, есть ли вариант упаковки из 32-битного сразу в 8-битный юнит).
- Самая быстрая реализация второго алгоритма дает результат за три операции, что на 1 условную операцию меньше, чем на
SO
. Кроме того, этот алгоритм дает правильный топологический порядок четвертей, что позволяет сравнивать и упорядочивать точки по квадрантам, минуя дополнительные трансформации, т.е. еще быстрее.
- При использовании типа
float
, оба алгоритма сохраняют скорости.
- При использовании типа
double
или long long
, число операций для обоих алгоритмов увеличивается на 1
(на SSE <= 4
).
- На
SSE5
, при использовании типа double
или long long
, оба алгоритма сохраняют скорости.
Для справки: скорость PCMPEQD/PCMPEQW/PCMPEQB (32/16/8-битное SIMD-сравнение двух векторов) составляет в среднем от 4 до 6 нс. (<10тыс. итераций), от 25 до 48 нс. (<10млн. итераций). на 1Ггц (тестировал в конце прошлого года на своем кластере: 8x E5450 @ 3.0GHz, SSE4.1). Т.е., в сравнении с 0 (для PCMPGTD), могу предположить скорость не менее указанной, а то и выше, ибо в этом случае производится чтение только одного вектора.
но могли бы Вы чётко расписать/пометить/представить, какова именно "Самая быстрая реализация первого алгоритма", которая "дает результат за одну операцию"
Я и не писал быстрые варианты, я написал 3 общих варианта без векторизации. Векторизация может быть реализована по разному, смотрите примечание Note. К примеру, реализацию первого алгоритма для двух точек может быть такой (MSVS):
static const __m128i zero = _mm_set_epi32( 0, 0, 0, 0 );
__m128i pair = _mm_set_epi32( y2, y1, x2, x1 );
__m128i p1 = _mm_cmpgt_epi32( pair, zero );
__m128i p2 = _mm_cmpgt_epi32( zero, pair );
Результат для первой точки находится в векторе
{ p2.m128i_i32[2], p2.m128i_i32[0], p1.m128i_i32[2], p1.m128i_i32[0] },
а для второй, соответственно - в векторе
{ p2.m128i_i32[3], p2.m128i_i32[1], p1.m128i_i32[3], p1.m128i_i32[1] },
Как Вы их дальше будете использовать, зависит от Вашего потокового алгоритма: сохраните результаты в дополнительном векторе, или определите четверть сразу - дальнейшая скорость зависит от Вашей реализации.
Что ж, тогда выходит, что Вы меня до смерти запугали разными очень низкоуровневыми и совершенно непонятными мне вещами, но при этом не предлагаете ни одного решения, которое стоило бы сравнить с уже имеющимися другими решениями
Ну что ж, сочувствую. Понимаю, что ответ скорее другого уровня. Тем не менее, пример я Вам дал, сожалею, если Вы не понимаете как его можно использовать.
--
Если хотите просто измерить скорость, разделите исходную выборку точек по парам и вычислите для каждой эти три операции. На самом деле, ничего сложного, Вам просто нужно поставить его в те же условия, что и остальные алгоритмы, которые Вы тестируете.
--
Кстати, было бы хорошо, если бы Вы старались использовать в своих несомненно весомых оценках русскую терминологию. Особенно, если Вы хотите, чтобы я мог Вам чем-то помочь.
Ну смотрите, вот прототип static inline int quadrant(double x, double y, double centerX, double centerY);
@Stanislaw Pankevich, Еще раз Вам повторяю Note:
при задействовании PCMPGTD (SSE2), такую скорость можно получить только в потоковом режиме, т.е. на наборе точек > 2 штук
а Вы мне предлагаете 1 точку.
показать работу с потоком, главное, чтобы решение было законченным и работало
Ждите, когда я найду для этого время, либо программиста, который прочитает мои инструкции и реализует в своем решении. Но Вы можете просто рассчитать, сколько времени работают остальные. ))
Реализация для обоих алгоритмов (MSVS):
#include <windows.h>
#include <stdio.h>
#include <tchar.h>
#include <time.h>
#include <emmintrin.h>
#include <vector>
//
static const __m128i zero = _mm_set_epi32( 0, 0, 0, 0 );
static const __m128i xor = _mm_set_epi32( 0, 0, -1, -1 );
//
inline void Quadrant1( int x1, int y1, int x2, int y2, __m128i*p1, __m128i*p2 ){
__m128i pair = _mm_set_epi32( y2, x2, y1, x1 );
__m128i gz = _mm_cmpgt_epi32( pair, zero );
__m128i lz = _mm_cmpgt_epi32( zero, pair );
*p1 = _mm_unpacklo_epi32( gz, lz );
*p2 = _mm_unpackhi_epi32( gz, lz );
}
//
inline void Quadrant2( int x1, int y1, int x2, int y2, __m128i*p1, __m128i*p2 ){
__m128i pair = _mm_set_epi32( y2, x2, y1, x1 );
__m128i gz = _mm_cmpgt_epi32( pair, zero );
__m128i lz = _mm_cmpgt_epi32( zero, pair );
*p1 = _mm_unpacklo_epi32( gz, lz );
*p2 = _mm_unpackhi_epi32( gz, lz );
//
if( lz.m128i_i32[1] ){
*p1 = _mm_xor_si128( *p1, xor );
}
if( lz.m128i_i32[3] ){
*p2 = _mm_xor_si128( *p2, xor );
}
}
//
int _tmain( int argc, LPCTSTR argv[] ){
static const DWORD Measure = 5000;
//
std::vector< POINT >Points;
std::vector< __m128i >Q;
ULONGLONG Start, T1, T2;
DWORD_PTR Mask;
double ms;
SIZE_T ic, p1, p2, Count;
//
if( ( argc < 2 ) || ( 1 != _stscanf_s( argv[1], TEXT("%Iu"), &Count ) ) ){
return 0;
}
_tprintf( TEXT("count=%Iu\r\n" ), Count );
Mask = ::SetThreadAffinityMask( ::GetCurrentThread(), 1 );
_tprintf( TEXT("measure...\r\n" ) );
Start = __rdtsc();
Sleep( Measure );
ms = double( __rdtsc() - Start ) / Measure;
_tprintf( TEXT("Freq=%lf GHz\r\n" ), ms / 1000000 );
//
srand( (unsigned)time( nullptr ) );
Points.resize( Count );
Q.resize( Count );
_tprintf( TEXT("filling...\r\n" ) );
for( ic = Count ; ic-- ; ){
Points[ic].x = rand();
Points[ic].y = rand();
}
_tprintf( TEXT("calc...\r\n" ) );
Start = __rdtsc();
for( p1 = 0, p2 = 1, ic = Count / 2 ; ic-- ; ++p1, ++p2 ){
Quadrant1( Points[p1].x, Points[p1].y, Points[p2].x, Points[p2].y, &Q[p1], &Q[p2] );
}
T1 = __rdtsc() - Start;
Start = __rdtsc();
for( p1 = 0, p2 = 1, ic = Count / 2 ; ic-- ; ++p1, ++p2 ){
Quadrant2( Points[p1].x, Points[p1].y, Points[p2].x, Points[p2].y, &Q[p1], &Q[p2] );
}
T2 = __rdtsc() - Start;
_tprintf(
TEXT("SIMD time Quadrant1 (Freq/1GHz): %lf/%lf ms.\r\n"), T1 / ms, T1 / 1000000.0
);
_tprintf(
TEXT("SIMD time Quadrant2 (Freq/1GHz): %lf/%lf ms.\r\n"), T2 / ms, T2 / 1000000.0
);
::SetThreadAffinityMask( ::GetCurrentThread(), Mask );
return 0;
}
Вывод:
count=100000000
measure...
Freq=3.002832 GHz
filling...
calc...
SIMD time Quadrant1 (Freq/1GHz): 534.183610/1604.063475 ms.
SIMD time Quadrant2 (Freq/1GHz): 596.293807/1790.569944 ms.
Press any key to continue . . .
Это расчеты на одном ядре. Первое значение - время, потраченное на расчет на частоте ~3GHz, второй - на условном 1GHz. Запускал много раз, расхождение в среднем не более 10мс. "Прогрев" ядра в начале 5 сек. Число точек то же самое, что и у @avp: 100млн. ОЗУ достаточно медленное (8x 1GB PC2-5300), так что скорость на таком количестве точек еще может увеличиться, если запустить на быстрой памяти: каждый квадрант записывается в массив.
DASM Quadrant1 (контекст):
mov dword ptr [rsp+3Ch],eax
mov dword ptr [rsp+38h],ecx
mov dword ptr [rsp+34h],edx
mov dword ptr [rsp+30h],r8d
movdqa xmm0,xmmword ptr [rsp+30h]
movdqa xmm2,xmm0
pcmpgtd xmm2,xmm4
movdqa xmm1,xmm4
pcmpgtd xmm1,xmm0
movdqa xmm0,xmm2
punpckldq xmm0,xmm1
movdqa xmmword ptr [r11],xmm0
add r11,10h
punpckhdq xmm2,xmm1
movdqa xmmword ptr [r11],xmm2
DASM Quadrant2 (контекст):
mov dword ptr [rsp+3Ch],ecx
mov dword ptr [rsp+38h],edx
mov dword ptr [rsp+34h],eax
mov dword ptr [rsp+30h],r8d
movdqa xmm0,xmmword ptr [rsp+30h]
movdqa xmm2,xmm0
pcmpgtd xmm2,xmm4
movdqa xmm1,xmm4
pcmpgtd xmm1,xmm0
movdqa xmmword ptr [Points],xmm1
movdqa xmm0,xmm2
punpckldq xmm0,xmm1
movdqa xmmword ptr [r9],xmm0
punpckhdq xmm2,xmm1
movdqa xmmword ptr [r9+10h],xmm2
cmp dword ptr [rsp+44h],0
je wmain+2DBh (13FA7133Bh)
movdqa xmm0,xmmword ptr [r9]
pxor xmm0,xmm3
movdqa xmmword ptr [r9],xmm0
cmp dword ptr [rsp+4Ch],0
je wmain+2F2h (13FA71352h)
movdqa xmm0,xmmword ptr [r9+10h]
pxor xmm0,xmm3
movdqa xmmword ptr [r9+10h],xmm0
Для произвольного центра результат такой:
SIMD time Quadrant1 (Freq/1GHz): 559.389193/1678.360437 ms.
SIMD time Quadrant2 (Freq/1GHz): 666.485334/1999.685781 ms.
Для произвольного центра, который учитывается через SIMD(PSUBD) результат такой:
SIMD time Quadrant1 (Freq/1GHz): 553.127711/1661.376060 ms.
SIMD time Quadrant2 (Freq/1GHz): 614.841356/1846.739349 ms.
На самом деле, быстрее всего расчет с произвольным центром будет в гомогенном потоке инструкций, т.е. предварительным циклом. Но в таком случае, эта "добавка" к расчетам будет ощущаться на уровне той же стартовой погрешности ))
--
То же самое касается и _mm_xor_si128
, который xor'ит
квадрант во второй версии алгоритма. Т.е. в конечном счете, оба алгоритма могут работать примерно на одной скорости, если их грамотно специализировать под задачу. Собственно xor вообще можно вызвать обычный, 64-битный
(на x64
).
--
Ну и учтите, что я вызываю еще две лишних команды для первой версии (punpckhdq/punpckldq
). Это нужно только для того, чтобы разделить результаты и записать их в память. Я мог бы просто записать в память пары квадрантов, это будет еще быстрее.
--
Основное время съедает, конечно, работа с памятью.
Немного доработал реализацию, чтобы упаковать 128-битные квадранты в 32-битные слова:
static const __m128i shfl = _mm_set_epi8(
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 12, 8, 4, 0
);
...
inline void Quadrant1( int x1, int y1, int x2, int y2, const __m128i&c, int*q1, int*q2 ){
__m128i pair = _mm_sub_epi32( _mm_set_epi32( y2, x2, y1, x1 ), c );
__m128i gz = _mm_cmpgt_epi32( pair, zero );
__m128i lz = _mm_cmpgt_epi32( zero, pair );
__m128i p1 = _mm_unpacklo_epi32( gz, lz );
__m128i p2 = _mm_unpackhi_epi32( gz, lz );
*q1 = _mm_shuffle_epi8( p1, shfl ).m128i_i32[0];
*q2 = _mm_shuffle_epi8( p2, shfl ).m128i_i32[0];
}
//
inline void Quadrant2( int x1, int y1, int x2, int y2, const __m128i&c, int*q1, int*q2 ){
__m128i pair = _mm_sub_epi32( _mm_set_epi32( y2, x2, y1, x1 ), c );
__m128i gz = _mm_cmpgt_epi32( pair, zero );
__m128i lz = _mm_cmpgt_epi32( zero, pair );
__m128i p1 = _mm_unpacklo_epi32( gz, lz );
__m128i p2 = _mm_unpackhi_epi32( gz, lz );
//
if( lz.m128i_i32[1] ){
p1 = _mm_xor_si128( p1, xor );
}
if( lz.m128i_i32[3] ){
p2 = _mm_xor_si128( p2, xor );
}
*q1 = _mm_shuffle_epi8( p1, shfl ).m128i_i32[0];
*q2 = _mm_shuffle_epi8( p2, shfl ).m128i_i32[0];
}
На этом этапе каждый квадрант 4-битный, с размером бита в 1 байт, поэтому, если нужен switch
, то он должен быть уже таким для первого варианта:
switch( a ){
case 0x00ff00ff: return 1;
case 0x00ffff00: return 2;
case 0xff0000ff: return 4;
case 0xff00ff00: return 3;
default: return -1;
}
и для второго:
switch( a ){
case 0x00ff00ff: return 1;
case 0x00ffff00: return 2;
case 0xff0000ff: return 3;
case 0xff00ff00: return 4;
default: return -1;
}
Скорость этого варианта:
count=100000000
measure...
Freq=3.002314 GHz
filling...
calc...
SIMD time Quadrant1 (Freq/1GHz): 557.680927/1674.333081 ms.
SIMD time Quadrant2 (Freq/1GHz): 648.186217/1946.058354 ms.
Press any key to continue . . .
Т.е., как ни крути, а скорость обоих моих вариантов больше по крайней мере в три раза.
Полностью нагрузил узел, на 500 млн., вышел такой результат:
count=500000000
measure...
Freq=3.002184 GHz
filling...
calc...
SIMD time Quadrant1 (Freq/1GHz): 2731.684208/8201.019582 ms.
SIMD time Quadrant2 (Freq/1GHz): 3532.774079/10606.039056 ms.
SIMD avg. time Quadrant1 (Freq/1GHz): 0.000005/0.000016 ms.
SIMD avg. time Quadrant2 (Freq/1GHz): 0.000007/0.000021 ms.
SIMD speed Quadrant1 (Freq/1GHz): 183037262.683/60968029.036 q./sec.
SIMD speed Quadrant2 (Freq/1GHz): 141531835.549/47142952.931 q./sec.
Press any key to continue . . .
3 секунды все таки более стабильный результат, можно точнее посчитать скорость: 140/47-180/60
млн. квадрантов в секунду (Freq/1GHz).
Изменил алгоритм генерации координат, как у @avp:
rand() - INT_MAX / 2
Запретил запись в вектор Q
(удалилил его, но прописал volatile int*q1, volatile int*q2
в параметрах к функциям Quadrant
). Увеличил число координат до 700 млн
. Выиграл 1нс.
на втором алгоритме, т.е. они сблизились:
count=700000000
measure...
Freq=3.000889 GHz
filling...
calc...
SIMD time Quadrant1 (Freq/1GHz): 3779.277685/11341.192149 ms.
SIMD time Quadrant2 (Freq/1GHz): 4058.144274/12178.039779 ms.
SIMD avg. time Quadrant1 (Freq/1GHz): 0.000005/0.000016 ms.
SIMD avg. time Quadrant2 (Freq/1GHz): 0.000006/0.000017 ms.
SIMD speed Quadrant1 (Freq/1GHz): 185220578.727715/61721906.374871 q./ms.
SIMD speed Quadrant2 (Freq/1GHz): 172492635.220518/57480515.148841 q./ms.
Press any key to continue . . .
Для @avp, после анализа его исходников я доработал реализацию с 2 добавлениями:
- Готовлю сразу SIMD-векторы, для подачи в Quadrant
- Немного разворачиваю циклы
Исходник реализации: http://pastebin.com/u63d8abp
Вывод:
count=100000000
measure...
Freq=3.000481 GHz
fill...
calc...
SIMD time Recenter (Freq/1GHz): 321.122363/963.521667 ms.
SIMD time Quadrant1 (Freq/1GHz): 314.661862/944.137053 ms.
SIMD time Quadrant2 (Freq/1GHz): 444.856547/1334.783781 ms.
SIMD avg. time Recenter (Freq/1GHz): 0.000003/0.000010 ms.
SIMD avg. time Quadrant1 (Freq/1GHz): 0.000003/0.000009 ms.
SIMD avg. time Quadrant2 (Freq/1GHz): 0.000004/0.000013 ms.
SIMD speed Recenter (Freq/1GHz): 311407773.002369/103785938.007349 p./s.
SIMD speed Quadrant1 (Freq/1GHz): 317801462.834867/105916826.039450 q./s.
SIMD speed Quadrant2 (Freq/1GHz): 224791566.118078/74918500.976302 q./s.
Press any key to continue . . .
А теперь делаем ошибку: p += Exp
меняем на ++p
. Т.е. объем вычислений не меняем, просто уменьшаем число обращений к ОЗУ
.
Вывод:
count=100000000
measure...
Freq=3.001383 GHz
fill...
calc...
SIMD time Recenter (Freq/1GHz): 26.564881/79.731387 ms.
SIMD time Quadrant1 (Freq/1GHz): 133.211034/399.817350 ms.
SIMD time Quadrant2 (Freq/1GHz): 195.823826/587.742327 ms.
SIMD avg. time Recenter (Freq/1GHz): 0.000000/0.000001 ms.
SIMD avg. time Quadrant1 (Freq/1GHz): 0.000001/0.000004 ms.
SIMD avg. time Quadrant2 (Freq/1GHz): 0.000002/0.000006 ms.
SIMD speed Recenter (Freq/1GHz): 3764368388.323660/1254211217.973670 p./s.
SIMD speed Quadrant1 (Freq/1GHz): 750688565.116046/250114208.400411 q./s.
SIMD speed Quadrant2 (Freq/1GHz): 510663089.915592/170142586.991187 q./s.
Press any key to continue . . .
Теперь, надеюсь, понятно, какую роль в вычислениях играет память и частота.
>>62
тоже получает знак, но оставляет его не в 0-ой позиции, а в 1-ой. То есть для положительного числа получается 0, а для отрицательного — 2. Там наверняка должен ещё быть&
, чтобы выключить остальные биты. – VladD 23 мар '14 в 13:54