/**
* Copyright: Copyright Auburn Sounds 2016-2019.
* License:   $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
* Authors:   Guillaume Piolat
*/
module inteli.xmmintrin;

public import inteli.types;

import inteli.internals;

import inteli.mmx;

import core.stdc.stdlib: malloc, free;
import core.exception: onOutOfMemoryError;

version(D_InlineAsm_X86)
    version = InlineX86Asm;
else version(D_InlineAsm_X86_64)
    version = InlineX86Asm;


// SSE1

nothrow @nogc:


enum int _MM_EXCEPT_INVALID    = 0x0001;
enum int _MM_EXCEPT_DENORM     = 0x0002;
enum int _MM_EXCEPT_DIV_ZERO   = 0x0004;
enum int _MM_EXCEPT_OVERFLOW   = 0x0008;
enum int _MM_EXCEPT_UNDERFLOW  = 0x0010;
enum int _MM_EXCEPT_INEXACT    = 0x0020;
enum int _MM_EXCEPT_MASK       = 0x003f;

enum int _MM_MASK_INVALID      = 0x0080;
enum int _MM_MASK_DENORM       = 0x0100;
enum int _MM_MASK_DIV_ZERO     = 0x0200;
enum int _MM_MASK_OVERFLOW     = 0x0400;
enum int _MM_MASK_UNDERFLOW    = 0x0800;
enum int _MM_MASK_INEXACT      = 0x1000;
enum int _MM_MASK_MASK         = 0x1f80;

enum int _MM_ROUND_NEAREST     = 0x0000;
enum int _MM_ROUND_DOWN        = 0x2000;
enum int _MM_ROUND_UP          = 0x4000;
enum int _MM_ROUND_TOWARD_ZERO = 0x6000;
enum int _MM_ROUND_MASK        = 0x6000;

enum int _MM_FLUSH_ZERO_MASK   = 0x8000;
enum int _MM_FLUSH_ZERO_ON     = 0x8000;
enum int _MM_FLUSH_ZERO_OFF    = 0x0000;

__m128 _mm_add_ps(__m128 a, __m128 b) pure @safe
{
    return a + b;
}

unittest
{
    __m128 a = [1, 2, 3, 4];
    a = _mm_add_ps(a, a);
    assert(a.array[0] == 2);
    assert(a.array[1] == 4);
    assert(a.array[2] == 6);
    assert(a.array[3] == 8);
}

__m128 _mm_add_ss(__m128 a, __m128 b) pure @safe
{
    a[0] += b[0];
    return a;
}
unittest
{
    __m128 a = [1, 2, 3, 4];
    a = _mm_add_ss(a, a);
    assert(a.array == [2.0f, 2, 3, 4]);
}

__m128 _mm_and_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128)(cast(__m128i)a & cast(__m128i)b);
}
unittest
{
    // Note: tested in emmintrin.d
}

__m128 _mm_andnot_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128)( (~cast(__m128i)a) & cast(__m128i)b );
}


// TODO: _mm_avg_pu16
// TODO: _mm_avg_pu8

__m128 _mm_cmpeq_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.oeq)(a, b);
}

__m128 _mm_cmpeq_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.oeq)(a, b);
}

__m128 _mm_cmpge_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.oge)(a, b);
}

__m128 _mm_cmpge_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.oge)(a, b);
}

__m128 _mm_cmpgt_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ogt)(a, b);
}

__m128 _mm_cmpgt_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ogt)(a, b);
}

__m128 _mm_cmple_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ole)(a, b);
}

__m128 _mm_cmple_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ole)(a, b);
}

__m128 _mm_cmplt_ps (__m128 a, __m128 b) pure @safe 
{
    return cast(__m128) cmpps!(FPComparison.olt)(a, b);
}

__m128 _mm_cmplt_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.olt)(a, b);
}

__m128 _mm_cmpneq_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.une)(a, b);
}

__m128 _mm_cmpneq_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.une)(a, b);
}

__m128 _mm_cmpnge_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ult)(a, b);
}

__m128 _mm_cmpnge_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ult)(a, b);
}

__m128 _mm_cmpngt_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ule)(a, b);
}

__m128 _mm_cmpngt_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ule)(a, b);
}

__m128 _mm_cmpnle_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ugt)(a, b);
}

__m128 _mm_cmpnle_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ugt)(a, b);
}

__m128 _mm_cmpnlt_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.uge)(a, b);
}

__m128 _mm_cmpnlt_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.uge)(a, b);
}

__m128 _mm_cmpord_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.ord)(a, b);
}

__m128 _mm_cmpord_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.ord)(a, b);
}

__m128 _mm_cmpunord_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpps!(FPComparison.uno)(a, b);
}

__m128 _mm_cmpunord_ss (__m128 a, __m128 b) pure @safe
{
    return cast(__m128) cmpss!(FPComparison.uno)(a, b);
}

// Note: we've reverted clang and GCC behaviour with regards to EFLAGS
// Some such comparisons yields true for NaNs, other don't.

