///////////////////////////////////////////////////////////////////////////

//

// Copyright (c) 2009-2014 DreamWorks Animation LLC. 

//

// All rights reserved.

//

// Redistribution and use in source and binary forms, with or without

// modification, are permitted provided that the following conditions are

// met:

// *       Redistributions of source code must retain the above copyright

// notice, this list of conditions and the following disclaimer.

// *       Redistributions in binary form must reproduce the above

// copyright notice, this list of conditions and the following disclaimer

// in the documentation and/or other materials provided with the

// distribution.

// *       Neither the name of DreamWorks Animation nor the names of

// its contributors may be used to endorse or promote products derived

// from this software without specific prior written permission.

//

// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT

// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,

// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT

// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,

// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY

// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT

// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

//

///////////////////////////////////////////////////////////////////////////

#ifndef IMF_DWACOMPRESSORSIMD_H_HAS_BEEN_INCLUDED

#define IMF_DWACOMPRESSORSIMD_H_HAS_BEEN_INCLUDED

//

// Various SSE accelerated functions, used by Imf::DwaCompressor. 

// These have been separated into a separate .h file, as the fast

// paths are done with template specialization.

//

// Unless otherwise noted, all pointers are assumed to be 32-byte 

// aligned. Unaligned pointers may risk seg-faulting.

//

#include "ImfNamespace.h"

#include "ImfSimd.h"

#include "ImfSystemSpecific.h"

#include "OpenEXRConfig.h"

#include <half.h>

#include <assert.h>

OPENEXR_IMF_INTERNAL_NAMESPACE_HEADER_ENTER

#define _SSE_ALIGNMENT        32

#define _SSE_ALIGNMENT_MASK 0x0F

#define _AVX_ALIGNMENT_MASK 0x1F

//

// Test if we should enable GCC inline asm paths for AVX

//

#ifdef OPENEXR_IMF_HAVE_GCC_INLINE_ASM_AVX 

    #define IMF_HAVE_GCC_INLINEASM

    #ifdef __LP64__

        #define IMF_HAVE_GCC_INLINEASM_64

    #endif /* __LP64__ */

#endif /* OPENEXR_IMF_HAVE_GCC_INLINE_ASM_AVX */

//

// A simple 64-element array, aligned properly for SIMD access. 

//

template <class T>

class SimdAlignedBuffer64

{

    public:

        SimdAlignedBuffer64(): _buffer (0), _handle (0)           

        {

            alloc();

        }

        SimdAlignedBuffer64(const SimdAlignedBuffer64 &rhs): _handle(0)

        {

            alloc();

            memcpy (_buffer, rhs._buffer, 64 * sizeof (T));

        }

        ~SimdAlignedBuffer64 ()

        {

            EXRFreeAligned (_handle);

            _handle = 0;

            _buffer = 0;

        }

        void alloc()

        {

            //

            // Try EXRAllocAligned first - but it might fallback to

            // unaligned allocs. If so, overalloc.

            //

            _handle = (char *) EXRAllocAligned

                (64 * sizeof(T), _SSE_ALIGNMENT);

            if (((size_t)_handle & (_SSE_ALIGNMENT - 1)) == 0)

            {

                _buffer = (T *)_handle;

                return;

            }

            EXRFreeAligned(_handle);

            _handle = (char *) EXRAllocAligned

                (64 * sizeof(T) + _SSE_ALIGNMENT, _SSE_ALIGNMENT);

            char *aligned = _handle;

            while ((size_t)aligned & (_SSE_ALIGNMENT - 1))

                aligned++;

            _buffer = (T *)aligned;    

        }

        T     *_buffer;

    private:

        char  *_handle;

};

typedef SimdAlignedBuffer64<float>          SimdAlignedBuffer64f;

typedef SimdAlignedBuffer64<unsigned short> SimdAlignedBuffer64us;

