#define GGML_COMMON_IMPL_C

#include "ggml-common.h"

#include "ggml-quants.h"

#include "ggml-impl.h"

#include "ggml-cpu.h"

#include "simd-mappings.h"

#include "../../quants.h"

#include "../../ggml-cpu-impl.h"

#include <math.h>

#include <string.h>

#include <assert.h>

#include <stdlib.h> // for qsort

#include <stdio.h>  // for GGML_ASSERT

#define GROUP_MAX_EPS 1e-15f

#define GROUP_MAX_EPS_IQ3_XXS 1e-8f

#define GROUP_MAX_EPS_IQ2_S 1e-8f

#define GROUP_MAX_EPS_IQ1_M 1e-7f

#define GROUP_MAX_EPS_IQ1_S 1e-12f

#define UNUSED GGML_UNUSED

// some compilers don't provide _mm256_set_m128i, e.g. gcc 7

#define MM256_SET_M128I(a, b) _mm256_insertf128_si256(_mm256_castsi128_si256(b), (a), 1)

#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__) || defined(__SSSE3__)

// multiply int8_t, add results pairwise twice

static inline __m128i mul_sum_i8_pairs(const __m128i x, const __m128i y) {

    // Get absolute values of x vectors

    const __m128i ax = _mm_sign_epi8(x, x);

    // Sign the values of the y vectors

    const __m128i sy = _mm_sign_epi8(y, x);

    // Perform multiplication and create 16-bit values

    const __m128i dot = _mm_maddubs_epi16(ax, sy);

    const __m128i ones = _mm_set1_epi16(1);

    return _mm_madd_epi16(ones, dot);

}

#if __AVX__ || __AVX2__ || __AVX512F__

// horizontally add 8 floats

static inline float hsum_float_8(const __m256 x) {

    __m128 res = _mm256_extractf128_ps(x, 1);

    res = _mm_add_ps(res, _mm256_castps256_ps128(x));

    res = _mm_add_ps(res, _mm_movehl_ps(res, res));

    res = _mm_add_ss(res, _mm_movehdup_ps(res));

    return _mm_cvtss_f32(res);

}

// horizontally add 8 int32_t

static inline int hsum_i32_8(const __m256i a) {

    const __m128i sum128 = _mm_add_epi32(_mm256_castsi256_si128(a), _mm256_extractf128_si256(a, 1));

    const __m128i hi64 = _mm_unpackhi_epi64(sum128, sum128);

    const __m128i sum64 = _mm_add_epi32(hi64, sum128);

    const __m128i hi32  = _mm_shuffle_epi32(sum64, _MM_SHUFFLE(2, 3, 0, 1));

    return _mm_cvtsi128_si32(_mm_add_epi32(sum64, hi32));

}

// horizontally add 4 int32_t

static inline int hsum_i32_4(const __m128i a) {

    const __m128i hi64 = _mm_unpackhi_epi64(a, a);

    const __m128i sum64 = _mm_add_epi32(hi64, a);

    const __m128i hi32  = _mm_shuffle_epi32(sum64, _MM_SHUFFLE(2, 3, 0, 1));

    return _mm_cvtsi128_si32(_mm_add_epi32(sum64, hi32));

}

#if defined(__AVX2__) || defined(__AVX512F__)

static inline __m256i mul_add_epi8(const __m256i x, const __m256i y) {

    const __m256i ax = _mm256_sign_epi8(x, x);

    const __m256i sy = _mm256_sign_epi8(y, x);

    return _mm256_maddubs_epi16(ax, sy);

}

// spread 32 bits to 32 bytes { 0x00, 0xFF }

static inline __m256i bytes_from_bits_32(const uint8_t * x) {

    uint32_t x32;

    memcpy(&x32, x, sizeof(uint32_t));

    const __m256i shuf_mask = _mm256_set_epi64x(

            0x0303030303030303, 0x0202020202020202,

            0x0101010101010101, 0x0000000000000000);

    __m256i bytes = _mm256_shuffle_epi8(_mm256_set1_epi32(x32), shuf_mask);

    const __m256i bit_mask = _mm256_set1_epi64x(0x7fbfdfeff7fbfdfe);

    bytes = _mm256_or_si256(bytes, bit_mask);

    return _mm256_cmpeq_epi8(bytes, _mm256_set1_epi64x(-1));

}

// Unpack 32 4-bit fields into 32 bytes

// The output vector contains 32 bytes, each one in [ 0 .. 15 ] interval

static inline __m256i bytes_from_nibbles_32(const uint8_t * rsi)

{

    const __m128i tmp = _mm_loadu_si128((const __m128i *)rsi);

    const __m256i bytes = MM256_SET_M128I(_mm_srli_epi16(tmp, 4), tmp);

    const __m256i lowMask = _mm256_set1_epi8( 0xF );

    return _mm256_and_si256(lowMask, bytes);

}

// add int16_t pairwise and return as float vector

static inline __m256 sum_i16_pairs_float(const __m256i x) {

    const __m256i ones = _mm256_set1_epi16(1);

    const __m256i summed_pairs = _mm256_madd_epi16(ones, x);

    return _mm256_cvtepi32_ps(summed_pairs);

}

static inline __m256 mul_sum_us8_pairs_float(const __m256i ax, const __m256i sy) {

#if defined(__AVX512VNNI__) && defined(__AVX512VL__)

    const __m256i zero = _mm256_setzero_si256();

    const __m256i summed_pairs = _mm256_dpbusd_epi32(zero, ax, sy);

    return _mm256_cvtepi32_ps(summed_pairs);

#elif defined(__AVXVNNI__)

    const __m256i zero = _mm256_setzero_si256();

    const __m256i summed_pairs = _mm256_dpbusd_avx_epi32(zero, ax, sy);

    return _mm256_cvtepi32_ps(summed_pairs);

#else

    // Perform multiplication and create 16-bit values

    const __m256i dot = _mm256_maddubs_epi16(ax, sy);

    return sum_i16_pairs_float(dot);

#endif

}

// multiply int8_t, add results pairwise twice and return as float vector

static inline __m256 mul_sum_i8_pairs_float(const __m256i x, const __m256i y) {

#if __AVXVNNIINT8__

    const __m256i zero = _mm256_setzero_si256();

    const __m256i summed_pairs = _mm256_dpbssd_epi32(zero, x, y);

    return _mm256_cvtepi32_ps(summed_pairs);

#else

    // Get absolute values of x vectors

    const __m256i ax = _mm256_sign_epi8(x, x);

    // Sign the values of the y vectors

    const __m256i sy = _mm256_sign_epi8(y, x);

    return mul_sum_us8_pairs_float(ax, sy);

#endif

}

