代码

/* References:
- saffah_cc_v41_agg1 https://duck.ac/submission/126801: reuse the accepted
  complete 2D rank engine and its safe staging-buffer alias. Inherited credits
  remain below; no separate license notice accompanies the public submission.
- saffah_cc_v41_agg1 https://duck.ac/submission/122619 (also preserved in
  https://duck.ac/submission/126794): adapt the idea of MADV_HUGEPAGE advice
  for scratch/output storage, with attribution to that public implementation.
*/
// ===== REFERENCES (r2d_ 席, 2026-10-15) =====
// [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#126306** <https://duck.ac/submission/126306>
//     （判题机 181214.045 µs = 本发入场时的当刻 T）—— **本文件正文（引擎部分）为该提交的逐字节复刻**
//     （`§2.18.349` 第一步"抄满"取回 + 剥注释核对；取回件 `problems/1008/ref2/rival_126306.cpp`）。
//     沿用内容：整套「x 分桶 MSD-11（`x>>13`）+ 桶内 13 位单趟计数排序（u16 计数 + 16 路向量前缀进位）
//     + y 直方图低/高字节拆分（`gHistLow`/`gHistHigh`，只进低字节、进位时才写高字节，再向量加宽重建）
//     + L1 段 9 位计数器树 `gQ`（含 `gInc` 增量表与 `cnt_q2/cnt_ins2` 两级特化）
//     + 融合 L2+L3（6 位数字 `r64` 掩码行树 / 9 位数字 `rank8`/`rank16` 掩码行树）
//     + 6 字节 staging（`gBOff` 按 idx 高位分桶）+ flush 字节指针遍历」引擎。
//     它自带的全部 References / Idea / Purpose 注释块**原样保留在本文件头下方**，未作删改。
// [2] 本账号（saffah_cc_v41_agg1）**#122394** <https://duck.ac/submission/122394>（183.221181 ms）
//     —— [1] 的继承链上游（[1] 的注释自述 "saffah_cc_v41_agg1, .../122394 : inherited compact
//     rank-loop fields"），即本席接手时我方在册最好件（`work/o08a_sub1.cpp`）；本发**未复制**其任何代码。
// [3] 本账号更早的自造刀（均在 [1] 的注释块里被具名继承）：#109641（计数器树单基址 `gQ`）·
//     #121277 / #122394（`cursor+3` / `cursor+4` 免掩码直接载入）。
// 许可证合规：duck.ac 公开代码按站规公开可读；[1] 的公开源**未附独立许可证声明**（其注释自述
//   "no separate license notice"），各件仅在同一竞赛站内复刻与再创作，不作他用。
// ======================
// ===== 思路 =====
// [r2d-P1]（r2d_ 席）把 60 MB 的 staging 缓冲，叠到 stage A 已触碰过、fused 段已死的 x 窗口区。
// 【口径】入场 `mine=183.221181`（#122394）· `T=181.214045`（#126306，**晚于**我方 ⇒ 宽支 1.005T=182.121115）。
//   本发落地后 T 的持有提交将早于本件 ⇒ 判据翻**严支** `0.99*T+1µs = 179.402905 ms`。
//
// 【取证 1 ★判题机 "Memory = 碰过的页" 的直接定价（本席用 pretouch A/B 做，非引用他人结论）】
//   同一份引擎、同一个 probe main，只在 `count_2d` **计时区之外**加一趟"把 PKA(80MB)+PKB(窗口区 144MB)+
//   g_STG(60MB) 每页各碰一次"的预触碰循环：
//     | 臂 | probe ticks | 说明 |
//     | pretouch=0 | 629,444,724 | 计时区内发生全部首触页错 |
//     | pretouch=1 | 593,057,498 | 首触页错全部移到计时区外 |
//   ⇒ **差 36.4 M ticks = 10.4 ms = 5.8%**（这批页 ≈ 55 k 页 ⇒ **≈ 663 ticks/页 ≈ 0.19 µs/页**，
//     与协调者给的 `mem/4 × 0.244 µs` 同量级、符号一致）。
//   旁证（同一 tail）：普通 probe tail（不在计时区外碰 `out[]`）636.49 M ↔ 加了 40 MB `memset(o)` 的
//   tail 629.44 M ⇒ **40 MB = 9765 页 = 7.05 M ticks**（722 ticks/页）✓ 模型自洽。
//
// 【本刀的机制】全题的三块"我方自有"大数组里，`g_STG`（6 B/记录 × n = 60 MB = 14648 页）**只在 fused
//   段被写、在 flush 段被读**；而 fused 段开始时，PKB 的 x 桶窗口区（stage A 已经全部触碰过、L1 段读完后
//   全部死亡）还完整驻留。⇒ 把 staging 基址从独立的 `g_STG` 换成 `PKB + 1 MB`：
//   * 这 60 MB 的页**已经被 stage A 触碰过** ⇒ 不再产生任何首触页错（只余窗口末尾 ~6% 的空隙页）；
//   * `PKB + 1 MB` 与 fused 段 L2 散射写的 `PKB[0 .. 最大 y 桶)` 不重叠（随机数据最大桶 ≈ 32768 条 = 256 KB；
//     守卫阈值取 `STGPKOFF/8 = 131072` 条 = 1 MB，超限**自动回落**到独立 `g_STG_std` ⇒ 语义与旧码逐位相同）；
//   * 排序暂存区在 `PKB + MSN*XBSTEP`（144 MB 处，只用前 65 KB），与 1 MB 处的 staging 互不干涉。
//   纯地址改动：**不改变任何算法、任何记录布局、任何计数顺序**。
//
// 【剂量（判题机"立即体验"，`work/r2d__tmp/r2d_ab.sh`，n=1e7 生产几何，交替各 2 轮）】
//     | 臂 | ticks（两轮） | 相对 |
//     | #126306 正文（`r2d_eng.cpp`） | 636,490,184 / 636,487,058 | — |
//     | **本件（`r2d_alias1.cpp`）** | **629,226,024 / 629,223,544** | **−1.14%**（−7.26 M ticks）|
//   * 两轮同名件差 ≤ 0.0004%，`h=25293157591` **逐位相同** ✓（n=1e7 生产输出未变）。
// 【闸门】`work/w8s_gate.py r2d_eng.cpp r2d_alias1.cpp`：**15 configs, 0 failing**（n=3000/20000 六形态与
//   O(n²) 暴力对拍 BRUTE-OK；n=1e5/1e6/1e7 × 三形态 FNV 逐位相同）✓ ·
//   `ref/gcc9`（**判题机同款 gcc 9.3**）`-O2 -std=c++17` 编译 0 错 ✓ · md5 与 `work/` 下所有 .cpp 均不同 ✓。
// 【目的（实验性提交声明）】"页错预算"轴在本题的第一次定价与第一次收割；若判题机判负，榜单只保留 #122394，无下行 ✓
// ================
/* Idea: extend huge-page advice to this engine's large static scratch arrays
and caller output. Restrict each request to whole pages contained in the given
buffer and keep all counting, allocation and staging behavior unchanged.
Purpose: experimental measurement on duck.ac of reduced first-touch page and
address-translation work for the existing static-buffer implementation.
*/
#include <sys/mman.h>
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <immintrin.h>