int _mm_comieq_ss (__m128 a, __m128 b) pure @safe // comiss + sete
{
    return comss!(FPComparison.ueq)(a, b); // yields true for NaN!
}

int _mm_comige_ss (__m128 a, __m128 b) pure @safe // comiss + setae
{
    return comss!(FPComparison.oge)(a, b); 
}

int _mm_comigt_ss (__m128 a, __m128 b) pure @safe // comiss + seta
{
    return comss!(FPComparison.ogt)(a, b);
}

int _mm_comile_ss (__m128 a, __m128 b) pure @safe // comiss + setbe
{
    return comss!(FPComparison.ule)(a, b); // yields true for NaN!
}

int _mm_comilt_ss (__m128 a, __m128 b) pure @safe // comiss + setb
{
    return comss!(FPComparison.ult)(a, b); // yields true for NaN!
}

int _mm_comineq_ss (__m128 a, __m128 b) pure @safe // comiss + setne
{
    return comss!(FPComparison.one)(a, b);
}


// TODO: __m128 _mm_cvt_pi2ps (__m128 a, __m64 b)
// TODO: __m64 _mm_cvt_ps2pi (__m128 a)


__m128 _mm_cvt_si2ss(__m128 v, int x) pure @safe
{
    v[0] = cast(float)x;
    return v;
}
unittest
{
    __m128 a = _mm_cvt_si2ss(_mm_set1_ps(0.0f), 42);
    assert(a.array == [42f, 0, 0, 0]);
}

// Note: is just another name for _mm_cvtss_si32
alias _mm_cvt_ss2si = _mm_cvtss_si32;


// TODO: __m128 _mm_cvtpi16_ps (__m64 a)
// TODO: __m128 _mm_cvtpi32_ps (__m128 a, __m64 b)
// TODO: __m128 _mm_cvtpi32x2_ps (__m64 a, __m64 b)
// TODO: __m128 _mm_cvtpi8_ps (__m64 a)
// TODO: __m64 _mm_cvtps_pi16 (__m128 a)
// TODO: __m64 _mm_cvtps_pi32 (__m128 a)
// TODO: __m64 _mm_cvtps_pi8 (__m128 a)
// TODO: __m128 _mm_cvtpu16_ps (__m64 a)
// TODO: __m128 _mm_cvtpu8_ps (__m64 a)

__m128 _mm_cvtsi32_ss(__m128 v, int x) pure @safe
{
    v[0] = cast(float)x;
    return v;
}
unittest
{
    __m128 a = _mm_cvtsi32_ss(_mm_set1_ps(0.0f), 42);
    assert(a.array == [42.0f, 0, 0, 0]);
}

// Note: on macOS, using "llvm.x86.sse.cvtsi642ss" was buggy
__m128 _mm_cvtsi64_ss(__m128 v, long x) pure @safe
{
    v[0] = cast(float)x;
    return v;
}
unittest
{
    __m128 a = _mm_cvtsi64_ss(_mm_set1_ps(0.0f), 42);
    assert(a.array == [42.0f, 0, 0, 0]);
}

float _mm_cvtss_f32(__m128 a) pure @safe
{
    return a[0];
}

version(LDC)
{
    alias _mm_cvtss_si32 = __builtin_ia32_cvtss2si;
}
else
{
    int _mm_cvtss_si32 (__m128 a) pure @safe
    {
        return convertFloatToInt32UsingMXCSR(a[0]);
    }
}
unittest
{
    assert(1 == _mm_cvtss_si32(_mm_setr_ps(1.0f, 2.0f, 3.0f, 4.0f)));
}

version(LDC)
{
    version(X86_64)
        alias _mm_cvtss_si64 = __builtin_ia32_cvtss2si64;
    else
    {
        // Note: __builtin_ia32_cvtss2si64 crashes LDC in 32-bit
        long _mm_cvtss_si64 (__m128 a) pure @safe
        {
            return convertFloatToInt64UsingMXCSR(a[0]);
        }
    }
}
else
{
    long _mm_cvtss_si64 (__m128 a) pure @safe
    {
        return convertFloatToInt64UsingMXCSR(a[0]);
    }
}
unittest
{
    assert(1 == _mm_cvtss_si64(_mm_setr_ps(1.0f, 2.0f, 3.0f, 4.0f)));

    uint savedRounding = _MM_GET_ROUNDING_MODE();

    _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
    assert(-86186 == _mm_cvtss_si64(_mm_set1_ps(-86186.5f)));

    _MM_SET_ROUNDING_MODE(_MM_ROUND_DOWN);
    assert(-86187 == _mm_cvtss_si64(_mm_set1_ps(-86186.1f)));

    _MM_SET_ROUNDING_MODE(_MM_ROUND_UP);
    assert(86187 == _mm_cvtss_si64(_mm_set1_ps(86186.1f)));

    _MM_SET_ROUNDING_MODE(_MM_ROUND_TOWARD_ZERO);
    assert(-86186 == _mm_cvtss_si64(_mm_set1_ps(-86186.9f)));

    _MM_SET_ROUNDING_MODE(savedRounding);
}


