Modulo 2*Pi using SSE/SSE2

Question

I'm still pretty new to using SSE and am trying to implement a modulo of 2*Pi for double-precision inputs of the order 1e8 (the result of which will be fed into some vectorised trig calculations).

My current attempt at the code is based around the idea that mod(x, 2*Pi) = x - floor(x/(2*Pi))*2*Pi and looks like:

#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);
}

Unfortunately, this is currently returning a lot of NaN with a couple of answers of 1.128e119 as well (some what outside the range of 0 -> 2*Pi that I was aiming for!). I suspect that where I'm going wrong is in the double-to-int-to-double conversion that I'm trying to use to do the floor.

Can anyone suggest where I've gone wrong and how to improve it?

P.S. sorry about the format of that code, it's the first time I've posted a question on here and can't seem to get it to give me empty lines within the code block to make it readable.

Perhaps I should make it a wee bit clearer - I'm not super-concerned by accuracy of result here (in fact if it returned with single-precision accuracy, that would be fine). What I'm most concerned with is a) getting it to run in SSE so I can get rid of an unnecessary store/load that I've got at the moment between two other SSE blocks and b) making it run quickly if possible. — BigA
First of all, use _mm_set1_pd(6.283185307179586) like a normal person; the compiler is already smart enough to turn that into the equivalent of static const double[2], but with duplicate elimination (like for duplicate string literals). Second, you can use _mm_cvttpd_epi32 (note the extra t for truncate) to convert to integer with truncation instead of the current rounding mode. If your numbers are positive, that's the same as floor. If you have SSE4.1, you can use _mm_floor_pd (roundpd) to round to integer directly instead of converting to signed int32_t and back. — Peter Cordes
It's not clear why you'd get NaN, though. Cvt outside the -2^31..2^31-1 range will give you 0x80000000 (i.e. -2^31, the most negative integer), which Intel uses as the "integer indefinite value". Converting that back to double should succeed. Those might be your 1.128e119 results? It's expected that you get catastrophic cancellation errors, especially since your 2Pi * recip_2Pi probably isn't exactly 1.0 because the rounding errors will be different. (Actually 1.0- 2Pi*recip_2Pi with gcc -O0 is only ~2.1E-15, so that's several orders of magnitude smaller than your inputs.) — Peter Cordes

janneb janneb · Accepted Answer · 2012-07-20T09:08:05

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)

Modulo 2*Pi using SSE/SSE2

2 Answers