#pragma GCC push_options
#pragma GCC optimize("O3,unroll-loops,rename-registers")
#pragma GCC target("avx2,popcnt,bmi2")

typedef uint64_t u64;
typedef uint32_t u32;

#define MSH 11
#define MSN (1 << MSH)
#define XSH (24 - MSH)
static const int NMAX = 10000002;

alignas(64) static u32 XA[NMAX], XB[NMAX];
#define XBCAP 9216u
#define XBSTEP (XBCAP + 8u)
alignas(64) static u64 PKA[NMAX], PKB[(size_t)MSN * XBSTEP + (size_t)NMAX + 64];
alignas(64) static u32 DUPA[NMAX];
static u32 cntb[MSN], xend[MSN];
static u32 cnb2[128];
alignas(32) static uint16_t local_count14[1u << XSH];
static void local_prefix14() {
  __m256i carry=_mm256_setzero_si256();
  const __m256i last128=_mm256_set1_epi16(0x0f0e);
  for(u32 i=0;i<(1u<<XSH);i+=16){
    const __m256i v=_mm256_load_si256((const __m256i*)(local_count14+i));
    __m256i p=_mm256_add_epi16(v,_mm256_slli_si256(v,2));
    p=_mm256_add_epi16(p,_mm256_slli_si256(p,4));
    p=_mm256_add_epi16(p,_mm256_slli_si256(p,8));
    const __m256i half_totals=_mm256_shuffle_epi8(p,last128);
    const __m256i cross=_mm256_permute2x128_si256(half_totals,half_totals,0x08);
    p=_mm256_add_epi16(p,cross);
    _mm256_store_si256((__m256i*)(local_count14+i),_mm256_add_epi16(_mm256_sub_epi16(p,v),carry));
    const __m256i total=_mm256_add_epi16(half_totals,_mm256_permute2x128_si256(half_totals,half_totals,0x01));
    carry=_mm256_add_epi16(carry,total);
  }
}