version(LDC)
{
    alias _mm_cvtt_ss2si = __builtin_ia32_cvttss2si;
}
else
{
    int _mm_cvtt_ss2si (__m128 a) pure @safe
    {
        return cast(int)(a[0]);
    }
}
unittest
{
    assert(1 == _mm_cvtt_ss2si(_mm_setr_ps(1.9f, 2.0f, 3.0f, 4.0f)));
}

// TODO: __m64 _mm_cvtt_ps2pi (__m128 a)

alias _mm_cvttss_si32 = _mm_cvtt_ss2si; // it's actually the same op

// Note: __builtin_ia32_cvttss2si64 crashes LDC when generating 32-bit x86 code.
long _mm_cvttss_si64 (__m128 a) pure @safe
{
    return cast(long)(a[0]); // Generates cvttss2si as expected
}
unittest
{
    assert(1 == _mm_cvttss_si64(_mm_setr_ps(1.9f, 2.0f, 3.0f, 4.0f)));
}

__m128 _mm_div_ps(__m128 a, __m128 b) pure @safe
{
    return a / b;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_div_ps(a, a);
    float[4] correct = [1.0f, 1.0f, 1.0f, 1.0f];
    assert(a.array == correct);
}

__m128 _mm_div_ss(__m128 a, __m128 b) pure @safe
{
    a[0] /= b[0];
    return a;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_div_ss(a, a);
    float[4] correct = [1.0f, -2.0, 3.0f, 1.0f];
    assert(a.array == correct);
}

// TODO: int _mm_extract_pi16 (__m64 a, int imm8)

/// Free aligned memory that was allocated with `_mm_malloc`.
void _mm_free(void * mem_addr) @trusted
{
    // support for free(NULL)
    if (mem_addr is null)
        return;

    // Technically we don't need to store size and alignement in the chunk, but we do in case we 
    // have to implement _mm_realloc

    size_t pointerSize = (void*).sizeof;
    void** rawLocation = cast(void**)(cast(char*)mem_addr - size_t.sizeof);
    size_t* alignmentLocation = cast(size_t*)(cast(char*)mem_addr - 3 * pointerSize);
    size_t alignment = *alignmentLocation;
    assert(alignment != 0);
    assert(isPointerAligned(mem_addr, alignment));
    free(*rawLocation);
}

uint _MM_GET_EXCEPTION_MASK() pure @safe
{
    return _mm_getcsr() & _MM_MASK_MASK;
}

uint _MM_GET_EXCEPTION_STATE() pure @safe
{
    return _mm_getcsr() & _MM_EXCEPT_MASK;
}

uint _MM_GET_FLUSH_ZERO_MODE() pure @safe
{
    return _mm_getcsr() & _MM_FLUSH_ZERO_MASK;
}

uint _MM_GET_ROUNDING_MODE() pure @safe
{
    return _mm_getcsr() & _MM_ROUND_MASK;
}

uint _mm_getcsr() pure @safe
{
    version (InlineX86Asm)
    {
        uint controlWord;
        asm nothrow @nogc pure @safe
        {
            stmxcsr controlWord;
        }
        return controlWord;
    }
    else
        static assert(0, "Not yet supported");
}

// TODO: __m64 _mm_insert_pi16 (__m64 a, int i, int imm8)

__m128 _mm_load_ps(const(float)*p) pure @trusted
{
    return *cast(__m128*)p;
}

__m128 _mm_load_ps1(const(float)*p) pure @trusted
{
    float[4] f = [ *p, *p, *p, *p ];
    return loadUnaligned!(float4)(f.ptr);
}

__m128 _mm_load_ss (const(float)* mem_addr) pure @trusted
{
    float[4] f = [ *mem_addr, 0.0f, 0.0f, 0.0f ];
    return loadUnaligned!(float4)(f.ptr);
}

alias _mm_load1_ps = _mm_load_ps1;

__m128 _mm_loadh_pi (__m128 a, const(__m64)* mem_addr) pure @safe
{
    long2 la = cast(long2)a;
    la[1] = (*mem_addr)[0];
    return cast(__m128)la;
}

__m128 _mm_loadl_pi (__m128 a, const(__m64)* mem_addr) pure @safe
{
    long2 la = cast(long2)a;
    la[0] = (*mem_addr)[0];
    return cast(__m128)la;
}

__m128 _mm_loadr_ps (const(float)* mem_addr) pure @trusted
{
    __m128* aligned = cast(__m128*)mem_addr;
    __m128 a = *aligned;
    return shufflevector!(__m128, 3, 2, 1, 0)(a, a);
}

__m128 _mm_loadu_ps(const(float)*p) pure @safe
{
    return loadUnaligned!(__m128)(p);
}

