Вычисление квадранта, в котором находится точка двумерной координатной плоскости.

#битовые_операции #c #геометрия

Здравствуйте, существует ли машинный способ определения квадранта, в котором находится
точка с координатами x и y?



Наивное решение делает простой перебор вариантов:

if (x > 0 && y > 0)

    printf("point (%d, %d) lies in the First quandrant\n");

else if (x < 0 && y > 0)

    printf("point (%d, %d) lies in the Second quandrant\n");

else if (x < 0 && y < 0)

    printf("point (%d, %d) lies in the Third quandrant\n");

else if (x > 0 && y < 0)

    printf("point (%d, %d) lies in the Fourth quandrant\n");

else if (x == 0 && y == 0)

    printf("point (%d, %d) lies at the origin\n");




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

if (point.x > centerPoint.x) {
    if (point.y > centerPoint.y) {
        return KPClusterDistributionQuadrantFour;
    } else {
        return KPClusterDistributionQuadrantOne;
    }
} else if (point.x < centerPoint.x) {
    if (point.y > centerPoint.y) {
        return KPClusterDistributionQuadrantThree;
    } else {
        return KPClusterDistributionQuadrantTwo;
    }
} else {
    if (point.y > centerPoint.y) {
        return KPClusterDistributionQuadrantThree;
    } else if (point.y < centerPoint.y) {
        return KPClusterDistributionQuadrantTwo;
    }
}

return KPClusterDistributionQuadrantNone;




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

квадрант I)   x > 0, y > 0  => ++
квадрант II)  x < 0, y > 0  => -+
квадрант III) x < 0, y < 0  => --
квадрант IV)  x > 0, y < 0  => +-


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

+ ? + = 1
- ? + = 2
- ? - = 3
+ ? - = 4

или в контексте моей координатной системы с иным распределением квадрантов:

+ ? - = 1
- ? - = 2
- ? + = 3
+ ? + = 4


Вопрос в том, что может стоять за этим "?".

Буду рад любым советам. Интересуюсь этим вопросом больше для образования, хотя практика
тоже обещает выиграть от этого, если такое решение возможно и существует.

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

Спасибо.



ОБНОВЛЕНО ПОЗЖЕ

Только что нашёл решение на SO, только теперь нужно попытаться его, как следует,
понять: Optimizing quadrant selection.


  The fastest way in C/C++ it would be
  (((unsigned int)x >> 30) & 2) | ((unsigned int)y >> 31)
  (30/31 or 62/63, depending on size of int). This will give the quadrants in order
0, 2, 3, 1.


И решение относительно данного центра:


  (((unsigned int)(x - center.x) >> 30) & 2) | ((unsigned int)(y-center.y) >> 31)


Буду признателен, если кто подскажет, что нужно изменить в этом решении, чтобы оно
стало работать для точек x, y с типом double.



Комментарий по поводу принятого ответа

Оказалось, что всё-таки обычное решение с двойным if-else ветвлением работает быстрее,
чем решение в принятом ответе (соотношение, которое я получаю достаточно стабильно:
160 к 180). 



Вот окончательные версии, которые я сравнивал между собой (система координат инвертирована
относительно оси y относительно декартовой):

Версия с обычным ветвлением:

static inline KPClusterDistributionQuadrant KPClusterDistributionQuadrantForPointInsideMapRectBranching(MKMapRect
mapRect, MKMapPoint point) {
    MKMapPoint centerPoint = (MKMapPoint){
        mapRect.origin.x + mapRect.size.width / 2,
        mapRect.origin.y + mapRect.size.height / 2
    };

    if (point.x >= centerPoint.x) {
        if (point.y >= centerPoint.y) {
            return KPClusterDistributionQuadrantFour;
        } else {
            return KPClusterDistributionQuadrantOne;
        }
    } else {
        if (point.y >= centerPoint.y) {
            return KPClusterDistributionQuadrantThree;
        } else {
            return KPClusterDistributionQuadrantTwo;
        }
    }
}


Вариант @paulgri:

NS_INLINE int KPClusterDistributionQuadrantForPointInsideMapRect(MKMapRect mapRect,
MKMapPoint point) {
    MKMapPoint centerPoint = (MKMapPoint){
        mapRect.origin.x + mapRect.size.width / 2,
        mapRect.origin.y + mapRect.size.height / 2
    };

    double dx = point.x - centerPoint.x;
    double dy = point.y - centerPoint.y;

    return ((*((long long *)&dy) >> 63) & 3) ^ ((*((long long *)&dx) >> 63) & 1);
}