static u32 gH512[512];
#define NBKT 512
#define BSH 15
#define STGPADE 32   /* [q8z-STG6] staging-row pad in RECORDS: 21*6 = 126 B ~= the old 16*8 = 128 B */
alignas(4096) static unsigned char g_STG_std[6 * (size_t)NMAX + (size_t)NBKT * STGPADE * 6 + 128];   /* [q8z-STG6] 6 B/record (fallback) */
/* [r2d-P1] The staging buffer is 6 B/record = 60 MB of first-touch page faults
   (~14.6 k pages ~ 2.8 ms).  It is only live in the fused phase, which starts
   AFTER the x-bucket window region of PKB is dead, so it aliases that region:
   those pages were faulted by stage A already.  STGPKOFF stays clear of the
   fused L2 scatter (which writes PKB[0..max bucket)). */
#define STGPKOFF ((size_t)1 << 20)
#define STG6EP(stg,k) ((unsigned char *)(stg) + (size_t)(k))
static u32 gBkt[NBKT + 1];
static u32 gBOff[NBKT + 1];

/* [w8x-Q1] The 9-bit counter tree AND its increment masks live in ONE array so the
   hot loops need a single base register.  Layout (u32 units) / byte offsets:
     gQ[  0..  7]  P1  level 1, 8 lanes           byte    0
     gQ[  8.. 71]  P2  level 2, 8 rows x 8        byte   32
     gQ[ 72..583]  P3  level 3, 64 rows x 8       byte  288
     gQ[584..647]  INC increment masks (INC[c][i]=(i>c))  byte 2336
   With d = 9-bit digit, a=d>>6, b=(d>>3)&7, c=d&7:
     gP1[a]       == gQ[a]            &gP1[0]      == gQ
     gP2[a][b]    == gQ[8  + (d>>3)]  &gP2[a][0]  == gQ + 8  + a*8
     gP3[a][b][c] == gQ[72 + d]       &gP3[a][b][0] == gQ + 72 + (d>>3)*8      */
#define QINCB 584u
alignas(32) static u32 gQ[648];
alignas(16) static uint16_t gR64Q[136];
#define gR64H (gR64Q)
#define gR64L ((uint16_t (*)[8])(gR64Q + 8))
#define gR64Mask ((uint16_t (*)[8])(gR64Q + 72))
static bool gR64MasksReady = false;
static inline void r64_init() {
  if (!gR64MasksReady) {
    for (int d=0; d<8; d++) for (int j=0; j<8; j++)
      gR64Mask[d][j] = (uint16_t)(j>d);
    gR64MasksReady = true;
  }
  memset(gR64H, 0, 16);
  memset(gR64L, 0, 128);
}
static inline u32 r64_q(u32 d) {
  return (u32)gR64H[d>>3] + gR64L[d>>3][d&7];
}
static inline void r64_ins(u32 d) {
  u32 hi=d>>3, lo=d&7;
  __m128i h=_mm_load_si128((const __m128i*)gR64H);
  __m128i hm=_mm_load_si128((const __m128i*)gR64Mask[hi]);
  _mm_store_si128((__m128i*)gR64H, _mm_add_epi16(h,hm));
  __m128i *row=(__m128i*)gR64L[hi];
  __m128i l=_mm_load_si128(row);
  __m128i lm=_mm_load_si128((const __m128i*)gR64Mask[lo]);
  _mm_store_si128(row, _mm_add_epi16(l,lm));
}

alignas(32) static u32 gMsk[8][8];
#define gInc(i) (gQ + QINCB + (i) * 8u)   /* [w8x-Q1] increment masks now live in gQ */
static u32 gCur[512];

static u32 gBnd[514];
static u32 gHis[513];
alignas(64) static uint16_t gHis16[512][64];
alignas(64) static uint8_t gHistLow[32768],gHistHigh[32768];
static u32 gHis32[512][64];

