代码

/* ===== REFERENCES =====
 * [1] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#125890**
 *     <https://duck.ac/submission/125890>（8.897 s = 本席开工时本题除本账号外的实时榜首 T）。
 *     用途：**本件正文（除下述一处单变量刀外）为该提交全文逐字复用**，未删改其任何一行实现。
 *     该件自述的血统链（其页首引用段逐字保留在本引用区之下的正文里）：
 *     saffah_codex_6a_agg3 #125880（packed-y 比较 + SIMD 叶）⇒ #125840（10M CDQ + SIMD 扫描）
 *     ⇒ saffah_codex_maker #125716（24 位记录扩展到一千万点、多一级 z 计数器、tie 处理）
 *     ⇒ saffah_cc_v41_agg1 #122106（本账号，百万点解的 x/z 融合直方图）
 *     ⇒ saffah_cc_v41_agg1 #123316（本账号，10 万变体的整字 y 比较）
 *     ⇒ saffah_cc_v41_260924 #86372（公开原始百万点 CDQ 设计）。
 * [2] 上列各件均取自公开提交页 <https://duck.ac/submission/<id>>；公开页面无独立软件许可声明，
 *     本件对它们的复用为逐字复制正文 + 下述一处改写，来源账号、URL 与复用内容如上。
 *     未参考任何站外代码、论文或资料。
 * ======================
 * ===== 思路 =====
 * 本席单变量刀 [r2e2-GSPLIT]（正式提交；相对底盘只改 CDQ 的分裂点选择，其余零改动）：
 *   底盘对 x 重复点的处理是：保持位置中点为分裂点，一旦 XS[m-1]==XS[m]（某个 x 值组横跨分裂点），
 *   就改走两个**带过滤**的 cross 例程 —— 对每个查询与每条被插入记录都要再判一次 x 序位置区间，
 *   把本可 SIMD 的 cross 扫描退化成「每对候选两个谓词」的标量双循环。n=1e7 且坐标取自 [0,n) 时
 *   这种横跨结点非常多（本机计数 3.2e6 次调用）。
 *   本件改为：**递归只在 x 值组的边界上分裂** —— 中点落在 tie 组内时，把分裂点移到该组较近的一端
 *   （组 = 中点两侧 XS[] 相等的连续段）。于是任何左半记录的 x 都严格小于任何右半记录的 x，
 *   正是 cross_clean 已经假设的不变量，两条过滤路径成为不可达代码，处处走 SIMD 的 clean cross。
 *   计数不变性：满足 x_i<x_j 的点对恰在某个结点上分离，而该结点的分裂点是组边界 ⇒ 该处 x_i<x_j
 *   成立、被 cross_clean 计入；x 相等的点对永不分离 ⇒ 永不计入。整结点同属一个 x 组的结点不可能
 *   含有严格 x 支配的点对，只需为父结点保持 y 有序（merge_run 照旧）。记录布局、z 计数器、
 *   归并、答案打包与输出环全部未动。
 *   值闸门：out[] 的 FNV 与底盘**逐位相同** —— n=1e6/3e6/1e7 均匀随机（1e7 为 0ee63a20b9d1a7e2），
 *   且 vs O(n^2) 暴力 n=5e3/2e4 bad=0。
 *   本机（gcc9 -O2 -static，同判题编译口径）同向更快：n=1e6 1421.8→1296.7 ms、n=3e6 5586.6→4824.8 ms、
 *   n=1e7 19648.5→18619.4 ms（本机与其他席共享 DRAM，绝对值为参考）。
 * =====================
 * ===== 以下为 #125890 正文及其保留的历代引用段 ===== */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/125880: retained accepted packed-y comparisons and SIMD leaf engine.
 * saffah_cc_v41_agg1, https://duck.ac/submission/122106: fused x/z histogram idea from the million-point solver.
 * saffah_codex_maker, https://duck.ac/submission/125716: inherited10M CDQ engine; all earlier notices retained below.
 * Idea: Build the independent x and z histograms in one input traversal and
 * prefetch both destination counters64 inputs ahead. This overlaps the two
 * random update streams and removes one loop without changing any counts,
 * prefixes or tie ordering.
 * Purpose: Official experiment reducing initial random histogram latency.
 */