namespace {

//

// Color space conversion, Inverse 709 CSC, Y'CbCr -> R'G'B'

//

void

csc709Inverse (float &comp0, float &comp1, float &comp2)

{

    float src[3];

    src[0] = comp0;

    src[1] = comp1;

    src[2] = comp2;

    comp0 = src[0]                    + 1.5747f * src[2];

    comp1 = src[0] - 0.1873f * src[1] - 0.4682f * src[2];

    comp2 = src[0] + 1.8556f * src[1];

}

#ifndef IMF_HAVE_SSE2

//

// Scalar color space conversion, based on 709 primiary chromaticies.

// No scaling or offsets, just the matrix

//

void

csc709Inverse64 (float *comp0, float *comp1, float *comp2)

{

    for (int i = 0; i < 64; ++i)

        csc709Inverse (comp0[i], comp1[i], comp2[i]);

}

#else /* IMF_HAVE_SSE2 */

//

// SSE2 color space conversion

//

void

csc709Inverse64 (float *comp0, float *comp1, float *comp2)

{

    __m128 c0 = { 1.5747f,  1.5747f,  1.5747f,  1.5747f};

    __m128 c1 = { 1.8556f,  1.8556f,  1.8556f,  1.8556f};

    __m128 c2 = {-0.1873f, -0.1873f, -0.1873f, -0.1873f};

    __m128 c3 = {-0.4682f, -0.4682f, -0.4682f, -0.4682f}; 

    __m128 *r = (__m128 *)comp0;

    __m128 *g = (__m128 *)comp1;

    __m128 *b = (__m128 *)comp2;

    __m128 src[3];

    #define CSC_INVERSE_709_SSE2_LOOP(i)                       \

            src[0] = r[i];                                     \

            src[1] = g[i];                                     \

            src[2] = b[i];                                     \

                                                               \

            r[i] = _mm_add_ps (r[i], _mm_mul_ps (src[2], c0)); \

                                                               \

            g[i]   = _mm_mul_ps (g[i], c2);                    \

            src[2] = _mm_mul_ps (src[2], c3);                  \

            g[i]   = _mm_add_ps (g[i], src[0]);                \

            g[i]   = _mm_add_ps (g[i], src[2]);                \

                                                               \

            b[i] = _mm_mul_ps (c1,   src[1]);                  \

            b[i] = _mm_add_ps (b[i], src[0]);

    CSC_INVERSE_709_SSE2_LOOP (0)

    CSC_INVERSE_709_SSE2_LOOP (1)

    CSC_INVERSE_709_SSE2_LOOP (2)

    CSC_INVERSE_709_SSE2_LOOP (3)

    CSC_INVERSE_709_SSE2_LOOP (4)

    CSC_INVERSE_709_SSE2_LOOP (5)

    CSC_INVERSE_709_SSE2_LOOP (6)

    CSC_INVERSE_709_SSE2_LOOP (7)

    CSC_INVERSE_709_SSE2_LOOP (8)

    CSC_INVERSE_709_SSE2_LOOP (9)

    CSC_INVERSE_709_SSE2_LOOP (10)

    CSC_INVERSE_709_SSE2_LOOP (11)

    CSC_INVERSE_709_SSE2_LOOP (12)

    CSC_INVERSE_709_SSE2_LOOP (13)

    CSC_INVERSE_709_SSE2_LOOP (14)

    CSC_INVERSE_709_SSE2_LOOP (15)

}

#endif /* IMF_HAVE_SSE2 */

//

// Color space conversion, Forward 709 CSC, R'G'B' -> Y'CbCr

//

// Simple FPU color space conversion. Based on the 709

// primary chromaticies, with no scaling or offsets.

//

void

csc709Forward64 (float *comp0, float *comp1, float *comp2)

{

    float src[3];

    for (int i = 0; i<64; ++i)

    {

        src[0] = comp0[i];    

        src[1] = comp1[i]; 

        src[2] = comp2[i];     

        comp0[i] =  0.2126f * src[0] + 0.7152f * src[1] + 0.0722f * src[2];

        comp1[i] = -0.1146f * src[0] - 0.3854f * src[1] + 0.5000f * src[2];

        comp2[i] =  0.5000f * src[0] - 0.4542f * src[1] - 0.0458f * src[2];

    }

}

//

// Byte interleaving of 2 byte arrays:

//    src0 = AAAA 

//    src1 = BBBB

//    dst  = ABABABAB

//

// numBytes is the size of each of the source buffers

//

#ifndef IMF_HAVE_SSE2 

// 

// Scalar default implementation 

//

void

interleaveByte2 (char *dst, char *src0, char *src1, int numBytes)

{

    for (int x = 0; x < numBytes; ++x)

    {

        dst[2 * x]     = src0[x];

        dst[2 * x + 1] = src1[x];

    }

}

#else  /* IMF_HAVE_SSE2 */

// 

// SSE2 byte interleaving

//

void

interleaveByte2 (char *dst, char *src0, char *src1, int numBytes)

{

    int dstAlignment  = (size_t)dst  % 16;

    int src0Alignment = (size_t)src0 % 16;

    int src1Alignment = (size_t)src1 % 16;

    __m128i *dst_epi8  = (__m128i*)dst;

    __m128i *src0_epi8 = (__m128i*)src0;

    __m128i *src1_epi8 = (__m128i*)src1;

    int sseWidth  =  numBytes / 16;

    if ((!dstAlignment) && (!src0Alignment) && (!src1Alignment))

    {

        __m128i tmp0, tmp1;

        //

        // Aligned loads and stores

        //

        for (int x = 0; x < sseWidth; ++x)

        {

            tmp0 = src0_epi8[x];

            tmp1 = src1_epi8[x];

            _mm_stream_si128 (&dst_epi8[2 * x],

                              _mm_unpacklo_epi8 (tmp0, tmp1));

            _mm_stream_si128 (&dst_epi8[2 * x + 1],

                              _mm_unpackhi_epi8 (tmp0, tmp1));

        }

        //

        // Then do run the leftovers one at a time

        //

        for (int x = 16 * sseWidth; x < numBytes; ++x)

        {

            dst[2 * x]     = src0[x];

            dst[2 * x + 1] = src1[x];

        }

    }

    else if ((!dstAlignment) && (src0Alignment == 8) && (src1Alignment == 8))

    {

        //

        // Aligned stores, but catch up a few values so we can 

        // use aligned loads

        //

        for (int x = 0; x < 8; ++x)

        {

            dst[2 * x]     = src0[x];

            dst[2 * x + 1] = src1[x];

        }

        dst_epi8  = (__m128i*)&dst[16];

        src0_epi8 = (__m128i*)&src0[8];

        src1_epi8 = (__m128i*)&src1[8];

        sseWidth  =  (numBytes - 8) / 16;

        for (int x=0; x

        {

            _mm_stream_si128 (&dst_epi8[2 * x],

                              _mm_unpacklo_epi8 (src0_epi8[x], src1_epi8[x]));

            _mm_stream_si128 (&dst_epi8[2 * x + 1],

                              _mm_unpackhi_epi8 (src0_epi8[x], src1_epi8[x]));

        }

        //

        // Then do run the leftovers one at a time

        //

        for (int x = 16 * sseWidth + 8; x < numBytes; ++x)

        {

            dst[2 * x]     = src0[x];

            dst[2 * x + 1] = src1[x];

        }

    }

    else

    {

        //

        // Unaligned everything

        //

        for (int x = 0; x < sseWidth; ++x)

        {

            __m128i tmpSrc0_epi8 = _mm_loadu_si128 (&src0_epi8[x]);

            __m128i tmpSrc1_epi8 = _mm_loadu_si128 (&src1_epi8[x]);

            _mm_storeu_si128 (&dst_epi8[2 * x],

                              _mm_unpacklo_epi8 (tmpSrc0_epi8, tmpSrc1_epi8));

            _mm_storeu_si128 (&dst_epi8[2 * x + 1],

                              _mm_unpackhi_epi8 (tmpSrc0_epi8, tmpSrc1_epi8));

        }

        //

        // Then do run the leftovers one at a time

        //

        for (int x = 16 * sseWidth; x < numBytes; ++x)

        {

            dst[2 * x]     = src0[x];

            dst[2 * x + 1] = src1[x];

        }

    }

}

#endif /* IMF_HAVE_SSE2 */

//

// Float -> half float conversion

//

// To enable F16C based conversion, we can't rely on compile-time

// detection, hence the multiple defined versions. Pick one based

// on runtime cpuid detection.

//

//

// Default boring conversion

//

void 

convertFloatToHalf64_scalar (unsigned short *dst, float *src)

{

    for (int i=0; i<64; ++i)

        dst[i] = ((half)src[i]).bits();

}

//

// F16C conversion - Assumes aligned src and dst

//

void

convertFloatToHalf64_f16c (unsigned short *dst, float *src)

{

    //

    // Ordinarly, I'd avoid using inline asm and prefer intrinsics. 

    // However, in order to get the intrinsics, we need to tell 

    // the compiler to generate VEX instructions.

    //

    // (On the GCC side, -mf16c goes ahead and activates -mavc,

    //  resulting in VEX code. Without -mf16c, no intrinsics..)

    //

    // Now, it's quite likely that we'll find ourselves in situations

    // where we want to build *without* VEX, in order to maintain

    // maximum compatability. But to get there with intrinsics,

    // we'd need to break out code into a separate file. Bleh.

    // I'll take the asm.

    //

    #if defined IMF_HAVE_GCC_INLINEASM

        __asm__

           ("vmovaps       (%0),     %%ymm0         \n"

            "vmovaps   0x20(%0),     %%ymm1         \n"

            "vmovaps   0x40(%0),     %%ymm2         \n"

            "vmovaps   0x60(%0),     %%ymm3         \n"

            "vcvtps2ph $0,           %%ymm0, %%xmm0 \n"

            "vcvtps2ph $0,           %%ymm1, %%xmm1 \n"

            "vcvtps2ph $0,           %%ymm2, %%xmm2 \n"

            "vcvtps2ph $0,           %%ymm3, %%xmm3 \n"

            "vmovdqa   %%xmm0,       0x00(%1)       \n"

            "vmovdqa   %%xmm1,       0x10(%1)       \n"

            "vmovdqa   %%xmm2,       0x20(%1)       \n"

            "vmovdqa   %%xmm3,       0x30(%1)       \n"

            "vmovaps   0x80(%0),     %%ymm0         \n"

            "vmovaps   0xa0(%0),     %%ymm1         \n"

            "vmovaps   0xc0(%0),     %%ymm2         \n"

            "vmovaps   0xe0(%0),     %%ymm3         \n"

            "vcvtps2ph $0,           %%ymm0, %%xmm0 \n"

            "vcvtps2ph $0,           %%ymm1, %%xmm1 \n"

            "vcvtps2ph $0,           %%ymm2, %%xmm2 \n"

            "vcvtps2ph $0,           %%ymm3, %%xmm3 \n"

            "vmovdqa   %%xmm0,       0x40(%1)       \n"

            "vmovdqa   %%xmm1,       0x50(%1)       \n"

            "vmovdqa   %%xmm2,       0x60(%1)       \n"

            "vmovdqa   %%xmm3,       0x70(%1)       \n"

        #ifndef __AVX__

            "vzeroupper                             \n"

        #endif /* __AVX__ */

            : /* Output  */                

            : /* Input   */ "r"(src), "r"(dst)

        #ifndef __AVX__

            : /* Clobber */ "%xmm0", "%xmm1", "%xmm2", "%xmm3", "memory"

        #else

            : /* Clobber */ "%ymm0", "%ymm1", "%ymm2", "%ymm3", "memory"

        #endif /* __AVX__ */

           );

    #else

        convertFloatToHalf64_scalar (dst, src);

    #endif /* IMF_HAVE_GCC_INLINEASM */

}

//

// Convert an 8x8 block of HALF from zig-zag order to

// FLOAT in normal order. The order we want is:

//

//          src                           dst 

//  0  1  2  3  4  5  6  7       0  1  5  6 14 15 27 28

//  8  9 10 11 12 13 14 15       2  4  7 13 16 26 29 42

// 16 17 18 19 20 21 22 23       3  8 12 17 25 30 41 43

// 24 25 26 27 28 29 30 31       9 11 18 24 31 40 44 53  

// 32 33 34 35 36 37 38 39      10 19 23 32 39 45 52 54

// 40 41 42 43 44 45 46 47      20 22 33 38 46 51 55 60

// 48 49 50 51 52 53 54 55      21 34 37 47 50 56 59 61

// 56 57 58 59 60 61 62 63      35 36 48 49 57 58 62 63

//

void

fromHalfZigZag_scalar (unsigned short *src, float *dst)

{

    half *srcHalf = (half *)src;

    dst[0] = (float)srcHalf[0];

    dst[1] = (float)srcHalf[1];

    dst[2] = (float)srcHalf[5];

    dst[3] = (float)srcHalf[6];

    dst[4] = (float)srcHalf[14];

    dst[5] = (float)srcHalf[15];

    dst[6] = (float)srcHalf[27];

    dst[7] = (float)srcHalf[28];

    dst[8] = (float)srcHalf[2];

    dst[9] = (float)srcHalf[4];

    dst[10] = (float)srcHalf[7];

    dst[11] = (float)srcHalf[13];

    dst[12] = (float)srcHalf[16];

    dst[13] = (float)srcHalf[26];

    dst[14] = (float)srcHalf[29];

    dst[15] = (float)srcHalf[42];

    dst[16] = (float)srcHalf[3];

    dst[17] = (float)srcHalf[8];

    dst[18] = (float)srcHalf[12];

    dst[19] = (float)srcHalf[17];

    dst[20] = (float)srcHalf[25];

    dst[21] = (float)srcHalf[30];

    dst[22] = (float)srcHalf[41];

    dst[23] = (float)srcHalf[43];

    dst[24] = (float)srcHalf[9];

    dst[25] = (float)srcHalf[11];

    dst[26] = (float)srcHalf[18];

    dst[27] = (float)srcHalf[24];

    dst[28] = (float)srcHalf[31];

    dst[29] = (float)srcHalf[40];

    dst[30] = (float)srcHalf[44];

    dst[31] = (float)srcHalf[53];

    dst[32] = (float)srcHalf[10];

    dst[33] = (float)srcHalf[19];

    dst[34] = (float)srcHalf[23];

    dst[35] = (float)srcHalf[32];

    dst[36] = (float)srcHalf[39];

    dst[37] = (float)srcHalf[45];

    dst[38] = (float)srcHalf[52];

    dst[39] = (float)srcHalf[54];

    dst[40] = (float)srcHalf[20];

    dst[41] = (float)srcHalf[22];

    dst[42] = (float)srcHalf[33];

    dst[43] = (float)srcHalf[38];

    dst[44] = (float)srcHalf[46];

    dst[45] = (float)srcHalf[51];

    dst[46] = (float)srcHalf[55];

    dst[47] = (float)srcHalf[60];

    dst[48] = (float)srcHalf[21];

    dst[49] = (float)srcHalf[34];

    dst[50] = (float)srcHalf[37];

    dst[51] = (float)srcHalf[47];

    dst[52] = (float)srcHalf[50];

    dst[53] = (float)srcHalf[56];

    dst[54] = (float)srcHalf[59];

    dst[55] = (float)srcHalf[61];

    dst[56] = (float)srcHalf[35];

    dst[57] = (float)srcHalf[36];

    dst[58] = (float)srcHalf[48];

    dst[59] = (float)srcHalf[49];

    dst[60] = (float)srcHalf[57];

    dst[61] = (float)srcHalf[58];

    dst[62] = (float)srcHalf[62];

    dst[63] = (float)srcHalf[63];

}

//

// If we can form the correct ordering in xmm registers,

// we can use F16C to convert from HALF -> FLOAT. However,

// making the correct order isn't trivial. 

// 

// We want to re-order a source 8x8 matrix from:

//

//  0  1  2  3  4  5  6  7       0  1  5  6 14 15 27 28

//  8  9 10 11 12 13 14 15       2  4  7 13 16 26 29 42

// 16 17 18 19 20 21 22 23       3  8 12 17 25 30 41 43

// 24 25 26 27 28 29 30 31       9 11 18 24 31 40 44 53   (A)

// 32 33 34 35 36 37 38 39  --> 10 19 23 32 39 45 52 54

// 40 41 42 43 44 45 46 47      20 22 33 38 46 51 55 60

// 48 49 50 51 52 53 54 55      21 34 37 47 50 56 59 61

// 56 57 58 59 60 61 62 63      35 36 48 49 57 58 62 63

//

// Which looks like a mess, right? 

//

// Now, check out the NE/SW diagonals of (A). Along those lines, 

// we have runs of contiguous values! If we rewrite (A) a bit, we get:

//

//  0

//  1  2

//  5  4  3

//  6  7  8  9

// 14 13 12 11 10

// 15 16 17 18 19 20

// 27 26 25 24 23 22 21            (B)

// 28 29 30 31 32 33 34 35

//    42 41 40 39 38 37 36

//       43 44 45 46 47 48

//          53 52 51 50 49

//             54 55 56 57

//                60 59 58

//                   61 62

//                      63

//

// In this ordering, the columns are the rows (A). If we can 'transpose' 

// (B), we'll achieve our goal. But we want this to fit nicely into 

// xmm registers and still be able to load large runs efficiently.  

// Also, notice that the odd rows are in ascending order, while 

// the even rows are in descending order. 

//

// If we 'fold' the bottom half up into the top, we can preserve ordered

// runs accross rows, and still keep all the correct values in columns. 

// After transposing, we'll need to rotate things back into place. 

// This gives us:

//

//  0 | 42   41   40   39   38   37   36

//  1    2 | 43   44   45   46   47   48

//  5    4    3 | 53   52   51   50   49

//  6    7    8    9 | 54   55   56   57      (C)

// 14   13   12   11   10 | 60   59   58

// 15   16   17   18   19   20 | 61   62

// 27   26   25   24   23   22   21 | 61

// 28   29   30   31   32   33   34   35

//

// But hang on. We still have the backwards descending rows to deal with.

// Lets reverse the even rows so that all values are in ascending order

//

//  36   37  38   39   40   41   42 | 0

//  1    2 | 43   44   45   46   47   48

//  49   50  51   52   53 |  3    4    5  

//  6    7    8    9 | 54   55   56   57      (D)

// 58   59   60 | 10   11   12   13   14  

// 15   16   17   18   19   20 | 61   62

// 61 | 21   22   23   24   25   26   27 

// 28   29   30   31   32   33   34   35

//

// If we can form (D),  we will then:

//   1) Reverse the even rows

//   2) Transpose

//   3) Rotate the rows 

//

// and we'll have (A).

//

void 

fromHalfZigZag_f16c (unsigned short *src, float *dst)

{

    #if defined IMF_HAVE_GCC_INLINEASM_64

        __asm__

           /* x3 <- 0                    

            * x8 <- [ 0- 7]              

            * x6 <- [56-63]              

            * x9 <- [21-28]              

            * x7 <- [28-35]              

            * x3 <- [ 6- 9] (lower half) */

          ("vpxor   %%xmm3,  %%xmm3, %%xmm3   \n"

           "vmovdqa    (%0), %%xmm8           \n"

           "vmovdqa 112(%0), %%xmm6           \n"

           "vmovdqu  42(%0), %%xmm9           \n"

           "vmovdqu  56(%0), %%xmm7           \n"

           "vmovq    12(%0), %%xmm3           \n"

           /* Setup rows 0-2 of A in xmm0-xmm2 

            * x1 <- x8 >> 16 (1 value)     

            * x2 <- x8 << 32 (2 values)    

            * x0 <- alignr([35-42], x8, 2) 

            * x1 <- blend(x1, [41-48])     

            * x2 <- blend(x2, [49-56])     */

           "vpsrldq      $2, %%xmm8, %%xmm1   \n"      

           "vpslldq      $4, %%xmm8, %%xmm2   \n"      

           "vpalignr     $2, 70(%0), %%xmm8, %%xmm0 \n"

           "vpblendw  $0xfc, 82(%0), %%xmm1, %%xmm1 \n"

           "vpblendw  $0x1f, 98(%0), %%xmm2, %%xmm2 \n"

           /* Setup rows 4-6 of A in xmm4-xmm6 

            * x4 <- x6 >> 32 (2 values)   

            * x5 <- x6 << 16 (1 value)    

            * x6 <- alignr(x6,x9,14)      

            * x4 <- blend(x4, [ 7-14])    

            * x5 <- blend(x5, [15-22])    */

           "vpsrldq      $4, %%xmm6, %%xmm4         \n"

           "vpslldq      $2, %%xmm6, %%xmm5         \n"

           "vpalignr    $14, %%xmm6, %%xmm9, %%xmm6 \n"

           "vpblendw  $0xf8, 14(%0), %%xmm4, %%xmm4 \n"

           "vpblendw  $0x3f, 30(%0), %%xmm5, %%xmm5 \n"

           /* Load the upper half of row 3 into xmm3 

            * x3 <- [54-57] (upper half) */

           "vpinsrq      $1, 108(%0), %%xmm3, %%xmm3\n"

           /* Reverse the even rows. We're not using PSHUFB as

            * that requires loading an extra constant all the time,

            * and we're alreadly pretty memory bound.

            */

           "vpshuflw $0x1b, %%xmm0, %%xmm0          \n" 

           "vpshuflw $0x1b, %%xmm2, %%xmm2          \n" 

           "vpshuflw $0x1b, %%xmm4, %%xmm4          \n" 

           "vpshuflw $0x1b, %%xmm6, %%xmm6          \n" 

           "vpshufhw $0x1b, %%xmm0, %%xmm0          \n" 

           "vpshufhw $0x1b, %%xmm2, %%xmm2          \n" 

           "vpshufhw $0x1b, %%xmm4, %%xmm4          \n" 

           "vpshufhw $0x1b, %%xmm6, %%xmm6          \n" 

           "vpshufd $0x4e, %%xmm0, %%xmm0          \n" 

           "vpshufd $0x4e, %%xmm2, %%xmm2          \n" 

           "vpshufd $0x4e, %%xmm4, %%xmm4          \n" 

           "vpshufd $0x4e, %%xmm6, %%xmm6          \n" 

           /* Transpose xmm0-xmm7 into xmm8-xmm15 */

           "vpunpcklwd %%xmm1, %%xmm0, %%xmm8       \n"

           "vpunpcklwd %%xmm3, %%xmm2, %%xmm9       \n"

           "vpunpcklwd %%xmm5, %%xmm4, %%xmm10      \n"

           "vpunpcklwd %%xmm7, %%xmm6, %%xmm11      \n"

           "vpunpckhwd %%xmm1, %%xmm0, %%xmm12      \n"

           "vpunpckhwd %%xmm3, %%xmm2, %%xmm13      \n"

           "vpunpckhwd %%xmm5, %%xmm4, %%xmm14      \n"

           "vpunpckhwd %%xmm7, %%xmm6, %%xmm15      \n"

           "vpunpckldq  %%xmm9,  %%xmm8, %%xmm0     \n"

           "vpunpckldq %%xmm11, %%xmm10, %%xmm1     \n"

           "vpunpckhdq  %%xmm9,  %%xmm8, %%xmm2     \n"

           "vpunpckhdq %%xmm11, %%xmm10, %%xmm3     \n"

           "vpunpckldq %%xmm13, %%xmm12, %%xmm4     \n"

           "vpunpckldq %%xmm15, %%xmm14, %%xmm5     \n"

           "vpunpckhdq %%xmm13, %%xmm12, %%xmm6     \n"

           "vpunpckhdq %%xmm15, %%xmm14, %%xmm7     \n"

           "vpunpcklqdq %%xmm1,  %%xmm0, %%xmm8     \n"

           "vpunpckhqdq %%xmm1,  %%xmm0, %%xmm9     \n"

           "vpunpcklqdq %%xmm3,  %%xmm2, %%xmm10    \n"

           "vpunpckhqdq %%xmm3,  %%xmm2, %%xmm11    \n"

           "vpunpcklqdq %%xmm4,  %%xmm5, %%xmm12    \n"

           "vpunpckhqdq %%xmm5,  %%xmm4, %%xmm13    \n"

           "vpunpcklqdq %%xmm7,  %%xmm6, %%xmm14    \n"

           "vpunpckhqdq %%xmm7,  %%xmm6, %%xmm15    \n"

           /* Rotate the rows to get the correct final order. 

            * Rotating xmm12 isn't needed, as we can handle

            * the rotation in the PUNPCKLQDQ above. Rotating

            * xmm8 isn't needed as it's already in the right order           

            */

           "vpalignr  $2,  %%xmm9,  %%xmm9,  %%xmm9 \n"

           "vpalignr  $4, %%xmm10, %%xmm10, %%xmm10 \n"

           "vpalignr  $6, %%xmm11, %%xmm11, %%xmm11 \n"

           "vpalignr $10, %%xmm13, %%xmm13, %%xmm13 \n"

           "vpalignr $12, %%xmm14, %%xmm14, %%xmm14 \n"

           "vpalignr $14, %%xmm15, %%xmm15, %%xmm15 \n"

            /* Convert from half -> float */

           "vcvtph2ps  %%xmm8, %%ymm8            \n"  

           "vcvtph2ps  %%xmm9, %%ymm9            \n"

           "vcvtph2ps %%xmm10, %%ymm10           \n"

           "vcvtph2ps %%xmm11, %%ymm11           \n"

           "vcvtph2ps %%xmm12, %%ymm12           \n"

           "vcvtph2ps %%xmm13, %%ymm13           \n"

           "vcvtph2ps %%xmm14, %%ymm14           \n"

           "vcvtph2ps %%xmm15, %%ymm15           \n"

           /* Move float values to dst */

           "vmovaps    %%ymm8,    (%1)           \n"

           "vmovaps    %%ymm9,  32(%1)           \n"

           "vmovaps   %%ymm10,  64(%1)           \n" 

           "vmovaps   %%ymm11,  96(%1)           \n" 

           "vmovaps   %%ymm12, 128(%1)           \n" 

           "vmovaps   %%ymm13, 160(%1)           \n" 

           "vmovaps   %%ymm14, 192(%1)           \n" 

           "vmovaps   %%ymm15, 224(%1)           \n"

        #ifndef __AVX__

            "vzeroupper                          \n"

        #endif /* __AVX__ */

            : /* Output  */                

            : /* Input   */ "r"(src), "r"(dst)

            : /* Clobber */ "memory",

        #ifndef __AVX__

                            "%xmm0",  "%xmm1",  "%xmm2",  "%xmm3", 

                            "%xmm4",  "%xmm5",  "%xmm6",  "%xmm7",

                            "%xmm8",  "%xmm9",  "%xmm10", "%xmm11",

                            "%xmm12", "%xmm13", "%xmm14", "%xmm15"

        #else

                            "%ymm0",  "%ymm1",  "%ymm2",  "%ymm3", 

                            "%ymm4",  "%ymm5",  "%ymm6",  "%ymm7",

                            "%ymm8",  "%ymm9",  "%ymm10", "%ymm11",

                            "%ymm12", "%ymm13", "%ymm14", "%ymm15"

        #endif /* __AVX__ */

        );

    #else

        fromHalfZigZag_scalar(src, dst);

    #endif /* defined IMF_HAVE_GCC_INLINEASM_64 */

}

//

// Inverse 8x8 DCT, only inverting the DC. This assumes that

// all AC frequencies are 0.

//

#ifndef IMF_HAVE_SSE2

void 

dctInverse8x8DcOnly (float *data)

{

    float val = data[0] * 3.535536e-01f * 3.535536e-01f;

    for (int i = 0; i < 64; ++i)

        data[i] = val;

}

#else  /* IMF_HAVE_SSE2 */

void

dctInverse8x8DcOnly (float *data)

{

    __m128 src = _mm_set1_ps (data[0] * 3.535536e-01f * 3.535536e-01f);

    __m128 *dst = (__m128 *)data;

    for (int i = 0; i < 16; ++i)

        dst[i] = src;

}

#endif /* IMF_HAVE_SSE2 */

//

// Full 8x8 Inverse DCT:

//

// Simple inverse DCT on an 8x8 block, with scalar ops only.

//  Operates on data in-place.

//

// This is based on the iDCT formuation (y = frequency domain,

//                                       x = spatial domain)

//

//    [x0]    [        ][y0]    [        ][y1] 

//    [x1] =  [  M1    ][y2]  + [  M2    ][y3] 

//    [x2]    [        ][y4]    [        ][y5] 

//    [x3]    [        ][y6]    [        ][y7]

//

//    [x7]    [        ][y0]    [        ][y1] 

//    [x6] =  [  M1    ][y2]  - [  M2    ][y3] 

//    [x5]    [        ][y4]    [        ][y5] 

//    [x4]    [        ][y6]    [        ][y7]

//

// where M1:             M2:

//

//   [a  c  a   f]     [b  d  e  g]

//   [a  f -a  -c]     [d -g -b -e]

//   [a -f -a   c]     [e -b  g  d]

//   [a -c  a  -f]     [g -e  d -b]

//

// and the constants are as defined below..

//

// If you know how many of the lower rows are zero, that can

// be passed in to help speed things up. If you don't know, 

// just set zeroedRows=0.

//

//

// Default implementation

//