将 % 与 SSE2 一起使用？

Question

这是我要转换为 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）。

有什么花哨的内在函数吗？在此 list 上找不到任何内容。除法 + 某种火箭移位？

Answer 1

正如评论所说，看起来你可以使它只是一个带有比较 + 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，转换为整数并返回更糟，因为打包转换指令也解码为一些随机微指令。更不用说 32 位整数不能覆盖 double 的全部范围，但如果您的输入那么大，您的范围缩小精度就完蛋了。

如果你有FMA，你可以避免乘法和减法的y * integer_quotient部分出现舍入错误。 _mm_fmsub_pd.

将 % 与 SSE2 一起使用？

Using % with SSE2?

c++

intrinsics

sse2