static inline __m128i packNibbles( __m256i bytes )

{

    // Move bits within 16-bit lanes from 0000_abcd_0000_efgh into 0000_0000_abcd_efgh

#if __AVX512F__

    const __m256i bytes_srli_4 = _mm256_srli_epi16(bytes, 4);   // 0000_0000_abcd_0000

    bytes = _mm256_or_si256(bytes, bytes_srli_4);               // 0000_abcd_abcd_efgh

    return _mm256_cvtepi16_epi8(bytes);                         // abcd_efgh

#else

    const __m256i lowByte = _mm256_set1_epi16( 0xFF );

    __m256i high = _mm256_andnot_si256( lowByte, bytes );

    __m256i low = _mm256_and_si256( lowByte, bytes );

    high = _mm256_srli_epi16( high, 4 );

    bytes = _mm256_or_si256( low, high );

    // Compress uint16_t lanes into bytes

    __m128i r0 = _mm256_castsi256_si128( bytes );

    __m128i r1 = _mm256_extracti128_si256( bytes, 1 );

    return _mm_packus_epi16( r0, r1 );

#endif

}

#elif defined(__AVX__)

static inline __m128i packNibbles( __m128i bytes1, __m128i bytes2 )

{

    // Move bits within 16-bit lanes from 0000_abcd_0000_efgh into 0000_0000_abcd_efgh

    const __m128i lowByte = _mm_set1_epi16( 0xFF );

    __m128i high = _mm_andnot_si128( lowByte, bytes1 );

    __m128i low = _mm_and_si128( lowByte, bytes1 );

    high = _mm_srli_epi16( high, 4 );

    bytes1 = _mm_or_si128( low, high );

    high = _mm_andnot_si128( lowByte, bytes2 );

    low = _mm_and_si128( lowByte, bytes2 );

    high = _mm_srli_epi16( high, 4 );

    bytes2 = _mm_or_si128( low, high );

    return _mm_packus_epi16( bytes1, bytes2);

}

static inline __m128i mul_add_epi8_sse(const __m128i x, const __m128i y) {

    const __m128i ax = _mm_sign_epi8(x, x);

    const __m128i sy = _mm_sign_epi8(y, x);

    return _mm_maddubs_epi16(ax, sy);

}

// spread 32 bits to 32 bytes { 0x00, 0xFF }

static inline __m256i bytes_from_bits_32(const uint8_t * x) {

    uint32_t x32;

    memcpy(&x32, x, sizeof(uint32_t));

    const __m128i shuf_maskl = _mm_set_epi64x(0x0101010101010101, 0x0000000000000000);

    const __m128i shuf_maskh = _mm_set_epi64x(0x0303030303030303, 0x0202020202020202);

    __m128i bytesl = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskl);

    __m128i bytesh = _mm_shuffle_epi8(_mm_set1_epi32(x32), shuf_maskh);

    const __m128i bit_mask = _mm_set1_epi64x(0x7fbfdfeff7fbfdfe);

    bytesl = _mm_or_si128(bytesl, bit_mask);

    bytesh = _mm_or_si128(bytesh, bit_mask);

    bytesl = _mm_cmpeq_epi8(bytesl, _mm_set1_epi64x(-1));

    bytesh = _mm_cmpeq_epi8(bytesh, _mm_set1_epi64x(-1));

    return MM256_SET_M128I(bytesh, bytesl);

}

// Unpack 32 4-bit fields into 32 bytes

// The output vector contains 32 bytes, each one in [ 0 .. 15 ] interval

static inline __m256i bytes_from_nibbles_32(const uint8_t * rsi)

{

    // Load 16 bytes from memory

    __m128i tmpl = _mm_loadu_si128((const __m128i *)rsi);

    __m128i tmph = _mm_srli_epi16(tmpl, 4);

    const __m128i lowMask = _mm_set1_epi8(0xF);

    tmpl = _mm_and_si128(lowMask, tmpl);

    tmph = _mm_and_si128(lowMask, tmph);

    return MM256_SET_M128I(tmph, tmpl);

}

// add int16_t pairwise and return as float vector

static inline __m256 sum_i16_pairs_float(const __m128i xh, const __m128i xl) {

    const __m128i ones = _mm_set1_epi16(1);

    const __m128i summed_pairsl = _mm_madd_epi16(ones, xl);

    const __m128i summed_pairsh = _mm_madd_epi16(ones, xh);

    const __m256i summed_pairs = MM256_SET_M128I(summed_pairsh, summed_pairsl);

    return _mm256_cvtepi32_ps(summed_pairs);

}

static inline __m256 mul_sum_us8_pairs_float(const __m256i ax, const __m256i sy) {

    const __m128i axl = _mm256_castsi256_si128(ax);

    const __m128i axh = _mm256_extractf128_si256(ax, 1);

    const __m128i syl = _mm256_castsi256_si128(sy);

    const __m128i syh = _mm256_extractf128_si256(sy, 1);

    // Perform multiplication and create 16-bit values

    const __m128i dotl = _mm_maddubs_epi16(axl, syl);

    const __m128i doth = _mm_maddubs_epi16(axh, syh);

    return sum_i16_pairs_float(doth, dotl);

}

// multiply int8_t, add results pairwise twice and return as float vector

static inline __m256 mul_sum_i8_pairs_float(const __m256i x, const __m256i y) {

    const __m128i xl = _mm256_castsi256_si128(x);

    const __m128i xh = _mm256_extractf128_si256(x, 1);

    const __m128i yl = _mm256_castsi256_si128(y);

    const __m128i yh = _mm256_extractf128_si256(y, 1);

    // Get absolute values of x vectors

    const __m128i axl = _mm_sign_epi8(xl, xl);

    const __m128i axh = _mm_sign_epi8(xh, xh);

    // Sign the values of the y vectors

    const __m128i syl = _mm_sign_epi8(yl, xl);

    const __m128i syh = _mm_sign_epi8(yh, xh);

    // Perform multiplication and create 16-bit values

    const __m128i dotl = _mm_maddubs_epi16(axl, syl);

    const __m128i doth = _mm_maddubs_epi16(axh, syh);

    return sum_i16_pairs_float(doth, dotl);

}