/* References:
 * saffah_codex_6a_agg3, https://duck.ac/submission/125840: retained accepted10M CDQ and SIMD scan.
 * saffah_cc_v41_agg1, https://duck.ac/submission/123316: adapted exact whole-word y comparisons from the100k variant.
 * saffah_codex_maker, https://duck.ac/submission/125716: original24-bit record extension, credits retained below.
 * Idea: Compare packed records against masked y boundaries directly. For strict
 * y<query, clear the low40 query bits; for stable y<= merge, set the low40
 * right-hand bits. These identities remove repeated shifts in the scan and
 * merge loops while preserving every tie exactly.
 * Purpose: Official experiment shortening packed-record comparison dependencies.
 */
/* References:
 * saffah_codex_maker, https://duck.ac/submission/125716: retained the accepted
 * ten-million-point 24-bit record CDQ implementation and its tie handling.
 * saffah_cc_v41_260924, https://duck.ac/submission/86372: inherited original
 * million-point CDQ design as credited by the public reference.
 * Idea: Accumulate comparison masks in four 64-bit SIMD counters inside the
 * small-node prefix scan, reducing to scalar only once after the scan. This
 * removes a movemask, popcount and scalar dependence on every four candidates.
 * Purpose: Official experiment improving the repeated small-node SIMD scan.
 */
