代码

// 1011: 5D strict dominance counting, n = 300000 (see sol/1011a for n = 100000).
//
// Idea (mine): answer(i) = |intersection_d { j : pos_d(j) < t_d(i) }| where
// pos_d is the position of a point in the order sorted by dim d and t_d(i) is the
// number of points strictly smaller than i in dim d.
//
// Let a = argmin_d t_d(i) ("leading dim").  Every dominated point has pos_a < t_a,
// so the intersection can be materialised as a bitset over *dim-a positions* and
// only the first ceil(t_a/64) words ever have to be touched.  Of the remaining four
// dimensions one is handled by a sweep (points activated in that dim's order, which
// gives an exact prefix for free) and the other three by prefix bitsets stored at
// block granularity B (block b covers positions [b*B,(b+1)*B)).  The block prefix is
// a subset of the true prefix, so the points of the residual range [floor(t/B)*B, t)
// that satisfy the other four conditions exactly are added back by scanning that
// range (<= B points) against a precomputed per-(a,d) rank record.
#pragma GCC target("avx2,popcnt")
#include <cstdlib>
#include <cstring>
#include <immintrin.h>

typedef unsigned u32;
typedef unsigned long long u64;

#ifndef NMAX
#define NMAX 300000
#endif
#ifndef B
#define B 128
#endif
#ifndef KD
#define KD 5
#endif
#ifndef PF
#define PF 0
#endif
#ifndef RPF
#define RPF 4
#endif

static u32 SD[KD][NMAX];     // SD[d][p] = point id at sorted position p of dim d
static u32 RP[KD][NMAX];     // RP[d][i] = sorted position of point i in dim d
static u32 LT[KD][NMAX];     // LT[d][i] = #{ j : x_d(j) < x_d(i) }
static u32 OWN[NMAX];        // leading dim of each query
static u32 QL[KD][NMAX];     // queries owned by a, sorted by LT[SW[a]][.]
static u32 QC[KD];           // number of queries owned by a
static u32 MARKV[NMAX];
static u32 cntx[2 * NMAX + 2];

static u64 *TT[KD][3];       // prefix bitsets: [a][k] -> (entries) x W words
static u64 *ACC;             // build accumulator (W words)
static u32 *CX[KD][3];       // rank records: [a][k] -> n * 4 u32
static int W_, NBLK_;
static int TD[KD][3];
static int SW[KD];

// Entry E of a prefix table (bitset of {pos_d < E*B} over dim-a positions) is only
// ever read in its first (E+1)*B/64 words: a query using entry E has leading threshold
// t_a <= t_d < (E+1)*B.  Each entry therefore stores just that prefix (triangular
// layout, 64-byte aligned), halving both the table memory and the build traffic.
#define LW (B / 64)
static size_t EOFFS[NMAX / B + 4];

// ------------------------------------------------------------------ helpers
static inline u64 popc64(u64 v) { return (u64)__builtin_popcountll(v); }

static inline __m256i vpc(__m256i v, __m256i lut, __m256i m4) {
    __m256i lo = _mm256_and_si256(v, m4);
    __m256i hi = _mm256_and_si256(_mm256_srli_epi16(v, 4), m4);
    __m256i c = _mm256_add_epi8(_mm256_shuffle_epi8(lut, lo), _mm256_shuffle_epi8(lut, hi));
    return _mm256_sad_epu8(c, _mm256_setzero_si256());
}

#define CSA(h, l, a, b, c)                                    \
    {                                                         \
        __m256i u_ = _mm256_xor_si256(a, b);                   \
        h = _mm256_or_si256(_mm256_and_si256(a, b),            \
                            _mm256_and_si256(u_, c));         \
        l = _mm256_xor_si256(u_, c);                          \
    }

// popcount of the 4-way AND of p0..p3 over `words` words (Harley-Seal CSA tree)
__attribute__((target("avx2")))
static u64 and4pc(const u64 *p0, const u64 *p1, const u64 *p2, const u64 *p3, u32 words) {
    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);
    __m256i tot = _mm256_setzero_si256();     // accumulates popcounts of the "16" planes
    __m256i extra = _mm256_setzero_si256();   // accumulates plain popcounts of the tail
    __m256i ones = tot, twos = tot, fours = tot, eights = tot;
    u32 i = 0;
