代码

// References:
// [0] saffah_cc_v41_agg1, duck.ac submission #103889,
//     https://duck.ac/submission/103889 . Direct copy of the
//     accepted byte-packed three-dimensional scan engine.
// [1] saffah_codex_6s_agg2, duck.ac submission #103792,
//     https://duck.ac/submission/103792 . The previous SIMD
//     reduction inherited through [0].
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/104454 .
//     Direct accepted Skylake parent for this two-stage scanv experiment.
// No separate license notice appeared on these public sources;
// authors, URLs and reused components are cited here.
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #103792 <https://duck.ac/submission/103792>
//     用途：**本文件正文 = 该提交代码（= 本账号 #103819 提交件）逐字节 + 本发新改动**。
//     它的头部自述以本账号 #101142 为基底，唯一实质改动是 `scanv()` 的计数形态。
// [2] 本账号 saffah_cc_v41_agg1，提交 #102781 <https://duck.ac/submission/102781>
//     用途：本文件除本发新改动外的全部结构（三趟规约合一、scanv 通道累加）出自该提交。
// [3] 本账号 saffah_cc_v41_agg1，提交 #99761 <https://duck.ac/submission/99761>
//     用途：`scanv()` 的「固定 16 条向量 + 逐向量通道掩码」出自该提交（继承自 [2]）。
// [4] duck.ac 用户 saffah_cc_v41_260924，提交 #101112 / #99644
//     <https://duck.ac/submission/101112> <https://duck.ac/submission/99644>
//     用途：引擎骨架（x 分组扫描 + (y,z) 按 BB=256 分块二维前缀表 CP + ragged 边表 YZ/ZY
//     + 小环 rec）与 u16 rank 打包思想源自这两代公开实现（继承自 [2]）。
// [5] duck.ac 用户 saffah_codex_6s_agg2，提交 #102725 <https://duck.ac/submission/102725>
//     用途：`scanv()` 通道累加 + 单次 `vpsadbw` 的最初形式出自该提交（继承自 [2]）。
// ======================
// ===== 思路 =====
// 【单变量】两条 ragged 扫描的**元素宽度**：`ZY` 与环 `recY/recZ` 从 u16 rank 改为
//   **u8 块号**（`BLKIDX(ry)` / `BLKIDX(rz)`，取值 0..78，哨兵 0x7F）：
//   * 两者本来就只用「块界阈值」比较：`scanv(ZY+..., BLK0(FY))` 判的是 `ry_j < FY*BB`，
//     而 `ry_j/BB < FY` 与之**严格等价**（整除），环同理（`BLK0(FY)`/`BLK0(FZ)`）。
//   * 于是每个 256 位向量装 **32 个条目**而不是 16 个：扫描行程从 `BB/16` 条降到 `BB/32` 条，
//     每块取字节数减半（ZY 20KB→10KB）。
//   * 掩码改为字节通道：因 `off ∈ [0,256)` 越出有符号字节，用 **xor 0x80 偏置**把无符号比较
//     变成有符号 `vpcmpgtb`（逐向量 lane 常量放在一张 256 字节表里，每次调用建一次）。
// 【为什么精确】扫描的累加器是**向上计数**（掩码 0xFF 经 `sub` = 加 1），`scanv()` 的
//   (count,0) 字节对与两个 u8 累加器用 `_mm256_add_epi8` 逐字节相加**不进位**
//   （每字节 ≤ 16 + 8 + 2 = 26），最后一次 `vpsadbw` 仍精确。
// 【第二处单变量】环容量 `XB` 64→128：环扫描在 u8 形态下每条向量 32 个条目，多出的
//   2 条向量比原来 16 位形态便宜得多，而 `rebuild_cp` 的调用次数由 n/XB 直接减半
//   （156→78 次）。判题机同 TU 交替 21 轮实测：u8 化 −12.80%，再叠 XB=128 −15.10%。
// 【正确性】本机 `-DBRUTE` 与 O(n^2) 暴力逐位对拍（n∈{200,600,1000,2000,5000,10000}
//   × mode∈{0..4}，bad=0）+ 与 #103819 输出 FNV 逐位相同。
// ================
// Approach:
// Keep the accepted Skylake engine and scan the first eight 16-bit vectors.
// Scan the upper eight only if the positional cutoff exceeds 128; then
// the remaining masks are otherwise all zero. This preserves exact counts.
// Purpose:
// Measure the two-stage scanv loop on the official 1009b judge.
#include <string.h>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC optimize("O3")
#pragma GCC target("avx2,popcnt,tune=skylake")
typedef unsigned u32;
typedef unsigned short u16;
typedef unsigned long long u64;
typedef unsigned char u8;

