a bitonic sorting network for 16 signed short in SSE2

We saw a sorting network for 8 float numbers.  These values fit nicely into two 128bit registers.
Twice as many signed shorts can be fit into the same space.  And SSE2 has intrinsics for maximum and minimum for 8 signed short numbers at a time using _mm_min_epi16 and _mm_max_epi16. 

While the bitonic sorting network for 16 values requires 10 min/max pairs rather than the 9 required in the network shown in Knuth, the shuffles are much easier for the plain vanilla bitonic sorting network.

A lot of shuffling has to happen to make this all work.  The unpack operations will happen in pairs:

void unpack_epi16_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
    __m128i aa = _mm_unpacklo_epi16(a, b);
    __m128i bb = _mm_unpackhi_epi16(a, b);
    aOut = aa;
    bOut = bb;
}
void unpack_epi32_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
    __m128i aa = _mm_unpacklo_epi32(a, b);
    __m128i bb = _mm_unpackhi_epi32(a, b);
    aOut = aa;
    bOut = bb;
}
void unpack_epi64_ab(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
    __m128i aa = _mm_unpacklo_epi64(a, b);
    __m128i bb = _mm_unpackhi_epi64(a, b);
    aOut = aa;
    bOut = bb;
}
Note that it is safe to call all of these with the same variables used as input and output.

void BitonicSortingNetwork_epi16_SSE2(const __m128i &a, const __m128i &b, __m128i &aOut, __m128i &bOut)
{
    __m128i a1 = a;
    __m128i b1 = b;
    __m128i a1Min = _mm_min_epi16(a1, b1);
    __m128i b1Max = _mm_max_epi16(a1, b1);
    __m128i a2 = a1Min;
    __m128i b2 = _mm_shufflelo_epi16(b1Max, _MM_SHUFFLE(2, 3, 0, 1));
    b2 = _mm_shufflehi_epi16(b2, _MM_SHUFFLE(2, 3, 0, 1));
    __m128i a2Min = _mm_min_epi16(a2, b2);
    __m128i b2Max = _mm_max_epi16(a2, b2);
    __m128i a3, b3;
    unpack_epi16_ab(a2Min, b2Max, a3, b3);
    unpack_epi16_ab(a3, b3, a3, b3);
    unpack_epi16_ab(a3, b3, a3, b3);
    __m128i a3Min = _mm_min_epi16(a3, b3);
    __m128i b3Max = _mm_max_epi16(a3, b3);
    __m128i a4, b4;
    unpack_epi32_ab(a3Min, b3Max, a4, b4);
    a4 = _mm_shufflelo_epi16(a4, _MM_SHUFFLE(0, 1, 2, 3));
    b4 = _mm_shufflehi_epi16(b4, _MM_SHUFFLE(0, 1, 2, 3));
    __m128i a4Min = _mm_min_epi16(a4, b4);
    __m128i b4Max = _mm_max_epi16(a4, b4);
    __m128i a5, b5;
    unpack_epi16_ab(a4Min, b4Max, a5, b5);
    unpack_epi64_ab(a5, b5, a5, b5);
    a5 = _mm_shuffle_epi32(a5, _MM_SHUFFLE(2, 3, 0, 1));
    __m128i a5Min = _mm_min_epi16(a5, b5);
    __m128i b5Max = _mm_max_epi16(a5, b5);
    __m128i a6, b6;
    unpack_epi16_ab(a5Min, b5Max, a6, b6);
    unpack_epi16_ab(a6, b6, a6, b6);
    __m128i a6Min = _mm_min_epi16(a6, b6);
    __m128i b6Max = _mm_max_epi16(a6, b6);
    __m128i a7, b7;
    b7 = _mm_shuffle_epi32(b6Max, _MM_SHUFFLE(1, 0, 3, 2));
    b7 = _mm_shufflehi_epi16(b7, _MM_SHUFFLE(0, 1, 2, 3));
    a7 = _mm_shufflelo_epi16(a6Min, _MM_SHUFFLE(0, 1, 2, 3));
    __m128i a7Min = _mm_min_epi16(a7, b7);
    __m128i b7Max = _mm_max_epi16(a7, b7);
    __m128i a8, b8;
    unpack_epi16_ab(a7Min, b7Max, b8, a8);
    a8 = _mm_shuffle_epi32(a8, _MM_SHUFFLE(0, 1, 2, 3));
    __m128i a8Min = _mm_min_epi16(a8, b8);
    __m128i b8Max = _mm_max_epi16(a8, b8);
    __m128i a9, b9;
    unpack_epi16_ab(a8Min, b8Max, a9, b9);
    __m128i a9Min = _mm_min_epi16(a9, b9);
    __m128i b9Max = _mm_max_epi16(a9, b9);
    __m128i a10, b10;
    unpack_epi16_ab(a9Min, b9Max, a10, b10);
    __m128i a10Min = _mm_min_epi16(a10, b10);
    __m128i b10Max = _mm_max_epi16(a10, b10);
    unpack_epi16_ab(a10Min, b10Max, aOut, bOut);
}
This is verified using the standard 0/1 technique described in Knuth.
The compiled code uses only five registers, so avoids register spills.

I am still not sure if this was worth the effort.  But it's done.