代码

// 5D strict dominance via prefix-bitset blocking + AVX2.
// Port of the 4D bitset engine (1010a): 5 prefix bitsets, brute force checks 4 "other dims".
#pragma GCC target("avx2")
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <sys/mman.h>
#include <immintrin.h>

typedef unsigned u32;
typedef unsigned long long u64;

static int order[5][300010];
static int first_pos[5][300010];
static int cnt[300010];
static int seen[300010];
static int ordrank[5][5][300010];
static int stamp = 0;
static u32 partial[300010];
static int proc[300010];

static const unsigned char pc_lut_bytes[32] = {
        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};
static inline __m256i pc_bytes(__m256i v){
    const __m256i lut = _mm256_loadu_si256((const __m256i*)pc_lut_bytes);
    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);
    return _mm256_add_epi8(lo, hi);
}
static inline u32 sum4(__m256i v){
    return (u32)(_mm256_extract_epi64(v,0) + _mm256_extract_epi64(v,1)
               + _mm256_extract_epi64(v,2) + _mm256_extract_epi64(v,3));
}

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

    for (int d = 0; d < 5; 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;
    }

    for (int d = 0; d < 5; d++)
        for (int a = 0; a < 5; a++)
            for (int t = 0; t < n; t++)
                ordrank[d][a][t] = first_pos[a][order[d][t]];

    size_t prefsz = (size_t)5 * (M + 1) * NW * sizeof(u64);
    const int S = 8192;
    {
        int p = 0;
        for (int s = 0; s < n; s += S) {
            int e = s + S; if (e > n) e = n;
            for (int t = s; t < e; t++) proc[p++] = order[0][t];
            std::sort(proc + p - (e - s), proc + p, [](int a, int b){ return first_pos[1][a] < first_pos[1][b]; });
        }
    }

    u64 *pref = (u64*)mmap(NULL, prefsz, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS|MAP_POPULATE, -1, 0);
    if (pref == (u64*)MAP_FAILED) pref = (u64*)malloc(prefsz);
    for (int d = 0; d < 5; d++) {
        u64 *p = pref + (size_t)d * (M + 1) * NW;
        memset(p, 0, (size_t)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;
            int w = 0;
            for (; w + 16 <= NW; w += 16) {
                _mm256_storeu_si256((__m256i*)(cur + w),     _mm256_loadu_si256((const __m256i*)(prev + w)));
                _mm256_storeu_si256((__m256i*)(cur + w + 4), _mm256_loadu_si256((const __m256i*)(prev + w + 4)));
                _mm256_storeu_si256((__m256i*)(cur + w + 8), _mm256_loadu_si256((const __m256i*)(prev + w + 8)));
                _mm256_storeu_si256((__m256i*)(cur + w + 12), _mm256_loadu_si256((const __m256i*)(prev + w + 12)));
            }
            for (; w < NW; w++) cur[w] = prev[w];
            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;
    u64 *p4 = pref + (size_t)4 * (M + 1) * NW;

    // ---- pass 1: AND + popcount, 2 points at a time ----
    for (int pos = 0; pos < n; pos += 2) {
        int j0 = proc[pos];
        int j1 = (pos + 1 < n) ? proc[pos + 1] : -1;
        {
            int pf = pos + 16;
            if (pf < n) {
                int jf = proc[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);
                _mm_prefetch((const char*)(p4 + (size_t)(first_pos[4][jf] / B) * NW), _MM_HINT_T0);
            }
        }
        const u64 *q0a = p0 + (size_t)(first_pos[0][j0] / B) * NW;
        const u64 *q1a = p1 + (size_t)(first_pos[1][j0] / B) * NW;
        const u64 *q2a = p2 + (size_t)(first_pos[2][j0] / B) * NW;
        const u64 *q3a = p3 + (size_t)(first_pos[3][j0] / B) * NW;
        const u64 *q4a = p4 + (size_t)(first_pos[4][j0] / B) * NW;
        const u64 *q0b = q0a, *q1b = q1a, *q2b = q2a, *q3b = q3a, *q4b = q4a;
        if (j1 >= 0) {
            q0b = p0 + (size_t)(first_pos[0][j1] / B) * NW;
            q1b = p1 + (size_t)(first_pos[1][j1] / B) * NW;
            q2b = p2 + (size_t)(first_pos[2][j1] / B) * NW;
            q3b = p3 + (size_t)(first_pos[3][j1] / B) * NW;
            q4b = p4 + (size_t)(first_pos[4][j1] / B) * NW;
        }
        __m256i acc0 = _mm256_setzero_si256();
        __m256i acc1 = _mm256_setzero_si256();
        int w = 0;
        for (; w + 4 <= NW; w += 4) {
            __m256i a0 = _mm256_loadu_si256((const __m256i*)(q0a + w));
            __m256i a1 = _mm256_loadu_si256((const __m256i*)(q1a + w));
            __m256i a2 = _mm256_loadu_si256((const __m256i*)(q2a + w));
            __m256i a3 = _mm256_loadu_si256((const __m256i*)(q3a + w));
            __m256i a4 = _mm256_loadu_si256((const __m256i*)(q4a + w));
            __m256i b0 = _mm256_loadu_si256((const __m256i*)(q0b + w));
            __m256i b1 = _mm256_loadu_si256((const __m256i*)(q1b + w));
            __m256i b2 = _mm256_loadu_si256((const __m256i*)(q2b + w));
            __m256i b3 = _mm256_loadu_si256((const __m256i*)(q3b + w));
            __m256i b4 = _mm256_loadu_si256((const __m256i*)(q4b + w));
            __m256i ra = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(a0,a1), _mm256_and_si256(a2,a3)), a4);
            __m256i rb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(b0,b1), _mm256_and_si256(b2,b3)), b4);
            acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(pc_bytes(ra), _mm256_setzero_si256()));
            acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(pc_bytes(rb), _mm256_setzero_si256()));
        }
        u32 c0 = sum4(acc0), c1 = sum4(acc1);
        for (; w < NW; w++) {
            c0 += (u32)__builtin_popcountll(q0a[w] & q1a[w] & q2a[w] & q3a[w] & q4a[w]);
            if (j1 >= 0) c1 += (u32)__builtin_popcountll(q0b[w] & q1b[w] & q2b[w] & q3b[w] & q4b[w]);
        }
        partial[j0] = c0;
        if (j1 >= 0) partial[j1] = c1;
    }

    // ---- pass 2: partial-block correction (16-wide SIMD, checks 4 other dims) ----
    for (int pos = 0; pos < n; pos++) {
        int j = proc[pos];
        u32 c = partial[j];
        stamp++;
        int f0 = first_pos[0][j], f1 = first_pos[1][j], f2 = first_pos[2][j], f3 = first_pos[3][j], f4 = first_pos[4][j];
        int t, mask, k;
        // dim 0: other 1,2,3,4
        {
            int* r1 = ordrank[0][1]; int* r2 = ordrank[0][2]; int* r3 = ordrank[0][3]; int* r4 = ordrank[0][4];
            __m256i vf1 = _mm256_set1_epi32(f1), vf2 = _mm256_set1_epi32(f2), vf3 = _mm256_set1_epi32(f3), vf4 = _mm256_set1_epi32(f4);
            for (t = (f0 / B) * B; t + 16 <= f0; t += 16) {
                __m256i r1a = _mm256_loadu_si256((const __m256i*)(r1+t));
                __m256i r1b = _mm256_loadu_si256((const __m256i*)(r1+t+8));
                __m256i r2a = _mm256_loadu_si256((const __m256i*)(r2+t));
                __m256i r2b = _mm256_loadu_si256((const __m256i*)(r2+t+8));
                __m256i r3a = _mm256_loadu_si256((const __m256i*)(r3+t));
                __m256i r3b = _mm256_loadu_si256((const __m256i*)(r3+t+8));
                __m256i r4a = _mm256_loadu_si256((const __m256i*)(r4+t));
                __m256i r4b = _mm256_loadu_si256((const __m256i*)(r4+t+8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1a), _mm256_cmpgt_epi32(vf2,r2a)), _mm256_cmpgt_epi32(vf3,r3a)), _mm256_cmpgt_epi32(vf4,r4a));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1b), _mm256_cmpgt_epi32(vf2,r2b)), _mm256_cmpgt_epi32(vf3,r3b)), _mm256_cmpgt_epi32(vf4,r4b));
                mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                while (mask) { k = __builtin_ctz(mask); int q = order[0][t+k]; if (seen[q] != stamp) { seen[q] = stamp; c++; } mask &= mask - 1; }
            }
            for (; t < f0; t++) if (r1[t] < f1 && r2[t] < f2 && r3[t] < f3 && r4[t] < f4) { int q = order[0][t]; if (seen[q] != stamp) { seen[q] = stamp; c++; } }
        }
        // dim 1: other 0,2,3,4
        {
            int* r1 = ordrank[1][0]; int* r2 = ordrank[1][2]; int* r3 = ordrank[1][3]; int* r4 = ordrank[1][4];
            __m256i vf1 = _mm256_set1_epi32(f0), vf2 = _mm256_set1_epi32(f2), vf3 = _mm256_set1_epi32(f3), vf4 = _mm256_set1_epi32(f4);
            for (t = (f1 / B) * B; t + 16 <= f1; t += 16) {
                __m256i r1a = _mm256_loadu_si256((const __m256i*)(r1+t));
                __m256i r1b = _mm256_loadu_si256((const __m256i*)(r1+t+8));
                __m256i r2a = _mm256_loadu_si256((const __m256i*)(r2+t));
                __m256i r2b = _mm256_loadu_si256((const __m256i*)(r2+t+8));
                __m256i r3a = _mm256_loadu_si256((const __m256i*)(r3+t));
                __m256i r3b = _mm256_loadu_si256((const __m256i*)(r3+t+8));
                __m256i r4a = _mm256_loadu_si256((const __m256i*)(r4+t));
                __m256i r4b = _mm256_loadu_si256((const __m256i*)(r4+t+8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1a), _mm256_cmpgt_epi32(vf2,r2a)), _mm256_cmpgt_epi32(vf3,r3a)), _mm256_cmpgt_epi32(vf4,r4a));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1b), _mm256_cmpgt_epi32(vf2,r2b)), _mm256_cmpgt_epi32(vf3,r3b)), _mm256_cmpgt_epi32(vf4,r4b));
                mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                while (mask) { k = __builtin_ctz(mask); int q = order[1][t+k]; if (seen[q] != stamp) { seen[q] = stamp; c++; } mask &= mask - 1; }
            }
            for (; t < f1; t++) if (r1[t] < f0 && r2[t] < f2 && r3[t] < f3 && r4[t] < f4) { int q = order[1][t]; if (seen[q] != stamp) { seen[q] = stamp; c++; } }
        }
        // dim 2: other 0,1,3,4
        {
            int* r1 = ordrank[2][0]; int* r2 = ordrank[2][1]; int* r3 = ordrank[2][3]; int* r4 = ordrank[2][4];
            __m256i vf1 = _mm256_set1_epi32(f0), vf2 = _mm256_set1_epi32(f1), vf3 = _mm256_set1_epi32(f3), vf4 = _mm256_set1_epi32(f4);
            for (t = (f2 / B) * B; t + 16 <= f2; t += 16) {
                __m256i r1a = _mm256_loadu_si256((const __m256i*)(r1+t));
                __m256i r1b = _mm256_loadu_si256((const __m256i*)(r1+t+8));
                __m256i r2a = _mm256_loadu_si256((const __m256i*)(r2+t));
                __m256i r2b = _mm256_loadu_si256((const __m256i*)(r2+t+8));
                __m256i r3a = _mm256_loadu_si256((const __m256i*)(r3+t));
                __m256i r3b = _mm256_loadu_si256((const __m256i*)(r3+t+8));
                __m256i r4a = _mm256_loadu_si256((const __m256i*)(r4+t));
                __m256i r4b = _mm256_loadu_si256((const __m256i*)(r4+t+8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1a), _mm256_cmpgt_epi32(vf2,r2a)), _mm256_cmpgt_epi32(vf3,r3a)), _mm256_cmpgt_epi32(vf4,r4a));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1b), _mm256_cmpgt_epi32(vf2,r2b)), _mm256_cmpgt_epi32(vf3,r3b)), _mm256_cmpgt_epi32(vf4,r4b));
                mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                while (mask) { k = __builtin_ctz(mask); int q = order[2][t+k]; if (seen[q] != stamp) { seen[q] = stamp; c++; } mask &= mask - 1; }
            }
            for (; t < f2; t++) if (r1[t] < f0 && r2[t] < f1 && r3[t] < f3 && r4[t] < f4) { int q = order[2][t]; if (seen[q] != stamp) { seen[q] = stamp; c++; } }
        }
        // dim 3: other 0,1,2,4
        {
            int* r1 = ordrank[3][0]; int* r2 = ordrank[3][1]; int* r3 = ordrank[3][2]; int* r4 = ordrank[3][4];
            __m256i vf1 = _mm256_set1_epi32(f0), vf2 = _mm256_set1_epi32(f1), vf3 = _mm256_set1_epi32(f2), vf4 = _mm256_set1_epi32(f4);
            for (t = (f3 / B) * B; t + 16 <= f3; t += 16) {
                __m256i r1a = _mm256_loadu_si256((const __m256i*)(r1+t));
                __m256i r1b = _mm256_loadu_si256((const __m256i*)(r1+t+8));
                __m256i r2a = _mm256_loadu_si256((const __m256i*)(r2+t));
                __m256i r2b = _mm256_loadu_si256((const __m256i*)(r2+t+8));
                __m256i r3a = _mm256_loadu_si256((const __m256i*)(r3+t));
                __m256i r3b = _mm256_loadu_si256((const __m256i*)(r3+t+8));
                __m256i r4a = _mm256_loadu_si256((const __m256i*)(r4+t));
                __m256i r4b = _mm256_loadu_si256((const __m256i*)(r4+t+8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1a), _mm256_cmpgt_epi32(vf2,r2a)), _mm256_cmpgt_epi32(vf3,r3a)), _mm256_cmpgt_epi32(vf4,r4a));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1b), _mm256_cmpgt_epi32(vf2,r2b)), _mm256_cmpgt_epi32(vf3,r3b)), _mm256_cmpgt_epi32(vf4,r4b));
                mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                while (mask) { k = __builtin_ctz(mask); int q = order[3][t+k]; if (seen[q] != stamp) { seen[q] = stamp; c++; } mask &= mask - 1; }
            }
            for (; t < f3; t++) if (r1[t] < f0 && r2[t] < f1 && r3[t] < f2 && r4[t] < f4) { int q = order[3][t]; if (seen[q] != stamp) { seen[q] = stamp; c++; } }
        }
        // dim 4: other 0,1,2,3
        {
            int* r1 = ordrank[4][0]; int* r2 = ordrank[4][1]; int* r3 = ordrank[4][2]; int* r4 = ordrank[4][3];
            __m256i vf1 = _mm256_set1_epi32(f0), vf2 = _mm256_set1_epi32(f1), vf3 = _mm256_set1_epi32(f2), vf4 = _mm256_set1_epi32(f3);
            for (t = (f4 / B) * B; t + 16 <= f4; t += 16) {
                __m256i r1a = _mm256_loadu_si256((const __m256i*)(r1+t));
                __m256i r1b = _mm256_loadu_si256((const __m256i*)(r1+t+8));
                __m256i r2a = _mm256_loadu_si256((const __m256i*)(r2+t));
                __m256i r2b = _mm256_loadu_si256((const __m256i*)(r2+t+8));
                __m256i r3a = _mm256_loadu_si256((const __m256i*)(r3+t));
                __m256i r3b = _mm256_loadu_si256((const __m256i*)(r3+t+8));
                __m256i r4a = _mm256_loadu_si256((const __m256i*)(r4+t));
                __m256i r4b = _mm256_loadu_si256((const __m256i*)(r4+t+8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1a), _mm256_cmpgt_epi32(vf2,r2a)), _mm256_cmpgt_epi32(vf3,r3a)), _mm256_cmpgt_epi32(vf4,r4a));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vf1,r1b), _mm256_cmpgt_epi32(vf2,r2b)), _mm256_cmpgt_epi32(vf3,r3b)), _mm256_cmpgt_epi32(vf4,r4b));
                mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                while (mask) { k = __builtin_ctz(mask); int q = order[4][t+k]; if (seen[q] != stamp) { seen[q] = stamp; c++; } mask &= mask - 1; }
            }
            for (; t < f4; t++) if (r1[t] < f0 && r2[t] < f1 && r3[t] < f2 && r4[t] < f3) { int q = order[4][t]; if (seen[q] != stamp) { seen[q] = stamp; c++; } }
        }
        out[j] = c;
    }
    munmap(pref, prefsz);
}

评测结果