__m128i _mm_loadu_si16(const(void)* mem_addr) pure @trusted
{
    short r = *cast(short*)(mem_addr);
    short8 result = [0, 0, 0, 0, 0, 0, 0, 0];
    result[0] = r;
    return cast(__m128i)result;
}
unittest
{
    short r = 13;
    short8 A = cast(short8) _mm_loadu_si16(&r);
    short[8] correct = [13, 0, 0, 0, 0, 0, 0, 0];
    assert(A.array == correct);
}

__m128i _mm_loadu_si64(const(void)* mem_addr) pure @trusted
{
    long r = *cast(long*)(mem_addr);
    long2 result = [0, 0];
    result[0] = r;
    return cast(__m128i)result;
}
unittest
{
    long r = 446446446446;
    long2 A = cast(long2) _mm_loadu_si64(&r);
    long[2] correct = [446446446446, 0];
    assert(A.array == correct);
}

/// Allocate size bytes of memory, aligned to the alignment specified in align, 
/// and return a pointer to the allocated memory. `_mm_free` should be used to free 
/// memory that is allocated with `_mm_malloc`.
void* _mm_malloc(size_t size, size_t alignment) @trusted
{
    assert(alignment != 0);
    size_t request = requestedSize(size, alignment);
    void* raw = malloc(request);
    if (request > 0 && raw == null) // malloc(0) can validly return anything
        onOutOfMemoryError();
    return storeRawPointerPlusInfo(raw, size, alignment); // PERF: no need to store size
}

// TODO: _mm_maskmove_si64
// TODO: _m_maskmovq

// TODO: _mm_max_pi16
version(LDC)
{
    alias _mm_max_ps = __builtin_ia32_maxps;
}
else
{
    __m128 _mm_max_ps(__m128 a, __m128 b) pure @safe
    {
        __m128 r;
        r[0] = (a[0] > b[0]) ? a[0] : b[0];
        r[1] = (a[1] > b[1]) ? a[1] : b[1];
        r[2] = (a[2] > b[2]) ? a[2] : b[2];
        r[3] = (a[3] > b[3]) ? a[3] : b[3];
        return r;
    }
}
unittest
{
    __m128 A = _mm_setr_ps(1, 2, float.nan, 4);
    __m128 B = _mm_setr_ps(4, 1, 4, float.nan);
    __m128 M = _mm_max_ps(A, B);
    assert(M[0] == 4);
    assert(M[1] == 2);
    assert(M[2] == 4);    // in case of NaN, second operand prevails (as it seems)
    assert(M[3] != M[3]); // in case of NaN, second operand prevails (as it seems)
}

// TODO: _mm_max_pu8

version(LDC)
{
    alias _mm_max_ss = __builtin_ia32_maxss;
}
else
{
    __m128 _mm_max_ss(__m128 a, __m128 b) pure @safe
    {
        __m128 r = a;
        r[0] = (a[0] > b[0]) ? a[0] : b[0];
        return r;
    }
}
unittest
{
    __m128 A = _mm_setr_ps(1, 2, 3, 4);
    __m128 B = _mm_setr_ps(4, 1, 4, 1);
    __m128 C = _mm_setr_ps(float.nan, 1, 4, 1);
    __m128 M = _mm_max_ss(A, B);
    assert(M[0] == 4);
    assert(M[1] == 2);
    assert(M[2] == 3);
    assert(M[3] == 4);
    M = _mm_max_ps(A, C); // in case of NaN, second operand prevails
    assert(M[0] != M[0]);
    M = _mm_max_ps(C, A); // in case of NaN, second operand prevails
    assert(M[0] == 1);
}

// TODO: _mm_min_pi16

version(LDC)
{
    alias _mm_min_ps = __builtin_ia32_minps;
}
else
{
    __m128 _mm_min_ps(__m128 a, __m128 b) pure @safe
    {
        __m128 r;
        r[0] = (a[0] < b[0]) ? a[0] : b[0];
        r[1] = (a[1] < b[1]) ? a[1] : b[1];
        r[2] = (a[2] < b[2]) ? a[2] : b[2];
        r[3] = (a[3] < b[3]) ? a[3] : b[3];
        return r;
    }
}
unittest
{
    __m128 A = _mm_setr_ps(1, 2, float.nan, 4);
    __m128 B = _mm_setr_ps(4, 1, 4, float.nan);
    __m128 M = _mm_min_ps(A, B);
    assert(M[0] == 1);
    assert(M[1] == 1);
    assert(M[2] == 4);    // in case of NaN, second operand prevails (as it seems)
    assert(M[3] != M[3]); // in case of NaN, second operand prevails (as it seems)
}

// TODO: _mm_min_pi8