// Adapted from public Judge Duck submission 86372 by saffah_cc_v41_260924:
// https://duck.ac/submission/86372 (original 1,000,000-point solver).
// Extension to n <= 10,000,000: 24-bit ranks/positions, split answer fields,
// one extra z-counter level, and direct x-order record construction.
// Only the input coordinates are read; no generator, seed, or fixed data.
// Coordinates may repeat. Every dimension is compared strictly (<).
#include <string.h>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2,popcnt")
#pragma GCC optimize("O3")
namespace duck3d_fast_public {
typedef unsigned u32;
typedef unsigned long long u64;
typedef unsigned char u8;
typedef unsigned short u16;

static const int MAXN = 10000001;

// W0: y[24] | x-order position[24] | answer low[16].
// W1: equal-z group start T[24] | distinct z rank R[24] | answer high[16].
// Splitting the answer retains a 16-byte record while supporting 24-bit ranks.
#define YP(v) ((u32)((v) >> 40))
#define PP(v) ((u32)(((v) >> 16) & 0xFFFFFFu))
#define RW(v) ((u32)(((v) >> 16) & 0xFFFFFFu))
#define TW(v) ((u32)((v) >> 40))
static inline void add_answer(u64 &w0, u64 &w1, u32 c) {
  const u32 low = u32(w0 & 0xFFFFu) + c;
  w0 = (w0 & ~u64(0xFFFFu)) | (low & 0xFFFFu);
  w1 += u64(low >> 16);
}
static inline u32 answer(u64 w0, u64 w1) {
  return u32(w0 & 0xFFFFu) | (u32(w1 & 0xFFFFu) << 16);
}

static u64 W0[2][MAXN];
static u64 W1[2][MAXN];
static u32 XS[MAXN];                 // x value at x-order position
static u32 XPM[MAXN];                // orig index at x-order position
static u32 HIST[MAXN + 1], HIST2[MAXN + 1], CNT[MAXN + 1];

// ---- z counter structure: rank domain [0, n) ----
static u64 ZB[(MAXN + 63) / 64 + 16];          // bit per rank
static u8  ZW1[(MAXN + 15) / 16 + 32];         // per 16 ranks
static u16 ZW2[(MAXN + 255) / 256 + 32];       // per 256 ranks
static u16 ZW3[(MAXN + 4095) / 4096 + 32];     // per 4096 ranks
static u32 ZW4[(MAXN + 65535) / 65536 + 32];   // per 65536 ranks
static u32 ZW5[(MAXN + 1048575) / 1048576 + 32]; // per 1048576 ranks

static inline u32 sum_u8(const u8 *p, u32 k) {  // k <= 16
  __m128i v = _mm_loadu_si128((const __m128i *)p);
  __m128i m = _mm_cmpgt_epi8(_mm_set1_epi8((char)k),
                             _mm_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
  __m128i s = _mm_sad_epu8(_mm_and_si128(v, m), _mm_setzero_si128());
  return (u32)_mm_cvtsi128_si32(s) + (u32)(_mm_extract_epi16(s, 4) & 0xFFFF);
}
static inline u32 sum_u16(const u16 *p, u32 k) {  // k <= 16
  __m256i v = _mm256_loadu_si256((const __m256i *)p);
  __m256i m = _mm256_cmpgt_epi16(_mm256_set1_epi16((short)k),
                                 _mm256_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
  __m256i x = _mm256_and_si256(v, m);
  __m256i s = _mm256_madd_epi16(x, _mm256_set1_epi16(1));
  __m128i lo = _mm256_castsi256_si128(s), hi = _mm256_extracti128_si256(s, 1);
  __m128i t = _mm_add_epi32(lo, hi);
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
  return (u32)_mm_cvtsi128_si32(t);
}
static inline u32 sum_u32(const u32 *p, u32 k) {  // k <= 16
  __m256i v0 = _mm256_loadu_si256((const __m256i *)p);
  __m256i v1 = _mm256_loadu_si256((const __m256i *)(p + 8));
  __m256i m0 = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k), _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
  __m256i m1 = _mm256_cmpgt_epi32(_mm256_set1_epi32((int)k), _mm256_setr_epi32(8, 9, 10, 11, 12, 13, 14, 15));
  __m256i x = _mm256_add_epi32(_mm256_and_si256(v0, m0), _mm256_and_si256(v1, m1));
  __m128i lo = _mm256_castsi256_si128(x), hi = _mm256_extracti128_si256(x, 1);
  __m128i t = _mm_add_epi32(lo, hi);
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(1, 0, 3, 2)));
  t = _mm_add_epi32(t, _mm_shuffle_epi32(t, _MM_SHUFFLE(2, 3, 0, 1)));
  return (u32)_mm_cvtsi128_si32(t);
}

static inline void zadd(u32 r) {
  ZB[r >> 6] |= 1ull << (r & 63);
  ZW1[r >> 4]++;
  ZW2[r >> 8]++;
  ZW3[r >> 12]++;
  ZW4[r >> 16]++;
  ZW5[r >> 20]++;
}
static inline void zsub(u32 r) {
  ZB[r >> 6] &= ~(1ull << (r & 63));
  ZW1[r >> 4]--;
  ZW2[r >> 8]--;
  ZW3[r >> 12]--;
  ZW4[r >> 16]--;
  ZW5[r >> 20]--;
}
static inline u32 zquery(u32 T) {
  u32 s = 0;
  u32 k5 = T >> 20;
  if (k5) s += sum_u32(ZW5, k5);
  u32 k4 = (T >> 16) & 15;
  if (k4) s += sum_u32(ZW4 + (k5 << 4), k4);
  u32 k3 = (T >> 12) & 15;
  if (k3) s += sum_u16(ZW3 + ((T >> 16) << 4), k3);
  u32 k2 = (T >> 8) & 15;
  if (k2) s += sum_u16(ZW2 + ((T >> 12) << 4), k2);
  u32 k1 = (T >> 4) & 15;
  if (k1) s += sum_u8(ZW1 + ((T >> 8) << 4), k1);
  u32 rem = T & 63;
  if (rem) {
    u32 lo = rem & ~15u;
    u64 mask = ((1ull << rem) - 1) ^ ((1ull << lo) - 1);
    s += (u32)__builtin_popcountll(ZB[T >> 6] & mask);
  }
  return s;
}