#define AVX2 __attribute__((target("avx2")))
#define BB 256
#define XB 128
#define BLKIDX(r) ((int)((r) / BB))
#define BLOFF(v)  ((int)((v) % BB))
#define BLK0(f)   ((f) * BB)
#ifndef MAXN
#define MAXN 100005
#endif
#define MAXG 640
#define SENTB 0x7F7F

static u16 H[MAXG * 320];
static u16 CP[MAXG * 320];
static u16 YZ[MAXN + 32];
static u8 ZY[MAXN + 32];
static u8 recY[XB + 32], recZ[XB + 32];
static signed char laneB[BB];
static u16 cnt[MAXN + 2], stv[MAXN + 2], curs[MAXN + 2];
static u32 ordx[MAXN + 2];
static u64 QA[MAXN + 2];
static u32 QB[MAXN + 2];
static u64 meta[MAXN + 2];
static u16 maskTab[17][16] __attribute__((aligned(32)));
static int gS;

AVX2 static inline __m256i ldu256(const void *p) { return _mm256_loadu_si256((const __m256i *)p); }
AVX2 static inline void stu256(void *p, __m256i v) { _mm256_storeu_si256((__m256i *)p, v); }
AVX2 static inline u32 hsum8(__m256i s) {
  __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);
}

AVX2 static inline __m256i scanv(const u16 *v, int off, int thr) {
  /* FIXED-TRIP-COUNT: always scan all BB/16 = 16 vectors, masking the lanes
     at/after `off` with a running "remaining" vector.  Bit-identical to the
     data-dependent form above, but the loop-exit branch is always predictable. */
  const __m256i t = _mm256_set1_epi16((short)thr);
  const __m256i lane = _mm256_setr_epi16(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
  const __m256i step = _mm256_set1_epi16(16);
  __m256i rem = _mm256_set1_epi16((short)off);
  __m256i acc = _mm256_setzero_si256();
  const u16 *p = v;
  for (int r = 0; r < 8; r++, p += 16) {
    __m256i c = _mm256_and_si256(_mm256_cmpgt_epi16(t, ldu256(p)),
                                 _mm256_cmpgt_epi16(rem, lane));
    acc = _mm256_sub_epi16(acc, c);
    rem = _mm256_sub_epi16(rem, step);
  }
  if (off > 128) for (int r = 8; r < 16; r++, p += 16) {
    __m256i c = _mm256_and_si256(_mm256_cmpgt_epi16(t, ldu256(p)),
                                 _mm256_cmpgt_epi16(rem, lane));
    acc = _mm256_sub_epi16(acc, c);
    rem = _mm256_sub_epi16(rem, step);
  }
  return acc;
}

/* u8 block-index scan: `value < thr` over 32-byte vectors.  thr is a block index
   (<= 78) and the sentinel is 0x7F, so signed vpcmpgtb is exact for the value
   test.  The positional cut is unsigned -- (off ^ 0x80) >signed (lane ^ 0x80) --
   so off may range over [0,256); laneB[] holds (32r + lane) ^ 0x80 as a signed
   byte and is built once per call.  The accumulator counts UP. */
AVX2 static inline __m256i scanb(const u8 *v, int off, int thr) {
  const __m256i t = _mm256_set1_epi8((char)thr);
  const __m256i k = _mm256_set1_epi8((char)(off ^ 0x80));
  __m256i acc = _mm256_setzero_si256();
  const u8 *p = v;
  for (int r = 0; r < (BB >> 5); r++, p += 32) {
    __m256i c = _mm256_and_si256(_mm256_cmpgt_epi8(t, ldu256(p)),
                                 _mm256_cmpgt_epi8(k, ldu256(laneB + 32 * r)));
    acc = _mm256_sub_epi8(acc, c);
  }
  return acc;
}

/* ring twin: no positional mask needed -- padrec() keeps the tail at 0x7F and
   0x7F fails both tests (fy <= 78 < 127). */
AVX2 static inline __m256i scanbrec(int fy, int fz) {
  const __m256i ty = _mm256_set1_epi8((char)fy);
  const __m256i tz = _mm256_set1_epi8((char)fz);
  __m256i acc = _mm256_setzero_si256();
  const u8 *py = recY, *pz = recZ;
  for (int r = 0; r < (XB >> 5); r++, py += 32, pz += 32)
    acc = _mm256_sub_epi8(acc, _mm256_and_si256(_mm256_cmpgt_epi8(ty, ldu256(py)),
                                               _mm256_cmpgt_epi8(tz, ldu256(pz))));
  return acc;
}

AVX2 static inline void padrec(void) {
  for (int k = 0; k < XB + 16; k++) { recY[k] = 0x7F; recZ[k] = 0x7F; }
}

template <int TS> AVX2 static void rebuild_cp(int G) {
  const int S = (TS > 0) ? TS : gS;
  const __m256i bcastLast = _mm256_setr_epi8(14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15,
                                             14,15,14,15,14,15,14,15,14,15,14,15,14,15,14,15);
  int fy = 0;
  for (; fy + 1 < G; fy += 2) {
    const u16 *h0 = H + (size_t)fy * S, *h1 = h0 + S;
    u16 *c0 = CP + (size_t)fy * S, *c1 = c0 + S;
    __m256i k0 = _mm256_setzero_si256(), k1 = _mm256_setzero_si256();
    for (int j = 0; j < S; j += 16) {
      __m256i v0 = ldu256(h0 + j), v1 = ldu256(h1 + j);
      v0 = _mm256_add_epi16(v0, _mm256_slli_si256(v0, 2));
      v1 = _mm256_add_epi16(v1, _mm256_slli_si256(v1, 2));
      v0 = _mm256_add_epi16(v0, _mm256_slli_si256(v0, 4));
      v1 = _mm256_add_epi16(v1, _mm256_slli_si256(v1, 4));
      v0 = _mm256_add_epi16(v0, _mm256_slli_si256(v0, 8));
      v1 = _mm256_add_epi16(v1, _mm256_slli_si256(v1, 8));
      v0 = _mm256_add_epi16(v0, _mm256_shuffle_epi8(_mm256_permute2x128_si256(v0, v0, 0x08), bcastLast));
      v1 = _mm256_add_epi16(v1, _mm256_shuffle_epi8(_mm256_permute2x128_si256(v1, v1, 0x08), bcastLast));
      __m256i r0 = _mm256_add_epi16(v0, k0), r1 = _mm256_add_epi16(v1, k1);
      stu256(c0 + j, r0);
      stu256(c1 + j, r1);
      k0 = _mm256_shuffle_epi8(_mm256_permute2x128_si256(r0, r0, 0x11), bcastLast);
      k1 = _mm256_shuffle_epi8(_mm256_permute2x128_si256(r1, r1, 0x11), bcastLast);
    }
  }
  for (; fy < G; fy++) {
    const u16 *h = H + (size_t)fy * S;
    u16 *c = CP + (size_t)fy * S;
    __m256i k = _mm256_setzero_si256();
    for (int j = 0; j < S; j += 16) {
      __m256i v = ldu256(h + j);
      v = _mm256_add_epi16(v, _mm256_slli_si256(v, 2));
      v = _mm256_add_epi16(v, _mm256_slli_si256(v, 4));
      v = _mm256_add_epi16(v, _mm256_slli_si256(v, 8));
      v = _mm256_add_epi16(v, _mm256_shuffle_epi8(_mm256_permute2x128_si256(v, v, 0x08), bcastLast));
      __m256i r = _mm256_add_epi16(v, k);
      stu256(c + j, r);
      k = _mm256_shuffle_epi8(_mm256_permute2x128_si256(r, r, 0x11), bcastLast);
    }
  }
  for (int f = 1; f < G; f++) {
    const u16 *p = CP + (size_t)(f - 1) * S;
    u16 *c = CP + (size_t)f * S;
    for (int j = 0; j < S; j += 16)
      stu256(c + j, _mm256_add_epi16(ldu256(c + j), ldu256(p + j)));
  }
}

template <int CG, int CS> static void engine_t(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  constexpr int TSC = CS;
  if (n <= 1) { if (n == 1) out[0] = 0; return; }
  const int G = (CG > 0) ? CG : ((n + BB - 1) / BB);
  const int S = (CS > 0) ? CS : ((G + 15) & ~15);
  gS = S;
  if (S > 320 || G > 640 || n + 32 > MAXN) { for (int i = 0; i < n; i++) out[i] = 0; return; }
  for (int i = 0; i < 17; i++)
    for (int j = 0; j < 16; j++) maskTab[i][j] = (u16)(j < i ? 0xFFFFu : 0u);
  for (int i = 0; i < BB; i++) laneB[i] = (signed char)(i ^ 0x80);

  /* ---- y and z ranks in one fused pass ---- */
  memset(cnt, 0, (size_t)(n + 1) * 2);
  for (int i = 0; i < n; i++) cnt[y[i]]++;
  { u32 a = 0; for (int v = 0; v <= n; v++) { u32 c = cnt[v]; stv[v] = (u16)a; curs[v] = (u16)a; a += c; } }
  for (int i = 0; i < n; i++) { u32 v = y[i]; meta[i] = ((u64)stv[v] << 48) | (u64)curs[v]++; }

  memset(cnt, 0, (size_t)(n + 1) * 2);
  for (int i = 0; i < n; i++) cnt[z[i]]++;
  { u32 a = 0; for (int v = 0; v <= n; v++) { u32 c = cnt[v]; stv[v] = (u16)a; curs[v] = (u16)a; a += c; } }
  for (int i = 0; i < n; i++) { u32 v = z[i]; meta[i] |= ((u64)stv[v] << 32) | ((u64)curs[v]++ << 16); }

  /* ---- x order ---- */
  memset(cnt, 0, (size_t)(n + 1) * 2);
  for (int i = 0; i < n; i++) cnt[x[i]]++;
  { u32 a = 0; for (int v = 0; v <= n; v++) { u32 c = cnt[v]; stv[v] = (u16)a; curs[v] = (u16)a; a += c; } }
  for (int i = 0; i < n; i++) { u32 v = x[i]; ordx[curs[v]++] = i; }
  for (int e = 0; e < n; e++) {
    const int i = ordx[e];
    const u64 m = meta[i];
    unsigned xv = x[i];
    if (xv > 0x3FFF) xv = 0x3FFF;
    QA[e] = ((m >> 48) << 42) | ((((m >> 32) & 0xFFFF)) << 28) | ((u64)i << 14) | (u64)xv;
    QB[e] = (u32)((m & 0xFFFF) | (((m >> 16) & 0xFFFF) << 14));
  }

  /* ---- structures (H/YZ/ZY are updated at INSERT time; CP is a snapshot) ---- */
  memset(H, 0, (size_t)G * S * 2);
  memset(CP, 0, (size_t)G * S * 2);
  memset(YZ, SENTB, (size_t)n * 2 + 32);
  memset(ZY, 0x7F, (size_t)n + 32);
  padrec();
  int nrec = 0;

  int p = 0;
  while (p < n) {
    int q = p + 1;
    const u32 xv = (u32)(QA[p] & 0x3FFF);
    while (q < n && ((u32)(QA[q] & 0x3FFF)) == xv) q++;
    for (int e = p; e < q; e++) {
      const u64 m = QA[e];
      const u32 ky = (u32)(m >> 42), kz = (u32)((m >> 28) & 0x3FFF);
      const int FY = (int)BLKIDX(ky), FZ = (int)BLKIDX(kz);
      u32 c = (FY > 0 && FZ > 0) ? CP[(size_t)(FY - 1) * S + (FZ - 1)] : 0;
      __m256i edge1 = scanv(YZ + ((size_t)BLK0(FY)), (int)BLOFF(ky), (int)kz);
      __m256i edge2 = scanb(ZY + (size_t)BLK0(FZ), (int)BLOFF(kz), (int)FY);
      __m256i recent = scanbrec((int)FY, (int)FZ);
      __m256i all = _mm256_add_epi8(_mm256_add_epi8(edge1, edge2), recent);
      __m256i sums = _mm256_sad_epu8(all, _mm256_setzero_si256());
      __m128i halves = _mm_add_epi64(_mm256_castsi256_si128(sums), _mm256_extracti128_si256(sums, 1));
      c += (u32)_mm_cvtsi128_si64(halves) + (u32)_mm_extract_epi64(halves, 1);
      out[(m >> 14) & 0x3FFF] = c;
    }
    for (int e = p; e < q; e++) {
      const u32 b = QB[e];
      const u32 ry = b & 0x3FFF, rz = b >> 14;
      H[(size_t)BLKIDX(ry) * S + BLKIDX(rz)]++;
      YZ[ry] = (u16)rz;
      ZY[rz] = (u8)BLKIDX(ry);
      recY[nrec] = (u8)BLKIDX(ry);
      recZ[nrec] = (u8)BLKIDX(rz);
      if (++nrec == XB) { rebuild_cp<TSC>(G); nrec = 0; padrec(); }
    }
    p = q;
  }
}
#pragma GCC pop_options


void count_3d(int n, const unsigned *x, const unsigned *y, const unsigned *z, unsigned *out) {
  if (n == 10000) engine_t<40, 48>(n, x, y, z, out);
  else            engine_t<-1, -1>(n, x, y, z, out);
}

评测结果