用于设计几何库的 C 结构或数组

Question

我需要编写一个库来管理 C 中二维 space 中点的几何变换。这些点将以形状聚合，我希望能够（自动）向量化处理通过 OpenMP 的完整形状。

我遇到的问题是继续声明要点的最佳方式：

typedef __attribute__((aligned(8))) float point_t[2];

或

typedef struct point_t
{
  float x, y;
} point_t;

知道了，以后我会用盒子类型的:

typedef __attribute__((aligned(64))) point_t box_t[4];

从编程的角度来看，访问 box[1].y 比 box[1][1]（矩形框第二个点的 y 坐标）更易读。现在，编译器是否会理解该结构只是一个很好的数组处理程序并相应地进行矢量化？

Answer 1

这将取决于编译器。唯一可以确定的方法是检查结果。

godbolt.org 上的编译器资源管理器是检查编译器输出内容的便捷方式。我写了一个简单的 translate 函数：

#ifdef USE_XY
#define X(p) ((p).x)
#define Y(p) ((p).y)
typedef struct point_t
{
  float x, y;
} point_t;
#else
#define X(p) ((p)[0])
#define Y(p) ((p)[1])
typedef __attribute__((aligned(8))) float point_t[2];
#endif

typedef __attribute__((aligned(64))) point_t box_t[4];

void translate(box_t* box, float dx, float dy) {
    #pragma omp simd
    for (int i=0; i<4; ++i) {
        X((*box)[i]) += dx;
        Y((*box)[i]) += dy;
    }
}

使用 ARM64 gcc 8.2 编译（结果为 https://gcc.godbolt.org/z/EEP17Yd7G），我们得到 -O2 -fopenmp 和 -O2 -fopenmp -DUSE_XY:

translate:
        dup     v0.4s, v0.s[0]
        ldr     q3, [x0]
        ldr     q2, [x0, 16]
        ins     v0.s[1], v1.s[0]
        ins     v0.s[3], v1.s[0]
        fadd    v3.4s, v3.4s, v0.4s
        fadd    v0.4s, v2.4s, v0.4s
        str     q3, [x0]
        str     q0, [x0, 16]
        ret

...-O2 和 -O2 -DUSE_XY:

translate:
        add     x1, x0, 32
.L2:
        ldp     s3, s2, [x0]
        fadd    s3, s3, s0
        fadd    s2, s2, s1
        stp     s3, s2, [x0]
        add     x0, x0, 8
        cmp     x0, x1
        bne     .L2
        ret

前者使用SIMD指令，后者不用。 -DUSE_XY 是否存在并没有什么区别。所以，我们知道对于这个确切的代码和这些确切的编译器标志，这个确切的编译器版本是能够做到这一点。当然，这并不能保证您的所有代码都能成功。

Answer 2

您可能要考虑的另一种替代方法是使用复数，例如

typedef _Complex float pointT;

当然，您会想要检查生成的程序集，但希望优化器在内置类型上做得很好似乎是合理的。

复数可以简化编码。

加减法可以直接写，不用写函数：如果a和b是pointT的那么 a+b 是点 T，它是 a 和 b 的总和。类似地，用标量（整数或浮点数）缩放点 T 可以直接写为点乘以标量。一个点的长度由 cabs() 给出，它的方向由 carg() 给出。

旋转和（各向同性）缩放特别好。如果p是一个点，m是一个长度为s，方向为d的复数，那么

m = s*cexp( I*d)
m*p is p rotated through d and scaled by s.

而且如果m和n表示这样的变换，那么

m*n represents n then m

（特别是 1/m 是由 m 表示的变换的逆）。

如果 p 和 q 是复数那么

p*conj(q) is dot(p,q) + I*cross(p,q)

最后一个很容易产生错误的烦人的事情是数学角度之间的约定差异，从 x 轴逆时针方向，和一个常见的几何约定，角度从 y 轴顺时针方向。如果你认为复数z代表y+Ix，所以y = creal(z) 和x = cimag(z) 那么上面的复数旋转与几何约定相同

一个缺点是，如果您需要实现其他变换（例如非各向同性缩放和剪切），那么您将需要对它们进行不同的表示，如果没有一种表示任何方式的方法，可能会很尴尬（线性）变换。

Answer 3

无论如何，您将传递给向量函数的是 float*。您唯一的问题是确保您的 x,y 结构正确映射到您的 2 元素数组。

尽管 C 标准不保证结构中不会有填充（开头除外），但我看不出编译器不执行您期望的操作的原因。我很确定 GNU 和 Microsoft 都会默认将两个浮点数打包成 8 个字节。

我想说一个简单的 typedef struct { float x,y; } point_t; 应该是安全的。

为了安全起见，您可以添加偏执检查，例如：

struct { float x,y; } coords; // expected to map to float[2]
float vector[2];
assert (offsetof(coords,x) == 0); // already guaranteed by the C standard
assert (offsetof(coords,y) == sizeof(float));
assert (sizeof(coords) == sizeof(vector));

如果该代码运行良好（或者更确切地说，编译器优化了该代码），我看不出它以后怎么会对你耍花招。

用于设计几何库的 C 结构或数组

C struct or array to design a geometric library

c

lib