static inline void cnt_init(void) { memset(gQ, 0, QINCB * sizeof(u32)); }
static inline u32 cnt_q(u32 d) {
  return gQ[d >> 6] + gQ[8u + (d >> 3)] + gQ[72u + d];
}
static inline void cnt_ins(u32 d) {
  u32 a = d >> 6, t = d >> 3, c = d & 7;
  __m256i *p3 = (__m256i *)(gQ + 72u + (t << 3));
  __m256i *p2 = (__m256i *)(gQ + 8u + (a << 3));
  __m256i *p1 = (__m256i *)gQ;
  __m256i i3 = _mm256_load_si256((const __m256i *)gInc(c));
  __m256i v3 = _mm256_load_si256(p3);
  _mm256_store_si256(p3, _mm256_add_epi32(v3, i3));
  __m256i i2 = _mm256_load_si256((const __m256i *)gInc(t & 7u));
  __m256i v2 = _mm256_load_si256(p2);
  _mm256_store_si256(p2, _mm256_add_epi32(v2, i2));
  __m256i i1 = _mm256_load_si256((const __m256i *)gInc(a));
  __m256i v1 = _mm256_load_si256(p1);
  _mm256_store_si256(p1, _mm256_add_epi32(v1, i1));
}

static inline u32 cnt_q2(u32 d) {
  u32 b = (d >> 3) & 7;
  return gQ[8u + b] + gQ[72u + (b << 3) + (d & 7u)];
}
static inline void cnt_ins2(u32 d) {
  u32 b = (d >> 3) & 7, c = d & 7;
  __m256i *p3 = (__m256i *)(gQ + 72u + (b << 3));
  __m256i *p2 = (__m256i *)(gQ + 8u);
  __m256i i3 = _mm256_load_si256((const __m256i *)gInc(c));
  __m256i v3 = _mm256_load_si256(p3);
  _mm256_store_si256(p3, _mm256_add_epi32(v3, i3));
  __m256i i2 = _mm256_load_si256((const __m256i *)gInc(b));
  __m256i v2 = _mm256_load_si256(p2);
  _mm256_store_si256(p2, _mm256_add_epi32(v2, i2));
}
static inline void cnt_init2(void) { memset(gQ + 8u, 0, 8u * sizeof(u32)); memset(gQ + 72u, 0, 64u * sizeof(u32)); }


alignas(32) static uint16_t gR16H[16], gR16L[16][32];
alignas(32) static uint16_t gR16MaskH[16][16], gR16MaskL[32][32];
static int gR16MasksReady = 0;

static inline void rank16_init() {
  if (!gR16MasksReady) {
    for (int a=0;a<16;a++) for (int i=0;i<16;i++) gR16MaskH[a][i]=(uint16_t)(i>a);
    for (int a=0;a<32;a++) for (int i=0;i<32;i++) gR16MaskL[a][i]=(uint16_t)(i>a);
    gR16MasksReady=1;
  }
  memset(gR16H,0,sizeof(gR16H));
  memset(gR16L,0,sizeof(gR16L));
}
static inline u32 rank16_q(u32 d) {
  return (u32)gR16H[d>>5] + gR16L[d>>5][d&31];
}
static inline void rank16_ins(u32 d) {
  u32 hi=d>>5,lo=d&31;
  __m256i h=_mm256_load_si256((const __m256i*)gR16H);
  __m256i hm=_mm256_load_si256((const __m256i*)gR16MaskH[hi]);
  _mm256_store_si256((__m256i*)gR16H,_mm256_add_epi16(h,hm));
  __m256i *row=(__m256i*)gR16L[hi];
  const __m256i *mask=(const __m256i*)gR16MaskL[lo];
  __m256i a=_mm256_load_si256(row), b=_mm256_load_si256(row+1);
  _mm256_store_si256(row,_mm256_add_epi16(a,_mm256_load_si256(mask)));
  _mm256_store_si256(row+1,_mm256_add_epi16(b,_mm256_load_si256(mask+1)));
}

