Build 2^n in double/simd

Question

I am trying to build 2^n using the double representation. The trick is (well) known

// tips to calculate 2^n using the exponent of the double IEEE representation
union ieee754{
    double d;
    uint32_t i[2];
};

// Converts an unsigned long long to a double
inline double uint642dp(uint64_t ll) {
  ieee754 tmp;
  tmp.ll=ll;
  return tmp.d;
}
-----------------------------------
// e^x = 2^n e^y, n = round(x/log(2)+0.5)

double x = 4.3; // start from there
double y = round(x/log(2)+0.5)
int n = (int)y; // convert in to double 
uint642dp(( ((uint64_t)n) +1023)<<52); // calculate 2^n fastly

if n = 4, it will return 16.

Presently I am looking for the same this but for SIMD calculation. Considering SSE2 double, after the round function I will get a register sse2 __m128d v = (4.0, 3.0); form this register how calculate 2^v ... I am blocked mainly due to the cast __m128d to __m128i, it does not exist (it still exists a cast but it does not move bit, just change the "interpretation" of the register double/int).

I do not want return the datas form simd register to normal register to make the transformation. It exists certainly a tips with SIMD, but I do not know it.

So Help ^_^'

Best,

This is probably a pointless optimisation, since int -> float/double conversion on SSE is already very fast and generating 2^n as an int is trivial of course. — Paul R
The "trick" (the union abuse) unfortunately is not valid C++ ... — Kerrek SB

Florent DUGUET Florent DUGUET · Accepted Answer · 2016-05-01T15:37:00

The trick you are looking for is actually the _mm256_castsi256_pd intrinsic which allows you to convert a SIMD-array of integers into a SIMD-array of doubles. This is only for C/C++ type checking and does not translate into any instruction.

Here is a snippet of code to perform this operation (only valid within some bounds of exponent):

#include <immintrin.h>

__m256d pow2n (__m256i n)
{
    const __m256i bias = _mm256_set1_epi64x( 1023 );
    __m256i t = _mm256_add_epi64 (n, bias);
    t = _mm256_slli_epi64 (t, 52);
    return _mm256_castsi256_pd (t) ;
}

#include <cstdio>

int main ()
{
        __m256i rn = _mm256_set_epi64x( 7, 9, 4, 2 );
        __m256d pn = pow2n (rn) ;
        double v [4] ;
        _mm256_storeu_pd (v, pn) ;
        printf ("v = %lf %lf %lf %lf\n", v[0], v[1], v[2], v[3]) ;
        return 0 ;
}

pow2n compiles to only 2 insns, as you can see on the Godbolt Compiler Explorer

Output is :

v = 4.000000 16.000000 512.000000 128.000000

This requires AVX2, but a 128bit version would only require SSE2. Note the exponents are provided as 64bit integers. If you have 32bit integers, use _mm256_cvtepu32_epi64.

If you're starting with a vector of double, like the OP, then use _mm256_cvtpd_epi32 and _mm256_cvtepu32_epi64. (double -> int64 directly isn't available until AVX-512).

Build 2^n in double/simd

1 Answers