使用SSE/SSE2的Modulo 2 * Pi [英] Modulo 2*Pi using SSE/SSE2

问题描述

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

我目前对代码的尝试基于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的double-to-int-double转换.

有人可以建议我哪里出了问题以及如何改进吗?

P.S.对于该代码的格式感到抱歉，这是我第一次在此处发布问题，而且似乎无法理解它，无法在代码块中给我空行以使其可读.

解决方案

如果您想要任何一种准确性，那么简单的算法就非常糟糕.有关准确的范围缩小算法，请参见 Ng等人，针对大型论证的论据减少:至善至美(现在可以通过Wayback Machine获得: 解决方案

If you want any kind of accuracy, the simple algorithm is terribly bad. For an accurate range reduction algorithm, see e.g. Ng et al., ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit (now available via the Wayback Machine: 2012-12-24)

这篇关于使用SSE/SSE2的Modulo 2 * Pi的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！