如何加速这个 LUT 查找的直方图?

This version takes advantage of those ideas and a few histogram tricks, and compiles to asm that looks good (Godbolt compiler explorer: gcc6.2 -O3 -march=haswell (AVX2)):

// untested
//#include <algorithm>
#include <stdint.h>

const int PLANES = 1000;
void use_bins(uint16_t bins[PLANES][32]);   // pass the result to an extern function so it doesn't optimize away

                       // 65536 or higher triggers the static_assert
alignas(64) static uint16_t a[PLANES][1000];  // static/global, I guess?

void lut_and_histogram(uint8_t __restrict__ lut[32768])
{

    alignas(16) uint16_t bins[PLANES][32];  // don't zero the whole thing up front: that would evict more data from cache than necessary
                                            // Better would be zeroing the relevant plane of each bin right before using.
                                            // you pay the rep stosq startup overhead more times, though.

    for (int i = 0; i < PLANES;i++) {
        alignas(16) uint16_t tmpbins[4][32] = {0};
        constexpr int a_elems = sizeof(a[0])/sizeof(uint16_t);
        static_assert(a_elems > 1, "someone changed a[] into a* and forgot to update this code");
        static_assert(a_elems <= UINT16_MAX, "bins could overflow");
        const uint16_t *ai = a[i];
        for (int j = 0 ; j<a_elems ; j+=4) { //hashing a[i] to 32 bins.
            // Unrolling to separate bin arrays reduces serial dependencies
            // to avoid bottlenecks when the same bin is used repeatedly.
            // This has to be balanced against using too much L1 cache for the bins.

            // TODO: load a vector of data from ai[j] and unpack it with pextrw.
            // even just loading a uint64_t and unpacking it to 4 uint16_t would help.
            tmpbins[0][ lut[ai[j+0]] ]++;
            tmpbins[1][ lut[ai[j+1]] ]++;
            tmpbins[2][ lut[ai[j+2]] ]++;
            tmpbins[3][ lut[ai[j+3]] ]++;
            static_assert(a_elems % 4 == 0, "unroll factor doesn't divide a element count");
        }

        // TODO: do multiple a[i] in parallel instead of slicing up a single run.
        for (int k = 0 ; k<32 ; k++) {
            // gcc does auto-vectorize this with a short fully-unrolled VMOVDQA / VPADDW x3
            bins[i][k] = tmpbins[0][k] + tmpbins[1][k] +
                         tmpbins[2][k] + tmpbins[3][k];
        }
    }
  
    // do something with bins.  An extern function stops it from optimizing away.
    use_bins(bins);
}

内循环 asm 如下所示:

.L2:
    movzx   ecx, WORD PTR [rdx]
        add     rdx, 8                    # pointer increment over ai[]
    movzx   ecx, BYTE PTR [rsi+rcx]
    add     WORD PTR [rbp-64272+rcx*2], 1     # memory-destination increment of a histogram element
    movzx   ecx, WORD PTR [rdx-6]
    movzx   ecx, BYTE PTR [rsi+rcx]
    add     WORD PTR [rbp-64208+rcx*2], 1
    ... repeated twice more

    const size_t nb = 32;  // number of buckets
    const int VEC_WIDTH = 16;  // sizeof(__m512i) / sizeof(uint32_t)
    alignas(__m512i) uint32_t counts[nb * VEC_WIDTH] = {0};

// then in your histo loop

    __m512i idx = ...;  // in this case from LUT lookups
    idx = _mm512_add_epi32(idx, _mm512_setr_epi32(
       0*nb, 1*nb,  2*nb,  3*nb,   4*nb,  5*nb,  6*nb,  7*nb,
       8*nb, 9*nb, 10*nb, 11*nb,  12*nb, 13*nb, 14*nb, 15*nb));
       // note these are C array indexes, not byte offsets
    __m512i vc = _mm512_i32gather_epi32(idx, counts, sizeof(counts[0]));
    vc = _mm512_add_epi32(vc, _mm512_set1_epi32(1));
    _mm512_i32scatter_epi32(counts, idx, vc, sizeof(counts[0]));

// Optionally with size_t nb as an arg
// also optionally use restrict if you never reduce in-place, into the bottom of the input.
void reduce_counts(int *output, const int *counts)
{
    for (int i = 0 ; i < nb - (VEC_WIDTH-1) ; i+=VEC_WIDTH) {
        __m512i v = _mm512_load_si512(&counts[i]);  // aligned load, full cache line
        // optional: unroll this and accumulate two vectors in parallel for better spatial locality and more ILP
        for (int offset = nb; offset < nb*VEC_WIDTH ; offset+=nb) {
            __m512i tmp = _mm512_loadu_si512(&counts[i + offset]);
            v = _mm512_add_epi32(v, tmp);
        }
        _mm512_storeu_si512(&output[i], v);
    }

   // if  nb isn't a multiple of the vector width, do some cleanup here
   // possibly using a masked store to write into a final odd-sized destination
}