version(LDC)
{
    alias _mm_min_ss = __builtin_ia32_minss;
}
else
{
    __m128 _mm_min_ss(__m128 a, __m128 b) pure @safe
    {
        __m128 r = a;
        r[0] = (a[0] < b[0]) ? a[0] : b[0];
        return r;
    }
}
unittest
{
    __m128 A = _mm_setr_ps(1, 2, 3, 4);
    __m128 B = _mm_setr_ps(4, 1, 4, 1);
    __m128 C = _mm_setr_ps(float.nan, 1, 4, 1);
    __m128 M = _mm_min_ss(A, B);
    assert(M[0] == 1);
    assert(M[1] == 2);
    assert(M[2] == 3);
    assert(M[3] == 4);
    M = _mm_min_ps(A, C); // in case of NaN, second operand prevails
    assert(M[0] != M[0]);
    M = _mm_min_ps(C, A); // in case of NaN, second operand prevails
    assert(M[0] == 1);
}

__m128 _mm_move_ss (__m128 a, __m128 b) pure @safe
{
    return shufflevector!(__m128, 4, 1, 2, 3)(a, b);
}

__m128 _mm_movehl_ps (__m128 a, __m128 b) pure @safe
{
    return shufflevector!(float4, 2, 3, 6, 7)(a, b);
}

__m128 _mm_movelh_ps (__m128 a, __m128 b) pure @safe
{
    return shufflevector!(float4, 0, 1, 4, 5)(a, b);
}


version(LDC)
{
    alias _mm_movemask_ps = __builtin_ia32_movmskps;
}
else
{
    int _mm_movemask_ps (__m128 a) pure @safe
    {
        int4 ai = cast(int4)a;
        int r = 0;
        if (ai[0] < 0) r += 1;
        if (ai[1] < 0) r += 2;
        if (ai[2] < 0) r += 4;
        if (ai[3] < 0) r += 8;
        return r;
    }
}
unittest
{
    int4 A = [-1, 0, -43, 0];
    assert(5 == _mm_movemask_ps(cast(float4)A));
}

__m128 _mm_mul_ps(__m128 a, __m128 b) pure @safe
{
    return a * b;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_mul_ps(a, a);
    float[4] correct = [2.25f, 4.0f, 9.0f, 1.0f];
    assert(a.array == correct);
}

__m128 _mm_mul_ss(__m128 a, __m128 b) pure @safe
{
    a[0] *= b[0];
    return a;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_mul_ss(a, a);
    float[4] correct = [2.25f, -2.0f, 3.0f, 1.0f];
    assert(a.array == correct);
}

// TODO: _mm_mulhi_pu16

__m128 _mm_or_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128)(cast(__m128i)a | cast(__m128i)b);
}

// TODO: __m64 _m_pavgb (__m64 a, __m64 b)
// TODO: __m64 _m_pavgw (__m64 a, __m64 b)
// TODO: int _m_pextrw (__m64 a, int imm8)
// TODO: __m64 _m_pinsrw (__m64 a, int i, int imm8)
// TODO: __m64 _m_pmaxsw (__m64 a, __m64 b)
// TODO: __m64 _m_pmaxub (__m64 a, __m64 b)
// TODO: __m64 _m_pminsw (__m64 a, __m64 b)
// TODO: __m64 _m_pminub (__m64 a, __m64 b)
// TODO: int _m_pmovmskb (__m64 a)

// TODO: __m64 _m_pmulhuw (__m64 a, __m64 b)

enum _MM_HINT_NTA = 0;
enum _MM_HINT_T0 = 1;
enum _MM_HINT_T1 = 2;
enum _MM_HINT_T2 = 3;

// Note: locality must be compile-time, unlike Intel Intrinsics API
void _mm_prefetch(int locality)(void* p) pure @safe
{
    llvm_prefetch(p, 0, locality, 1);
}

// TODO: __m64 _m_psadbw (__m64 a, __m64 b)
// TODO: __m64 _m_pshufw (__m64 a, int imm8)

version(LDC)
{
    alias _mm_rcp_ps = __builtin_ia32_rcpps;
}
else
{
    __m128 _mm_rcp_ps (__m128 a) pure @safe
    {
        a[0] = 1.0f / a[0];
        a[1] = 1.0f / a[1];
        a[2] = 1.0f / a[2];
        a[3] = 1.0f / a[3];
        return a;
    }
}

version(LDC)
{
    alias _mm_rcp_ss = __builtin_ia32_rcpss;
}
else
{
    __m128 _mm_rcp_ss (__m128 a) pure @safe
    {
        a[0] = 1.0f / a[0];
        return a;
    }
}

version(LDC)
{
    alias _mm_rsqrt_ps = __builtin_ia32_rsqrtps;
}
else
{
    __m128 _mm_rsqrt_ps (__m128 a) pure @safe
    {
        a[0] = 1.0f / sqrt(a[0]);
        a[1] = 1.0f / sqrt(a[1]);
        a[2] = 1.0f / sqrt(a[2]);
        a[3] = 1.0f / sqrt(a[3]);
        return a;
    }
}

version(LDC)
{
    alias _mm_rsqrt_ss = __builtin_ia32_rsqrtss;
}
else
{
    __m128 _mm_rsqrt_ss (__m128 a) pure @safe
    {
        a[0] = 1.0f / sqrt(a[0]);
        return a;
    }
}