#define ANDV(off)                                                                        \
    _mm256_and_si256(                                                                    \
        _mm256_loadu_si256((const __m256i *)(p0 + i + (off))),                           \
        _mm256_and_si256(                                                                \
            _mm256_loadu_si256((const __m256i *)(p1 + i + (off))),                       \
            _mm256_and_si256(_mm256_loadu_si256((const __m256i *)(p2 + i + (off))),      \
                             _mm256_loadu_si256((const __m256i *)(p3 + i + (off))))))
    for (; i + 64 <= words; i += 64) {
        __m256i tA, tB, fA, fB, eA, eB;
        CSA(tA, ones, ones, ANDV(0), ANDV(4));
        CSA(tB, ones, ones, ANDV(8), ANDV(12));
        CSA(fA, twos, twos, tA, tB);
        CSA(tA, ones, ones, ANDV(16), ANDV(20));
        CSA(tB, ones, ones, ANDV(24), ANDV(28));
        CSA(fB, twos, twos, tA, tB);
        CSA(eA, fours, fours, fA, fB);
        CSA(tA, ones, ones, ANDV(32), ANDV(36));
        CSA(tB, ones, ones, ANDV(40), ANDV(44));
        CSA(fA, twos, twos, tA, tB);
        CSA(tA, ones, ones, ANDV(48), ANDV(52));
        CSA(tB, ones, ones, ANDV(56), ANDV(60));
        CSA(fB, twos, twos, tA, tB);
        CSA(eB, fours, fours, fA, fB);
        {
            __m256i sixt;
            CSA(sixt, eights, eights, eA, eB);
            tot = _mm256_add_epi64(tot, vpc(sixt, lut, m4));
        }
    }
    for (; i + 4 <= words; i += 4) extra = _mm256_add_epi64(extra, vpc(ANDV(0), lut, m4));
#undef ANDV
    tot = _mm256_slli_epi64(tot, 4);
    tot = _mm256_add_epi64(tot, _mm256_slli_epi64(vpc(eights, lut, m4), 3));
    tot = _mm256_add_epi64(tot, _mm256_slli_epi64(vpc(fours, lut, m4), 2));
    tot = _mm256_add_epi64(tot, _mm256_slli_epi64(vpc(twos, lut, m4), 1));
    tot = _mm256_add_epi64(tot, vpc(ones, lut, m4));
    tot = _mm256_add_epi64(tot, extra);
    u64 t[4];
    _mm256_storeu_si256((__m256i *)t, tot);
    u64 s = t[0] + t[1] + t[2] + t[3];
    for (; i < words; i++) s += popc64(p0[i] & p1[i] & p2[i] & p3[i]);
    return s;
}

// Residual range [lo, lo+cntp) of one table dim: count the points whose four
// recorded ranks are all below the four thresholds, each point at most once per query.
__attribute__((target("avx2,popcnt"), always_inline))
static inline u64 residual(const u32 *rec, const u32 *sd, u32 cntp,
                           u32 t0, u32 t1, u32 t2, u32 t3, u32 token, u32 *mark) {
    u64 res = 0;
    __m256i vt = _mm256_setr_epi32((int)t0, (int)t1, (int)t2, (int)t3,
                                   (int)t0, (int)t1, (int)t2, (int)t3);
    u32 q = 0;
    for (; q + 2 <= cntp; q += 2) {
        __m256i r = _mm256_loadu_si256((const __m256i *)(rec + (size_t)q * 4));
        u32 msk = (u32)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(vt, r)));
        if ((msk & 0x0F) == 0x0F) {
            u32 j = sd[q];
            if (mark[j] != token) { mark[j] = token; res++; }
        }
        if ((msk & 0xF0) == 0xF0) {
            u32 j = sd[q + 1];
            if (mark[j] != token) { mark[j] = token; res++; }
        }
    }
    for (; q < cntp; q++) {
        const u32 *r = rec + (size_t)q * 4;
        if (r[0] < t0 && r[1] < t1 && r[2] < t2 && r[3] < t3) {
            u32 j = sd[q];
            if (mark[j] != token) { mark[j] = token; res++; }
        }
    }
    return res;
}