// larger version of mul_sum_i8_pairs_float where x and y are each represented by four 128-bit vectors

static inline __m256 mul_sum_i8_quad_float(const __m128i x_1_0, const __m128i x_1_1, const __m128i x_2_0, const __m128i x_2_1,

                                           const __m128i y_1_0, const __m128i y_1_1, const __m128i y_2_0, const __m128i y_2_1) {

    const __m128i mone = _mm_set1_epi16(1);

    const __m128i p16_1_0 = mul_add_epi8_sse(x_1_0, y_1_0);

    const __m128i p16_1_1 = mul_add_epi8_sse(x_1_1, y_1_1);

    const __m128i p16_2_0 = mul_add_epi8_sse(x_2_0, y_2_0);

    const __m128i p16_2_1 = mul_add_epi8_sse(x_2_1, y_2_1);

    const __m128i p_1_0 = _mm_madd_epi16(p16_1_0, mone);

    const __m128i p_1_1 = _mm_madd_epi16(p16_1_1, mone);

    const __m128i p_2_0 = _mm_madd_epi16(p16_2_0, mone);

    const __m128i p_2_1 = _mm_madd_epi16(p16_2_1, mone);

    const __m128i p_1 = _mm_add_epi32(p_1_0, p_1_1);

    const __m128i p_2 = _mm_add_epi32(p_2_0, p_2_1);

    return _mm256_cvtepi32_ps(MM256_SET_M128I(p_2, p_1));

}

// quad fp16 delta calculation

static inline __m256 quad_fp16_delta_float(const float x0, const float y0, const float x1, const float y1) {

    // GGML_CPU_FP16_TO_FP32 is faster than Intel F16C

    return _mm256_set_m128(_mm_set1_ps(GGML_CPU_FP16_TO_FP32(x1) * GGML_CPU_FP16_TO_FP32(y1)),

                           _mm_set1_ps(GGML_CPU_FP16_TO_FP32(x0) * GGML_CPU_FP16_TO_FP32(y0)));

}

static inline __m256 quad_mx_delta_float(const int8_t x0, const float y0, const int8_t x1, const float y1) {

    return _mm256_set_m128(_mm_set1_ps(GGML_E8M0_TO_FP32_HALF(x1) * GGML_CPU_FP16_TO_FP32(y1)),

                           _mm_set1_ps(GGML_E8M0_TO_FP32_HALF(x0) * GGML_CPU_FP16_TO_FP32(y0)));

}

#endif

#elif defined(__SSSE3__)

// horizontally add 4x4 floats

static inline float hsum_float_4x4(const __m128 a, const __m128 b, const __m128 c, const __m128 d) {

    __m128 res_0 =_mm_hadd_ps(a, b);

    __m128 res_1 =_mm_hadd_ps(c, d);

    __m128 res =_mm_hadd_ps(res_0, res_1);

    res =_mm_hadd_ps(res, res);

    res =_mm_hadd_ps(res, res);

    return _mm_cvtss_f32(res);

}

#endif // __AVX__ || __AVX2__ || __AVX512F__

#endif // defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__) || defined(__SSSE3__)

void quantize_row_q8_0(const float * GGML_RESTRICT x, void * GGML_RESTRICT vy, int64_t k) {

    assert(QK8_0 == 32);

    assert(k % QK8_0 == 0);

    const int nb = k / QK8_0;

    block_q8_0 * GGML_RESTRICT y = vy;

#if defined(__AVX2__) || defined(__AVX__)

    for (int i = 0; i < nb; i++) {

        // Load elements into 4 AVX vectors

        __m256 v0 = _mm256_loadu_ps( x );

        __m256 v1 = _mm256_loadu_ps( x + 8 );

        __m256 v2 = _mm256_loadu_ps( x + 16 );

        __m256 v3 = _mm256_loadu_ps( x + 24 );

        x += 32;

        // Compute max(abs(e)) for the block

        const __m256 signBit = _mm256_set1_ps( -0.0f );

        __m256 maxAbs = _mm256_andnot_ps( signBit, v0 );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v1 ) );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v2 ) );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v3 ) );

        __m128 max4 = _mm_max_ps( _mm256_extractf128_ps( maxAbs, 1 ), _mm256_castps256_ps128( maxAbs ) );

        max4 = _mm_max_ps( max4, _mm_movehl_ps( max4, max4 ) );

        max4 = _mm_max_ss( max4, _mm_movehdup_ps( max4 ) );

        const float maxScalar = _mm_cvtss_f32( max4 );

        // Quantize these floats

        const float d = maxScalar / 127.f;

        y[i].d = GGML_CPU_FP32_TO_FP16(d);

        const float id = ( maxScalar != 0.0f ) ? 127.f / maxScalar : 0.0f;

        const __m256 mul = _mm256_set1_ps( id );

        // Apply the multiplier

        v0 = _mm256_mul_ps( v0, mul );

        v1 = _mm256_mul_ps( v1, mul );

        v2 = _mm256_mul_ps( v2, mul );

        v3 = _mm256_mul_ps( v3, mul );

        // Round to nearest integer

        v0 = _mm256_round_ps( v0, _MM_ROUND_NEAREST );

        v1 = _mm256_round_ps( v1, _MM_ROUND_NEAREST );

        v2 = _mm256_round_ps( v2, _MM_ROUND_NEAREST );

        v3 = _mm256_round_ps( v3, _MM_ROUND_NEAREST );

        // Convert floats to integers

        __m256i i0 = _mm256_cvtps_epi32( v0 );

        __m256i i1 = _mm256_cvtps_epi32( v1 );

        __m256i i2 = _mm256_cvtps_epi32( v2 );

        __m256i i3 = _mm256_cvtps_epi32( v3 );

#if defined(__AVX2__)

        // Convert int32 to int16

        i0 = _mm256_packs_epi32( i0, i1 );  // 0, 1, 2, 3,  8, 9, 10, 11,  4, 5, 6, 7, 12, 13, 14, 15

        i2 = _mm256_packs_epi32( i2, i3 );  // 16, 17, 18, 19,  24, 25, 26, 27,  20, 21, 22, 23, 28, 29, 30, 31

                                            // Convert int16 to int8

        i0 = _mm256_packs_epi16( i0, i2 );  // 0, 1, 2, 3,  8, 9, 10, 11,  16, 17, 18, 19,  24, 25, 26, 27,  4, 5, 6, 7, 12, 13, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31

        // We got our precious signed bytes, but the order is now wrong

        // These AVX2 pack instructions process 16-byte pieces independently

        // The following instruction is fixing the order

        const __m256i perm = _mm256_setr_epi32( 0, 4, 1, 5, 2, 6, 3, 7 );

        i0 = _mm256_permutevar8x32_epi32( i0, perm );

        _mm256_storeu_si256((__m256i *)y[i].qs, i0);

