使用SSE/SSE2的模2*Pi

Question

我对使用SSE还很陌生，并且正在尝试为1e8级的双精度输入实现一个2*Pi模数(其结果将被输入到一些矢量化的trig计算中)。

我目前的代码尝试是基于这样的想法，即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。

有没有人能建议我哪里出了错，以及如何改进？

另外，关于代码的格式很抱歉，这是我第一次在这里发布问题，似乎不能让它在代码块中给我空行来使其可读。

Answer 1

如果你想要任何精度，简单的算法是非常糟糕的。有关精确的距离缩减算法，请参见Ng et al., ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit (现可通过Wayback Machine：2012-12-24获得)

Answer 2

对于大型参数，通常使用Hayne-Panek algorithm。然而，海恩-帕内克的论文很难阅读，我建议看看Chapter 11 in the Handbook of Floating-Point Arithmetic，以获得更容易理解的解释。