Yesterday, we looked at optimising an SDI CRC computation. After a bunch of mathematical optimisations, we reached this point:

uint64_t pack60(const uint16_t* data) {
    return (((uint64_t)data[ 0]) <<  4) ^
           (((uint64_t)data[ 2]) << 14) ^
           (((uint64_t)data[ 4]) << 24) ^
           (((uint64_t)data[ 6]) << 34) ^
           (((uint64_t)data[ 8]) << 44) ^
           (((uint64_t)data[10]) << 54);
}

__m128i pack120(const uint16_t* data) {
    return _mm_set_epi64x(pack60(data + 12) >> 4, pack60(data));
}

__m128i xor_clmul(__m128i a, __m128i b) {
    return _mm_xor_si128(_mm_clmulepi64_si128(a, b, 0x00),
                         _mm_clmulepi64_si128(a, b, 0x11));
}

void crc_sdi(uint32_t* crcs, const uint16_t* data, size_t n) {
    __m128i c = _mm_cvtsi32_si128(crcs[0]);
    __m128i y = _mm_cvtsi32_si128(crcs[1]);
    { // *= x^-14 semi-mod P
        __m128i k = _mm_cvtsi32_si128(
            0x9f380000 /* x^-14-(64-18)-32-1 mod P */);
        c = _mm_clmulepi64_si128(c, k, 0x00);
        y = _mm_clmulepi64_si128(y, k, 0x00);
    }
    for (size_t i = 0; i < n; i += 24) {
        { // *= x^120 semi-mod P
            __m128i k = _mm_setr_epi32(
                0, 0x4b334000 /* x^120+64-1 mod P */,
                0, 0x96d30000 /* x^120-1    mod P */);
            c = xor_clmul(c, k);
            y = xor_clmul(y, k);
        }
        { // +=
            c = _mm_xor_si128(c, pack120(data + i));
            y = _mm_xor_si128(y, pack120(data + i + 1));
        }
    }
    { // *= x^14 semi-mod P
        __m128i k = _mm_setr_epi32(
            0, 0x14980000 /* x^14+64-1 mod P */,
            0, 0x00040000 /* x^14-1    mod P */);
        c = xor_clmul(c, k);
        y = xor_clmul(y, k);
    }
    { // mod P
        __m128i k = _mm_setr_epi32( /* x^128-1 div P */
            0x14980559, 0x4c9bb5d5,
            0x80040000, 0x5e405011);
        c = _mm_xor_si128(_mm_srli_si128(xor_clmul(c, k), 8),
                          _mm_clmulepi64_si128(c, k, 0x01));
        y = _mm_xor_si128(_mm_srli_si128(xor_clmul(y, k), 8),
                          _mm_clmulepi64_si128(y, k, 0x01));
        __m128i P = _mm_cvtsi32_si128(0x46001 /* P */);
        c = _mm_clmulepi64_si128(c, P, 0x00);
        y = _mm_clmulepi64_si128(y, P, 0x00);
    }
    crcs[0] = _mm_cvtsi128_si32(c);
    crcs[1] = _mm_cvtsi128_si32(y);
}

This ran at approximately 11 bits per cycle:

VersionBroadwellSkylakeIce Lake
_mm_clmulepi64_si12810.811.012.8

From here, the question is how to best implement pack120. Compared to yesterday, there is a very different approach based around _mm_shuffle_epi8:

void crc_sdi_pack(uint32_t* crcs, const uint16_t* data, size_t n) {
...
        { // +=
            __m128i d1 = _mm_loadu_si128((__m128i*)(data + i));
            __m128i d2 = _mm_loadu_si128((__m128i*)(data + i + 8));
            __m128i d3 = _mm_loadu_si128((__m128i*)(data + i + 16));
            __m128i k1 = _mm_setr_epi16(16, 16, 64, 64, 1, 1, 4, 4);
            __m128i k2 = _mm_setr_epi16(16,  0, 64,  0, 1, 0, 4, 0);
            __m128i k3 = _mm_setr_epi16( 0, 16,  0, 64, 0, 1, 0, 4);
            d1 = _mm_mullo_epi16(k1, d1);
            __m128i k4=_mm_setr_epi8( 0, 1, -1,  8,  9, -1, -1, -1,
                                      2, 3, -1, 10, 11, -1, -1, -1);
            __m128i k5=_mm_setr_epi8(-1, 4,  5, -1, 12, 13, -1, -1,
                                     -1, 6,  7, -1, 14, 15, -1, -1);
            d1 = _mm_shuffle_epi8(d1, k4) ^ _mm_shuffle_epi8(d1, k5);
            __m128i d1m; // Force a store to memory
            asm("vmovdqa %1, %0" : "=m"(d1m) : "x"(d1) : "memory");
            __m128i cd = _mm_packus_epi32(_mm_madd_epi16(d2, k2),
                                          _mm_madd_epi16(d3, k2));
            __m128i yd = _mm_packus_epi32(_mm_madd_epi16(d2, k3),
                                          _mm_madd_epi16(d3, k3));
            __m128i k6=_mm_setr_epi8(-1, -1, -1, -1, -1,  0,  1, -1,
                                      4,  5,  8,  9, -1, 12, 13, -1);
            __m128i k7=_mm_setr_epi8(-1, -1, -1, -1, -1, -1,  2,  3,
                                     -1,  6,  7, 10, 11, -1, 14, 15);
            cd = _mm_shuffle_epi8(cd, k6) ^ _mm_shuffle_epi8(cd, k7)
               ^ _mm_loadu_si64(&d1m);
            yd = _mm_shuffle_epi8(yd, k6) ^ _mm_shuffle_epi8(yd, k7)
               ^ _mm_loadu_si64(8 + (char*)&d1m);
            c = _mm_xor_si128(c, cd);
            y = _mm_xor_si128(y, yd);
        }
...
}

This is very impressive, getting us to the region of 20 bits per cycle:

VersionBroadwellSkylakeIce Lake
_mm_shuffle_epi820.219.723.4

If targetting AVX2, then some pairs of operations can be collapsed:

__m256i broadcast128(const uint16_t* data) {
    __m256i result;
    asm("vbroadcasti128 %1, %0" : "=x"(result) :
                                   "m"(*(const __m128i*)data));
    return result;
}

void crc_sdi_pack(uint32_t* crcs, const uint16_t* data, size_t n) {
...
        { // +=
            __m128i d1 = _mm_loadu_si128((__m128i*)(data + i));
            __m256i d2 = broadcast128(data + i + 8);
            __m256i d3 = broadcast128(data + i + 16);
            __m128i k1 = _mm_setr_epi16(
                16, 16, 64, 64, 1, 1, 4, 4);
            __m256i k23 = _mm256_setr_epi16(
                16,  0, 64,  0, 1, 0, 4, 0,
                 0, 16,  0, 64, 0, 1, 0, 4);
            d1 = _mm_mullo_epi16(k1, d1);
            __m128i k4 = _mm_setr_epi8(
                0,  1, -1,  8,  9, -1, -1, -1,
                2,  3, -1, 10, 11, -1, -1, -1);
            __m128i k5 = _mm_setr_epi8(
                -1, 4,  5, -1, 12, 13, -1, -1,
                -1, 6,  7, -1, 14, 15, -1, -1);
            d1 = _mm_shuffle_epi8(d1, k4) ^ _mm_shuffle_epi8(d1, k5);
            __m128i d1m; // Force a store to memory
            asm("vmovdqa %1, %0" : "=m"(d1m) : "x"(d1) : "memory");
            __m256i cdyd = _mm256_packus_epi32(
                _mm256_madd_epi16(d2, k23),
                _mm256_madd_epi16(d3, k23));
            __m256i k6 = _mm256_setr_epi8(
                -1, -1, -1, -1, -1,  0,  1, -1,
                 4,  5,  8,  9, -1, 12, 13, -1,
                -1, -1, -1, -1, -1,  0,  1, -1, 
                 4,  5,  8,  9, -1, 12, 13, -1);
            __m256i k7 = _mm256_setr_epi8(
                -1, -1, -1, -1, -1, -1,  2,  3,
                -1,  6,  7, 10, 11, -1, 14, 15,
                -1, -1, -1, -1, -1, -1,  2,  3,
                -1,  6,  7, 10, 11, -1, 14, 15);
            cdyd = _mm256_shuffle_epi8(cdyd, k6)
                 ^ _mm256_shuffle_epi8(cdyd, k7);
            __m256i m; // Force a store to memory
            asm("vmovdqa %1, %0" : "=m"(m) : "x"(cdyd) : "memory");
            __m128i cd = _mm256_castsi256_si128(cdyd)
                       ^ _mm_loadu_si64(&d1m);
            __m128i yd = _mm_loadu_si128(1 + (__m128i*)&m)
                       ^ _mm_loadu_si64(8 + (char*)&d1m);
            c = _mm_xor_si128(c, cd);
            y = _mm_xor_si128(y, yd);
        }
...
}

