代码

/* 3D dominance counting: out[i] = #{ j : x[j]<x[i] && y[j]<y[i] && z[j]<z[i] }.
 * CDQ over the x-sorted (time) axis split at x-group boundaries; merge by y;
 * Fenwick tree over z.  Base blocks (<=128) solved with an SSE2 brute force and
 * sorted by y with a bottom-up merge. */
#include <emmintrin.h>

typedef unsigned u32;

static int N;
static const u32 *PX;
static u32 *AY, *AZ;
static int *AO;
static u32 *TMPY, *TMPZ;
static int *TMPO;
static u32 *BIT;
static u32 *ANS;
static int *GSTART, *GNEXT;

static inline void bit_add(int i) { for (; i <= N; i += i & (-i)) BIT[i]++; }
static inline u32 bit_sum(int i) { u32 s = 0; for (; i > 0; i -= i & (-i)) s += BIT[i]; return s; }

static inline void sort_y(int l, int r) {
    for (int i = l + 1; i < r; i++) {
        u32 vy = AY[i], vz = AZ[i]; int vo = AO[i];
        int j = i - 1;
        while (j >= l && AY[j] > vy) { AY[j+1]=AY[j]; AZ[j+1]=AZ[j]; AO[j+1]=AO[j]; j--; }
        AY[j+1]=vy; AZ[j+1]=vz; AO[j+1]=vo;
    }
}

static void merge_sort_y(int l, int r) {
    int len = r - l;
    if (len < 2) return;
    if (len <= 16) { sort_y(l, r); return; }
    for (int w = 1; w < len; w <<= 1) {
        for (int s = l; s < r; s += (w << 1)) {
            int m = s + w; if (m > r) m = r;
            int e = s + (w << 1); if (e > r) e = r;
            if (m >= e) continue;
            int a = s, b = m, k = s;
            while (a < m && b < e) {
                if (AY[a] <= AY[b]) { TMPY[k]=AY[a]; TMPZ[k]=AZ[a]; TMPO[k]=AO[a]; a++; }
                else { TMPY[k]=AY[b]; TMPZ[k]=AZ[b]; TMPO[k]=AO[b]; b++; }
                k++;
            }
            while (a < m) { TMPY[k]=AY[a]; TMPZ[k]=AZ[a]; TMPO[k]=AO[a]; a++; k++; }
            while (b < e) { TMPY[k]=AY[b]; TMPZ[k]=AZ[b]; TMPO[k]=AO[b]; b++; k++; }
            for (k = s; k < e; k++) { AY[k]=TMPY[k]; AZ[k]=TMPZ[k]; AO[k]=TMPO[k]; }
        }
    }
}

#define BASE 128
#define PAD (BASE + 4)

