代码

// 1008e8 -- 2D dominance counting, n = 1e8:  out[i] = #{j: x[j]<x[i] and y[j]<y[i]}
//
// v6 = v5's memory-fitted design with the random y gathers removed.  On the judge
// v5 was AC at 14.263 s, and a variant whose gather was replaced by a cheap hash
// of idx ran 6.207 s -- so ~8 s of v5 was 1e8 random reads into the 400 MB y
// array (4 KB pages -> TLB miss + page walk on almost every one).  Here y is
// permuted in lockstep with the x radix sort (an extra 4-byte stream per pass,
// ~2.4 GB of extra traffic in total, far cheaper than the gathers).
//
// The problem is split by the y median so each half fits:
//   y_i <  K : out[i] = #{j: y_j<y_i, x_j<x_i}            (all such j are low)
//   y_i >= K : out[i] = #{j low: x_j < x_i} + #{j high: y_j<y_i, x_j<x_i}
// The cross term is one sequential merge of the two x-sorted key arrays.
//
// Record encoding
//   during split/sort : [spare:11 | idx:27 | x:27]        (x in the low bits)
//   after attach      : [cnt:11  | idx:27 | y:26]         (y shifted into [0,2^26))
// The per-element count is 27 bits: 11 in the record + the top 16 in a u16 array,
// so every count travels with its element and out[] never receives a random RMW.
//
// Buffers (M = n/2 + slack): RA,RB,RC u64[M] (~400 MB each), Y1,Y2 u32[M]
// (~200 MB each), U1 u16[M] (100 MB) = ~1.66 GB of our own, plus the tasklib's
// out (400 MB) which holds the answers (and nothing else).
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <sys/mman.h>
#include <immintrin.h>

#ifndef NO_AVX2
#pragma GCC push_options
#pragma GCC target("avx2")
#endif

typedef uint16_t u16;
typedef uint32_t u32;
typedef uint64_t u64;

static const u64 M27 = 0x7FFFFFFull;
static const u64 LOW53 = 0x1FFFFFFFFFFFFFull;   // idx:27 | y:26

static u64 *RA, *RB, *RC;
static u32 *Y1, *Y2;
static u16 *U1;
static size_t gCap;
static const unsigned *gX, *gY;
static unsigned *gOut;

static u32 gCnt[4096];
static size_t gOutLim = 0;
static u32 gCur[4096];

// ---- branchless 8x8x8 cumulative digit counter: query = 3 loads, insert = 3
//      masked SIMD suffix adds over 8-lane u32; 2 KB of state, all L1 resident.
static u32 gP1[8];
static u32 gP2[8][8];
static u32 gP3[8][8][8];
static u32 gMsk[8][8] __attribute__((aligned(32)));
static int gMskInit = 0;

static inline void cnt_init(void) {
  if (!gMskInit) {
    for (u32 m = 0; m < 8; m++)
      for (u32 i = 0; i < 8; i++) gMsk[m][i] = (i > m);
    gMskInit = 1;
  }
  memset(gP1, 0, sizeof(gP1));
  memset(gP2, 0, sizeof(gP2));
  memset(gP3, 0, sizeof(gP3));
}
static inline u32 cnt_q(u32 d) {
  u32 a = d >> 6, b = (d >> 3) & 7, c = d & 7;
  return gP1[a] + gP2[a][b] + gP3[a][b][c];
}
static inline void cnt_ins(u32 d) {
  u32 a = d >> 6, b = (d >> 3) & 7, c = d & 7;
  __m256i one = _mm256_set1_epi32(1);
  __m256i m3 = _mm256_load_si256((const __m256i *)gMsk[c]);
  __m256i v3 = _mm256_loadu_si256((const __m256i *)&gP3[a][b][0]);
  _mm256_storeu_si256((__m256i *)&gP3[a][b][0],
                      _mm256_add_epi32(v3, _mm256_and_si256(m3, one)));
  __m256i m2 = _mm256_load_si256((const __m256i *)gMsk[b]);
  __m256i v2 = _mm256_loadu_si256((const __m256i *)&gP2[a][0]);
  _mm256_storeu_si256((__m256i *)&gP2[a][0],
                      _mm256_add_epi32(v2, _mm256_and_si256(m2, one)));
  __m256i m1 = _mm256_load_si256((const __m256i *)gMsk[a]);
  __m256i v1 = _mm256_loadu_si256((const __m256i *)&gP1[0]);
  _mm256_storeu_si256((__m256i *)&gP1[0],
                      _mm256_add_epi32(v1, _mm256_and_si256(m1, one)));
}

