代码

/* 2D dominance counting: out[i] = #{ j : x[j]<x[i] && y[j]<y[i] }, n up to 1e7,
 * coordinates in [0,n).
 *
 * 1. stable counting sort by y, then by x  ->  points in (x,y) order.
 * 2. sweep that order maintaining a 4-level hierarchical counter over the y
 *    values seen so far (per-value u16 counts, plus u32 sums per 64 / 4096 /
 *    262144 values) so that rank(y) = #seen values < y costs ~60 cycles in
 *    L1/L2 with SSE2.
 * 3. within an x-group y is non-decreasing, so the same-x pairs counted by the
 *    sweep are exactly (index of first equal y in the group - group start). */
#include <stdlib.h>
#include <string.h>
#include <emmintrin.h>

typedef unsigned u32;
typedef unsigned short u16;
typedef unsigned long long u64;

static inline u32 hsum(__m128i acc) {
    __m128i h = _mm_add_epi32(acc, _mm_shuffle_epi32(acc, 0x4E));
    h = _mm_add_epi32(h, _mm_shuffle_epi32(h, 0xB1));
    return (u32)_mm_cvtsi128_si32(h);
}
static inline u32 sum_u32(const u32 *p, u32 cnt) {
    __m128i acc = _mm_setzero_si128();
    u32 k = 0;
    for (; k + 4 <= cnt; k += 4) acc = _mm_add_epi32(acc, _mm_loadu_si128((const __m128i *)(p + k)));
    u32 s = hsum(acc);
    for (; k < cnt; k++) s += p[k];
    return s;
}
static inline u32 sum_u16(const u16 *p, u32 cnt) {
    __m128i acc = _mm_setzero_si128();
    u32 k = 0;
    for (; k + 8 <= cnt; k += 8) {
        __m128i v = _mm_loadu_si128((const __m128i *)(p + k));
        __m128i z = _mm_setzero_si128();
        acc = _mm_add_epi32(acc, _mm_unpacklo_epi16(v, z));
        acc = _mm_add_epi32(acc, _mm_unpackhi_epi16(v, z));
    }
    u32 s = hsum(acc);
    for (; k < cnt; k++) s += p[k];
    return s;
}

void count_2d(int n, const unsigned *x, const unsigned *y, unsigned *out) {
    if (n <= 0) return;
    u32 W1 = (u32)(((u32)n + 63u) >> 6);
    u32 W2 = (u32)(((u32)n + 4095u) >> 12);
    u32 W3 = (u32)(((u32)n + 262143u) >> 18);

    u32 *cnt = (u32 *)calloc((size_t)n + 1, sizeof(u32));
    u32 *Sx = (u32 *)malloc((size_t)n * 4);
    u32 *Sy = (u32 *)malloc((size_t)n * 4);
    u32 *So = (u32 *)malloc((size_t)n * 4);
    u32 *Tx = (u32 *)malloc((size_t)n * 4);
    u32 *Ty = (u32 *)malloc((size_t)n * 4);
    u32 *To = (u32 *)malloc((size_t)n * 4);
    u32 *ans = (u32 *)malloc((size_t)n * 4);
    u16 *cval = (u16 *)calloc((size_t)n + 64, sizeof(u16));
    u32 *pop1 = (u32 *)calloc(W1 + 1, sizeof(u32));
    u32 *pop2 = (u32 *)calloc(W2 + 1, sizeof(u32));
    u32 *pop3 = (u32 *)calloc(W3 + 1, sizeof(u32));
    if (!cnt || !Sx || !Sy || !So || !Tx || !Ty || !To || !ans || !cval || !pop1 || !pop2 || !pop3) return;

    /* counting sort by y (stable) */
    for (int i = 0; i < n; i++) cnt[y[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++) { u32 p = cnt[y[i]]++; Sx[p] = x[i]; Sy[p] = y[i]; So[p] = (u32)i; }

    /* counting sort by x (stable) */
    memset(cnt, 0, ((size_t)n + 1) * 4);
    for (int i = 0; i < n; i++) cnt[Sx[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++) { u32 p = cnt[Sx[i]]++; Tx[p] = Sx[i]; Ty[p] = Sy[i]; To[p] = So[i]; }
    free(Sx); free(Sy); free(So); Sx = Sy = So = 0;

    /* sweep */
    for (int i = 0; i < n; i++) {
        u32 v = Ty[i];
        u32 w = v >> 6, g = v >> 12, G = v >> 18;
        u32 s = G ? sum_u32(pop3, G) : 0;
        u32 m = g - (G << 6);
        if (m) s += sum_u32(pop2 + (G << 6), m);
        u32 mg = g << 6;
        u32 mm = w - mg;
        if (mm) s += sum_u32(pop1 + mg, mm);
        u32 b = v & 63;
        if (b) s += sum_u16(cval + (w << 6), b);
        ans[i] = s;
        cval[v]++;
        pop1[w]++; pop2[g]++; pop3[G]++;
    }

    /* subtract same-x-group pairs */
    {
        int i = 0;
        while (i < n) {
            u32 xv = Tx[i];
            int e = i;
            while (e < n && Tx[e] == xv) e++;
            int k = i;
            while (k < e) {
                u32 yv = Ty[k];
                int k2 = k;
                while (k2 < e && Ty[k2] == yv) k2++;
                u32 sub = (u32)(k - i);
                for (int t = k; t < k2; t++) ans[t] -= sub;
                k = k2;
            }
            i = e;
        }
    }

    for (int i = 0; i < n; i++) out[To[i]] = ans[i];
    free(cnt); free(Tx); free(Ty); free(To); free(ans);
    free(cval); free(pop1); free(pop2); free(pop3);
}

评测结果