根据二阶导数计算的曲线的SIMD优化

因此，您只是在寻找一种生成4个phase_current值的向量的方法，您可以将它们作为arg传递给任意函数。

TL：DR：设置增量和步长的初始向量，以便每个向量元素跨4步跨序列，从而phase_current[i+0..i+3]仅使用两个向量ADD操作（垂直而不是水平）给的向量。您可以使用代数/数学排除此串行依赖性。

这有点像一个前缀和（您可以使用log2(vector_width)带有vector_width元素的向量的shuffle + add操作对SIMD进行SIMD。）您还可以使用两步计算来并行化多个线程的前缀和，其中每个线程对一个对象的区域求和。数组，然后合并结果，并使每个线程将目标数组的区域偏移一个常数（该区域的第一个元素的总和。有关多线程的信息，请参见链接的问题。）

但是，您拥有一个巨大的简化，即phase_increment_step（您想要的值的二阶导数）为常数。我假设这USEFUL_FUNC(phase_current);是通过值而不是通过非const引用获取arg，因此对它的唯一修改phase_current是+=在循环中进行。而且，这useful_func无法以某种方式改变增量或增量。

实现此目的的一种选择是仅在4个独立的SIMD向量元素中运行标量算法，每次偏移1次迭代。使用整数加法时，尤其是在矢量整型加法延迟仅为1个周期的英特尔CPU上，运行运行总计的4次迭代比较便宜，我们可以在调用之间进行USEFUL_FUNC。那将是一种生成USEFUL_FUNC向量输入的方法，其完成的工作量与标量代码一样多（假设SIMD整数加法与标量整数加法一样便宜，如果我们将数据依赖性限制为每个时钟2个加法，则这是正确的）。

上面的方法有些通用，可以解决这个问题，因为存在真正的序列依赖，而我们无法通过简单的数学方法廉价地消除它。

如果我们很聪明，我们可以做得比一次每次执行4个序列的前缀和或强力更好。理想情况下，我们可以得出一种闭合形式的方法，以按值的顺序（或SIMD向量宽度是乘以多个累加器所需的展开因子）乘以4 USEFUL_FUNC。

总结序列step，step*2，step*3，...会给我们一个恒定次高斯的封闭形式的公式整数最多的总和n：sum(1..n) = n*(n+1)/2。该序列为0、1、3、6、10、15、21、28，...（https://oeis.org/A000217）。（我已经排除了首字母缩写phase_increment）。

诀窍是按此顺序前进4。(n+4)*(n+5)/2 - n*(n+1)/2 简化为4*n + 10。再次采用该导数，我们得到4。但是在第二个积分中走4步，我们就有了4*4 = 16。因此，我们可以维持一个向量phase_increment，该向量可以通过将SIMD添加向量来增加16*phase_increment_step。

我不能完全确定我是否拥有正确的步数推理功能（将4乘以16的额外系数有点令人惊讶）。制定正确的公式并在向量序列中求一阶和二阶差分，可以很清楚地看出：

 // design notes, working through the first couple vectors
 // to prove this works correctly.

S = increment_step (constant)
inc0 = increment initial value
p0 = phase_current initial value

// first 8 step-increases:
[ 0*S,  1*S,   2*S,  3*S ]
[ 4*S,  5*S,   6*S,  7*S ]

// first vector of 4 values:
[ p0,  p0+(inc0+S),  p0+(inc0+S)+(inc0+2*S),  p0+(inc0+S)+(inc0+2*S)+(inc0+3*S) ]
[ p0,  p0+inc0+S,  p0+2*inc0+3*S,  p0+3*inc0+6*S ]  // simplified

// next 4 values:
[ p0+4*inc0+10*S,  p0+5*inc0+15*S,  p0+6*inc0+21*S,  p0+7*inc0+28*S ]

使用这个和以前的4*n + 10公式：

// first 4 vectors of of phase_current
[ p0,              p0+1*inc0+ 1*S,  p0+2*inc0+3*S,   p0+ 3*inc0+ 6*S ]
[ p0+4*inc0+10*S,  p0+5*inc0+15*S,  p0+6*inc0+21*S,  p0+ 7*inc0+28*S ]
[ p0+8*inc0+36*S,  p0+9*inc0+45*S,  p0+10*inc0+55*S, p0+11*inc0+66*S ]
[ p0+12*inc0+78*S,  p0+13*inc0+91*S,  p0+14*inc0+105*S, p0+15*inc0+120*S ]

 first 3 vectors of phase_increment (subtract consecutive phase_current vectors):
[ 4*inc0+10*S,     4*inc0 + 14*S,   4*inc0 + 18*S,   4*inc0 + 22*S  ]
[ 4*inc0+26*S,     4*inc0 + 30*S,   4*inc0 + 34*S,   4*inc0 + 38*S  ]
[ 4*inc0+42*S,     4*inc0 + 46*S,   4*inc0 + 50*S,   4*inc0 + 54*S  ]

 first 2 vectors of phase_increment_step:
[        16*S,              16*S,            16*S,            16*S  ]
[        16*S,              16*S,            16*S,            16*S  ]
Yes, as expected, a constant vector works for phase_increment_step

因此，我们可以使用Intel的SSE / AVX内部函数编写如下代码：

#include <stdint.h>
#include <immintrin.h>

void USEFUL_FUNC(__m128i);

// TODO: more efficient generation of initial vector values
void double_integral(uint32_t phase_start, uint32_t phase_increment_start, uint32_t phase_increment_step, unsigned blockSize)
{

    __m128i pstep1 = _mm_set1_epi32(phase_increment_step);

    // each vector element steps by 4
    uint32_t inc0=phase_increment_start, S=phase_increment_step;
    __m128i pincr  = _mm_setr_epi32(4*inc0 + 10*S,  4*inc0 + 14*S,   4*inc0 + 18*S,   4*inc0 + 22*S);

    __m128i phase = _mm_setr_epi32(phase_start,  phase_start+1*inc0+ 1*S,  phase_start+2*inc0+3*S,   phase_start + 3*inc0+ 6*S );
     //_mm_set1_epi32(phase_start); and add.
     // shuffle to do a prefix-sum initializer for the first vector?  Or SSE4.1 pmullo by a vector constant?

    __m128i pstep_stride = _mm_slli_epi32(pstep1, 4);  // stride by pstep * 16
    for (unsigned i = 0; i < blockSize; ++i)  {
        USEFUL_FUNC(phase);
        pincr = _mm_add_epi32(pincr, pstep_stride);
        phase = _mm_add_epi32(phase, pincr);
    }
}

进一步阅读：有关SIMD的更多信息，但主要是x86 SSE / AVX，请参阅https://stackoverflow.com/tags/sse/info，尤其是Insomniac Games（GDC 2015）上来自SIMD的幻灯片，其中有一些有关如何通常考虑一下SIMD，以及如何布置数据以便可以使用它。