#else

        // Since we don't have in AVX some necessary functions,

        // we split the registers in half and call AVX2 analogs from SSE

        __m128i ni0 = _mm256_castsi256_si128( i0 );

        __m128i ni1 = _mm256_extractf128_si256( i0, 1);

        __m128i ni2 = _mm256_castsi256_si128( i1 );

        __m128i ni3 = _mm256_extractf128_si256( i1, 1);

        __m128i ni4 = _mm256_castsi256_si128( i2 );

        __m128i ni5 = _mm256_extractf128_si256( i2, 1);

        __m128i ni6 = _mm256_castsi256_si128( i3 );

        __m128i ni7 = _mm256_extractf128_si256( i3, 1);

        // Convert int32 to int16

        ni0 = _mm_packs_epi32( ni0, ni1 );

        ni2 = _mm_packs_epi32( ni2, ni3 );

        ni4 = _mm_packs_epi32( ni4, ni5 );

        ni6 = _mm_packs_epi32( ni6, ni7 );

        // Convert int16 to int8

        ni0 = _mm_packs_epi16( ni0, ni2 );

        ni4 = _mm_packs_epi16( ni4, ni6 );

        _mm_storeu_si128((__m128i *)(y[i].qs +  0), ni0);

        _mm_storeu_si128((__m128i *)(y[i].qs + 16), ni4);

#endif

    }

#else

    GGML_UNUSED(nb);

    // scalar

    quantize_row_q8_0_ref(x, y, k);

#endif

}

void quantize_row_q8_1(const float * GGML_RESTRICT x, void * GGML_RESTRICT vy, int64_t k) {

    assert(k % QK8_1 == 0);

    const int nb = k / QK8_1;

    block_q8_1 * GGML_RESTRICT y = vy;

#if defined(__AVX2__) || defined(__AVX__)

    for (int i = 0; i < nb; i++) {

        // Load elements into 4 AVX vectors

        __m256 v0 = _mm256_loadu_ps( x );

        __m256 v1 = _mm256_loadu_ps( x + 8 );

        __m256 v2 = _mm256_loadu_ps( x + 16 );

        __m256 v3 = _mm256_loadu_ps( x + 24 );

        x += 32;

        // Compute max(abs(e)) for the block

        const __m256 signBit = _mm256_set1_ps( -0.0f );

        __m256 maxAbs = _mm256_andnot_ps( signBit, v0 );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v1 ) );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v2 ) );

        maxAbs = _mm256_max_ps( maxAbs, _mm256_andnot_ps( signBit, v3 ) );

        __m128 max4 = _mm_max_ps( _mm256_extractf128_ps( maxAbs, 1 ), _mm256_castps256_ps128( maxAbs ) );

        max4 = _mm_max_ps( max4, _mm_movehl_ps( max4, max4 ) );

        max4 = _mm_max_ss( max4, _mm_movehdup_ps( max4 ) );

        const float max_scalar = _mm_cvtss_f32( max4 );

        // Quantize these floats

        const float d = max_scalar / 127.f;

        y[i].d = GGML_CPU_FP32_TO_FP16(d);

        const float id = ( max_scalar != 0.0f ) ? 127.f / max_scalar : 0.0f;

        const __m256 mul = _mm256_set1_ps( id );

        // Apply the multiplier

        v0 = _mm256_mul_ps( v0, mul );

        v1 = _mm256_mul_ps( v1, mul );

        v2 = _mm256_mul_ps( v2, mul );

        v3 = _mm256_mul_ps( v3, mul );

        // Round to nearest integer

        v0 = _mm256_round_ps( v0, _MM_ROUND_NEAREST );

        v1 = _mm256_round_ps( v1, _MM_ROUND_NEAREST );

        v2 = _mm256_round_ps( v2, _MM_ROUND_NEAREST );

        v3 = _mm256_round_ps( v3, _MM_ROUND_NEAREST );

        // Convert floats to integers

        __m256i i0 = _mm256_cvtps_epi32( v0 );

        __m256i i1 = _mm256_cvtps_epi32( v1 );

        __m256i i2 = _mm256_cvtps_epi32( v2 );

        __m256i i3 = _mm256_cvtps_epi32( v3 );

#if defined(__AVX2__)

        // Compute the sum of the quants and set y[i].s

        y[i].s = GGML_CPU_FP32_TO_FP16(d * hsum_i32_8(_mm256_add_epi32(_mm256_add_epi32(i0, i1), _mm256_add_epi32(i2, i3))));

        // Convert int32 to int16

        i0 = _mm256_packs_epi32( i0, i1 );  // 0, 1, 2, 3,  8, 9, 10, 11,  4, 5, 6, 7, 12, 13, 14, 15

        i2 = _mm256_packs_epi32( i2, i3 );  // 16, 17, 18, 19,  24, 25, 26, 27,  20, 21, 22, 23, 28, 29, 30, 31

                                            // Convert int16 to int8

        i0 = _mm256_packs_epi16( i0, i2 );  // 0, 1, 2, 3,  8, 9, 10, 11,  16, 17, 18, 19,  24, 25, 26, 27,  4, 5, 6, 7, 12, 13, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31

        // We got our precious signed bytes, but the order is now wrong

        // These AVX2 pack instructions process 16-byte pieces independently

        // The following instruction is fixing the order

        const __m256i perm = _mm256_setr_epi32( 0, 4, 1, 5, 2, 6, 3, 7 );

        i0 = _mm256_permutevar8x32_epi32( i0, perm );

        _mm256_storeu_si256((__m256i *)y[i].qs, i0);

