C 结构体或数组来设计几何库

如何解决C 结构体或数组来设计几何库

我需要编写一个库来管理 C 中 2D 空间中点的几何变换。这些点将聚合成形状，我希望能够（自动）矢量化通过 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 坐标）更清晰。现在，编译器是否会理解该结构只是一个很好的数组处理程序并相应地进行向量化？

解决方法

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

godbolt.org 上的 Compiler Explorer 是检查编译器输出内容的便捷方式。我写了一个简单的 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,16]
        ret

...这对于 -O2 和 -O2 -DUSE_XY：

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

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

无论如何，您将传递给向量函数的是 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));

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

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

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) 那么上面的复数旋转与几何约定相同。

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