Решение со StackOverflow (адаптировано @avp для работы с double):

NS_INLINE int KPClusterDistributionQuadrantForPointInsideMapRect_avp(MKMapRect mapRect,
MKMapPoint point) {
    MKMapPoint centerPoint = (MKMapPoint){
        mapRect.origin.x + mapRect.size.width / 2,
        mapRect.origin.y + mapRect.size.height / 2
    };

    double dx = point.x - centerPoint.x;
    double dy = point.y - centerPoint.y;

    return ((*((uint64_t *)&dx) >> 62) & 2) | (*((uint64_t *)&dy) >> 63);
}


Еще одно решение, предложенное моим коллегой - оно более медленное, но все же рабочее,
поэтому я добавляю его сюда для коллекции:

uint32_t x = ...;
uint32_t y = ...;

char xsign = x >> 31 & 1;
char ysign = y >> 31 & 1;

int quadrant = 1 + (xsign ^ ysign) + 2 * ysign;

    

Ответы

Ответ 1

Если квадрант отыскивается относительно точки О(0,0), то можно обойтись битовыми
операциями и манипулировать знаками чисел, зная, что сдвиг вправо является арифметическим
(надеюсь, переменныне типа int?), например так (к размерам типов и пр. прошу не придираться,
это набросок кода - только идея):
n = ((y>>31)&3 ^ (x>>31)&1) + 1;

Здесь +1 надо только для того, чтобы квадратны начинались с 1.
Если относительно некоторой точки (x0,y0), то достаточно вместо x и y поставить разности.
Как вариант - используйте тернарную операцию и обычные сравнения.
UPD: одновременно написали ))
Для типов с плавающей запятой битовые операции не применимы, используйте тернарную.
UPD2: обратите внимание, это не тот же вариант, в моем варианте квадратны нумеруются
правильно, по порядку 1,2,3,4.
UPD3:
Если x и y переменные типа double, то можно указатели на них привести к указателям
на целочисленный тип так, чтобы знаковые разряды мантиссы оказались на месте знаковых
разрядов целого и дальше все по той же схеме. Только посмотреть надо, где они у double.
Приведение указателей на примитивные типы - фактически фиктивная операция, время не
тратится ни на какие преобразования.


Ответ 2

Под машинным способом я главным образом подразумеваю использование битовых операций

Если Вы под этим завуалировали "как получить быстрое решение?", то мой вариант: слепите
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 
#include 
#include 
#include 
#include 
#include 
//
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 . . .

Теперь, надеюсь, понятно, какую роль в вычислениях играет память и частота.


Ответ 3

@Stanislaw Pankevich, решение SO (только для несдвинутого нуля) для double
  for (double a = 1; a > -2; a -= 2) 
    for (double b = 1; b > -2; b -= 2) {
      uint64_t *pa = (void *)&a, *pb = (void *)&b;
      int qn = ((*pa >> 62) & 2) | (*pb >> 63);
      printf("[%f , %f] qn = %d\n", a, b, qn);
    }

основано на том, что знак double хранится в старшем разряде long long.
--
Кстати, странно, что оно оказалось медленнее вложенных if-ов. Попробуйте один-в-один
в "заинлайненном варианте" (вообще без функций, структур и т.п.).
UPDATE
Я измерил время вот так
#include 
#include 
#include 
#include 
#include 
#include 

long long
mtime()
{
  struct timeval t;

  gettimeofday(&t, NULL);
  long long mt = (long long)t.tv_sec * 1000 + t.tv_usec / 1000;
  return mt;
}

int 
main (int ac, char *av[])
{
  double a = 10, b = -10;
  uint64_t *pa = (void *)&a, *pb = (void *)&b;

  int n = atoi(av[1] ? av[1] : "10000000");
  int qn, cnt[4] = {0};
  srand(0);
  long long start = mtime();
  for (int i = 0; i < n; i++) {
    a = rand() - INT_MAX / 2;
    b = rand() - INT_MAX / 2;
    qn = ((*pa >> 62) & 2) | (*pb >> 63);
    cnt[qn]++;
  }
  printf ("shift: %lld (msec) n: %d [%d %d %d %d]\n",
          mtime() - start, n, cnt[0], cnt[1], cnt[2], cnt[3]);

  memset(cnt, 0, sizeof(cnt));
  srand(0);
  start = mtime();
  for (int i = 0; i < n; i++) {
    a = rand() - INT_MAX / 2;
    b = rand() - INT_MAX / 2;
    if ( a < 0) 
      if (b < 0)
        qn = 3;
      else
        qn = 2;
    else
      if (b < 0)
        qn = 1;
      else
        qn = 0;
    cnt[qn]++;
  }
  printf ("if: %lld (msec) n: %d [%d %d %d %d]\n",
          mtime() - start, n, cnt[0], cnt[1], cnt[2], cnt[3]);

}