#else

        // Since we don't have in AVX some necessary functions,

        // we split the registers in half and call AVX2 analogs from SSE

        __m128i ni0 = _mm256_castsi256_si128( i0 );

        __m128i ni1 = _mm256_extractf128_si256( i0, 1);

        __m128i ni2 = _mm256_castsi256_si128( i1 );

        __m128i ni3 = _mm256_extractf128_si256( i1, 1);

        __m128i ni4 = _mm256_castsi256_si128( i2 );

        __m128i ni5 = _mm256_extractf128_si256( i2, 1);

        __m128i ni6 = _mm256_castsi256_si128( i3 );

        __m128i ni7 = _mm256_extractf128_si256( i3, 1);

        // Compute the sum of the quants and set y[i].s

        const __m128i s0 = _mm_add_epi32(_mm_add_epi32(ni0, ni1), _mm_add_epi32(ni2, ni3));

        const __m128i s1 = _mm_add_epi32(_mm_add_epi32(ni4, ni5), _mm_add_epi32(ni6, ni7));

        y[i].s = GGML_CPU_FP32_TO_FP16(d * hsum_i32_4(_mm_add_epi32(s0, s1)));

        // Convert int32 to int16

        ni0 = _mm_packs_epi32( ni0, ni1 );

        ni2 = _mm_packs_epi32( ni2, ni3 );

        ni4 = _mm_packs_epi32( ni4, ni5 );

        ni6 = _mm_packs_epi32( ni6, ni7 );

        // Convert int16 to int8

        ni0 = _mm_packs_epi16( ni0, ni2 );

        ni4 = _mm_packs_epi16( ni4, ni6 );

        _mm_storeu_si128((__m128i *)(y[i].qs +  0), ni0);

        _mm_storeu_si128((__m128i *)(y[i].qs + 16), ni4);

#endif

    }

#else

    GGML_UNUSED(nb);

    // scalar

    quantize_row_q8_1_ref(x, y, k);

#endif

}

// placeholder implementation for Apple targets

void quantize_row_q8_K(const float * GGML_RESTRICT x, void * GGML_RESTRICT y, int64_t k) {

    quantize_row_q8_K_ref(x, y, k);

}

//===================================== Dot products =================================

//

// Helper functions

//

#if __AVX__ || __AVX2__ || __AVX512F__

// shuffles to pick the required scales in dot products

static inline __m256i get_scale_shuffle_q3k(int i) {

    static const uint8_t k_shuffle[128] = {

         0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,     2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3,

         4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5,     6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7,

         8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9,    10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,

        12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,    14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,

    };

    return _mm256_loadu_si256((const __m256i*)k_shuffle + i);

}

static inline __m256i get_scale_shuffle_k4(int i) {

    static const uint8_t k_shuffle[256] = {

         0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,

         2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3,

         4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5,

         6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7,

         8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9,

        10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,10,11,

        12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,12,13,

        14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15

    };

    return _mm256_loadu_si256((const __m256i*)k_shuffle + i);

}

static inline __m128i get_scale_shuffle(int i) {

    static const uint8_t k_shuffle[128] = {

         0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,

         2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,

         4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5,

         6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7,

         8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9,

        10,10,10,10,10,10,10,10, 11,11,11,11,11,11,11,11,

        12,12,12,12,12,12,12,12, 13,13,13,13,13,13,13,13,

        14,14,14,14,14,14,14,14, 15,15,15,15,15,15,15,15

    };

    return _mm_loadu_si128((const __m128i*)k_shuffle + i);

}

#endif

void ggml_vec_dot_q4_0_q8_0(int n, float * GGML_RESTRICT s, size_t bs, const void * GGML_RESTRICT vx, size_t bx, const void * GGML_RESTRICT vy, size_t by, int nrc) {

    const int qk = QK8_0;

    const int nb = n / qk;

    assert(n % qk == 0);

    assert(nrc == 1);

    UNUSED(nrc);

    UNUSED(bx);

    UNUSED(by);

    UNUSED(bs);

    const block_q4_0 * GGML_RESTRICT x = vx;

    const block_q8_0 * GGML_RESTRICT y = vy;

    int ib = 0;

    float sumf = 0;

#if defined(__AVX2__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    // Main loop

    for (; ib < nb; ++ib) {

        /* Compute combined scale for the block */

        const __m256 d = _mm256_set1_ps( GGML_CPU_FP16_TO_FP32(x[ib].d) * GGML_CPU_FP16_TO_FP32(y[ib].d) );

        __m256i qx = bytes_from_nibbles_32(x[ib].qs);

        // Now we have a vector with bytes in [ 0 .. 15 ] interval. Offset them into [ -8 .. +7 ] interval.

        const __m256i off = _mm256_set1_epi8( 8 );

        qx = _mm256_sub_epi8( qx, off );

        __m256i qy = _mm256_loadu_si256((const __m256i *)y[ib].qs);

        const __m256 q = mul_sum_i8_pairs_float(qx, qy);

        /* Multiply q with scale and accumulate */

        acc = _mm256_fmadd_ps( d, q, acc );

    }

    sumf = hsum_float_8(acc);

#elif defined(__AVX__)

    __m256 accum = _mm256_setzero_ps();

    for (; ib + 1 < nb; ib += 2) {

        const __m128i q4bits_1 = _mm_loadu_si128((const __m128i *)x[ib + 0].qs);

        const __m128i q4bits_2 = _mm_loadu_si128((const __m128i *)x[ib + 1].qs);

        const __m128i q8b_1_0 = _mm_loadu_si128((const __m128i *)y[ib + 0].qs);

        const __m128i q8b_1_1 = _mm_loadu_si128((const __m128i *)y[ib + 0].qs + 1);

        const __m128i q8b_2_0 = _mm_loadu_si128((const __m128i *)y[ib + 1].qs);

        const __m128i q8b_2_1 = _mm_loadu_si128((const __m128i *)y[ib + 1].qs + 1);

        const __m128i q4b_1_0 = _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), q4bits_1), _mm_set1_epi8(8));

        const __m128i q4b_1_1 = _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), _mm_srli_epi16(q4bits_1, 4)), _mm_set1_epi8(8));

        const __m128i q4b_2_0 = _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), q4bits_2), _mm_set1_epi8(8));

        const __m128i q4b_2_1 = _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), _mm_srli_epi16(q4bits_2, 4)), _mm_set1_epi8(8));

        const __m128i p16_1_0 = mul_add_epi8_sse(q4b_1_0, q8b_1_0);

        const __m128i p16_1_1 = mul_add_epi8_sse(q4b_1_1, q8b_1_1);

        const __m128i p16_2_0 = mul_add_epi8_sse(q4b_2_0, q8b_2_0);

        const __m128i p16_2_1 = mul_add_epi8_sse(q4b_2_1, q8b_2_1);

        const __m128i p_1 = _mm_add_epi16(p16_1_0, p16_1_1);

        const __m128i p_2 = _mm_add_epi16(p16_2_0, p16_2_1);

        const __m256 p =  sum_i16_pairs_float(p_2, p_1);

        const __m256 deltas = quad_fp16_delta_float(x[ib].d, y[ib].d, x[ib + 1].d, y[ib + 1].d);

        accum = _mm256_add_ps(_mm256_mul_ps(deltas, p), accum);

    }

    sumf = hsum_float_8(accum);