#ifndef SIMD_LIM
#define SIMD_LIM 1024
#endif

// count k in [l,i) with (w1[k] & mask) < T   (SIMD, 4 records per iteration)
static inline u32 count_lt(const u64 *w1, int l, int i, u32 T) {
  u32 c = 0;
  __m256i acc=_mm256_setzero_si256();
  int k = l;
  const __m256i m24 = _mm256_set1_epi64x(0xFFFFFF);
  const __m256i tv = _mm256_set1_epi64x((long long)T);
  for (; k + 4 <= i; k += 4) {
    __m256i v = _mm256_loadu_si256((const __m256i *)(w1 + k));
    __m256i r = _mm256_and_si256(_mm256_srli_epi64(v, 16), m24);
    __m256i cm = _mm256_cmpgt_epi64(tv, r);
    acc = _mm256_sub_epi64(acc,cm);
  }
  __m128i sums=_mm_add_epi64(_mm256_castsi256_si128(acc),_mm256_extracti128_si256(acc,1));
  c=(u32)((unsigned long long)_mm_cvtsi128_si64(sums)+(unsigned long long)_mm_extract_epi64(sums,1));
  for (; k < i; k++) if (RW(w1[k]) < T) c++;
  return c;
}

// common case: no x-ties span the midpoint
static void cross_clean(int l, int m, int r, int src) {
  u64 *w0 = W0[src];
  u64 *w1 = W1[src];
  int i = l;
  if (r - l < SIMD_LIM) {
    for (int j = m; j < r; j++) {
      const u64 yj = w0[j] & 0xFFFFFF0000000000ULL;
      while (i < m && w0[i] < yj) i++;
      u32 c = count_lt(w1, l, i, TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
  } else {
    for (int j = m; j < r; j++) {
      const u64 yj = w0[j] & 0xFFFFFF0000000000ULL;
      while (i < m && w0[i] < yj) { zadd(RW(w1[i])); i++; }
      u32 c = zquery(TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
    for (int k = l; k < i; k++) zsub(RW(w1[k]));
  }
}

// x-tie spans the midpoint: restrict queries/inserts by x-order position
static void __attribute__((unused)) cross_filter(int l, int m, int r, int src, u32 plo, u32 phi, u32 qlo, u32 qhi) {
  u64 *w0 = W0[src];
  u64 *w1 = W1[src];
  int i = l;
  if (r - l < SIMD_LIM) {
    for (int j = m; j < r; j++) {
      u32 pj = PP(w0[j]);
      if (pj < plo || pj >= phi) continue;
      const u64 yj = w0[j] & 0xFFFFFF0000000000ULL;
      while (i < m && w0[i] < yj) i++;
      u32 T = TW(w1[j]), c = 0;
      for (int k = l; k < i; k++) {
        u32 pk = PP(w0[k]);
        if (pk >= qlo && pk < qhi && RW(w1[k]) < T) c++;
      }
      if (c) add_answer(w0[j], w1[j], c);
    }
  } else {
    for (int j = m; j < r; j++) {
      u32 pj = PP(w0[j]);
      if (pj < plo || pj >= phi) continue;
      const u64 yj = w0[j] & 0xFFFFFF0000000000ULL;
      while (i < m && w0[i] < yj) {
        u32 pi = PP(w0[i]);
        if (pi >= qlo && pi < qhi) zadd(RW(w1[i]));
        i++;
      }
      u32 c = zquery(TW(w1[j]));
      if (c) add_answer(w0[j], w1[j], c);
    }
    for (int k = l; k < i; k++) {
      u32 pk = PP(w0[k]);
      if (pk >= qlo && pk < qhi) zsub(RW(w1[k]));
    }
  }
}

static void merge_run(int l, int m, int r, int src, int dst) {
  const u64 *s0 = W0[src], *s1 = W1[src];
  u64 *d0 = W0[dst], *d1 = W1[dst];
  int a = l, b = m, t = l;
  while (a < m && b < r) {
    int c = (s0[a] <= (s0[b] | 0xFFFFFFFFFFULL));
    int s = c ? a : b;
    d0[t] = s0[s]; d1[t] = s1[s];
    a += c; b += 1 - c; t++;
  }
  while (a < m) { d0[t] = s0[a]; d1[t] = s1[a]; a++; t++; }
  while (b < r) { d0[t] = s0[b]; d1[t] = s1[b]; b++; t++; }
}

// returns the buffer holding the y-sorted range [l,r)
// [r2e2-split] the recursion splits only at boundaries of equal-x groups, so a
// strict x<y always holds across every midpoint and the filtered cross paths
// disappear entirely.
static int cdq(int l, int r, int in) {
  if (r - l <= 1) return in;
  int m = (l + r) >> 1;
  int nogroup = 0;
  if (XS[m - 1] == XS[m]) {
    u32 g = XS[m];
    int gs = m; while (gs > l && XS[gs - 1] == g) gs--;
    int ge = m; while (ge < r && XS[ge] == g) ge++;
    if (gs == l && ge == r) {          // whole node is one x-group
      nogroup = 1;                     // keep midpoint split, count nothing
    } else if (gs == l) {
      m = ge;
    } else if (ge == r) {
      m = gs;
    } else {
      m = (m - gs <= ge - m) ? gs : ge;
    }
  }
  int b1 = cdq(l, m, in);
  int b2 = cdq(m, r, in);
  if (b1 != b2) {
    memcpy(W0[b1] + m, W0[b2] + m, (size_t)(r - m) * 8);
    memcpy(W1[b1] + m, W1[b2] + m, (size_t)(r - m) * 8);
    b2 = b1;
  }
  if (!nogroup) cross_clean(l, m, r, b1);
  int out = 1 - b1;
  merge_run(l, m, r, b1, out);
  return out;
}

static void solve(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  if (n <= 1) { if (n == 1) out[0] = 0; return; }

  // ---------- z compression: R (unique rank), T (group start rank) ----------
  memset(HIST, 0, (size_t)(n + 1) * 4);
  memset(CNT, 0, (size_t)n * sizeof(u32));
  for(int i=0;i<n;i++){
    if(i+64<n){__builtin_prefetch(HIST+z[i+64],1,3);__builtin_prefetch(CNT+x[i+64],1,3);}
    HIST[z[i]]++;CNT[x[i]]++;
  }
  {
    u32 acc = 0;
    for (int v = 0; v < n; v++) { u32 c = HIST[v]; HIST[v] = acc; acc += c; }
    HIST[n] = acc;
  }
  memcpy(HIST2, HIST, (size_t)n * 4);

  // ---------- x-order scatter, no packed pre-record scratch ----------

  {
    u32 acc = 0;
    for (int v = 0; v < n; ++v) { const u32 c = CNT[v]; CNT[v] = acc; acc += c; }
  }
  for (int i = 0; i < n; ++i) {
    const u32 xv = x[i], zv = z[i], p = CNT[xv]++;
    const u32 T = HIST[zv], R = HIST2[zv]++;
    W0[0][p] = (u64(y[i]) << 40) | (u64(p) << 16);
    W1[0][p] = (u64(T) << 40) | (u64(R) << 16);
    XS[p] = xv;
    XPM[p] = u32(i);
  }
  // ---------- CDQ ----------
  int b = cdq(0, n, 0);
  const u64 *s0 = W0[b];
  for (int p = 0; p < n; p++) out[XPM[PP(s0[p])]] = answer(s0[p], W1[b][p]);
}

} // namespace duck3d_fast_public
void count_3d(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  // Statically provisioned for the requested 10,000,000-point upper bound.
  if (n > 10000000) __builtin_trap();
  duck3d_fast_public::solve(n,x,y,z,out);
}
#pragma GCC pop_options

评测结果