unittest
{
    double maxRelativeError = 0.000245; // -72 dB
    void testInvSqrt(float number)
    {
        __m128 A = _mm_set1_ps(number);

        // test _mm_rcp_ps
        __m128 B = _mm_rcp_ps(A);
        foreach(i; 0..4)
        {
            double exact = 1.0f / A[i]; 
            double ratio = cast(double)(B[i]) / cast(double)(exact);
            assert(fabs(ratio - 1) <= maxRelativeError);
        }

        // test _mm_rcp_ss
        {
            B = _mm_rcp_ss(A);
            double exact = 1.0f / A[0];
            double ratio = cast(double)(B[0]) / cast(double)(exact);
            assert(fabs(ratio - 1) <= maxRelativeError);
        }

        // test _mm_rsqrt_ps
        B = _mm_rsqrt_ps(A);
        foreach(i; 0..4)
        {
            double exact = 1.0f / sqrt(A[i]);
            double ratio = cast(double)(B[i]) / cast(double)(exact);
            assert(fabs(ratio - 1) <= maxRelativeError);
        }

        // test _mm_rsqrt_ss
        {
            B = _mm_rsqrt_ss(A);
            double exact = 1.0f / sqrt(A[0]);
            double ratio = cast(double)(B[0]) / cast(double)(exact);
            assert(fabs(ratio - 1) <= maxRelativeError);
        }
    }

    testInvSqrt(1.1f);
    testInvSqrt(2.45674864151f);
    testInvSqrt(27841456468.0f);
}

// TODO: _mm_sad_pu8

void _MM_SET_EXCEPTION_MASK(int _MM_MASK_xxxx) pure @safe
{
    _mm_setcsr((_mm_getcsr() & ~_MM_MASK_MASK) | _MM_MASK_xxxx);
}

void _MM_SET_EXCEPTION_STATE(int _MM_EXCEPT_xxxx) pure @safe
{
    _mm_setcsr((_mm_getcsr() & ~_MM_EXCEPT_MASK) | _MM_EXCEPT_xxxx);
}

void _MM_SET_FLUSH_ZERO_MODE(int _MM_FLUSH_xxxx) pure @safe
{
    _mm_setcsr((_mm_getcsr() & ~_MM_FLUSH_ZERO_MASK) | _MM_FLUSH_xxxx);
}

__m128 _mm_set_ps (float e3, float e2, float e1, float e0) pure @trusted
{
    // Note: despite appearances, generates sensible code,
    //       inlines correctly and is constant folded
    float[4] result = [e0, e1, e2, e3];
    return loadUnaligned!(float4)(result.ptr); 
}

alias _mm_set_ps1 = _mm_set1_ps;

void _MM_SET_ROUNDING_MODE(int _MM_ROUND_xxxx) pure @safe
{
    _mm_setcsr((_mm_getcsr() & ~_MM_ROUND_MASK) | _MM_ROUND_xxxx);
}

__m128 _mm_set_ss (float a) pure @trusted
{
    __m128 r = _mm_setzero_ps();
    r[0] = a;
    return r;
}

__m128 _mm_set1_ps (float a) pure @trusted
{
    // Note: despite appearances, generates sensible code,
    //       inlines correctly and is constant folded
    float[4] result = [a, a, a, a];
    return loadUnaligned!(float4)(result.ptr);
}

void _mm_setcsr(uint controlWord) pure @safe
{
    version (InlineX86Asm)
    {
        asm pure nothrow @nogc @safe
        { 
            ldmxcsr controlWord;
        }
    }
    else
        static assert(0, "Not yet supported");
}

__m128 _mm_setr_ps (float e3, float e2, float e1, float e0) pure @trusted
{
    float[4] result = [e3, e2, e1, e0];
    return loadUnaligned!(float4)(result.ptr);
}

__m128 _mm_setzero_ps() pure @trusted
{
    // Compiles to xorps without problems
    float[4] result = [0.0f, 0.0f, 0.0f, 0.0f];
    return loadUnaligned!(float4)(result.ptr);
}

version(LDC)
{
    alias _mm_sfence = __builtin_ia32_sfence;
}
else
{
    void _mm_sfence() pure @safe
    {
        asm nothrow @nogc pure @safe
        {
            sfence;
        }
    }
}
unittest
{
    _mm_sfence();
}

// TODO: mm_shuffle_pi16

// Note: the immediate shuffle value is given at compile-time instead of runtime.
__m128 _mm_shuffle_ps(ubyte imm)(__m128 a, __m128 b) pure @safe
{
    return shufflevector!(__m128, imm & 3, (imm>>2) & 3, 4 + ((imm>>4) & 3), 4 + ((imm>>6) & 3) )(a, b);
}