#elif defined(__SSSE3__)

    // set constants

    const __m128i lowMask = _mm_set1_epi8(0xF);

    const __m128i off = _mm_set1_epi8(8);

    // Initialize accumulator with zeros

    __m128 acc_0 = _mm_setzero_ps();

    __m128 acc_1 = _mm_setzero_ps();

    __m128 acc_2 = _mm_setzero_ps();

    __m128 acc_3 = _mm_setzero_ps();

    for (; ib + 1 < nb; ib += 2) {

        _mm_prefetch(&x[ib] + sizeof(block_q4_0), _MM_HINT_T0);

        _mm_prefetch(&y[ib] + sizeof(block_q8_0), _MM_HINT_T0);

        // Compute combined scale for the block 0 and 1

        const __m128 d_0_1 = _mm_set1_ps( GGML_CPU_FP16_TO_FP32(x[ib].d) * GGML_CPU_FP16_TO_FP32(y[ib].d) );

        const __m128i tmp_0_1 = _mm_loadu_si128((const __m128i *)x[ib].qs);

        __m128i bx_0 = _mm_and_si128(lowMask, tmp_0_1);

        __m128i by_0 = _mm_loadu_si128((const __m128i *)y[ib].qs);

        bx_0 = _mm_sub_epi8(bx_0, off);

        const __m128i i32_0 = mul_sum_i8_pairs(bx_0, by_0);

        __m128i bx_1 = _mm_and_si128(lowMask, _mm_srli_epi64(tmp_0_1, 4));

        __m128i by_1 = _mm_loadu_si128((const __m128i *)(y[ib].qs + 16));

        bx_1 = _mm_sub_epi8(bx_1, off);

        const __m128i i32_1 = mul_sum_i8_pairs(bx_1, by_1);

        _mm_prefetch(&x[ib] + 2 * sizeof(block_q4_0), _MM_HINT_T0);

        _mm_prefetch(&y[ib] + 2 * sizeof(block_q8_0), _MM_HINT_T0);

        // Compute combined scale for the block 2 and 3

        const __m128 d_2_3 = _mm_set1_ps( GGML_CPU_FP16_TO_FP32(x[ib + 1].d) * GGML_CPU_FP16_TO_FP32(y[ib + 1].d) );

        const __m128i tmp_2_3 = _mm_loadu_si128((const __m128i *)x[ib + 1].qs);

        __m128i bx_2 = _mm_and_si128(lowMask, tmp_2_3);

        __m128i by_2 = _mm_loadu_si128((const __m128i *)y[ib + 1].qs);

        bx_2 = _mm_sub_epi8(bx_2, off);

        const __m128i i32_2 = mul_sum_i8_pairs(bx_2, by_2);

        __m128i bx_3 = _mm_and_si128(lowMask, _mm_srli_epi64(tmp_2_3, 4));

        __m128i by_3 = _mm_loadu_si128((const __m128i *)(y[ib + 1].qs + 16));

        bx_3 = _mm_sub_epi8(bx_3, off);

        const __m128i i32_3 = mul_sum_i8_pairs(bx_3, by_3);

        // Convert int32_t to float

        __m128 p0 = _mm_cvtepi32_ps(i32_0);

        __m128 p1 = _mm_cvtepi32_ps(i32_1);

        __m128 p2 = _mm_cvtepi32_ps(i32_2);

        __m128 p3 = _mm_cvtepi32_ps(i32_3);

        // Apply the scale

        __m128 p0_d = _mm_mul_ps( d_0_1, p0 );

        __m128 p1_d = _mm_mul_ps( d_0_1, p1 );

        __m128 p2_d = _mm_mul_ps( d_2_3, p2 );

        __m128 p3_d = _mm_mul_ps( d_2_3, p3 );

        // Acummulate

        acc_0 = _mm_add_ps(p0_d, acc_0);

        acc_1 = _mm_add_ps(p1_d, acc_1);

        acc_2 = _mm_add_ps(p2_d, acc_2);

        acc_3 = _mm_add_ps(p3_d, acc_3);

    }

    sumf = hsum_float_4x4(acc_0, acc_1, acc_2, acc_3);

#endif

    for (; ib < nb; ++ib) {

        int sumi0 = 0;

        int sumi1 = 0;

        for (int j = 0; j < qk/2; ++j) {

            const int v0 = (x[ib].qs[j] & 0x0F) - 8;

            const int v1 = (x[ib].qs[j] >>   4) - 8;

            sumi0 += (v0 * y[ib].qs[j]);

            sumi1 += (v1 * y[ib].qs[j + qk/2]);

        }

        int sumi = sumi0 + sumi1;

        sumf += sumi*GGML_CPU_FP16_TO_FP32(x[ib].d)*GGML_CPU_FP16_TO_FP32(y[ib].d);

    }

    *s = sumf;

}

void ggml_vec_dot_q4_1_q8_1(int n, float * GGML_RESTRICT s, size_t bs, const void * GGML_RESTRICT vx, size_t bx, const void * GGML_RESTRICT vy, size_t by, int nrc) {

    const int qk = QK8_1;

    const int nb = n / qk;

    assert(n % qk == 0);

    assert(nrc == 1);

    UNUSED(nrc);

    UNUSED(bx);

    UNUSED(by);

    UNUSED(bs);

    const block_q4_1 * GGML_RESTRICT x = vx;

    const block_q8_1 * GGML_RESTRICT y = vy;

    int ib = 0;

#if defined(__AVX2__) || defined(__AVX__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    float summs = 0;

    // Main loop

    for (; ib < nb; ++ib) {

        const float d0 = GGML_CPU_FP16_TO_FP32(x[ib].d);

        const float d1 = GGML_CPU_FP16_TO_FP32(y[ib].d);

        summs += GGML_CPU_FP16_TO_FP32(x[ib].m) * GGML_CPU_FP16_TO_FP32(y[ib].s);

        const __m256 d0v = _mm256_set1_ps( d0 );

        const __m256 d1v = _mm256_set1_ps( d1 );

        // Compute combined scales

        const __m256 d0d1 = _mm256_mul_ps( d0v, d1v );

        // Load 16 bytes, and unpack 4 bit fields into bytes, making 32 bytes

        const __m256i qx = bytes_from_nibbles_32(x[ib].qs);

        const __m256i qy = _mm256_loadu_si256( (const __m256i *)y[ib].qs );

        const __m256 xy = mul_sum_us8_pairs_float(qx, qy);

        // Accumulate d0*d1*x*y

#if defined(__AVX2__)

        acc = _mm256_fmadd_ps( d0d1, xy, acc );

#else

        acc = _mm256_add_ps( _mm256_mul_ps( d0d1, xy ), acc );

#endif

    }

    *s = hsum_float_8(acc) + summs;

#else

    UNUSED(nb);

    UNUSED(x);

    UNUSED(y);

    UNUSED(ib);

    ggml_vec_dot_q4_1_q8_1_generic(n, s, bs, vx, bx, vy, by, nrc);

#endif

}