avp@avp-ubu1:hashcode$ gcc tttx.c -std=gnu99 -O3
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 1874 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
if: 2409 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 1838 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
if: 2285 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 1848 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
if: 2310 (msec) n: 100000000 [25003913 24997412 24993577 25005098]
avp@avp-ubu1:hashcode$ grep CPU /proc/cpuinfo | head -1
model name  : Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz
avp@avp-ubu1:hashcode$

Получил, что if-ами медленнее.
UPDATE-2
Добавил double dx, dy и перемерил
  ....
  double dx = 1000, dy = 2000;
  ....
  for (int i = 0; i < n; i++) {
    a = rand() - INT_MAX / 2;
    b = rand() - INT_MAX / 2;
    a += dx; b += dy;  // вот это добавил
    qn = ((*pa >> 62) & 2) | (*pb >> 63);
    cnt[qn]++;
  }
  ......
    if ( a < -dx)   // вот так изменил тут
      if (b < -dy)
    ...
      if (b < -dy)

а вот результат
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 2077 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
if: 2487 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 1962 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
if: 2286 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
avp@avp-ubu1:hashcode$ ./a.out 100000000
shift: 2006 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
if: 2307 (msec) n: 100000000 [25003993 24997380 24993591 25005036]
avp@avp-ubu1:hashcode$

все равно показывает, что if-ами медленнее.


Ответ 4

У меня везде закончились комментарии, поэтому задам свой последний вопрос по поводу
обоих битовых вариантов здесь в виде ответа:
вариант @paulgri:
double dx = point.x - center.x;
double dy = point.y - center.y;

return ((*((long long *)&dy) >> 63) & 3) ^ ((*((long long *)&dx) >> 63) & 1);

вариант SO и @avp:
double dx = point.x - center.x;
double dy = point.y - center.y;

return ((*((uint64_t *)&dx) >> 62) & 2) | (*((uint64_t *)&dy) >> 63);

Что я не понимаю, так это то, почему, хотя мы и делаем очень похожие сдвиги в обоих
случаях, в первом варианте каст происходит в long long , а потом берётся pdx и *pdy,
а во втором варианте каст происходит в uint64t , а потом тоже берётся pdx и *pdy? Оба
варианта перестают работать, если я заменяю в первом long long * на uint64t * и наоборот:
во втором заменяю uint64t * на long long *. Прошу прокомментировать эту разницу использования
промежуточного указателя на знаковый 64-битный тип в одном варианте и использования
промежуточного указателя на беззнаковый 64-битный тип в другом варианте. Спасибо.
P.S. ХэшКод съедает своей разметкой символы _ и *, но смысл вроде всё равно ясен!


Ответ 5

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

function quadrant($x,$y) { 
    return ($x==0 || $y==0) ? 0 : ($x>0 ? ($y>0?1:4) : ($y>0?2:3));
} 
printf("x = %5.2f   y= %5.2f   quadrant= %d
", $x= 0,    $y= 0,    quadrant($x,$y));
printf("x = %5.2f   y= %5.2f   quadrant= %d
", $x= 1e-1, $y= 1e-1, quadrant($x,$y));
printf("x = %5.2f   y= %5.2f   quadrant= %d
", $x=-1e-1, $y= 1e-1, quadrant($x,$y));
printf("x = %5.2f   y= %5.2f   quadrant= %d
", $x=-1e-1, $y=-1e-1, quadrant($x,$y));
printf("x = %5.2f   y= %5.2f   quadrant= %d
", $x= 1e-1, $y=-1e-1, quadrant($x,$y));


Результаты:

x = 0.00   y= 0.00   quadrant= 0
x = 0.10   y= 0.10   quadrant= 1
x = -0.10   y= 0.10   quadrant= 2
x = -0.10   y= -0.10   quadrant= 3
x = 0.10   y= -0.10   quadrant= 4