首先，让我们详细地检查这段代码。我们有：
complex float M[rows][cols];
complex float v[cols];
float delta[rows];
complex float p = 1.0f;
float s = 1.0f;
delta[]在OP的代码中是
complex float类型，但它只包含
-1.0f或者
+1.0f。（此外，如果它是
-2.0f或
+2.0f，计算会更简单。）因此，我将其定义为实数而不是复数。
同样地，OP将s定义为int，但实际上只使用了-1.0f和+1.0f（在计算中）。这就是我将其显式定义为float的原因。
我省略了f数组，因为完全有一种简单的方法可以避免使用它。
代码中有趣部分的第一个循环：
for (i = 0; i < n; i++) {
v[i] = 0;
for (j = 0; j <n; j++) {
v[i] += M[j][i];
}
p *= v[i];
delta[i] = 1;
}
执行了几件事情。它将delta []数组中的所有元素初始化为1；它可以（并且可能应该）拆分为单独的循环。
由于外部循环增加了i，因此p将是v中元素的乘积；它也可以拆分成单独的循环。
由于内部循环将第i列中的所有元素相加到v [i]，因此外部和内部循环仅将每行作为向量相加到向量v。
因此，我们可以将上述内容重写为伪代码：
Copy first row of matrix M to vector v
For r = 1 .. rows-1:
Add complex values in row r of matrix M to vector v
p = product of complex elements in vector v
delta = 1.0f, 1.0f, 1.0f, .., 1.0f, 1.0f
接下来我们看第二个嵌套循环：
j = 0;
while (j < n-1) {
prod = 1.;
for (i = 0; i < n; i++) {
v[i] -= 2.*delta[j]*M[j][i];
prod *= v[i];
}
delta[j] = -delta[j];
s = -s;
p += s*prod;
f[0] = 0;
f[j] = f[j+1];
f[j+1] = j+1;
j = f[0];
}
只有在观察循环进展的 j 的值时才能看到，但是外层循环主体的最后四行实现了OEISA007814整数序列在j（0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,...）中。该循环迭代次数为2rows-1-1。这部分序列是对称的，并实现了高度为rows-1的二叉树：
4
3 3
2 2 2 2 (Read horizontally)
1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
原来，如果我们循环遍历
i = 1 .. 2
rows-1，那么
r 就是
i 中零位的数量。GCC 提供了一个
__builtin_ctz() 内置函数来计算这个值。（注意，
__builtin_ctz(0) 的值是未定义的；所以即使它在您的电脑上产生了特定的值，也不要这样做。）
内部循环从矩阵的第
j 行中减去按比例缩放的复数值，并从向量
v[] 中计算出复数项的乘积（在减法之后），并将其存储在变量
prod 中。
在内部循环之后，将反转
delta[j] 和比例因子
s 的值。将变量
prod 的值乘以
s 的值后加到
p 上。
在循环之后，最终（复数）结果是将
p 除以 2
rows-1 得到的。这可以使用 C99 函数
ldexp()（分别在实部和虚部上）更好地完成。
因此，我们可以将第二个嵌套循环重写为伪代码。
s = 1.0f
For k = 1 .. rows-1, inclusive:
r = __builtin_ctz(k), i.e. number of least
significant bits that
are zero in k
Subtract the complex values on row r of matrix M,
scaled by delta[r], from vector v[]
prod = the product of (complex) elements in vector v[]
Negate scale factor s (changing its sign)
Add prod times s to result p
在我的经验中，将实部和虚部分别拆分为不同的向量和矩阵是最好的。请考虑以下定义：
typedef struct {
float *real;
float *imag;
size_t floats_per_row;
size_t rows;
size_t cols;
} complex_float_matrix;
void float_set(float *, float, size_t);
void float_copy(float *, const float *, size_t);
float *float_malloc(size_t);
void float_add(float *, const float *, size_t);
void float_add_scaled(float *, const float *, float, size_t);
complex float complex_float_product(const float *, const float *, size_t);
所有上述内容都可以轻松进行矢量化，只要
float_malloc()返回的指针具有足够的对齐度（并且通过例如GCC函数属性
__attribute__((__assume_aligned__(BYTES_IN_VECTOR)));告知编译器），并且矩阵中的
floats_per_row成员是向量中浮点数数量的倍数。
（我不知道GCC是否可以自动矢量化上述函数，但我知道可以使用GCC矢量扩展手动进行矢量化。）
有了上述内容，整个代码的有趣部分，以伪C语言表示，变成了：
complex float permanent(const complex_float_matrix *m)
{
float *v_real, *v_imag;
float *scale;
complex float result;
complex float product;
float coeff = 1.0f;
size_t i = 1 << (m->rows - 1);
size_t r;
if (!m || m->cols < 1 || m->rows < 1 || !i) {
}
v_real = float_malloc(m->cols);
v_imag = float_malloc(m->cols);
scale = float_malloc(m->rows - 1);
if (!v_real || !v_imag || !scale) {
free(scale);
free(v_imag);
free(v_real);
}
float_set(scale, 2.0f, m->rows - 1);
float_copy(v_real, m->real, m->cols);
for (r = 1; r < m->rows; r++)
float_add(v_real, m->real + r * m->floats_per_row, m->cols);
float_copy(v_imag, m->imag, m->cols);
for (r = 1; r < m->rows; r++)
float_add(v_imag, m->imag + r * m->floats_per_row, m->cols);
result = complex_float_product(v_real, v_imag, m->cols);
while (--i) {
r = __builtin_ctz(i);
scale[r] = -scale[r];
float_add_scaled(v_real, m->real + r * m->floats_per_row, m->cols);
float_add_scaled(v_imag, m->imag + r * m->floats_per_row, m->cols);
product = complex_float_product(v_real, v_imag, m->cols);
coeff = -coeff;
result += coeff * product;
}
free(scale);
free(v_imag);
free(v_real);
return result;
}
目前，我个人会先不使用向量化实现上述内容，然后进行广泛测试，直到确定其正确性。
然后，我会检查GCC汇编输出（-S）以查看它是否可以足够地对各个操作（我之前列出的函数）进行向量化。
使用GCC向量扩展手动对函数进行向量化非常简单。声明一个浮点向量很容易：
typedef float vec2f __attribute__((vector_size (8), aligned (8)));
typedef float vec4f __attribute__((vector_size (16), aligned (16)));
typedef float vec8f __attribute__((vector_size (32), aligned (32)));
typedef float vec16f __attribute__((vector_size (64), aligned (64)));
每个向量中的单独组件可以使用数组表示法（对于
vec2f v;，使用
v [0]和
v [1]）。GCC可以对整个向量逐元素执行基本操作；在这里我们只需要加法和乘法。应避免水平操作（操作在同一向量中的元素），而应重新排序元素。
即使在没有此类向量化的体系结构上，GCC也会为上述向量大小生成可用的代码，但生成的代码可能会很慢。（至少到GCC版本5.4，在堆栈和返回之间会生成许多不必要的移动指令。）
为向量分配的内存需要具有足够的对齐性。
malloc()并非在所有情况下都提供足够对齐的内存；您应该改用
posix_memalign()。在本地或静态分配向量类型时，可以使用
aligned属性来增加GCC使用的对齐方式。在矩阵中，您需要确保行从足够对齐的边界开始；这就是我在结构中使用
floats_per_row变量的原因。
在向量（或行）中的元素数量很大，但又不是向量中浮点数数量的倍数的情况下，您应使用“惰性”值填充向量--这些值不会影响结果，例如加法和减法用
0.0f，乘法用
1.0f。
至少在x86和x86-64上，GCC将为仅使用指针的循环生成更好的代码。例如：
void float_set(float *array, float value, size_t count)
{
float *const limit = array + count;
while (array < limit)
*(array++) = value;
}
比起其他写法，这种写法能够生成更好的代码。
void float_set(float *array, float value, size_t count)
{
size_t i;
for (i = 0; i < count; i++)
array[i] = value;
}
or
void float_set(float *array, float value, size_t count)
{
while (count
*(array++) = value;
}
(这些通常会产生相似的代码)。个人而言，我会将其实现为
void float_set(float *array, float value, size_t count)
{
if (!((uintptr_t)array & 7) && !(count & 1)) {
uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8);
uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)array, 8);
uint64_t val;
__builtin_memcpy(&val, &value, 4);
__builtin_memcpy(4 + (char *)&val, &value, 4);
while (ptr < end)
*(ptr++) = val;
} else {
uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4);
uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4);
uint32_t val;
__builtin_memcpy(&val, &value, 4);
while (ptr < end)
*(ptr++) = val;
}
}
还有float_copy()
void float_copy(float *target, const float *source, size_t count)
{
if (!((uintptr_t)source & 7) &&
!((uintptr_t)target & 7) && !(count & 1)) {
uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8);
uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)target, 8);
uint64_t *src = (uint64_t *)__builtin_assume_aligned((void *)source, 8);
while (ptr < end)
*(ptr++) = *(src++);
} else {
uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4);
uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4);
uint32_t *src = (uint32_t *)__builtin_assume_aligned((void *)source, 4);
while (ptr < end)
*(ptr++) = *(src++);
}
}
或者类似于这样的东西。
最难向量化的是complex_float_product()函数。如果你用1.0f填充最终向量中未使用的元素作为实部，用0.0f填充虚部，就可以轻松地计算每个向量的复数积。请记住：
(a + b i) × (c + d i) = (a c - b d) + (a d + b c) i
这里的难点在于高效地计算向量中元素的复数积。幸运的是，这部分对于整体性能并不是很关键（除了非常短的向量或具有非常少列的矩阵），因此在实践中它应该并不重要。
（简而言之，“困难”的部分是找到重新排序元素以最大限度地使用打包向量乘法的方法，而不需要太多的洗牌/移动来减慢速度。）
对于float_add_scaled()函数，您应该创建一个填充有比例因子的向量，类似于以下内容：
void float_add_scaled(float *array, const float *source, float scale, size_t count)
{
const vec4f coeff = { scale, scale, scale, scale };
vec4f *ptr = (vec4f *)__builtin_assume_aligned((void *)array, 16);
vec4f *const end = (vec4f *)__builtin_assume_aligned((void *)(array + count), 16);
const vec4f *src = (vec4f *)__builtin_assume_aligned((void *)source, 16);
while (ptr < end)
*(ptr++) += *(src++) * coeff;
}
如果我们忽略对齐和大小的检查，以及备用实现。
malloc并不保证这一点。如果您切换到posix_memalign，会发生什么？ - Sean Brightmalloc将始终返回16字节对齐的指针，所以这可能是一个误导。 - Sean Brightcomplex类型，因为它是AOS打包。你需要重新排列数据到SOA。 - Mysticial