static void cdq(int l, int r) {
    int len = r - l;
    if (len <= 1) return;
    if (len <= BASE) {
        u32 lx[PAD], ly[PAD], lz[PAD], lacc[PAD];
        int lo[PAD];
        int p4 = (len + 3) & ~3;
        for (int k = 0; k < p4; k++) {
            if (k < len) { int t = l + k; lx[k] = PX[AO[t]]; ly[k] = AY[t]; lz[k] = AZ[t]; lo[k] = AO[t]; }
            else { lx[k] = 0x7fffffff; ly[k] = 0x7fffffff; lz[k] = 0x7fffffff; lo[k] = 0; }
            lacc[k] = 0;
        }
        const __m128i one = _mm_set1_epi32(1);
        for (int a = 0; a < len; a++) {
            __m128i xa = _mm_set1_epi32((int)lx[a]);
            __m128i ya = _mm_set1_epi32((int)ly[a]);
            __m128i za = _mm_set1_epi32((int)lz[a]);
            for (int b = a + 1; b < p4; b += 4) {
                __m128i bx = _mm_loadu_si128((const __m128i *)(lx + b));
                __m128i by = _mm_loadu_si128((const __m128i *)(ly + b));
                __m128i bz = _mm_loadu_si128((const __m128i *)(lz + b));
                __m128i m = _mm_and_si128(_mm_cmplt_epi32(xa, bx), _mm_cmplt_epi32(ya, by));
                m = _mm_and_si128(m, _mm_cmplt_epi32(za, bz));
                m = _mm_and_si128(m, one);
                __m128i acc = _mm_loadu_si128((const __m128i *)(lacc + b));
                _mm_storeu_si128((__m128i *)(lacc + b), _mm_add_epi32(acc, m));
            }
        }
        for (int k = 0; k < len; k++) ANS[lo[k]] += lacc[k];
        merge_sort_y(l, r);
        return;
    }
    if (GSTART[l] == GSTART[r - 1]) { merge_sort_y(l, r); return; }
    int mid = (l + r) >> 1;
    {
        int gs = GSTART[mid];
        int c1 = gs, c2 = GNEXT[gs];
        int h = (l + r) >> 1;
        int d1 = (c1 > l && c1 < r) ? (c1 > h ? c1 - h : h - c1) : 0x7fffffff;
        int d2 = (c2 > l && c2 < r) ? (c2 > h ? c2 - h : h - c2) : 0x7fffffff;
        mid = (d1 <= d2) ? c1 : c2;
    }
    if (mid <= l || mid >= r) { merge_sort_y(l, r); return; }

    cdq(l, mid);
    cdq(mid, r);

    int a = l;
    for (int j = mid; j < r; j++) {
        u32 yj = AY[j];
        while (a < mid && AY[a] < yj) { bit_add((int)AZ[a] + 1); a++; }
        ANS[AO[j]] += bit_sum((int)AZ[j]);
    }
    for (int k = l; k < a; k++) { int i = (int)AZ[k] + 1; for (; i <= N; i += i & (-i)) BIT[i]--; }

    int i = l, k = l, m = mid;
    while (i < mid && m < r) {
        if (AY[i] <= AY[m]) { TMPY[k]=AY[i]; TMPZ[k]=AZ[i]; TMPO[k]=AO[i]; i++; }
        else { TMPY[k]=AY[m]; TMPZ[k]=AZ[m]; TMPO[k]=AO[m]; m++; }
        k++;
    }
    while (i < mid) { TMPY[k]=AY[i]; TMPZ[k]=AZ[i]; TMPO[k]=AO[i]; i++; k++; }
    while (m < r) { TMPY[k]=AY[m]; TMPZ[k]=AZ[m]; TMPO[k]=AO[m]; m++; k++; }
    for (k = l; k < r; k++) { AY[k]=TMPY[k]; AZ[k]=TMPZ[k]; AO[k]=TMPO[k]; }
}

void count_3d(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
    N = n; PX = x;
    if (n <= 0) return;
    static u32 *buf = 0;
    if (!buf || (int)buf[0] < n) {
        delete[] buf;
        size_t words = (size_t)n * 8 + 4 * n + 16;
        buf = new u32[words];
        buf[0] = n;
        AY = buf + 1;
        AZ = AY + n;
        TMPY = AZ + n;
        TMPZ = TMPY + n;
        BIT = TMPZ + n;
        ANS = BIT + n + 1;
        AO = (int *)(ANS + n);
        TMPO = AO + n;
        GSTART = TMPO + n;
        GNEXT = GSTART + n;
    }
    u32 *cnt = (u32 *)(GNEXT + n + 1);
    for (int i = 0; i <= n; i++) cnt[i] = 0;
    for (int i = 0; i < n; i++) cnt[x[i]]++;
    u32 s = 0;
    for (int v = 0; v <= n; v++) { u32 c = cnt[v]; cnt[v] = s; s += c; }
    for (int i = 0; i < n; i++) { int p = (int)cnt[x[i]]++; AY[p] = y[i]; AZ[p] = z[i]; AO[p] = i; }
    for (int i = 0; i <= n; i++) BIT[i] = 0;
    for (int i = 0; i < n; i++) ANS[i] = 0;
    {
        int t = 0;
        while (t < n) {
            u32 xv = PX[AO[t]];
            int e = t;
            while (e < n && PX[AO[e]] == xv) e++;
            for (int k = t; k < e; k++) GSTART[k] = t;
            GNEXT[t] = e;
            t = e;
        }
        GNEXT[n] = n;
    }
    cdq(0, n);
    for (int t = 0; t < n; t++) out[AO[t]] = ANS[AO[t]];
}

评测结果