模 2*Pi 使用 SSE/SSE2答案

【问题标题】：Modulo 2*Pi using SSE/SSE2模 2*Pi 使用 SSE/SSE2
【发布时间】：2012-07-19 13:03:58
【问题描述】：

我对使用 SSE 还是很陌生，我正在尝试为 1e8 顺序的双精度输入实现 2*Pi 的模（其结果将被输入到一些向量化的三角计算中）。

我目前对代码的尝试是基于mod(x, 2*Pi) = x - floor(x/(2*Pi))*2*Pi 的想法，看起来像：

#define _PD_CONST(Name, Val)                                            \
static const double _pd_##Name[2] __attribute__((aligned(16))) = { Val, Val }  

_PD_CONST(2Pi, 6.283185307179586);  /* = 2*pi  */  
_PD_CONST(recip_2Pi, 0.159154943091895); /* = 1/(2*pi)  */

void vec_mod_2pi(const double * vec, int Size, double * modAns)
{
    __m128d sse_a, sse_b, sse_c;
    int i;
    int k = 0;
    double t = 0;

    unsigned int initial_mode;
    initial_mode = _MM_GET_ROUNDING_MODE();

    _MM_SET_ROUNDING_MODE(_MM_ROUND_DOWN);

    for (i = 0; i < Size; i += 2)
    {
        sse_a = _mm_loadu_pd(vec+i);
        sse_b = _mm_mul_pd( _mm_cvtepi32_pd( _mm_cvtpd_epi32( _mm_mul_pd(sse_a, *(__m128d*)_pd_recip_2Pi) ) ), *(__m128d*)_pd_2Pi);
        sse_c = _mm_sub_pd(sse_a, sse_b);
        _mm_storeu_pd(modAns+i,sse_c);
    }

    k = i-2;
    for (i = 0; i < Size%2; i++)
    {
        t = (double)((int)(vec[k+i] * 0.159154943091895)) * 6.283185307179586;
        modAns[k+i] = vec[k+i] - t;
    }

    _MM_SET_ROUNDING_MODE(initial_mode);
}

不幸的是，这目前返回了很多 NaN 以及几个 1.128e119 的答案（有些超出了我的目标 0 -> 2*Pi 的范围！）。我怀疑我出错的地方是我试图用来做floor的双到整数到双转换。

谁能建议我哪里出了问题以及如何改进它？

附：抱歉该代码的格式，这是我第一次在这里发布问题，似乎无法让它在代码块中给我空行以使其可读。

【问题讨论】：

也许我应该让它更清楚一点 - 我并不担心这里结果的准确性（事实上，如果它以单精度精度返回，那很好）。我最关心的是 a) 让它在 SSE 中运行，这样我就可以摆脱我目前在其他两个 SSE 块之间获得的不必要的存储/加载；b) 如果可能的话，让它快速运行。
首先，像正常人一样使用_mm_set1_pd(6.283185307179586)；编译器已经足够聪明，可以将其转换为 static const double[2] 的等价物，但具有重复消除功能（例如重复的字符串文字）。其次，您可以使用_mm_cvttpd_epi32（注意额外的 t 用于截断）转换为带有截断的整数，而不是当前的舍入模式。如果您的数字是正数，则与floor 相同。如果您有 SSE4.1，则可以使用 _mm_floor_pd (roundpd) 直接舍入为整数，而不是转换为有符号 int32_t 并返回。
目前还不清楚为什么你会得到 NaN。 -2^31..2^31-1 范围之外的 Cvt 将为您提供 0x80000000（即 -2^31，最大负整数），英特尔将其用作“整数不定值”。将其转换回double 应该会成功。这些可能是您的1.128e119 结果？预计您会遇到灾难性的取消错误，特别是因为您的 2Pi * recip_2Pi 可能不完全是 1.0，因为舍入错误会有所不同。（实际上 1.0- 2Pi*recip_2Pi 和 gcc -O0 仅约为 2.1E-15，因此比您的输入小几个数量级。）

标签： c optimization sse simd

【解决方案1】：

如果您想要任何类型的准确性，那么简单的算法是非常糟糕的。有关准确的范围缩小算法，请参见例如Ng et al., ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit（现在可以通过 Wayback Machine 获得：2012-12-24）

【讨论】：

我需要再次阅读那篇论文以确保我已经理解了它，但是这种方法不需要比 XMM 寄存器中可用的更多位吗？这种方法在 SSE 中真的可行吗？
@BigA：确实，您不能一次将所有位都放入一个 XMM reg 中。你必须做一些更复杂的事情。
一个不应忽视的简单“fmod 2pi”算法的优点是由此引入的舍入误差通常会抵消角度计算中较早的舍入误差。例如，如果 Sin 执行参数归约 mod Math.Pi 表达式 Math.Sin(x*(2*Math.Pi)) 的结果通常比执行归约 mod π 时更准确 [对于某些 x 值，归约 mod π 会稍微好一些，但对某些人来说效果更糟]。

【解决方案2】：

对于较大的参数，通常使用Hayne-Panek algorithm。但是，Hayne-Panek 论文很难阅读，我建议您查看Chapter 11 in the Handbook of Floating-Point Arithmetic 以获得更易于理解的解释。

【讨论】：