2次元空間内の点を回転させるための古典的な式が与えられると、2D回転

Question

の最適化：NPOINTS考える

cv::Point pt[NPOINTS]; 
cv::Point rotated[NPOINTS]; 
float angle = WHATEVER; 
float cosine = cos(angle); 
float sine = sin(angle); 

for (int i = 0; i < NPOINTS; i++) 
{ 
    rotated[i].x = pt[i].x * cosine - pt[i].y * sine; 
    rotated[i].y = pt[i].x * sine + pt[i].y * cosine; 
} 

が32であり、配列が整列され、どのようにしてSSEまたはAVXのコードを最適化について行きますか？この辺りの検索や他の場所で有益な何かを上げていなかった、と私は約ここに失われてしまった：

__m128i onePoint = _mm_set_epi32(pt[i].x, pt[i].y, pt[i].x, pt[i].y); 
__m128 onefPoint = _m128_cvtepi32_ps(onePoint); 
__m128 sinCos = _mm_set_ps(cosine, -sine, sine, cosine); 
__m128 rotated = _mm_mul_ps(onefPoint, sinCos); 

をしかし、どのように[y*cosine, -x*sine, x*sine, y*cosine]から[y*cosine + -x*sine, x*sine + y*cosine]に行くために？これが最善のアプローチですか？ __m512に簡単に拡大できますか？

UPDATE：私は少しより多くの研究を行なったし、私は今、約持っている：

__m128i onePoint = _mm_set_epi32(pt[i].x, pt[i].y, pt[i].x, pt[i].y); 
__m128 onefPoint = _m128_cvtepi32_ps(onePoint); 
__m128i twoPoint = _mm_set_epi32(pt[i+1].x, pt[i+1].y, pt[i+1].x, pt[i+1].y); 
__m128 twofPoint = _m128_cvtepi32_ps(twoPoint); 
__m128 sinCos = _mm_set_ps(cosine, -sine, sine, cosine); 
__m128 rotated1 = _mm_mul_ps(onefPoint, sinCos); 
__m128 rotated2 = _mm_mul_ps(twofPoint, sinCos); 
__m128 added = _mm_hadd_ps(rotated1, rotated2); 
__m128i intResult = _mm_cvtps_epi32(added); 
int results[4]; 
_mm_storeu_si128((__m128i*)results, intResult); 

これは、約6％のプロセッサ時間の11％から50％のスピードアップを提供します。 __m256に拡大し、一度に4つのポイントを実行すると、別のスピードアップが得られます。これはかなりひどいコードに見えますが、私は正しい方向に向かっていますか？

Answer 1

アレイの構造体（AoSoA）の配列を使用し、一度に8ポイントを処理します。以下のコードでは、point8は8点を含む配列の構造体です。関数rotate_point8は8つの点を回転し、1つの点を回転させる関数rotate_pointと同じ代数構造を持ちます。機能rotate_all8は、AoSoA point8*を使用して32ポイントを回転します。

シングルポイントローテーションコードは、4回の乗算、1回の加算、および1回の減算を行います。

the assembly for rotate_point8を見ると、GCCがループをアンロールし、アンロールごとに4つのSIMD乗算、1つのSIMD加算、1つのSIMD減算を行うことがわかります。それはあなたができる最高のものです：1つの価格のための8。

#include <x86intrin.h> 
#include <stdio.h> 
#include <math.h> 

struct point8 { 
    __m256 x; 
    __m256 y; 
}; 

struct point { 
    float x; 
    float y; 
}; 

static point rotate_point(point p, float a, float b) { 
    point r; 
    r.x = p.x*a - p.y*b; 
    r.y = p.x*b + p.y*a; 
    return r; 
} 

static point8 rotate_point8(point8 p, float a, float b) { 
    __m256 va = _mm256_set1_ps(a), vb = _mm256_set1_ps(b); 
    point8 r; 
    r.x = _mm256_sub_ps(_mm256_mul_ps(p.x,va), _mm256_mul_ps(p.y,vb)); 
    r.y = _mm256_add_ps(_mm256_mul_ps(p.x,vb), _mm256_mul_ps(p.y,va)); 
    return r; 
} 

void rotate_all(point* points, point* r, float angle) { 
    float a = cos(angle), b = sin(angle); 
    for(int i=0; i<32; i++) r[i] = rotate_point(points[i], a, b); 
} 

void rotate_all8(point8* points, point8* r8, float angle) { 
    float a = cos(angle), b = sin(angle); 
    for(int i=0; i<4; i++) r8[i] = rotate_point8(points[i], a, b); 
} 

int main(void) { 
    float x[32], y[32]; 
    point p[32], r[32]; 
    point8 p8[4], r8[4]; 
    float angle = 3.14159f/4; 

    for(int i=0; i<32; i++) y[i] = 1.0*i/31, x[i] = sqrt(1-y[i]*y[i]); 
    for(int i=0; i<32; i++) p[i].x = x[i], p[i].y = y[i]; 
    for(int i=0; i<4; i++) p8[i].x = _mm256_load_ps(&x[8*i]), p8[i].y = _mm256_load_ps(&y[8*i]); 

    for(int i=0; i<32; i++) printf("%f %f\n", p[i].x, p[i].y); puts(""); 

    rotate_all(p, r, angle); 
    for(int i=0; i<32; i++) printf("%f %f\n", r[i].x, r[i].y); puts(""); 

    rotate_all8(p8, r8, angle); 
    for(int i=0; i<4; i++) { 
    _mm256_storeu_ps(x, r8[i].x), _mm256_storeu_ps(y, r8[i].y); 
    for(int j=0; j<8; j++) printf("%f %f\n", x[j], y[j]); 
    } 
}