void ggml_vec_dot_mxfp4_q8_0(int n, float * GGML_RESTRICT s, size_t bs, const void * GGML_RESTRICT vx, size_t bx, const void * GGML_RESTRICT vy, size_t by, int nrc) {

    assert(nrc == 1);

    UNUSED(nrc);

    UNUSED(bx);

    UNUSED(by);

    UNUSED(bs);

    assert(n % QK_MXFP4 == 0);

    static_assert(QK_MXFP4 == QK8_0, "QK_MXFP4 and QK8_0 must be the same");

    const block_mxfp4 * GGML_RESTRICT x = vx;

    const block_q8_0 * GGML_RESTRICT y = vy;

    const int nb = n / QK_MXFP4;

    int ib = 0;

    float sumf = 0;

#if defined __AVX2__

    const __m128i values128 = _mm_loadu_si128((const __m128i*)kvalues_mxfp4);

    const __m128i m4b  = _mm_set1_epi8(0x0f);

    const __m256i mone = _mm256_set1_epi16(1);

    __m256 accum1 = _mm256_setzero_ps();

    __m256 accum2 = _mm256_setzero_ps();

    for (; ib + 1 < nb; ib += 2) {

        const __m128i q4bits_1 = _mm_loadu_si128((const __m128i*)x[ib + 0].qs);

        const __m128i q4bits_2 = _mm_loadu_si128((const __m128i*)x[ib + 1].qs);

        const __m256i q8b_1 = _mm256_loadu_si256((const __m256i *)y[ib + 0].qs);

        const __m256i q8b_2 = _mm256_loadu_si256((const __m256i *)y[ib + 1].qs);

        const __m256i q4b_1 = MM256_SET_M128I(_mm_shuffle_epi8(values128, _mm_and_si128(_mm_srli_epi16(q4bits_1, 4), m4b)),

                                              _mm_shuffle_epi8(values128, _mm_and_si128(q4bits_1, m4b)));

        const __m256i q4b_2 = MM256_SET_M128I(_mm_shuffle_epi8(values128, _mm_and_si128(_mm_srli_epi16(q4bits_2, 4), m4b)),

                                              _mm_shuffle_epi8(values128, _mm_and_si128(q4bits_2, m4b)));

        const __m256i p16_1 = mul_add_epi8(q4b_1, q8b_1);

        const __m256i p16_2 = mul_add_epi8(q4b_2, q8b_2);

        const __m256i p_1 = _mm256_madd_epi16(p16_1, mone);

        const __m256i p_2 = _mm256_madd_epi16(p16_2, mone);

        accum1 = _mm256_fmadd_ps(_mm256_set1_ps(GGML_CPU_FP16_TO_FP32(y[ib + 0].d)*GGML_E8M0_TO_FP32_HALF(x[ib + 0].e)),

                _mm256_cvtepi32_ps(p_1), accum1);

        accum2 = _mm256_fmadd_ps(_mm256_set1_ps(GGML_CPU_FP16_TO_FP32(y[ib + 1].d)*GGML_E8M0_TO_FP32_HALF(x[ib + 1].e)),

                _mm256_cvtepi32_ps(p_2), accum2);

    }

    sumf = hsum_float_8(_mm256_add_ps(accum1, accum2));

#elif defined __AVX__

    const __m128i values128 = _mm_loadu_si128((const __m128i*)kvalues_mxfp4);

    const __m128i m4b  = _mm_set1_epi8(0x0f);

    __m256 accum = _mm256_setzero_ps();

    for (; ib + 1 < nb; ib += 2) {

        const __m128i q4bits_1 = _mm_loadu_si128((const __m128i *)x[ib + 0].qs);

        const __m128i q4bits_2 = _mm_loadu_si128((const __m128i *)x[ib + 1].qs);

        const __m128i q8b_1_0 = _mm_loadu_si128((const __m128i *)y[ib + 0].qs);

        const __m128i q8b_1_1 = _mm_loadu_si128((const __m128i *)y[ib + 0].qs + 1);

        const __m128i q8b_2_0 = _mm_loadu_si128((const __m128i *)y[ib + 1].qs);

        const __m128i q8b_2_1 = _mm_loadu_si128((const __m128i *)y[ib + 1].qs + 1);

        const __m128i q4b_1_0 = _mm_shuffle_epi8(values128, _mm_and_si128(q4bits_1, m4b));

        const __m128i q4b_1_1 = _mm_shuffle_epi8(values128, _mm_and_si128(_mm_srli_epi16(q4bits_1, 4), m4b));

        const __m128i q4b_2_0 = _mm_shuffle_epi8(values128, _mm_and_si128(q4bits_2, m4b));

        const __m128i q4b_2_1 = _mm_shuffle_epi8(values128, _mm_and_si128(_mm_srli_epi16(q4bits_2, 4), m4b));

        const __m256 p = mul_sum_i8_quad_float(q4b_1_0, q4b_1_1, q4b_2_0, q4b_2_1, q8b_1_0, q8b_1_1, q8b_2_0, q8b_2_1);

        const __m256 deltas = quad_mx_delta_float(x[ib].e, y[ib].d, x[ib + 1].e, y[ib + 1].d);

        accum = _mm256_add_ps(_mm256_mul_ps(deltas, p), accum);

    }

    sumf = hsum_float_8(accum);

#endif

    for (; ib < nb; ++ib) {

        const float d = GGML_CPU_FP16_TO_FP32(y[ib].d)*GGML_E8M0_TO_FP32_HALF(x[ib].e);

        int sumi1 = 0;

        int sumi2 = 0;

        for (int j = 0; j < QK_MXFP4/2; ++j) {

            sumi1 += y[ib].qs[j +          0] * kvalues_mxfp4[x[ib].qs[j] & 0xf];

            sumi2 += y[ib].qs[j + QK_MXFP4/2] * kvalues_mxfp4[x[ib].qs[j] >>  4];

        }

        sumf += d * (sumi1 + sumi2);

    }

    *s = sumf;

}

