代码

// 4D strict dominance via prefix-bitset blocking + AVX2 (optimized).
#pragma GCC target("avx2")
#include <cstring>
#include <cstdlib>
#include <immintrin.h>

typedef unsigned u32;
typedef unsigned long long u64;

static int order[4][100010];
static int first_pos[4][100010];
static int cnt[100010];
static int seen[100010];
static int stamp = 0;

static inline u32 pc256(__m256i v){
    const __m256i lut = _mm256_setr_epi8(
        0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
        0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4);
    const __m256i m4 = _mm256_set1_epi8(0x0F);
    __m256i lo = _mm256_and_si256(v, m4);
    __m256i hi = _mm256_and_si256(_mm256_srli_epi16(v, 4), m4);
    lo = _mm256_shuffle_epi8(lut, lo);
    hi = _mm256_shuffle_epi8(lut, hi);
    __m256i s = _mm256_add_epi8(lo, hi);
    __m256i sad = _mm256_sad_epu8(s, _mm256_setzero_si256());
    return (u32)(_mm256_extract_epi64(sad,0) + _mm256_extract_epi64(sad,1)
               + _mm256_extract_epi64(sad,2) + _mm256_extract_epi64(sad,3));
}

void count_4d(int n, const unsigned *x[4], unsigned *out) {
    int NW = (n + 63) >> 6;
    const int B = 32;
    int M = (n + B - 1) / B;

    for (int d = 0; d < 4; d++) {
        const unsigned *xd = x[d];
        memset(cnt, 0, sizeof(int) * (n + 1));
        for (int i = 0; i < n; i++) cnt[xd[i]]++;
        int s = 0;
        for (int v = 0; v <= (int)n; v++) { int c = cnt[v]; cnt[v] = s; s += c; }
        for (int i = 0; i < n; i++) first_pos[d][i] = cnt[xd[i]];
        memset(cnt, 0, sizeof(int) * (n + 1));
        for (int i = 0; i < n; i++) cnt[xd[i]]++;
        for (int v = 1; v <= (int)n; v++) cnt[v] += cnt[v-1];
        for (int i = n - 1; i >= 0; i--) order[d][--cnt[xd[i]]] = i;
    }

    u64 *pref = (u64*)malloc((size_t)4 * (M + 1) * NW * sizeof(u64));
    for (int d = 0; d < 4; d++) {
        u64 *p = pref + (size_t)d * (M + 1) * NW;
        memset(p, 0, (size_t)(M + 1) * NW * sizeof(u64));
        for (int b = 0; b < M; b++) {
            u64 *cur = p + (size_t)(b + 1) * NW;
            u64 *prev = p + (size_t)b * NW;
            memcpy(cur, prev, (size_t)NW * sizeof(u64));
            int e = (b + 1) * B; if (e > n) e = n;
            for (int t = b * B; t < e; t++) {
                int pt = order[d][t];
                cur[pt >> 6] |= 1ULL << (pt & 63);
            }
        }
    }

    u64 *p0 = pref;
    u64 *p1 = pref + (size_t)(M + 1) * NW;
    u64 *p2 = pref + (size_t)2 * (M + 1) * NW;
    u64 *p3 = pref + (size_t)3 * (M + 1) * NW;

    // process points in x1-sorted order (cache-friendly for dim 0)
    for (int pos = 0; pos < n; pos++) {
        int j = order[0][pos];
        {
            int pf = pos + 16;
            if (pf < n) {
                int jf = order[0][pf];
                _mm_prefetch((const char*)(p0 + (size_t)(first_pos[0][jf] / B) * NW), _MM_HINT_T0);
                _mm_prefetch((const char*)(p1 + (size_t)(first_pos[1][jf] / B) * NW), _MM_HINT_T0);
                _mm_prefetch((const char*)(p2 + (size_t)(first_pos[2][jf] / B) * NW), _MM_HINT_T0);
                _mm_prefetch((const char*)(p3 + (size_t)(first_pos[3][jf] / B) * NW), _MM_HINT_T0);
            }
        }
        const u64 *q0 = p0 + (size_t)(first_pos[0][j] / B) * NW;
        const u64 *q1 = p1 + (size_t)(first_pos[1][j] / B) * NW;
        const u64 *q2 = p2 + (size_t)(first_pos[2][j] / B) * NW;
        const u64 *q3 = p3 + (size_t)(first_pos[3][j] / B) * NW;
        u32 c = 0;
        int w = 0;
        for (; w + 8 <= NW; w += 8) {
            __m256i a0 = _mm256_loadu_si256((const __m256i*)(q0 + w));
            __m256i a1 = _mm256_loadu_si256((const __m256i*)(q0 + w + 4));
            __m256i b0 = _mm256_loadu_si256((const __m256i*)(q1 + w));
            __m256i b1 = _mm256_loadu_si256((const __m256i*)(q1 + w + 4));
            __m256i c0 = _mm256_loadu_si256((const __m256i*)(q2 + w));
            __m256i c1 = _mm256_loadu_si256((const __m256i*)(q2 + w + 4));
            __m256i d0 = _mm256_loadu_si256((const __m256i*)(q3 + w));
            __m256i d1 = _mm256_loadu_si256((const __m256i*)(q3 + w + 4));
            __m256i r0 = _mm256_and_si256(_mm256_and_si256(a0, b0), _mm256_and_si256(c0, d0));
            __m256i r1 = _mm256_and_si256(_mm256_and_si256(a1, b1), _mm256_and_si256(c1, d1));
            c += pc256(r0) + pc256(r1);
        }
        for (; w < NW; w++) c += (u32)__builtin_popcountll(q0[w] & q1[w] & q2[w] & q3[w]);

        stamp++;
        int fp0 = first_pos[0][j];
        for (int t = (fp0 / B) * B; t < fp0; t++) {
            int q = order[0][t];
            if (x[1][q] < x[1][j] && x[2][q] < x[2][j] && x[3][q] < x[3][j])
                if (seen[q] != stamp) { seen[q] = stamp; c++; }
        }
        int fp1 = first_pos[1][j];
        for (int t = (fp1 / B) * B; t < fp1; t++) {
            int q = order[1][t];
            if (x[0][q] < x[0][j] && x[2][q] < x[2][j] && x[3][q] < x[3][j])
                if (seen[q] != stamp) { seen[q] = stamp; c++; }
        }
        int fp2 = first_pos[2][j];
        for (int t = (fp2 / B) * B; t < fp2; t++) {
            int q = order[2][t];
            if (x[0][q] < x[0][j] && x[1][q] < x[1][j] && x[3][q] < x[3][j])
                if (seen[q] != stamp) { seen[q] = stamp; c++; }
        }
        int fp3 = first_pos[3][j];
        for (int t = (fp3 / B) * B; t < fp3; t++) {
            int q = order[3][t];
            if (x[0][q] < x[0][j] && x[1][q] < x[1][j] && x[2][q] < x[2][j])
                if (seen[q] != stamp) { seen[q] = stamp; c++; }
        }
        out[j] = c;
    }
    free(pref);
}

评测结果