代码

// 4D strict dominance counting, n = 100000.
// Block-signature decomposition: K1-sized blocks per dimension; the empty-signature
// points come from a KD-dim prefix grid over block coordinates, one-signature points
// from ANDing coarse in-block prefix bitsets over the exact dim-d bit range, and
// >=2-signature points from small precomputed block-intersection cells.
#pragma GCC target("avx2")
#define KD 4
#define K1V 2048
#define NMAX 100000
#define QSORT 1
#define PF 8
#include <cstdlib>
#include <cstring>
#include <immintrin.h>

typedef unsigned u32;
typedef unsigned long long u64;

#ifndef KD
#define KD 4
#endif
#ifndef K1V
#define K1V 2048
#endif
#ifndef NMAX
#define NMAX 300000
#endif

static u32 ORD[KD][NMAX];
static u32 POS[KD][NMAX];
static u32 RK[KD][NMAX];
static u32 cntbuf[NMAX + 1];

static u32 *GRID;
static u64 *TB;
static u64 *TP[KD][KD];      // [m1][m1+1][W1]
static size_t strides_[KD];
static size_t TBwords_;
static u32 *COFF[KD][KD];    // [m1*m1+1]
static u32 *CDAT[KD][KD];    // n * CST
static int m1_, W1_, CST_;

static void buildAll(int n, const unsigned *x[KD]) {
    for (int d = 0; d < KD; d++) {
        memset(cntbuf, 0, (size_t)(n + 1) * sizeof(int));
        const unsigned *xd = x[d];
        for (int i = 0; i < n; i++) cntbuf[xd[i] + 1]++;
        for (int v = 1; v <= n; v++) cntbuf[v] += cntbuf[v - 1];
        for (int i = 0; i < n; i++) RK[d][i] = (u32)cntbuf[xd[i]];
        for (int i = 0; i < n; i++) ORD[d][cntbuf[xd[i]]++] = (u32)i;
        u32 *od = ORD[d], *pd = POS[d];
        for (int p = 0; p < n; p++) pd[od[p]] = (u32)p;
    }
    m1_ = (n + K1V - 1) / K1V;
    W1_ = K1V / 64;
    CST_ = ((KD + 3) / 4) * 4;

    // ---- grid over block coordinates ----
    size_t total = 1;
    for (int d = 0; d < KD; d++) { strides_[d] = total; total *= (size_t)m1_; }
    memset(GRID, 0, total * sizeof(u32));
    for (int j = 0; j < n; j++) {
        size_t idx = 0;
        for (int d = 0; d < KD; d++) idx += (size_t)(POS[d][j] / K1V) * strides_[d];
        GRID[idx]++;
    }
    for (int ax = KD - 1; ax >= 0; ax--) {
        size_t st = strides_[ax], blk = st * (size_t)m1_;
        for (size_t base = 0; base < total; base += blk)
            for (size_t k = 1; k < (size_t)m1_; k++) {
                u32 *dst = GRID + base + k * st, *src = GRID + base + (k - 1) * st;
                for (size_t r = 0; r < st; r++) dst[r] += src[r];
            }
    }

    // ---- coarse prefix bitsets inside each block ----
    // TP[d][e] : [block of dim d][entry b][W1] , bit q = point at ORD[d][block*K1+q]
    memset(TB, 0, TBwords_ * 8);
    u64 *run = (u64 *)calloc((size_t)m1_ * W1_, 8);
    size_t stride = (size_t)(m1_ + 1) * W1_;
    for (int d = 0; d < KD; d++)
        for (int e = 0; e < KD; e++) {
            if (e == d) continue;
            memset(run, 0, (size_t)m1_ * W1_ * 8);
            u64 *tbl = TP[d][e];
            for (int p = 0; p <= n; p++) {
                if (p % K1V == 0) {
                    int b = p / K1V;
                    if (b <= m1_)
                        for (int blk = 0; blk < m1_; blk++)
                            memcpy(tbl + (size_t)blk * stride + (size_t)b * W1_,
                                   run + (size_t)blk * W1_, (size_t)W1_ * 8);
                }
                if (p == n) break;
                u32 j = ORD[e][p];
                int bp = (int)(POS[d][j] / K1V);
                int q = (int)(POS[d][j] - (u32)(bp * K1V));
                run[(size_t)bp * W1_ + (q >> 6)] |= 1ULL << (q & 63);
            }
        }
    free(run);

    // ---- pair cells: intersection of two K1-blocks (sorted by pos of e1) ----
    u32 *hist = (u32 *)malloc(((size_t)m1_ * m1_ + 1) * 4);
    for (int e1 = 0; e1 < KD; e1++)
        for (int e2 = e1 + 1; e2 < KD; e2++) {
            memset(hist, 0, (size_t)m1_ * m1_ * 4);
            for (int p = 0; p < n; p++) {
                u32 j = ORD[e1][p];
                hist[(size_t)(p / K1V) * m1_ + (POS[e2][j] / K1V)]++;
            }
            u32 *co = COFF[e1][e2];
            u32 s = 0;
            for (size_t k = 0; k < (size_t)m1_ * m1_; k++) { co[k] = s; s += hist[k]; }
            co[(size_t)m1_ * m1_] = s;
            for (size_t k = 0; k < (size_t)m1_ * m1_; k++) hist[k] = co[k];
            u32 *cd = CDAT[e1][e2];
            for (int p = 0; p < n; p++) {
                u32 j = ORD[e1][p];
                u32 *dst = cd + (size_t)hist[(size_t)(p / K1V) * m1_ + (POS[e2][j] / K1V)]++ * CST_;
                for (int e = 0; e < KD; e++) dst[e] = POS[e][j];
                for (int k = KD; k < CST_; k++) dst[k] = 0;
            }
        }
    free(hist);
}