This gains us a few bits per cycle on all three processors:

VersionBroadwellSkylakeIce Lake
AVX222.323.225.4

If targetting AVX512, then the xor3 trick from yesterday can be applied on top:

__m128i xor3(__m128i a, __m128i b, __m128i c) {
    return _mm_ternarylogic_epi64(a, b, c, 0x96);
}

void crc_sdi_pack(uint32_t* crcs, const uint16_t* data, size_t n) {
...
    for (size_t i = 0; i < n; i += 24) {
        // *= x^120 semi-mod P
        // +=
        __m128i k = _mm_setr_epi32(
            0, 0x4b334000 /* x^120+64-1 mod P */,
            0, 0x96d30000 /* x^120-1    mod P */);
        __m128i d1 = _mm_loadu_si128((__m128i*)(data + i));
        __m256i d2 = broadcast128(data + i + 8);
        __m256i d3 = broadcast128(data + i + 16);
        __m128i k1 = _mm_setr_epi16(
            16, 16, 64, 64, 1, 1, 4, 4);
        __m256i k23 = _mm256_setr_epi16(
            16,  0, 64,  0, 1, 0, 4, 0,
             0, 16,  0, 64, 0, 1, 0, 4);
        d1 = _mm_mullo_epi16(k1, d1);
        __m128i k4 = _mm_setr_epi8(
            0,  1, -1,  8,  9, -1, -1, -1,
            2,  3, -1, 10, 11, -1, -1, -1);
        __m128i k5 = _mm_setr_epi8(
            -1, 4,  5, -1, 12, 13, -1, -1,
            -1, 6,  7, -1, 14, 15, -1, -1);
        d1 = _mm_shuffle_epi8(d1, k4) ^ _mm_shuffle_epi8(d1, k5);
        __m128i d1m; // Force a store to memory
        asm("vmovdqa %1, %0" : "=m"(d1m) : "x"(d1) : "memory");
        __m256i cdyd = _mm256_packus_epi32(
            _mm256_madd_epi16(d2, k23),
            _mm256_madd_epi16(d3, k23));
        __m256i k6 = _mm256_setr_epi8(
            -1, -1, -1, -1, -1,  0,  1, -1,
             4,  5,  8,  9, -1, 12, 13, -1,
            -1, -1, -1, -1, -1,  0,  1, -1, 
             4,  5,  8,  9, -1, 12, 13, -1);
        __m256i k7 = _mm256_setr_epi8(
            -1, -1, -1, -1, -1, -1,  2,  3,
            -1,  6,  7, 10, 11, -1, 14, 15,
            -1, -1, -1, -1, -1, -1,  2,  3,
            -1,  6,  7, 10, 11, -1, 14, 15);
        cdyd = _mm256_shuffle_epi8(cdyd, k6)
             ^ _mm256_shuffle_epi8(cdyd, k7);
        __m256i m; // Force a store to memory
        asm("vmovdqa %1, %0" : "=m"(m) : "x"(cdyd) : "memory");
        __m128i cd = _mm256_castsi256_si128(cdyd)
                   ^ _mm_loadu_si64(&d1m);
        __m128i yd = _mm_loadu_si128(1 + (__m128i*)&m)
                   ^ _mm_loadu_si64(8 + (char*)&d1m);
        c = xor3(_mm_clmulepi64_si128(c, k, 0x00),
                 _mm_clmulepi64_si128(c, k, 0x11), cd);
        y = xor3(_mm_clmulepi64_si128(y, k, 0x00),
                 _mm_clmulepi64_si128(y, k, 0x11), yd);
    }
...
}

On applicable processors, this gains another few bits per cycle:

VersionBroadwellSkylakeIce Lake
AVX512N/A25.828.9