alignas(32) static uint8_t gR8L[16][32], gR8MaskL[32][32];
static int gR8MaskReady=0;
static inline void rank8_init() {
  if (!gR16MasksReady) {
    for (int a=0;a<16;a++) for (int i=0;i<16;i++) gR16MaskH[a][i]=(uint16_t)(i>a);
    for (int a=0;a<32;a++) for (int i=0;i<32;i++) gR16MaskL[a][i]=(uint16_t)(i>a);
    gR16MasksReady=1;
  }
  if (!gR8MaskReady) {
    for (int a=0;a<32;a++) for (int i=0;i<32;i++) gR8MaskL[a][i]=(uint8_t)(i>a);
    gR8MaskReady=1;
  }
  const __m256i z=_mm256_setzero_si256();
  _mm256_store_si256((__m256i*)gR16H,z);
  #define CLEAR_R8(k) _mm256_store_si256((__m256i*)((u32*)gR8L + 8*(k)),z)
  CLEAR_R8(0); CLEAR_R8(1); CLEAR_R8(2); CLEAR_R8(3);
  CLEAR_R8(4); CLEAR_R8(5); CLEAR_R8(6); CLEAR_R8(7);
  CLEAR_R8(8); CLEAR_R8(9); CLEAR_R8(10); CLEAR_R8(11);
  CLEAR_R8(12); CLEAR_R8(13); CLEAR_R8(14); CLEAR_R8(15);
  #undef CLEAR_R8
}
static inline u32 rank8_q(u32 d) {
  return (u32)gR16H[d>>5] + gR8L[d>>5][d&31];
}
static inline void rank8_ins(u32 d) {
  u32 hi=d>>5,lo=d&31;
  __m256i h=_mm256_load_si256((const __m256i*)gR16H);
  __m256i hm=_mm256_load_si256((const __m256i*)gR16MaskH[hi]);
  _mm256_store_si256((__m256i*)gR16H,_mm256_add_epi16(h,hm));
  __m256i *row=(__m256i*)gR8L[hi];
  const __m256i *mask=(const __m256i*)gR8MaskL[lo];
  __m256i v=_mm256_load_si256(row);
  _mm256_store_si256(row,_mm256_add_epi8(v,_mm256_load_si256(mask)));
}
static inline void rank8_upgrade() {
  for (int a=0;a<16;a++)
    for (int b=0;b<32;b++) gR16L[a][b]=(uint16_t)gR8L[a][b];
}
static inline bool rank8_fits(const u64 *pk,u32 p,u32 qq) {
  uint16_t groups[16]={};
  for (u32 a=p;a<qq;a++) groups[(u32)(pk[a]>>29)&15]++;
  for (int h=0;h<16;h++) if (groups[h]>255) return false;
  return true;
}
template<int MODE>
static inline void rank_l3_run(u64 *pk,u32 p,u32 qq,unsigned char *stg) {
  if constexpr (MODE==2) rank8_init();
  else if constexpr (MODE==1) rank16_init();
  else cnt_init();
  for (u32 a=p;a<qq;a++) {
    u64 rec=pk[a];
    u32 da=(u32)(rec>>24)&511;
    u32 inc;
    if constexpr (MODE==2) { inc=rank8_q(da); rank8_ins(da); }
    else if constexpr (MODE==1) { inc=rank16_q(da); rank16_ins(da); }
    else { inc=cnt_q(da); cnt_ins(da); }
    u32 ix=(u32)(rec&0xFFFFFFu);
    { u32 _c=(u32)(rec>>40)+inc; unsigned char *_p=STG6EP(stg,((gBOff[ix>>BSH] += 6u) - 6u));
      *(u64*)_p=(((u64)_c)<<24)|ix; }
  }
}

__attribute__((aligned(64)))
__attribute__((noinline,cold)) static u32 dup_long(const u64 *pk, u32 t, u32 a, u32 ya0) {
  u32 dup = 0;
  for (u32 b = t - a; b < t; b++) {
    u64 r = pk[b];
    __asm__("" : "+r"(r));
    dup += (u32)(((u32)(r >> 24)) < ya0);
  }
  return dup;
}
static inline __attribute__((always_inline)) void flush_rank_output(
    u32 *__restrict out, const unsigned char *__restrict stage) {
  for(u32 b=0;b<NBKT;++b){
    const unsigned char *const end=stage+gBOff[b];
    const unsigned char *record=end-6u*gBkt[b];
    if(gBkt[b]>48u){
      const unsigned char *const prefetched_end=end-6u*48u;
      for(;record<prefetched_end;record+=6){
        u64 _pf; __builtin_memcpy(&_pf,record+6u*48u,8);
        __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(&out[(u32)(_pf & 0xFFFFFFu)])));
        u32 index_bytes,count_bytes;
        __builtin_memcpy(&index_bytes,record,4);
        __builtin_memcpy(&count_bytes,record+2,4);
        out[index_bytes&0xFFFFFFu]=count_bytes>>8;
      }
    }
    for(;record<end;record+=6){
      u32 index_bytes,count_bytes;
      __builtin_memcpy(&index_bytes,record,4);
      __builtin_memcpy(&count_bytes,record+2,4);
      out[index_bytes&0xFFFFFFu]=count_bytes>>8;
    }
  }
}