#if KD == 4
void count_4d(int n, const unsigned *x[4], unsigned *out)
#else
void count_5d(int n, const unsigned *x[5], unsigned *out)
#endif
{
    static int inited = 0;
    if (!inited) {
        inited = 1;
        int m1 = (NMAX + K1V - 1) / K1V;
        W1_ = K1V / 64;
        CST_ = ((KD + 3) / 4) * 4;
        size_t total = 1;
        for (int d = 0; d < KD; d++) total *= (size_t)m1;
        if (posix_memalign((void **)&GRID, 64, total * sizeof(u32))) return;
        size_t np = 0;
        for (int d = 0; d < KD; d++) np += KD - 1;
        TBwords_ = (size_t)np * m1 * (size_t)(m1 + 1) * W1_;
        if (posix_memalign((void **)&TB, 64, TBwords_ * 8)) return;
        size_t off = 0;
        for (int d = 0; d < KD; d++)
            for (int e = 0; e < KD; e++) {
                if (e == d) continue;
                TP[d][e] = TB + off;
                off += (size_t)m1 * (size_t)(m1 + 1) * W1_;
            }
        for (int e1 = 0; e1 < KD; e1++) {
            for (int e2 = e1 + 1; e2 < KD; e2++) {
                if (posix_memalign((void **)&COFF[e1][e2], 64, ((size_t)m1 * m1 + 1) * 4)) return;
                if (posix_memalign((void **)&CDAT[e1][e2], 64, (size_t)NMAX * CST_ * 4)) return;
            }
        }
    }
    buildAll(n, x);

#ifdef QSORT
    static u32 key_[NMAX], ord_[NMAX], tmp_[NMAX];
    for (int i = 0; i < n; i++) {
        u32 k = 0;
        for (int e = 0; e < KD; e++) k = k * (u32)m1_ + (u32)(RK[e][i] / K1V);
        key_[i] = k; ord_[i] = (u32)i;
    }
    {
        static u32 cnt[256];
        for (int sh = 0; sh < 32; sh += 8) {
            memset(cnt, 0, sizeof(cnt));
            for (int i = 0; i < n; i++) cnt[(key_[ord_[i]] >> sh) & 255]++;
            u32 s = 0;
            for (int b = 0; b < 256; b++) { u32 c = cnt[b]; cnt[b] = s; s += c; }
            for (int i = 0; i < n; i++) tmp_[cnt[(key_[ord_[i]] >> sh) & 255]++] = ord_[i];
            memcpy(ord_, tmp_, (size_t)n * 4);
        }
    }
#endif

    for (int qi = 0; qi < n; qi++) {
#ifdef QSORT
        int i = (int)ord_[qi];
#else
        int i = qi;
#endif
#ifdef PF
        {
            int ip = (PF < 0) ? 0 : (qi + (PF) < n ? (int)ord_[qi + (PF)] : (int)ord_[n - 1]);
            if (ip >= 0) {
                for (int e1 = 0; e1 < KD; e1++)
                    for (int e2 = e1 + 1; e2 < KD; e2++) {
                        size_t cidx = (size_t)(RK[e1][ip] / K1V) * m1_ + (RK[e2][ip] / K1V);
                        const u32 *cp = CDAT[e1][e2] + (size_t)COFF[e1][e2][cidx] * CST_;
                        for (int z = 0; z < 8; z++)
                            _mm_prefetch((const char *)cp + 64 * z, _MM_HINT_T0);
                    }
            }
        }
#endif
        u32 t[KD], B[KD];
        int b[KD];
        for (int e = 0; e < KD; e++) {
            t[e] = RK[e][i];
            b[e] = (int)(t[e] / K1V);
            B[e] = (u32)b[e] * K1V;
        }
        u64 ans = 0;
        {
            int ok = 1;
            size_t idx = 0;
            for (int e = 0; e < KD; e++) {
                if (b[e] == 0) { ok = 0; break; }
                idx += (size_t)(b[e] - 1) * strides_[e];
            }
            if (ok) ans = GRID[idx];
        }
        // ---- |S| = 1 terms ----
        for (int d = 0; d < KD; d++) {
            int r = (int)(t[d] - B[d]);
            if (r == 0) continue;
            int beta = b[d];
            int fullw = r >> 6, rem = r & 63;
            const u64 *p0, *p1, *p2, *p3 = 0;
            {
                int k = 0; const u64 *tmp[KD];
                for (int e = 0; e < KD; e++) {
                    if (e == d) continue;
                    tmp[k++] = TP[d][e] + ((size_t)beta * (m1_ + 1) + b[e]) * W1_;
                }
                p0 = tmp[0];
#if KD >= 3
                p1 = tmp[1];
#else
                p1 = 0;
#endif
#if KD >= 4
                p2 = tmp[2];
#else
                p2 = 0;
#endif
#if KD >= 5
                p3 = tmp[3];
#endif
            }
            __m256i acc = _mm256_setzero_si256();
            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);
            int w = 0;
            for (; w + 4 <= fullw; w += 4) {
                __m256i a = _mm256_loadu_si256((const __m256i *)(p0 + w));
                a = _mm256_and_si256(a, _mm256_loadu_si256((const __m256i *)(p1 + w)));
#if KD >= 4
                a = _mm256_and_si256(a, _mm256_loadu_si256((const __m256i *)(p2 + w)));
#endif
#if KD >= 5
                a = _mm256_and_si256(a, _mm256_loadu_si256((const __m256i *)(p3 + w)));
#endif
                __m256i lo = _mm256_and_si256(a, m4);
                __m256i hi = _mm256_and_si256(_mm256_srli_epi16(a, 4), m4);
                __m256i c = _mm256_add_epi8(_mm256_shuffle_epi8(lut, lo), _mm256_shuffle_epi8(lut, hi));
                acc = _mm256_add_epi64(acc, _mm256_sad_epu8(c, _mm256_setzero_si256()));
            }
            u64 part[4];
            _mm256_storeu_si256((__m256i *)part, acc);
            u32 cnt = (u32)(part[0] + part[1] + part[2] + part[3]);
            for (; w < fullw; w++) {
                u64 a = p0[w] & p1[w];
#if KD >= 4
                a &= p2[w];
#endif
#if KD >= 5
                a &= p3[w];
#endif
                cnt += (u32)__builtin_popcountll(a);
            }
            if (rem) {
                u64 a = p0[fullw] & p1[fullw];
#if KD >= 4
                a &= p2[fullw];
#endif
#if KD >= 5
                a &= p3[fullw];
#endif
                cnt += (u32)__builtin_popcountll(a & ((1ULL << rem) - 1));
            }
            ans += cnt;
        }
        // ---- |S| >= 2 terms ----
        for (int e1 = 0; e1 < KD; e1++) {
            for (int e2 = e1 + 1; e2 < KD; e2++) {
                size_t cidx = (size_t)b[e1] * m1_ + b[e2];
                u32 c0 = COFF[e1][e2][cidx], c1 = COFF[e1][e2][cidx + 1];
                if (c0 == c1) continue;
                const u32 *q = CDAT[e1][e2] + (size_t)c0 * CST_;
                int cnt = (int)(c1 - c0);
                u32 TT[8] __attribute__((aligned(32)));
                for (int e = 0; e < CST_; e++) {
                    if (e >= KD) { TT[e] = 0x7fffffffu; continue; }
                    u32 v = t[e];
                    if (e < e2 && e != e1 && B[e] < v) v = B[e];
                    TT[e] = v;
                }
                for (int k = 4; k < CST_; k++) TT[k] = TT[k & 3];
                // cell data is sorted by pos[e1]; only the prefix with pos[e1] < TT[e1] can pass
                int lo = 0, hi = cnt;
                while (lo < hi) {
                    int md = (lo + hi) >> 1;
                    if (q[(size_t)md * CST_ + e1] < TT[e1]) lo = md + 1; else hi = md;
                }
                cnt = lo;
                if (CST_ == 4) {
                    __m128i vt = _mm_load_si128((const __m128i *)TT);
                    int c = 0;
                    for (int k = 0; k < cnt; k++, q += 4) {
                        __m128i rr = _mm_loadu_si128((const __m128i *)q);
                        __m128i m1 = _mm_cmpgt_epi32(vt, rr);
                        c += _mm_movemask_ps(_mm_castsi128_ps(m1)) == 0x0f;
                    }
                    ans += (u32)c;
                } else {
                    __m256i vt = _mm256_load_si256((const __m256i *)TT);
                    int c = 0;
                    for (int k = 0; k < cnt; k++, q += CST_) {
                        __m256i rr = _mm256_loadu_si256((const __m256i *)q);
                        __m256i m1 = _mm256_cmpgt_epi32(vt, rr);
                        c += _mm256_movemask_ps(_mm256_castsi256_ps(m1)) == 0xff;
                    }
                    ans += (u32)c;
                }
            }
        }
        out[i] = (unsigned)ans;
    }
}

评测结果