version(LDC)
{
    // Disappeared with LDC 1.11
    static if (__VERSION__ < 2081)
        alias _mm_sqrt_ps = __builtin_ia32_sqrtps;
    else
    {
        __m128 _mm_sqrt_ps(__m128 vec) pure @safe
        {
            vec.array[0] = llvm_sqrt(vec.array[0]);
            vec.array[1] = llvm_sqrt(vec.array[1]);
            vec.array[2] = llvm_sqrt(vec.array[2]);
            vec.array[3] = llvm_sqrt(vec.array[3]);
            return vec;
        }
    }
}
else
{
    __m128 _mm_sqrt_ps(__m128 vec) pure @safe
    {
        vec.array[0] = sqrt(vec.array[0]);
        vec.array[1] = sqrt(vec.array[1]);
        vec.array[2] = sqrt(vec.array[2]);
        vec.array[3] = sqrt(vec.array[3]);
        return vec;
    }
}
unittest
{
    __m128 A = _mm_sqrt_ps(_mm_set1_ps(4.0f));
    assert(A.array[0] == 2.0f);
    assert(A.array[1] == 2.0f);
    assert(A.array[2] == 2.0f);
    assert(A.array[3] == 2.0f);
}

version(LDC)
{
    // Disappeared with LDC 1.11
    static if (__VERSION__ < 2081)
        alias _mm_sqrt_ss = __builtin_ia32_sqrtss;
    else
    {
        __m128 _mm_sqrt_ss(__m128 vec) pure @safe
        {
            vec.array[0] = llvm_sqrt(vec.array[0]);
            vec.array[1] = vec.array[1];
            vec.array[2] = vec.array[2];
            vec.array[3] = vec.array[3];
            return vec;
        }
    }
}
else
{
    __m128 _mm_sqrt_ss(__m128 vec) pure @safe
    {
        vec.array[0] = sqrt(vec.array[0]);
        return vec;
    }
}
unittest
{
    __m128 A = _mm_sqrt_ss(_mm_set1_ps(4.0f));
    assert(A.array[0] == 2.0f);
    assert(A.array[1] == 4.0f);
    assert(A.array[2] == 4.0f);
    assert(A.array[3] == 4.0f);
}

void _mm_store_ps (float* mem_addr, __m128 a) pure // not safe since nothing guarantees alignment
{
    __m128* aligned = cast(__m128*)mem_addr;
    *aligned = a;
}

alias _mm_store_ps1 = _mm_store1_ps;

void _mm_store_ss (float* mem_addr, __m128 a) pure @safe
{
    *mem_addr = a[0];
}

void _mm_store1_ps (float* mem_addr, __m128 a) pure // not safe since nothing guarantees alignment
{
    __m128* aligned = cast(__m128*)mem_addr;
    *aligned = shufflevector!(__m128, 0, 0, 0, 0)(a, a);
}

void _mm_storeh_pi(__m64* p, __m128 a) pure @safe
{
    long2 la = cast(long2)a;
    (*p)[0] = la[1];
}
unittest
{
    __m64 R = _mm_setzero_si64();
    long2 A = [13, 25];
    _mm_storeh_pi(&R, cast(__m128)A);
    assert(R[0] == 25);
}

void _mm_storel_pi(__m64* p, __m128 a) pure @safe
{
    long2 la = cast(long2)a;
    (*p)[0] = la[0];
}
unittest
{
    __m64 R = _mm_setzero_si64();
    long2 A = [13, 25];
    _mm_storel_pi(&R, cast(__m128)A);
    assert(R[0] == 13);
}

void _mm_storer_ps(float* mem_addr, __m128 a) pure // not safe since nothing guarantees alignment
{
    __m128* aligned = cast(__m128*)mem_addr;
    *aligned = shufflevector!(__m128, 3, 2, 1, 0)(a, a);
}

void _mm_storeu_ps(float* mem_addr, __m128 a) pure @safe
{
    storeUnaligned!(float4)(a, mem_addr);
}

// TODO: _mm_stream_pi, does not seem possible

// BUG: can't implement non-temporal store with LDC inlineIR since !nontemporal
// needs some IR outside this function that would say: 
//
//  !0 = !{ i32 1 }
//
// It's a LLVM IR metadata description.
// Regardless, non-temporal moves are really dangerous for performance...
void _mm_stream_ps (float* mem_addr, __m128 a)
{
    __m128* dest = cast(__m128*)mem_addr;
    *dest = a; // it's a regular move instead
}
unittest
{
    align(16) float[4] A;
    _mm_stream_ps(A.ptr, _mm_set1_ps(78.0f));
    assert(A[0] == 78.0f && A[1] == 78.0f && A[2] == 78.0f && A[3] == 78.0f);
}

__m128 _mm_sub_ps(__m128 a, __m128 b) pure @safe
{
    return a - b;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_sub_ps(a, a);
    float[4] correct = [0.0f, 0.0f, 0.0f, 0.0f];
    assert(a.array == correct);
}