static inline void advise_rank_pages(void *p,size_t bytes) {
  uintptr_t lo=((uintptr_t)p+4095u)&~(uintptr_t)4095u;
  uintptr_t hi=((uintptr_t)p+bytes)&~(uintptr_t)4095u;
  if(hi>lo) (void)madvise((void*)lo,hi-lo,MADV_HUGEPAGE);
}
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;
  }
  if (n > 16000000) {           // record fields: idx < 2^24, y < 2^24
    memset(out, 0, (size_t)n * 4);
    return;
  }
  static int inited = 0;
  if (!inited) {
    for (int m = 0; m < 8; m++)
      for (int i2 = 0; i2 < 8; i2++) { gMsk[m][i2] = (i2 > m) ? 0xFFFFFFFFu : 0u;
                                       gInc(m)[i2] = (i2 > m) ? 1u : 0u; }

    inited = 1;
  }
  u32 np = (u32)n;
  advise_rank_pages(PKA,sizeof(PKA));
  advise_rank_pages(PKB,sizeof(PKB));
  advise_rank_pages(out,(size_t)n*4);


  // ------------- sort by x: MSD 8 bits (x>>16) then two BUCKET-LOCAL LSD ----
  // The old form was 3 whole-array LSD passes = 680 MB of DRAM traffic.  Here the
  // two low passes run inside one 256-bucket window (~470 KB), so the second pass
  // reads the bucket from L3 and the write-back hits the line already resident:
  // DRAM drops to stage A (200 MB) + one read + one write-back per bucket (240 MB).
  // Stage A scatters straight from the caller's x/y, so no staging copy is needed.
  // XA / PKA are FREE during stage B and are used as the per-bucket temp (len <= np).
  static u32 bnd[MSN + 1];
  memset(cntb, 0, sizeof(cntb));
  memset(gHistLow,0,sizeof(gHistLow));memset(gHistHigh,0,sizeof(gHistHigh));
  { u32 _full = np >> BSH, _rem = np & ((1u << BSH) - 1u);
    for (u32 _b = 0; _b < NBKT; _b++)
      gBkt[_b] = (_b < _full) ? (1u << BSH) : (_b == _full ? _rem : 0u); }
  for (int d = 0; d < MSN; d++) { bnd[d] = (u32)d * XBSTEP; cntb[d] = bnd[d]; }
  bnd[MSN] = (u32)MSN * XBSTEP;
  u32 i = 0;
  for (; i + 10u < np; i++) {
    u32 future_dg = x[i + 10u] >> XSH;
    __asm__ __volatile__("prefetchw %0" :: "m"(*(const char *)(PKB + cntb[future_dg])));
    u32 xv = x[i];
    u32 dg = (xv >> XSH);
    u32 qq = cntb[dg]++;
    PKB[qq] = ((u64)(xv & 0xFFFFu) << 48) | ((u64)y[i] << 24) | i;
    { u32 h = (y[i] >> 9)   /* [w8s-K10] mask redundant: n <= 1.6e7 => y>>9 < 2^15 */;
      if(__builtin_expect(++gHistLow[h]==0,0))++gHistHigh[h]; }
  }
  for (; i < np; i++) {
    u32 xv = x[i];
    u32 dg = (xv >> XSH);
    u32 qq = cntb[dg]++;
    PKB[qq] = ((u64)(xv & 0xFFFFu) << 48) | ((u64)y[i] << 24) | i;
    { u32 h = (y[i] >> 9)   /* [w8s-K10] mask redundant: n <= 1.6e7 => y>>9 < 2^15 */;
      if(__builtin_expect(++gHistLow[h]==0,0))++gHistHigh[h]; }
  }
  for (int d = 0; d < MSN; d++) xend[d] = cntb[d];
  { u32 _bad = 0;
    for (int d = 0; d < MSN; d++) if (xend[d] - bnd[d] > XBCAP) _bad = 1;
    if (_bad) {
      memset(cntb, 0, sizeof(u32) * MSN);
      for (u32 i = 0; i < np; i++) cntb[(x[i] >> XSH)]++;
      { u32 acc = 0;
        for (int d = 0; d < MSN; d++) { u32 c = cntb[d]; bnd[d] = acc; cntb[d] = acc; xend[d] = acc + c; acc += c; } }
      bnd[MSN] = np;
      memset(gHistLow,0,sizeof(gHistLow));memset(gHistHigh,0,sizeof(gHistHigh));
      for (u32 i = 0; i < np; i++) {
        u32 dg = (x[i] >> XSH); u32 qq = cntb[dg]++;
        PKB[qq] = ((u64)(x[i] & 0xFFFFu) << 48) | ((u64)y[i] << 24) | i;
        { u32 h = (y[i] >> 9)   /* [w8s-K10] mask redundant: n <= 1.6e7 => y>>9 < 2^15 */; if(__builtin_expect(++gHistLow[h]==0,0))++gHistHigh[h]; }
      }
    } }
  for(u32 h=0;h<32768;h+=16){
    const __m256i low=_mm256_cvtepu8_epi16(_mm_load_si128((const __m128i*)(gHistLow+h)));
    const __m256i high=_mm256_cvtepu8_epi16(_mm_load_si128((const __m128i*)(gHistHigh+h)));
    _mm256_store_si256((__m256i*)((uint16_t*)gHis16+h),_mm256_or_si256(low,_mm256_slli_epi16(high,8)));
  }
  u32 histogram_total = 0;
  for (u32 d1 = 0; d1 < 512; d1++) {
    u32 total = 0;
    for (u32 d2 = 0; d2 < 64; d2++) total += gHis16[d1][d2];
    gH512[d1] = total;
    histogram_total += total;
  }
  const bool need32 = (histogram_total != np);
  if (need32) {
    memset(gHis32, 0, sizeof(gHis32));
    for (u32 i = 0; i < np; i++) ++gHis32[(y[i] >> 15) & 511][(y[i] >> 9) & 63];
    for (u32 d1 = 0; d1 < 512; d1++) {
      u32 total = 0;
      for (u32 d2 = 0; d2 < 64; d2++) total += gHis32[d1][d2];
      gH512[d1] = total;
    }
  }
  u64 *pk = PKB;
  u64 *alt = PKA;
  { u32 acc = 0; for (int b = 0; b < NBKT; b++) { gBOff[b] = acc; acc += 6u * (gBkt[b] + STGPADE); } gBOff[NBKT] = acc; }

  /* [r2d-P1] pick the staging base: PKB window region when the largest y-bucket
     fits below STGPKOFF, else the standalone array. */
  unsigned char *gSTG = g_STG_std;
  { u32 acc = 0;
    for (int d = 0; d < 512; d++) { u32 c = gH512[d]; gBnd[d] = acc; gCur[d] = acc; acc += c; }
    acc = 0;
    for (int d = 0; d < 512; d++) { u32 c = gH512[d]; if (c > acc) acc = c; }
    if (acc <= (u32)(STGPKOFF / 8u)) gSTG = (unsigned char *)PKB + STGPKOFF;
  }
  // ---------------- L1 bucket offsets (histogram filled in stage A) -------
  gBnd[512] = np;

  // ---------------- L1: MSD 8 bits (y bits 23..16) -> alt ----------------
  cnt_init();
  {
    for (u32 bb = 0; bb < MSN; bb++) {
    const u32 local_s=bnd[bb], local_e=xend[bb], local_len=local_e-local_s;
    u64 *sorted_pk=PKB;
    if (local_len>1 && local_len<=65535u) {
      memset(local_count14,0,sizeof(local_count14));
      for (u32 t=local_s;t<local_e;++t)
        ++local_count14[(u32)(PKB[t]>>48)&((1u<<XSH)-1u)];
      local_prefix14();
      u64 *const sorted=PKB+(size_t)MSN*XBSTEP;
      for (u32 t=local_s;t<local_e;++t) {
        const u64 r=PKB[t];
        const u32 d=(u32)(r>>48)&((1u<<XSH)-1u);
        sorted[local_count14[d]++]=r;
      }
      sorted_pk=sorted-local_s;
    } else if (local_len>1)
  {
    const u32 b=bb;
    u64 *const sort_tmp = PKB+(size_t)MSN*XBSTEP;
    u32 s = bnd[b], e = xend[b], len = e - s;
    if (len > 1) {
    memset(cntb, 0, 256 * sizeof(u32));
    memset(cnb2, 0, sizeof(cnb2));
    for (u32 t = s; t < e; t++) cntb[(u32)(PKB[t] >> 48) & 127]++;
    { u32 a2 = 0;
      for (int d = 0; d < 128; d++) { u32 c = cntb[d]; cntb[d] = a2; a2 += c; } }
    for (u32 t = s; t < e; t++) {
      u64 r = PKB[t];
      u32 dg = (u32)(r >> 48) & 127;
      u32 qq = cntb[dg]++;
      sort_tmp[qq] = r;
      cnb2[(u32)(r >> 55) & 127]++;     /* == the histogram pass 3 used to rebuild */
    }
    { u32 a2 = s;
      for (int d = 0; d < 128; d++) { u32 c = cnb2[d]; cnb2[d] = a2; a2 += c; } }
    for (u32 t = 0; t < len; t++) {
      u32 dg = (u32)(sort_tmp[t] >> 55) & 127;
      PKB[cnb2[dg]++] = sort_tmp[t];
    }
  }
  }

    u32 t = bnd[bb], bnde = xend[bb];
    u32 a=0;
    __m128i prevxv=_mm_set1_epi32(-1);
    __m256i yhist = _mm256_setzero_si256();
    u64 *cursor=sorted_pk+t,*const cursor_end=sorted_pk+bnde;
    for (; cursor < cursor_end; ++cursor) {
      u64 rec0 = *cursor;
      const __m128i xv=_mm_srli_epi64(_mm_cvtsi64_si128((long long)rec0),48);
      const u32 equal_mask=(u32)_mm_cvtsi128_si32(_mm_cmpeq_epi32(xv,prevxv));
      prevxv=xv;
      a=(a+1u)&equal_mask;
      u32 ya0; __builtin_memcpy(&ya0, (const char *)cursor + 3, 4);
      __m256i yb = _mm256_set1_epi32((int)ya0);
      u32 dup0;
      if (a < 8) {
        u32 mm = (u32)_mm256_movemask_ps(_mm256_castsi256_ps(_mm256_cmpgt_epi32(yb, yhist)));
        dup0 = (u32)__builtin_popcount(_bzhi_u32(mm, a));
      } else dup0 = dup_long(sorted_pk, (u32)(cursor-sorted_pk), a, ya0);
      yhist = _mm256_alignr_epi8(yhist,
          _mm256_permute2x128_si256(yb,yhist,0x20),12);
      {
        u64 rec = rec0;
        u32 dup = dup0;
        uint16_t _hi; __builtin_memcpy(&_hi, (const char*)cursor + 4, 2);
        u32 d = (u32)_hi >> 7;
        u32 inc = cnt_q(d);
        cnt_ins(d);
        u32 q = gCur[d]++;
        u32 cnt = (inc - dup) & 0xFFFFFFu;
        alt[q] = ((u64)cnt << 40) | _bzhi_u64(rec,39);
      }
    }
    }
  }

  // -------- fused L2 (y bits 15..8) + L3 (y bits 7..0) + store -----------
  for (u32 b0 = 0; b0 < 512; b0++) {
    u32 s = gBnd[b0], e = gBnd[b0 + 1];
    u32 m = e - s;
    if (m == 0) continue;
    u64 *pk = PKB - s;
    if (m == 1) {
      { u32 _ix = (u32)(alt[s] & 0xFFFFFFu); u32 _c=(u32)((alt[s]>>40)&0xFFFFFFu); unsigned char *_p=STG6EP(gSTG,((gBOff[_ix >> BSH] += 6u) - 6u));
        *(u64*)_p=(((u64)_c)<<24)|_ix; }
      continue;
    }
    // histogram + out prefetch while the bucket is read the first time
    { u32 acc = 0;
      for (int d = 0; d < 64; d++) {
        u32 c = need32 ? gHis32[b0][d] : (u32)gHis16[b0][d];
        gHis[d] = acc; gCur[d] = s + acc; acc += c;   /* [w8s-K5] absolute cursors */
      }
      gHis[64] = acc; }
    if (m <= 65535u) {
      r64_init();
      for (u32 t = s; t < e; t++) {
        u64 rec = alt[t];
        u32 d = ((const unsigned char*)(alt+t))[4] >> 1;
        u32 inc = r64_q(d);
        r64_ins(d);
        u32 q = gCur[d]++;
        pk[q] = rec + ((u64)inc << 40);
      }
    } else {
      cnt_init2();
      for (u32 t = s; t < e; t++) {
        u64 rec = alt[t];
        u32 d = ((const unsigned char*)(alt+t))[4] >> 1;
        u32 inc = cnt_q2(d);
        cnt_ins2(d);
        u32 q = gCur[d]++;
        pk[q] = rec + ((u64)inc << 40);
      }
    }
    for (u32 d1 = 0; d1 < 64; d1++) {
      const u32 p = s + gHis[d1], qq = s + gHis[d1 + 1];
      if (p == qq) continue;
      if (qq - p == 1) {
        { u32 _ix = (u32)(pk[p] & 0xFFFFFFu); u32 _c=(u32)((pk[p]>>40)&0xFFFFFFu); unsigned char *_p=STG6EP(gSTG,((gBOff[_ix >> BSH] += 6u) - 6u));
          *(u64*)_p=(((u64)_c)<<24)|_ix; }
      } else {
        if (qq-p > 65535u) rank_l3_run<0>(pk,p,qq,gSTG);
        else if (qq-p <= 512u) rank_l3_run<2>(pk,p,qq,gSTG);
        else rank_l3_run<1>(pk,p,qq,gSTG);
      }
    }
  }
  flush_rank_output(out,gSTG);

}

#pragma GCC pop_options

评测结果