优化二维旋转

Question

鉴于在二维空间中旋转的点的经典公式：

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; 
} 

鉴于NPOINTS是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？

更新：我做了一些更多的研究，我现在有大约：

__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); 

这使从处理器时间的11％一50％的增速约6％。扩展到__m256并一次完成四个点可以提高速度。这看起来很糟糕的代码，但我正朝着正确的方向吗？

Answer 1

使用数组结构数组（AoSoA）并一次处理8个点。在下面的代码中，point8是包含8个点的数组的结构。函数rotate_point8旋转8个点，并具有与旋转单个点的函数rotate_point相同的代数结构。功能rotate_all8使用AoSoA point8*旋转32个点。

单点旋转代码执行4次乘法，一次加法和一次减法。

如果我们看一下the assembly for rotate_point8，我们看到GCC展开了循环并进行了4次SIMD乘法，一次SIMD加法和一次SIMD减法。这是你可以做的最好的事情：八个是一个的价格。

#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]); 
    } 
}