__m128 _mm_sub_ss(__m128 a, __m128 b) pure @safe
{
    a[0] -= b[0];
    return a;
}
unittest
{
    __m128 a = [1.5f, -2.0f, 3.0f, 1.0f];
    a = _mm_sub_ss(a, a);
    float[4] correct = [0.0f, -2.0, 3.0f, 1.0f];
    assert(a.array == correct);
}


void _MM_TRANSPOSE4_PS (ref __m128 row0, ref __m128 row1, ref __m128 row2, ref __m128 row3) pure @safe
{
    __m128 tmp3, tmp2, tmp1, tmp0;
    tmp0 = _mm_unpacklo_ps(row0, row1);
    tmp2 = _mm_unpacklo_ps(row2, row3);
    tmp1 = _mm_unpackhi_ps(row0, row1);
    tmp3 = _mm_unpackhi_ps(row2, row3);
    row0 = _mm_movelh_ps(tmp0, tmp2);
    row1 = _mm_movehl_ps(tmp2, tmp0);
    row2 = _mm_movelh_ps(tmp1, tmp3);
    row3 = _mm_movehl_ps(tmp3, tmp1);
}

// Note: the only difference between these intrinsics is the signalling 
//       behaviour of quiet NaNs. This is incorrect but the case where
//       you would want to differentiate between qNaN and sNaN and then 
//       treat them differently on purpose seems extremely rare.
alias _mm_ucomieq_ss = _mm_comieq_ss;
alias _mm_ucomige_ss = _mm_comige_ss;
alias _mm_ucomigt_ss = _mm_comigt_ss;
alias _mm_ucomile_ss = _mm_comile_ss;
alias _mm_ucomilt_ss = _mm_comilt_ss;
alias _mm_ucomineq_ss = _mm_comineq_ss;


__m128 _mm_undefined_ps() pure @safe
{
    __m128 undef = void;
    return undef;
}

__m128 _mm_unpackhi_ps (__m128 a, __m128 b) pure @safe
{
    return shufflevector!(float4, 2, 6, 3, 7)(a, b);
}

__m128 _mm_unpacklo_ps (__m128 a, __m128 b) pure @safe
{
    return shufflevector!(float4, 0, 4, 1, 5)(a, b);
}

__m128 _mm_xor_ps (__m128 a, __m128 b) pure @safe
{
    return cast(__m128)(cast(__m128i)a ^ cast(__m128i)b);
}


private
{
    /// Returns: `true` if the pointer is suitably aligned.
    bool isPointerAligned(void* p, size_t alignment) pure
    {
        assert(alignment != 0);
        return ( cast(size_t)p & (alignment - 1) ) == 0;
    }

    /// Returns: next pointer aligned with alignment bytes.
    void* nextAlignedPointer(void* start, size_t alignment) pure
    {
        return cast(void*)nextMultipleOf(cast(size_t)(start), alignment);
    }

    // Returns number of bytes to actually allocate when asking
    // for a particular alignment
    @nogc size_t requestedSize(size_t askedSize, size_t alignment) pure
    {
        enum size_t pointerSize = size_t.sizeof;
        return askedSize + alignment - 1 + pointerSize * 3;
    }

    // Store pointer given my malloc, size and alignment
    @nogc void* storeRawPointerPlusInfo(void* raw, size_t size, size_t alignment) pure
    {
        enum size_t pointerSize = size_t.sizeof;
        char* start = cast(char*)raw + pointerSize * 3;
        void* aligned = nextAlignedPointer(start, alignment);
        void** rawLocation = cast(void**)(cast(char*)aligned - pointerSize);
        *rawLocation = raw;
        size_t* sizeLocation = cast(size_t*)(cast(char*)aligned - 2 * pointerSize);
        *sizeLocation = size;
        size_t* alignmentLocation = cast(size_t*)(cast(char*)aligned - 3 * pointerSize);
        *alignmentLocation = alignment;
        assert( isPointerAligned(aligned, alignment) );
        return aligned;
    }

    // Returns: x, multiple of powerOfTwo, so that x >= n.
    @nogc size_t nextMultipleOf(size_t n, size_t powerOfTwo) pure nothrow
    {
        // check power-of-two
        assert( (powerOfTwo != 0) && ((powerOfTwo & (powerOfTwo - 1)) == 0));

        size_t mask = ~(powerOfTwo - 1);
        return (n + powerOfTwo - 1) & mask;
    }
}

unittest
{
    assert(nextMultipleOf(0, 4) == 0);
    assert(nextMultipleOf(1, 4) == 4);
    assert(nextMultipleOf(2, 4) == 4);
    assert(nextMultipleOf(3, 4) == 4);
    assert(nextMultipleOf(4, 4) == 4);
    assert(nextMultipleOf(5, 4) == 8);

    {
        void* p = _mm_malloc(23, 16);
        assert(p !is null);
        assert(((cast(size_t)p) & 0xf) == 0);
        _mm_free(p);
    }

    void* nullAlloc = _mm_malloc(0, 32);
    assert(nullAlloc != null);
    _mm_free(nullAlloc);
}