// ------------------------------------------------------------------ build
static void PREP(int n, const unsigned *x[KD]) {
    W_ = (n + 63) >> 6;
    NBLK_ = (n - 1) / B + 1;

    for (int d = 0; d < KD; d++) {
        int maxv = 0;
        const unsigned *xd = x[d];
        for (int i = 0; i < n; i++) if ((int)xd[i] > maxv) maxv = (int)xd[i];
        memset(cntx, 0, (size_t)(maxv + 2) * 4);
        for (int i = 0; i < n; i++) cntx[xd[i] + 1]++;
        for (int v = 1; v <= maxv + 1; v++) cntx[v] += cntx[v - 1];
        for (int i = 0; i < n; i++) LT[d][i] = cntx[xd[i]];
        for (int i = 0; i < n; i++) {
            u32 v = xd[i];
            u32 p = cntx[v]++;
            RP[d][i] = p;
            SD[d][p] = (u32)i;
        }
    }

    // leading dim = smallest t, sweep dim = a+1, table dims = the rest
    for (int a = 0; a < KD; a++) {
        SW[a] = (a + 1) % KD;
        int k = 0;
        for (int d = 0; d < KD; d++) if (d != a && d != SW[a]) TD[a][k++] = d;
    }
    for (int i = 0; i < n; i++) {
        int best = 0;
        for (int d = 1; d < KD; d++) if (LT[d][i] < LT[best][i]) best = d;
        OWN[i] = (u32)best;
    }
    // per-owner query lists ordered by the sweep-dim threshold
    {
        static u32 cnt2[NMAX + 1];
        for (int a = 0; a < KD; a++) {
            memset(cnt2, 0, (size_t)(n + 1) * 4);
            int s = SW[a];
            for (int i = 0; i < n; i++) if ((int)OWN[i] == a) cnt2[LT[s][i]]++;
            u32 sum = 0;
            for (int v = 0; v <= n; v++) { u32 c = cnt2[v]; cnt2[v] = sum; sum += c; }
            for (int i = 0; i < n; i++) if ((int)OWN[i] == a) QL[a][cnt2[LT[s][i]]++] = (u32)i;
            QC[a] = sum;
        }
    }

    // tables + rank records
    for (int a = 0; a < KD; a++) {
        int s = SW[a];
        for (int k = 0; k < 3; k++) {
            int d = TD[a][k];
            int y = TD[a][(k + 1) % 3], z = TD[a][(k + 2) % 3];
            // rank record for every position p of dim d
            u32 *rec = CX[a][k];
            for (int p = 0; p < n; p++) {
                u32 j = SD[d][p];
                rec[4 * p + 0] = RP[a][j];
                rec[4 * p + 1] = RP[s][j];
                rec[4 * p + 2] = RP[y][j];
                rec[4 * p + 3] = RP[z][j];
            }
            // prefix bitsets over dim-a positions; the running accumulator stays in
            // L2 so the build only streams the table out once (no read-back).
            u64 *T = TT[a][k];
            u64 *acc = ACC;
            memset(acc, 0, (size_t)W_ * 8);
            memset(T, 0, (size_t)EOFFS[1] * 8);
            for (int b = 1; b <= NBLK_; b++) {
                int lo = (b - 1) * B, hi = b * B;
                if (hi > n) hi = n;
                for (int p = lo; p < hi; p++) {
                    u32 q = RP[a][SD[d][p]];
                    acc[q >> 6] |= 1ULL << (q & 63);
                }
                memcpy(T + EOFFS[b], acc, (size_t)(b + 1) * LW * 8);
            }
        }
    }
}

