3

In float, it seems pretty easy to floor() and than int(), such as:

float z = floor(LOG2EF * x + 0.5f);
const int32_t n = int32_t(z);

become:

__m128 z = _mm_add_ps(_mm_mul_ps(log2ef, x), half);
__m128 t = _mm_cvtepi32_ps(_mm_cvttps_epi32(z));
z = _mm_sub_ps(t, _mm_and_ps(_mm_cmplt_ps(z, t), one));

__m128i n = _mm_cvtps_epi32(z);

But how would you achieve this in double using only SSE2?

This is the double version I'd like to convert:

double z = floor(LOG2E * x + 0.5);
const int32_t n = int32_t(z);
markzzz
  • 47,390
  • 120
  • 299
  • 507
  • Any reason you can't just convert the double to an int? You don't need floor at all. (consider `int x = 5/2.0`) – UKMonkey Jan 28 '19 at 16:21
  • @UKMonkey: https://github.com/dpiparo/vdt/blob/master/include/exp.h#L73 – markzzz Jan 28 '19 at 16:23
  • 1
    @markzzz Did you take a look [here](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=4966&techs=SSE2) already? `_mm_cvttpd_epi32` appears to be just what you are looking for. – Max Langhof Jan 28 '19 at 16:37
  • Question why do you need use intrinsic? Compiler doesn't do it for you if proper flags are set? – Marek R Jan 28 '19 at 17:00
  • @MarekR no, it doesn't! And I'm learning... – markzzz Jan 28 '19 at 17:58
  • Do you actually want `round(x)` (that is almost what `floor(x+0.5)` would give you)? Do you care about border cases? (Is `round(1.5)==round(2.5)` ok?) And do you want to have the result as `int32`, `int64` or `double`? – chtz Jan 29 '19 at 07:53
  • Also, if you want an integer result, what shall happen with overflows? – chtz Jan 29 '19 at 07:59

1 Answers1

5

Just use the double precision equivalent (...pd...) of your single precision (...ps...) intrinsic:

__m128i n = _mm_cvtpd_epi32(z);

According to the Intel Intrinsics Guide, that intrinsic is indeed available for SSE2: https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=4966,1917&techs=SSE2

__m128i _mm_cvtpd_epi32 (__m128d a)

Convert packed double-precision (64-bit) floating-point elements in a to packed 32-bit integers, and store the results in dst.

FOR j := 0 to 1
  i := 32*j
  k := 64*j
  dst[i+31:i] := Convert_FP64_To_Int32(a[k+63:k])
ENDFOR
Max Langhof
  • 23,383
  • 5
  • 39
  • 72
  • Uhm, but I'm in double. 64bit, not 32. It will loose precision, no?. _mm_cvtpd_epi64 doesn't exist for SSE2 – markzzz Jan 28 '19 at 17:56
  • 4
    Your question explicitly wants a `const int32_t`. That's what the `32` in the intrinsic name stands for. The **d**ouble **p**recision (64 bit floating point) is denoted by `pd`, as opposed to `ps` for **s**ingle **p**recision. Check out the intrinsic guide, it notes the involved data types very clearly. I also added a quote of the most relevant information to the answer for your convenience. – Max Langhof Jan 28 '19 at 18:00
  • 2
    @markzzz and Max: `cvtpd_epi32` uses the current rounding mode. `cvttpd_epi32` uses truncation (toward 0). Neither of these is `floor` (towards -Infinity). So the question title is misleading, because you're actually asking about `(int)floor(d + 0.5)` being used as a bad emulation of `(int)nearbyint(d)` or `lrint(d)`. – Peter Cordes Jan 28 '19 at 23:24
  • Also, @Mark's actual problem (in [How to convert scalar code of the double version of VDT's Pade Exp fast\_ex() approx into SSE2?](//stackoverflow.com/q/54364694)) wants the integers lined up with the doubles, so they can be shifted and stuffed into the exponent field of a `double` as part of implementing `exp()`. So you'd need a `_mm_shuffle_epi32()` after this to create a vector of `[ 0,n2, 0,n1 ]`. Unfortunately there's no way to get SSE2/AVX/anything to not shuffle the int32 elements to packed in the bottom of the vector. – Peter Cordes Jan 28 '19 at 23:28
  • 2
    More importantly, @mark's actual problem includes rounding a `double` to nearest integer, with the result as a `double`. (Actually to get the fractional part, I think, by subtracting the integer part from the original (after some scaling?)). Converting to integer and back will work, though, because in this case it only needs to work for small-magnitude doubles. It's fine that doubles outside the -2^31 + 2^31-1 range will convert to integer values of `0x80000000` (INT_MIN). But it's not fine for the general case, so this question should ask for everything you actually need + details. – Peter Cordes Jan 28 '19 at 23:32
  • @PeterCordes: you are actually right (as always :)). `(int)floor(d + 0.5)` seems what I'm after! Why bad emulation? – markzzz Jan 29 '19 at 08:47