static void *bigalloc(size_t bytes) {
  void *p = mmap(0, bytes, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
  if (p != MAP_FAILED) return p;          // mmap fails on JudgeDuck -> malloc
  return malloc(bytes);
}

// ------------------------------------------------------------------
// 3-pass LSD radix sort by x of [spare:11|idx:27|x:27] records, permuting the
// parallel y array with the identical permutation; the result ends in X1 / Yb.
static void sort_x(u64 *X0, u64 *X1, u32 *Ya, u32 *Yb, size_t m) {
  memset(gCnt, 0, 512 * sizeof(u32));
  for (size_t i = 0; i < m; i++) gCnt[(u32)X0[i] & 511u]++;
  gCur[0] = 0;
  for (u32 b = 1; b < 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
  for (size_t i = 0; i < m; i++) {
    u64 rec = X0[i];
    u32 q = gCur[(u32)rec & 511u]++;
    X1[q] = rec;
    Yb[q] = Ya[i];
  }
  memset(gCnt, 0, 512 * sizeof(u32));
  for (size_t i = 0; i < m; i++) gCnt[(u32)(X1[i] >> 9) & 511u]++;
  gCur[0] = 0;
  for (u32 b = 1; b < 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
  for (size_t i = 0; i < m; i++) {
    u64 rec = X1[i];
    u32 q = gCur[(u32)(rec >> 9) & 511u]++;
    X0[q] = rec;
    Ya[q] = Yb[i];
  }
  memset(gCnt, 0, 512 * sizeof(u32));
  for (size_t i = 0; i < m; i++) gCnt[(u32)(X0[i] >> 18) & 511u]++;
  gCur[0] = 0;
  for (u32 b = 1; b < 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
  for (size_t i = 0; i < m; i++) {
    u64 rec = X0[i];
    u32 q = gCur[(u32)(rec >> 18) & 511u]++;
    X1[q] = rec;
    Yb[q] = Ya[i];
  }
}

// ------------------------------------------------------------------
// y is already in sequence order (Yv); fold the x-tie correction into the 27-bit
// count, rewrite the record as [cnt:11|idx:27|y:26], histogram y's top digit.
// `preseed` = the count already holds the half's cross term (high half).
static void attach(u64 *X, const u32 *Yv, u16 *Uv, size_t m, int preseed) {
  memset(gCnt, 0, 512 * sizeof(u32));
  size_t p = 0;
  while (p < m) {
    u64 xv = X[p] & M27;
    size_t q = p + 1;
    while (q < m && (X[q] & M27) == xv) q++;
    for (size_t a = p; a < q; a++) {
      u32 ia = (u32)((X[a] >> 27) & M27);
      u32 ya = Yv[a];
      u32 dup = 0;
      for (size_t b = p; b < a; b++) dup += (Yv[b] < ya);
      u32 cnt = preseed ? ((u32)(X[a] >> 54) | ((u32)Uv[a] << 10)) : 0u;
      cnt = (cnt - dup) & 0x7FFFFFFu;
      X[a] = ((u64)(cnt & 2047u) << 53) | ((u64)ia << 26) | (u64)ya;
      Uv[a] = (u16)(cnt >> 11);
      gCnt[(ya >> 17) & 511u]++;
    }
    p = q;
  }
}

// ------------------------------------------------------------------
// three MSD radix rank-accumulation levels over the 26-bit y field; X0 holds the
// x-sorted records on entry, Ca is the (u16) high part of the count, Cb a free
// u16 buffer.  Ends with the final store out[idx] = count.
static void levels(u64 *X0, u64 *X1, u16 *Ca, u16 *Cb, size_t m, unsigned *out) {
  // L1: y bits 17..25   (X0->X1, Ca->Cb)
  gCur[0] = 0;
  for (u32 b = 1; b <= 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
  for (u32 b = 0; b <= 512; b++) gCnt[1024 + b] = gCur[b];   // keep the offsets
  cnt_init();
  for (size_t i = 0; i < m; i++) {
    u64 rec = X0[i];
    u32 d = (u32)(rec >> 17) & 511u;
    u32 inc = cnt_q(d);
    cnt_ins(d);
    u32 q = gCur[d]++;
    u32 cnt = ((u32)(rec >> 53) | ((u32)Ca[i] << 11)) + inc;
    X1[q] = (rec & LOW53) | ((u64)(cnt & 2047u) << 53);
    Cb[q] = (u16)((cnt >> 11) & 0xFFFFu);
  }
  // L2: per L1 bucket, y bits 8..16   (X1->X0, Cb->Ca)
  for (u32 b0 = 0; b0 < 512; b0++) {
    size_t s = gCnt[1024 + b0], e = gCnt[1025 + b0];
    if (e - s <= 1) {
      if (e > s) { X0[s] = X1[s]; Ca[s] = Cb[s]; }
      continue;
    }
    memset(gCnt, 0, 512 * sizeof(u32));
    for (size_t t = s; t < e; t++) gCnt[(u32)(X1[t] >> 8) & 511u]++;
    gCur[0] = s;
    for (u32 b = 1; b < 512; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
    cnt_init();
    for (size_t t = s; t < e; t++) {
      u64 rec = X1[t];
      u32 d = (u32)(rec >> 8) & 511u;
      u32 inc = cnt_q(d);
      cnt_ins(d);
      u32 q = gCur[d]++;
      u32 cnt = ((u32)(rec >> 53) | ((u32)Cb[t] << 11)) + inc;
      X0[q] = (rec & LOW53) | ((u64)(cnt & 2047u) << 53);
      Ca[q] = (u16)((cnt >> 11) & 0xFFFFu);
    }
  }
  // L3: per L2 bucket, low 8 bits, in place; also histogram idx>>15 for L4
  memset(gCnt, 0, 4096 * sizeof(u32));
  {
    size_t i = 0;
    while (i < m) {
      size_t s = i;
      u32 pref = (u32)(X0[i] & 0x3FFFFFFull) >> 8;      // y bits 8..25
      while (i < m && (u32)(X0[i] & 0x3FFFFFFull) >> 8 == pref) i++;
      size_t e = i;
      for (size_t t = s; t < e; t++) gCnt[(u32)(X0[t] >> 41) & 0xFFFu]++;
      if (e - s <= 1) continue;
      cnt_init();
      for (size_t t = s; t < e; t++) {
        u64 rec = X0[t];
        u32 d = (u32)rec & 255u;
        u32 inc = cnt_q(d);
        cnt_ins(d);
        u32 cnt = ((u32)(rec >> 53) | ((u32)Ca[t] << 11)) + inc;
        X0[t] = (rec & LOW53) | ((u64)(cnt & 2047u) << 53);
        Ca[t] = (u16)((cnt >> 11) & 0xFFFFu);
      }
    }
  }
  // L4: partition by idx bits 15..26 (X0->X1, Ca->Cb), then store out[idx]
  gCur[0] = 0;
  for (u32 b = 1; b < 4096; b++) gCur[b] = gCur[b - 1] + gCnt[b - 1];
  for (size_t i = 0; i < m; i++) {
    u64 rec = X0[i];
    u32 q = gCur[(u32)(rec >> 41) & 0xFFFu]++;
    X1[q] = rec;
    Cb[q] = Ca[i];
  }
  for (size_t i = 0; i < m; i++) {
#ifdef DEBUG_BOUNDS
    u32 ix = (u32)((X1[i] >> 26) & M27);
    if (ix >= (u32)gOutLim) { printf("BAD idx=%u rec=%llx i=%zu m=%zu\n", ix, (unsigned long long)X1[i], i, m); exit(3); }
#endif
    out[(u32)((X1[i] >> 26) & M27)] = (u32)(X1[i] >> 53) | ((u32)Cb[i] << 11);
  }
}

// ------------------------------------------------------------------
void count_2d(int n, const unsigned *x, const unsigned *y, unsigned *out) {
  if (n <= 0) return;
  if (n == 1) { out[0] = 0; return; }
  if (n <= 2000) {
    memset(out, 0, (size_t)n * 4);
    for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++)
        if (x[j] < x[i] && y[j] < y[i]) out[i]++;
    return;
  }
  size_t N = (size_t)n;
  size_t M = N / 2 + 70000;                     // slack for the cut granularity
  if (gCap < M) {
    if (RA) { munmap(RA, gCap * 8); munmap(RB, gCap * 8); munmap(RC, gCap * 8); }
    if (Y1) { munmap(Y1, gCap * 4); munmap(Y2, gCap * 4); }
    if (U1) munmap(U1, gCap * 2);
    RA = (u64 *)bigalloc(M * 8);
    RB = (u64 *)bigalloc(M * 8);
    RC = (u64 *)bigalloc(M * 8);
    Y1 = (u32 *)bigalloc(M * 4);
    Y2 = (u32 *)bigalloc(M * 4);
    U1 = (u16 *)bigalloc(M * 2);
    gCap = M;
    if (!RA || !RB || !RC || !Y1 || !Y2 || !U1) { gCap = 0; return; }
  }
  gX = x; gY = y; gOut = out; gOutLim = N;

  // ---- choose the y cut K: exact bucket boundaries over the top 11 bits of y
  static u32 hist[2048];
  memset(hist, 0, sizeof(hist));
  for (size_t i = 0; i < N; i++) hist[y[i] >> 16]++;
  size_t want = N / 2, acc = 0;
  u32 K = 0;
  for (u32 b = 0; b < 2048; b++) {
    acc += hist[b];
    if (acc >= want) { K = (b + 1) << 16; break; }
  }
  if (K == 0) K = 1u << 16;
  size_t lo = acc, hi = N - acc;
  // constraints: both halves fit their buffers, the high half fits out (8 bytes
  // per record in 4 bytes per element), and both y ranges fit the 26-bit field
  if (hi == 0 || lo > M || hi > N / 2 || hi > M || (u64)K > (1u << 26) ||
      (u64)(N - K) > (1u << 26)) {
    memset(out, 0, N * 4);                      // pathological input: bail out
    for (size_t i = 0; i < N; i++)
      for (size_t j = 0; j < N; j++)
        if (x[j] < x[i] && y[j] < y[i]) out[i]++;
    return;
  }

  // ---- split A: low half -> RA (records) / Y1 (y, in the same order)
  {
    size_t a = 0;
    for (size_t i = 0; i < N; i++)
      if (y[i] < K) { RA[a] = ((u64)i << 27) | (u64)x[i]; Y1[a] = y[i]; a++; }
  }
  // ---- low half: x sort (records -> RC, y -> Y2)
  sort_x(RA, RC, Y1, Y2, lo);
  // keep the low half's sorted x values: the cross term needs them after the
  // low half's own pipeline has permuted its records away (Y1 is free again)
  for (size_t i = 0; i < lo; i++) Y1[i] = (u32)RC[i] & (u32)M27;
  attach(RC, Y2, U1, lo, 0);
  levels(RC, RA, U1, (u16 *)Y2, lo, out);       // low half answers -> out

  // ---- split B: high half -> RB (records) / Y2 (y shifted into [0,2^26))
  {
    size_t b = 0;
    for (size_t i = 0; i < N; i++)
      if (y[i] >= K) { RB[b] = ((u64)i << 27) | (u64)x[i]; Y2[b] = y[i] - K; b++; }
  }
  // ---- high half: x sort (records -> RC, y -> RA reused as a u32 buffer)
  sort_x(RB, RC, Y2, (u32 *)RA, hi);
  // ---- cross term: #{low: x_j < x_i} for every high element, one merge sweep
  {
    size_t a = 0; u32 run = 0;
    for (size_t b = 0; b < hi; b++) {
      u64 rec = RC[b];
      u32 xb = (u32)rec & (u32)M27;
      while (a < lo && Y1[a] < xb) { a++; run++; }
      RC[b] = rec | ((u64)(run & 1023u) << 54);   // low 10 bits into the record
      U1[b] = (u16)(run >> 10);                   // the top 16 into U1
    }
  }
  attach(RC, (u32 *)RA, U1, hi, 1);
  levels(RC, RB, U1, (u16 *)Y2, hi, out);        // high half answers -> out
}

评测结果