代码

// 4D dominance n=1e5: prefix-bitset B=512 + 3D blocking + 4-point AND + 4-pass exact brute force.
// AND loop caches dims 0,1,2 (L3), streams dim3 (NT). Partial-block correction is a 4-pass
// "responsibility" scan (q counted in smallest dim d where q is in j's partial block): no dedup,
// sequential ordrank/ordpos reads. e<d uses ordpos (position), e>d uses ordrank (first_pos).
#pragma GCC target("avx2,bmi")
#pragma GCC optimize("O3")
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <sys/mman.h>
#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 ordrank[4][4][100010];
static int ordpos[4][4][100010];
static int pos[4][100010];
static u32 partial[100010];
static u32 corr[100010];
static int proc[100010];
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_4d(int n, const unsigned *x[4], unsigned *out) {
    int NW = (n + 63) >> 6;
    const int B = 512;
    int M = (n + B - 1) / B;
    int NWp = (NW + 7) & ~7;
    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;
        for (int t = 0; t < n; t++) pos[d][order[d][t]] = t;
    }
    for (int d = 0; d < 4; d++)
        for (int a = 0; a < 4; a++)
            for (int t = 0; t < n; t++) {
                ordrank[d][a][t] = first_pos[a][order[d][t]];
                ordpos[d][a][t] = pos[a][order[d][t]];
            }
    size_t prefsz = (size_t)4 * (M + 1) * NWp * 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){
                int a1 = first_pos[1][a], b1 = first_pos[1][b];
                if (a1 != b1) return a1 < b1;
                return first_pos[2][a] < first_pos[2][b];
            });
        }
    }    void *rawp = malloc(prefsz + 64);
    u64 *pref = (u64*)(((size_t)rawp + 63) & ~(size_t)63);

    for (int d = 0; d < 4; d++) {
        u64 *p = pref + (size_t)d * (M + 1) * NWp;
        memset(p, 0, (size_t)NWp * sizeof(u64));
        for (int b = 0; b < M; b++) {
            u64 *cur = p + (size_t)(b + 1) * NWp;
            u64 *prev = p + (size_t)b * NWp;
            memcpy(cur, prev, (size_t)NWp * 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) * NWp;
    u64 *p2 = pref + (size_t)2 * (M + 1) * NWp;
    u64 *p3 = pref + (size_t)3 * (M + 1) * NWp;
    for (int pos = 0; pos < n; pos += 4) {
        int j[4];
        const u64 *q[4][4];
        for (int k = 0; k < 4; k++) {
            j[k] = (pos + k < n) ? proc[pos + k] : -1;
            if (j[k] >= 0) {
                q[k][0] = p0 + (size_t)(first_pos[0][j[k]] / B) * NWp;
                q[k][1] = p1 + (size_t)(first_pos[1][j[k]] / B) * NWp;
                q[k][2] = p2 + (size_t)(first_pos[2][j[k]] / B) * NWp;
                q[k][3] = p3 + (size_t)(first_pos[3][j[k]] / B) * NWp;
            } else {
                q[k][0] = q[0][0]; q[k][1] = q[0][1]; q[k][2] = q[0][2]; q[k][3] = q[0][3];
            }
        }
        __m256i acc0 = _mm256_setzero_si256();
        __m256i acc1 = _mm256_setzero_si256();
        __m256i acc2 = _mm256_setzero_si256();
        __m256i acc3 = _mm256_setzero_si256();
        int w = 0;
        for (; w + 4 <= NW; w += 4) {
            __m256i r0, r1, r2, r3;
            {
                __m256i a0 = _mm256_load_si256((const __m256i*)(q[0][0] + w));
                __m256i a1 = _mm256_load_si256((const __m256i*)(q[0][1] + w));
                __m256i a2 = _mm256_load_si256((const __m256i*)(q[0][2] + w));
                __m256i a3 = _mm256_stream_load_si256((__m256i*)(q[0][3] + w));
                r0 = _mm256_and_si256(_mm256_and_si256(a0,a1), _mm256_and_si256(a2,a3));
            }
            {
                __m256i a0 = _mm256_load_si256((const __m256i*)(q[1][0] + w));
                __m256i a1 = _mm256_load_si256((const __m256i*)(q[1][1] + w));
                __m256i a2 = _mm256_load_si256((const __m256i*)(q[1][2] + w));
                __m256i a3 = _mm256_stream_load_si256((__m256i*)(q[1][3] + w));
                r1 = _mm256_and_si256(_mm256_and_si256(a0,a1), _mm256_and_si256(a2,a3));
            }
            {
                __m256i a0 = _mm256_load_si256((const __m256i*)(q[2][0] + w));
                __m256i a1 = _mm256_load_si256((const __m256i*)(q[2][1] + w));
                __m256i a2 = _mm256_load_si256((const __m256i*)(q[2][2] + w));
                __m256i a3 = _mm256_stream_load_si256((__m256i*)(q[2][3] + w));
                r2 = _mm256_and_si256(_mm256_and_si256(a0,a1), _mm256_and_si256(a2,a3));
            }
            {
                __m256i a0 = _mm256_load_si256((const __m256i*)(q[3][0] + w));
                __m256i a1 = _mm256_load_si256((const __m256i*)(q[3][1] + w));
                __m256i a2 = _mm256_load_si256((const __m256i*)(q[3][2] + w));
                __m256i a3 = _mm256_stream_load_si256((__m256i*)(q[3][3] + w));
                r3 = _mm256_and_si256(_mm256_and_si256(a0,a1), _mm256_and_si256(a2,a3));
            }
            acc0 = _mm256_add_epi64(acc0, _mm256_sad_epu8(pc_bytes(r0), _mm256_setzero_si256()));
            acc1 = _mm256_add_epi64(acc1, _mm256_sad_epu8(pc_bytes(r1), _mm256_setzero_si256()));
            acc2 = _mm256_add_epi64(acc2, _mm256_sad_epu8(pc_bytes(r2), _mm256_setzero_si256()));
            acc3 = _mm256_add_epi64(acc3, _mm256_sad_epu8(pc_bytes(r3), _mm256_setzero_si256()));
        }
        u32 c0 = sum4(acc0), c1 = sum4(acc1), c2 = sum4(acc2), c3 = sum4(acc3);
        for (; w < NW; w++) {
            if (j[0] >= 0) c0 += (u32)__builtin_popcountll(q[0][0][w] & q[0][1][w] & q[0][2][w] & q[0][3][w]);
            if (j[1] >= 0) c1 += (u32)__builtin_popcountll(q[1][0][w] & q[1][1][w] & q[1][2][w] & q[1][3][w]);
            if (j[2] >= 0) c2 += (u32)__builtin_popcountll(q[2][0][w] & q[2][1][w] & q[2][2][w] & q[2][3][w]);
            if (j[3] >= 0) c3 += (u32)__builtin_popcountll(q[3][0][w] & q[3][1][w] & q[3][2][w] & q[3][3][w]);
        }
        if (j[0] >= 0) partial[j[0]] = c0;
        if (j[1] >= 0) partial[j[1]] = c1;
        if (j[2] >= 0) partial[j[2]] = c2;
        if (j[3] >= 0) partial[j[3]] = c3;
    }
    // 4-pass responsibility brute force
    memset(corr, 0, sizeof(u32) * n);
    for (int d = 0; d < 4; d++) {
        int e0 = (d == 0) ? 1 : 0;
        int e1 = (d <= 1) ? 2 : 1;
        int e2 = 3 - d;
        if (d == 0) { e0 = 1; e1 = 2; e2 = 3; }
        else if (d == 1) { e0 = 0; e1 = 2; e2 = 3; }
        else if (d == 2) { e0 = 0; e1 = 1; e2 = 3; }
        else { e0 = 0; e1 = 1; e2 = 2; }
        // select candidate arrays: e<d -> ordpos, e>d -> ordrank
        int* A = (e0 < d) ? ordpos[d][e0] : ordrank[d][e0];
        int* B_ = (e1 < d) ? ordpos[d][e1] : ordrank[d][e1];
        int* C = (e2 < d) ? ordpos[d][e2] : ordrank[d][e2];
        int* Dsel = ordrank[d][d];
        for (int t = 0; t < n; t++) {
            int j = order[d][t];
            int f = Dsel[t];
            int bs = (f / B) * B;
            int ta = (e0 < d) ? ((ordrank[d][e0][t] / B) * B) : ordrank[d][e0][t];
            int tb = (e1 < d) ? ((ordrank[d][e1][t] / B) * B) : ordrank[d][e1][t];
            int tc = (e2 < d) ? ((ordrank[d][e2][t] / B) * B) : ordrank[d][e2][t];
            __m256i vta = _mm256_set1_epi32(ta);
            __m256i vtb = _mm256_set1_epi32(tb);
            __m256i vtc = _mm256_set1_epi32(tc);
            int p = bs;
            int cntm = 0;
            for (; p + 16 <= f; p += 16) {
                __m256i a1 = _mm256_loadu_si256((const __m256i*)(A + p));
                __m256i a2 = _mm256_loadu_si256((const __m256i*)(A + p + 8));
                __m256i b1 = _mm256_loadu_si256((const __m256i*)(B_ + p));
                __m256i b2 = _mm256_loadu_si256((const __m256i*)(B_ + p + 8));
                __m256i c1 = _mm256_loadu_si256((const __m256i*)(C + p));
                __m256i c2 = _mm256_loadu_si256((const __m256i*)(C + p + 8));
                __m256i ma = _mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vta, a1), _mm256_cmpgt_epi32(vtb, b1)), _mm256_cmpgt_epi32(vtc, c1));
                __m256i mb = _mm256_and_si256(_mm256_and_si256(_mm256_cmpgt_epi32(vta, a2), _mm256_cmpgt_epi32(vtb, b2)), _mm256_cmpgt_epi32(vtc, c2));
                int mask = _mm256_movemask_ps(_mm256_castsi256_ps(ma)) | (_mm256_movemask_ps(_mm256_castsi256_ps(mb)) << 8);
                cntm += __builtin_popcount(mask);
            }
            for (; p < f; p++) {
                if (A[p] < ta && B_[p] < tb && C[p] < tc) cntm++;
            }
            corr[j] += cntm;
        }
    }
    for (int i = 0; i < n; i++) out[i] = partial[i] + corr[i];
    free(rawp);
}

评测结果