代码

#pragma GCC target("avx2")
// 4D strict dominance via bitset blocking + AVX2 (for n ~ 100k where n^2/64 fits).
#include <cstring>
#include <algorithm>
#include <cstdlib>
#include <immintrin.h>

typedef unsigned u32;
typedef unsigned long long u64;

#define B 64

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

void count_4d(int n, const unsigned *x[4], unsigned *out) {
    int NW = (n + 63) >> 6;
    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;

    for (int j = 0; j < n; j++) {
        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 + 4 <= NW; w += 4) {
            __m256i a = _mm256_loadu_si256((const __m256i*)(q0 + w));
            __m256i b = _mm256_loadu_si256((const __m256i*)(q1 + w));
            __m256i cc = _mm256_loadu_si256((const __m256i*)(q2 + w));
            __m256i d = _mm256_loadu_si256((const __m256i*)(q3 + w));
            __m256i r = _mm256_and_si256(_mm256_and_si256(a, b), _mm256_and_si256(cc, d));
            u64 r0 = _mm256_extract_epi64(r, 0);
            u64 r1 = _mm256_extract_epi64(r, 1);
            u64 r2 = _mm256_extract_epi64(r, 2);
            u64 r3 = _mm256_extract_epi64(r, 3);
            c += __builtin_popcountll(r0) + __builtin_popcountll(r1)
               + __builtin_popcountll(r2) + __builtin_popcountll(r3);
        }
        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);
}

评测结果