void ggml_vec_dot_q5_0_q8_0(int n, float * GGML_RESTRICT s, size_t bs, const void * GGML_RESTRICT vx, size_t bx, const void * GGML_RESTRICT vy, size_t by, int nrc) {

    const int qk = QK8_0;

    const int nb = n / qk;

    int ib = 0;

    assert(n % qk == 0);

    assert(qk == QK5_0);

    assert(nrc == 1);

    UNUSED(nrc);

    UNUSED(bx);

    UNUSED(by);

    UNUSED(bs);

    const block_q5_0 * GGML_RESTRICT x = vx;

    const block_q8_0 * GGML_RESTRICT y = vy;

#if defined(__AVX2__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    // Main loop

    for (; ib < nb; ++ib) {

        /* Compute combined scale for the block */

        const __m256 d = _mm256_set1_ps(GGML_CPU_FP16_TO_FP32(x[ib].d) * GGML_CPU_FP16_TO_FP32(y[ib].d));

        __m256i qx = bytes_from_nibbles_32(x[ib].qs);

        __m256i bxhi = bytes_from_bits_32(x[ib].qh);

        bxhi = _mm256_andnot_si256(bxhi, _mm256_set1_epi8((char)0xF0));

        qx = _mm256_or_si256(qx, bxhi);

        __m256i qy = _mm256_loadu_si256((const __m256i *)y[ib].qs);

        const __m256 q = mul_sum_i8_pairs_float(qx, qy);

        /* Multiply q with scale and accumulate */

        acc = _mm256_fmadd_ps(d, q, acc);

    }

    *s = hsum_float_8(acc);

#elif defined(__AVX__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    __m128i mask = _mm_set1_epi8((char)0xF0);

    // Main loop

    for (; ib < nb; ++ib) {

        /* Compute combined scale for the block */

        const __m256 d = _mm256_set1_ps(GGML_CPU_FP16_TO_FP32(x[ib].d) * GGML_CPU_FP16_TO_FP32(y[ib].d));

        __m256i bx_0 = bytes_from_nibbles_32(x[ib].qs);

        const __m256i bxhi = bytes_from_bits_32(x[ib].qh);

        __m128i bxhil = _mm256_castsi256_si128(bxhi);

        __m128i bxhih = _mm256_extractf128_si256(bxhi, 1);

        bxhil = _mm_andnot_si128(bxhil, mask);

        bxhih = _mm_andnot_si128(bxhih, mask);

        __m128i bxl = _mm256_castsi256_si128(bx_0);

        __m128i bxh = _mm256_extractf128_si256(bx_0, 1);

        bxl = _mm_or_si128(bxl, bxhil);

        bxh = _mm_or_si128(bxh, bxhih);

        bx_0 = MM256_SET_M128I(bxh, bxl);

        const __m256i by_0 = _mm256_loadu_si256((const __m256i *)y[ib].qs);

        const __m256 q = mul_sum_i8_pairs_float(bx_0, by_0);

        /* Multiply q with scale and accumulate */

        acc = _mm256_add_ps(_mm256_mul_ps(d, q), acc);

    }

    *s = hsum_float_8(acc);

#else

    UNUSED(nb);

    UNUSED(ib);

    UNUSED(x);

    UNUSED(y);

    ggml_vec_dot_q5_0_q8_0_generic(n, s, bs, vx, bx, vy, by, nrc);

#endif

}

void ggml_vec_dot_q5_1_q8_1(int n, float * GGML_RESTRICT s, size_t bs, const void * GGML_RESTRICT vx, size_t bx, const void * GGML_RESTRICT vy, size_t by, int nrc) {

    const int qk = QK8_1;

    const int nb = n / qk;

    int ib = 0;

    assert(n % qk == 0);

    assert(qk == QK5_1);

    assert(nrc == 1);

    UNUSED(nrc);

    UNUSED(bx);

    UNUSED(by);

    UNUSED(bs);

    const block_q5_1 * GGML_RESTRICT x = vx;

    const block_q8_1 * GGML_RESTRICT y = vy;

#if defined(__AVX2__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    float summs = 0.0f;

    // Main loop

    for (; ib < nb; ++ib) {

        const __m256 dx = _mm256_set1_ps(GGML_CPU_FP16_TO_FP32(x[ib].d));

        summs += GGML_CPU_FP16_TO_FP32(x[ib].m) * GGML_CPU_FP16_TO_FP32(y[ib].s);

        __m256i qx = bytes_from_nibbles_32(x[ib].qs);

        __m256i bxhi = bytes_from_bits_32(x[ib].qh);

        bxhi = _mm256_and_si256(bxhi, _mm256_set1_epi8(0x10));

        qx = _mm256_or_si256(qx, bxhi);

        const __m256 dy = _mm256_set1_ps(GGML_CPU_FP16_TO_FP32(y[ib].d));

        const __m256i qy = _mm256_loadu_si256((const __m256i *)y[ib].qs);

        const __m256 q = mul_sum_us8_pairs_float(qx, qy);

        acc = _mm256_fmadd_ps(q, _mm256_mul_ps(dx, dy), acc);

    }

    *s = hsum_float_8(acc) + summs;

#elif defined(__AVX__)

    // Initialize accumulator with zeros

    __m256 acc = _mm256_setzero_ps();

    __m128i mask = _mm_set1_epi8(0x10);

    float summs = 0.0f;

    // Main loop

    for (; ib < nb; ++ib) {

        const __m256 dx = _mm256_set1_ps(GGML_CPU_FP16_TO_FP32(x[ib].d));

        summs += GGML_CPU_FP16_TO_FP32(x[ib].m) * GGML_CPU_FP16_TO_FP32(y[ib].s);

        __m256i bx_0 = bytes_from_nibbles_32(x[ib].qs);

        const __m256i bxhi = bytes_from_bits_32(x[ib].qh);

        __m128i bxhil = _mm256_castsi256_si128(bxhi);