将 % 与 SSE2 一起使用?

Using % with SSE2?

本文关键字：一起 SSE2      更新时间：2023-10-16

这是我尝试转换为SSE2的代码：

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double *left = audioLeft;
double *right = audioRight;
double phase = 0.0;
double bp0 = mNoteFrequency * mHostPitch;
for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
// some other code (that will use phase)
phase += std::clamp(mRadiansPerSample * (bp0 * pB[sampleIndex] + pC[sampleIndex]), 0.0, PI);
while (phase >= TWOPI) { phase -= TWOPI; }
}

以下是我所取得的成就：

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];
double *left = audioLeft;
double *right = audioRight;
double phase = 0.0;
double bp0 = mNoteFrequency * mHostPitch;
__m128d v_boundLower = _mm_set1_pd(0.0);
__m128d v_boundUpper = _mm_set1_pd(PI);
__m128d v_bp0 = _mm_set1_pd(bp0);
__m128d v_radiansPerSample = _mm_set1_pd(mRadiansPerSample);
__m128d v_phase = _mm_set1_pd(phase);
__m128d v_pB = _mm_load_pd(pB);
__m128d v_pC = _mm_load_pd(pC);
__m128d v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);
for (int sampleIndex = 0; sampleIndex < roundintup8(blockSize); sampleIndex += 8, pB += 8, pC += 8) {
// some other code (that will use v_phase)
v_phase = _mm_add_pd(v_phase, v_result);
v_pB = _mm_load_pd(pB + 2);
v_pC = _mm_load_pd(pC + 2);
v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);
v_phase = _mm_add_pd(v_phase, v_result);
v_pB = _mm_load_pd(pB + 4);
v_pC = _mm_load_pd(pC + 4);
v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);
v_phase = _mm_add_pd(v_phase, v_result);
v_pB = _mm_load_pd(pB + 6);
v_pC = _mm_load_pd(pC + 6);
v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);
v_phase = _mm_add_pd(v_phase, v_result);
v_pB = _mm_load_pd(pB + 8);
v_pC = _mm_load_pd(pC + 8);
v_result = _mm_mul_pd(v_bp0, v_pB);
v_result = _mm_add_pd(v_result, v_pC);
v_result = _mm_mul_pd(v_result, v_radiansPerSample);
v_result = _mm_max_pd(v_result, v_boundLower);
v_result = _mm_min_pd(v_result, v_boundUpper);
// ... fmod?
}

但我真的不确定如何替换while (phase >= TWOPI) { phase -= TWOPI; }(这基本上是C++中的经典fmod)。

有什么花哨的内在吗？在此列表中找不到任何内容。 分裂+某种火箭位移位？

正如评论所说，看起来在这里你可以用比较 +andpd来做一个被屏蔽的减法。 只要您永远不能减去超过一个减法回到所需范围，这就可以了。

喜欢

const __m128d v2pi = _mm_set1_pd(TWOPI);

__m128d needs_range_reduction = _mm_cmpge_pd(vphase, v2pi);
__m128d offset = _mm_and_pd(needs_range_reduction, v2pi);  // 0.0 or 2*Pi
vphase = _mm_sub_pd(vphase, offset);

要实现实际(缓慢)fmod而不过多担心有效数的最后几位，您需要执行integer_quotient = floor(x/y)(或者可能是rint(x/y)或ceil)，然后x - y * integer_quotient。floor/rint/ceil使用 SSE4.1_mm_round_pd或_mm_floor_pd()都很便宜。 这将为您提供余数，它可以是负数，就像整数除法一样。

我敢肯定，有一些数字技术可以更好地避免在灾难性取消减去两个附近的数字之前四舍五入误差。 如果您关心精度，请去检查。 (当你不太关心精度时，使用double向量有点愚蠢;不妨使用float，每个向量完成两倍的工作)。 如果输入比模数大得多，则不可避免地会损失精度，并且最小化临时舍入误差可能非常重要。 但除此之外，精度只会是一个问题，除非您关心非常接近零的结果中的相对误差，而x几乎是y的精确倍数。 (接近零的结果，只有有效数的底部几位留给精度。

如果没有SSE4.1，有一些技巧，比如加减一个足够大的数字。 转换为整数并返回对pd来说更糟糕，因为打包转换指令也会解码为一些随机 uop。 更不用说 32 位整数并不能覆盖整个double范围，但如果您的输入如此之大，那么您的降程精度就会被搞砸。

如果您有 FMA，则可以避免乘法和 sub 的y * integer_quotient部分的舍入误差。_mm_fmsub_pd.

相关文章：