代码

// KD-dimensional strict dominance counting, n <= 100000, coordinates in [0,n).
// For each dimension d, points are ordered by rank; PB[d][b] is the bitset of the
// first b*K points of that order (a subset of the strict-smaller set of any point
// whose rank lies beyond b*K).  The answer is |AND_d PB[d][rank_d(i)/K]| plus the
// points inside the K-sized partial ranges of each dimension (deduplicated), which
// is exact because every point with a smaller rank either is in all the prefixes or
// falls into one of those partial ranges.
#include <cstdlib>
#include <cstring>
typedef unsigned u32;
typedef unsigned long long u64;

#define KD 4
#define NMAX 100000
#define K 512
#define WMAX ((NMAX + 63) / 64)
#define NB (NMAX / K + 2)

static int ORD[KD][NMAX];
static u32 RK[KD][NMAX];
static int STAMP[NMAX];
static u64 *PB;                     // [KD][NB][WMAX]
static u64 acc[WMAX], tmp[WMAX];

#define PBat(d, b) (PB + (((size_t)(d) * NB + (b)) * WMAX))

static void buildAll(int n, const unsigned *x[KD]) {
    int *cnt = (int *)malloc((size_t)(n + 1) * sizeof(int));
    for (int d = 0; d < KD; d++) {
        memset(cnt, 0, (size_t)(n + 1) * sizeof(int));
        for (int i = 0; i < n; i++) cnt[x[d][i] + 1]++;
        for (int v = 1; v <= n; v++) cnt[v] += cnt[v - 1];
        for (int i = 0; i < n; i++) ORD[d][cnt[x[d][i]]++] = i;
        // rank of point j = number of points with a strictly smaller coordinate
        for (int p = 0; p < n; ) {
            int q = p;
            unsigned val = x[d][ORD[d][p]];
            while (q < n && x[d][ORD[d][q]] == val) q++;
            for (int t = p; t < q; t++) RK[d][ORD[d][t]] = (u32)p;
            p = q;
        }
    }
    free(cnt);
    // prefix bitsets: PB[d][b] = points at positions [0, b*K) of order d
    int words = (n + 63) >> 6;
    u64 *cur = (u64 *)calloc((size_t)words, 8);
    for (int d = 0; d < KD; d++) {
        memset(cur, 0, (size_t)words * 8);
        memset(PBat(d, 0), 0, (size_t)words * 8);
        for (int p = 0; p < n; p++) {
            if (p && (p % K) == 0) memcpy(PBat(d, p / K), cur, (size_t)words * 8);
            int j = ORD[d][p];
            cur[j >> 6] |= 1ULL << (j & 63);
        }
    }
    free(cur);
}

static void solve(int n, const unsigned *x[KD], unsigned *out) {
    static u64 *allocated = 0;
    if (!allocated) {
        allocated = (u64 *)calloc((size_t)KD * NB * WMAX, 8);
        PB = allocated;
    }
    buildAll(n, x);
    int words = (n + 63) >> 6;
    memset(STAMP, 0, (size_t)n * sizeof(int));
    int stampCur = 0;
    for (int i = 0; i < n; i++) {
        int b0 = (int)(RK[0][i] / K);
        memcpy(acc, PBat(0, b0), (size_t)words * 8);
        for (int d = 1; d < KD; d++) {
            u64 *pd = PBat(d, RK[d][i] / K);
            for (int w = 0; w < words; w++) acc[w] &= pd[w];
        }
        u32 ans = 0;
        for (int w = 0; w < words; w++) ans += (u32)__builtin_popcountll(acc[w]);
        stampCur++;
        for (int d = 0; d < KD; d++) {
            int lo = (int)(RK[d][i] / K) * K, hi = (int)RK[d][i];
            for (int p = lo; p < hi; p++) {
                int j = ORD[d][p];
                if (STAMP[j] == stampCur) continue;
                STAMP[j] = stampCur;
                int ok = 1;
                for (int e = 0; e < KD; e++) if (RK[e][j] >= RK[e][i]) { ok = 0; break; }
                if (ok) ans++;
            }
        }
        out[i] = ans;
    }
}

#if KD == 4
void count_4d(int n, const unsigned *x[4], unsigned *out) { solve(n, x, out); }
#else
void count_5d(int n, const unsigned *x[5], unsigned *out) { solve(n, x, out); }
#endif

评测结果