// ------------------------------------------------------------------ query
// For each table dim the block boundary closest to the exact threshold is used
// (floor -> the residual range is added, ceil -> it is subtracted); the scan is
// then at most B/2 points instead of B.
__attribute__((target("avx2,popcnt")))
static void RUN(int n, unsigned *out) {
    static u64 active[(NMAX + 63) / 64];
    u32 token = 0;
    memset(MARKV, 0, (size_t)n * 4);
    for (int a = 0; a < KD; a++) {
        int s = SW[a];
        int d0 = TD[a][0], d1 = TD[a][1], d2 = TD[a][2];
        const u64 *T0 = TT[a][0], *T1 = TT[a][1], *T2 = TT[a][2];
        const u32 *C0 = CX[a][0], *C1 = CX[a][1], *C2 = CX[a][2];
        const u32 *sd0 = SD[d0], *sd1 = SD[d1], *sd2 = SD[d2];
        memset(active, 0, (size_t)W_ * 8);
        u32 ptr = 0;
        const u32 *ql = QL[a];
        u32 qn = QC[a];
        for (u32 q = 0; q < qn; q++) {
            u32 i = ql[q];
            u32 ts = LT[s][i];
            while (ptr < ts) {
                u32 j = SD[s][ptr++];
                u32 r = RP[a][j];
                active[r >> 6] |= 1ULL << (r & 63);
            }
            u32 ta = LT[a][i];
            u32 tt[3] = { LT[d0][i], LT[d1][i], LT[d2][i] };
            u32 bb[3], cnt[3];
            u32 md = 0;                       // bit k set <=> dim k uses the ceil boundary
            for (int k = 0; k < 3; k++) {
                u32 t_ = tt[k], blk = t_ / B, r_ = t_ - blk * B;
                if (r_ * 2 <= B) {
                    bb[k] = blk * B;
                    cnt[k] = r_;
                } else {
                    md |= 1u << k;
                    bb[k] = (blk + 1) * B;
                    u32 hi = bb[k] > (u32)n ? (u32)n : bb[k];
                    cnt[k] = hi - t_;
                }
            }
            u32 words = ta >> 6, rem = ta & 63;
#ifdef RPF
            {
                for (int k = 0; k < 3; k++) {
                    u32 st = (md >> k) & 1 ? tt[k] : bb[k];
                    const char *y = (const char *)((k == 0 ? C0 : k == 1 ? C1 : C2) + (size_t)st * 4);
                    for (int z = 0; z < RPF; z++) _mm_prefetch(y + 64 * z, _MM_HINT_T0);
                }
            }
#endif
            const u64 *E0 = T0 + EOFFS[bb[0] / B], *E1 = T1 + EOFFS[bb[1] / B], *E2 = T2 + EOFFS[bb[2] / B];
            u64 res = and4pc(active, E0, E1, E2, words);
            if (rem) {
                u64 m = (1ULL << rem) - 1;
                res += popc64(active[words] & E0[words] & E1[words] & E2[words] & m);
            }
            token += 2;
            u32 tkA = token - 1, tkC = token;
            {   // dim 0
                u32 st = (md & 1) ? tt[0] : bb[0];
                u32 add = residual(C0 + (size_t)st * 4, sd0 + st, cnt[0], ta, ts,
                                   (md & 1) ? bb[1] : tt[1], (md & 1) ? bb[2] : tt[2],
                                   (md & 1) ? tkC : tkA, MARKV);
                if (md & 1) res -= add; else res += add;
            }
            {   // dim 1
                u32 st = (md & 2) ? tt[1] : bb[1];
                u32 add = residual(C1 + (size_t)st * 4, sd1 + st, cnt[1], ta, ts,
                                   (md & 2) ? bb[2] : tt[2], (md & 2) ? bb[0] : tt[0],
                                   (md & 2) ? tkC : tkA, MARKV);
                if (md & 2) res -= add; else res += add;
            }
            {   // dim 2
                u32 st = (md & 4) ? tt[2] : bb[2];
                u32 add = residual(C2 + (size_t)st * 4, sd2 + st, cnt[2], ta, ts,
                                   (md & 4) ? bb[0] : tt[0], (md & 4) ? bb[1] : tt[1],
                                   (md & 4) ? tkC : tkA, MARKV);
                if (md & 4) res -= add; else res += add;
            }
            out[i] = (unsigned)res;
        }
    }
}

static int inited_ = 0;
static void INIT() {
    if (inited_) return;
    inited_ = 1;
    int nbm = (NMAX - 1) / B + 2;
    EOFFS[0] = 0;
    for (int e = 0; e < nbm + 1; e++) {
        size_t need = (size_t)(e + 1) * LW;
        size_t padded = (need + 7) & ~(size_t)7;
        EOFFS[e + 1] = EOFFS[e] + padded;
    }
    posix_memalign((void **)&ACC, 64, ((size_t)NMAX / 64 + 1) * 8);
    for (int a = 0; a < KD; a++)
        for (int k = 0; k < 3; k++) {
            posix_memalign((void **)&TT[a][k], 64, (size_t)EOFFS[nbm + 1] * 8);
            posix_memalign((void **)&CX[a][k], 64, (size_t)NMAX * 4 * 4);
        }
}

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

评测结果