代码

// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #102741 <https://duck.ac/submission/102741>
//     用途：**直接复制了该提交的全部正文**（本文件 = 其正文逐字节 + 一处 [AMP] 放大诊断改动）。
//     #102741 自述的上游引用（pdoom、iMMIQ、saffah_cc_v41_260924、Elo、Bert Dobbelaere 等，
//     以及本账号在 1001/1001b 上的历代单变量调优）一并在其头部注释里保留。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103978 / #104001 / #104157 / #104233
//     <https://duck.ac/submission/103978> <https://duck.ac/submission/104001>
//     <https://duck.ac/submission/104157> <https://duck.ac/submission/104233>
//     用途：参考了其"**放大法**"仪器（让**幂等**相位多跑 N−1 遍取 Δ/(N−1)；末态逐位不变 ⇒ 提交自带
//     正确性闸门）。原文把它用在 `1001` 与 `1001c` 上；本发首次搬到 `1001b`。
// ======================
// ===== 思路 =====
// 本次为**试验性提交（诊断，只读用时）**，目的：取得 `1001b` 的 **partition（byte3 分桶 + 暂存 + 整行 NT 刷出） 相位的判题机口径边际成本**，
// 以便给 `kbc_` 一张**本题自己的**相位地图。
//
// 为什么必须上判题机（BRIEF §2.18.787，本族两题实测）：**本机相位表是否可信，取决于"该引擎的关键
// 工作集是否随 n 超 L3"**。`1001b` 是 n=1<<27、workspace 400+ MB 的流式引擎 ⇒ **几乎肯定超 L3**
// ⇒ 本机读数不可信（对照组：`1001c` 的引擎把访存全压在 512 B 热表上，本机 22.7% vs 判题机 23.9%，
// 可信；而 `1001` 本机 24% vs 判题机 34.1%，不可信）。
// 另：`1001b` 的**数据分布与 `1001` 不同**（列均值 8 / 慢块率 87.9% vs 5.96 / 23.6%，BRIEF §2.18.356 的
// 原始现场）⇒ **相位比例很可能不同**，不能照抄 `1001` 的表。
//
// 方法：**放大法**。`partition<false>(a,n)` 只读 `a`、确定性地重写 workspace 的同一批字节，且 `byte_cursor[]` 每次进入时重新初始化；**每遍前把 `pos[k]` 复位到 `starts[k]`** ⇒ 再跑一遍 = 输出逐位不变。Δ/(N−1) = partition 的判题机边际。
//   本机闸门（n=2^25、均匀随机 = 判题形状）：基线 / AMPP=4 / AMPM=3 / AMPT=3 四臂 hash
//   **全部 = `ce6d38f7e329e707`（sorted=1）** ✓ ⇒ 末态逐位不变 ⇒ Accepted。
//   本机读数（仅量级参考，不可外推）：baseline 397.7 ms；partition 58.3 ms(14.7%)、byte2 散列 40.9 ms(10.3%)、
//   byte1 转置 92.1 ms(23.2%) —— **与 `1001` 的分布明显不同**（1001：34.1/26.8/19.8%）。
//
// 单变量：`AMPP` 1 → 4。改动只在 [AMP] 标记的循环包裹与 `pos[]` 复位处，正文其余逐字节未动。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #102462 <https://duck.ac/submission/102462>
//     用途：**本文件正文 = #102462 的完整源码**，只改 `middle` 的**页大小承载方式**（一处）。
//     1.1 <https://duck.ac/submission/101817> duck.ac 用户 codex_astra —— 原榜首（475.142 ms），
//         `partition` / `sort_leaf` / 16 输入网络 / emit 逐字节来自它。
//     1.2 <https://duck.ac/submission/101606> saffah_cc_v41_agg1 —— "Output prefetch"。
//     1.3 <https://duck.ac/submission/98318> / 1.4 <https://duck.ac/submission/93444> pdoom 系骨架；
//     1.5 <https://duck.ac/submission/48133> iMMIQ；1.6 <https://duck.ac/submission/48186> Elo。
// [2] 16 输入 / 60 比较器网络：Bert Dobbelaere,
//     <https://bertdobbelaere.github.io/sorting_networks.html#N16L60D10>
// [3] Judge Duck FAQ（i3-8100 @3.6GHz，L1D 32KB 8-way，L2 256KB 4-way，L3 6MB）。
// ======================
// ===== 思路 =====
// 单变量改动（一处结构 + 它依赖的那一句 madvise）：`middle` 由 `alignas(64)` 改为
// `alignas(2097152)`，并在 `sort()` 开头（**首次触碰 `middle` 之前**）加
// `madvise(middle, 4 MiB, MADV_HUGEPAGE=14)`。其余逐字节 = #102462。
//
// 机理（**确定性、与机器无关、本机可复现**：work/b2k_pageconc.py）：
//   `small()` 的 byte2 散射要同时向 **256 个 byte2 子桶**写 u16；子桶间距 = `2*cap` = 4800 B，
//   而 `cap >= mu + 7sigma = 2365`（安全性硬约束）⇒ **子桶间距永远 > 4096 B**。
//   ⇒ `middle` 的**活跃地址集 = 256 条并发写流，每条各占一个不同的 4 KB 页 ⇒ ~300 页**；
//     实测"连续 B 次访问的互异 4 KB 页数 / B" = 1.000 / 0.998 / 0.994 / 0.987 / 0.971
//     （B = 1/2/4/8/16）⇒ **并发度是最大值，不存在"有驻留"**。
//   ⇒ L1 DTLB 只有 64 项 ⇒ 该 store 流每元素几乎必 miss L1 DTLB（STLB ~1536 装得下
//     ⇒ 无页走，但每次仍要走一次 STLB 查找的延迟与吞吐）。
//   ★ 关键：`middle` 的**实际使用区间只有 1.23 MiB**（256*cap*2 B）⇒ 若它落在一个 **2 MiB
//     巨页**里，这 256 条流的翻译就只需要 **1 个 TLB 项**，整条 DTLB 压力消失。
//   ★ 零风险：若判题机 THP = `never`（或 madvise 被 EPERM），本改动是空操作，行为逐字节
//     等价于 #102462；`middle` 基址变了，但那只涉及"数组间差值 mod 4096"这一类效应，
//     而该类效应已由 work/b2k_loadcensus.py 的 Δ 普查判定为**净**（改动后重新核对：
//     middle 与 columns/staging/byte_cursor/hist/starts/pos/workspace 的 Δ mod 4096
//     = 3392/3328/2112/1088/64/3136/3328，全部 ≠ 0）。
//
// 判据：本发**只**定价这一处（不做其它任何改动）。若 ≥ +0.5% ⇒ DTLB 就是 `sscat` 那 1.37×
//   未解释因子的来源（则同法可推广到 `workspace` 的 256 条 NT 流）；若 ≈0 ⇒ THP 在本机
//   不可用，或该延迟被完全隐藏 ⇒ 这条机制被定价为 0。
// ★ 声明（不做"不可达"论断）：这只测"页大小承载"**这一个**机制；它失败只关掉这一条路。
// ================
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#pragma GCC target("avx2,bmi,bmi2,popcnt")
#include <immintrin.h>
#include <algorithm>
#include <cstring>
#include <cstdint>
#include <utility>
#include <unistd.h>
#if defined(__APPLE__)
#define syscall(...) ((long)-1)
#endif
#ifndef BUFFER_BITS
#define BUFFER_BITS 6
#endif
#ifndef TOP_BITS
#define TOP_BITS 8
#endif
#ifndef UNROLL
#define UNROLL 1
#endif
static constexpr unsigned TB=TOP_BITS, K=1u<<TB, BS=1u<<BUFFER_BITS;
static constexpr unsigned NPORT = 134217728u;
alignas(64) static unsigned char workspace[3u*(256u*557056u + 256u*64u) + 4096u + 16u];
alignas(256) static unsigned char staging_bytes[K*320];
static unsigned byte_cursor[K];
alignas(64) static unsigned hist[K], starts[K], pos[K];
static inline unsigned read24(const unsigned char *s){unsigned x;memcpy(&x,s,4);return x&0xffffff;}
static void small_original(unsigned char*s,unsigned*d,unsigned n,unsigned prefix,unsigned*){
for(unsigned i=0;i<n;++i)d[i]=prefix|read24(s+3*i);
std::sort(d,d+n);
}
alignas(2097152) static uint16_t middle[1<<21];
namespace fastsort {
constexpr unsigned kRadix=256;
constexpr unsigned kLeafLimit = 8192;
using U32 = unsigned; using Byte = uint8_t; using U16 = uint16_t; using U64 = uint64_t;
using Vec256 = __m256i; using Vec128 = __m128i;
template <class T> inline T load(const void *src) { T value; std::memcpy(&value, src, sizeof value); return value; }
template <class F, size_t... I> [[gnu::always_inline]] inline void each(F &&f, std::index_sequence<I...>) {
(f(std::integral_constant<size_t, I>{}), ...);
}
template <size_t N, class F> [[gnu::always_inline]] inline void repeat(F &&f) {
each(std::forward<F>(f), std::make_index_sequence<N>{});
}
}
namespace byte_sort {
using namespace fastsort;
[[gnu::noinline]] inline void counting_sort(const U16 *src, U32 n, U32 *dst, U32 base) {
alignas(64) U32 counts[65536] = {};
for (U32 i = 0; i < n; ++i) ++counts[load<U16>(src + i)];
for (U32 v = 0; v < 65536; ++v)
for (U32 c = counts[v]; c; --c) *dst++ = base | v;
}
[[gnu::always_inline]] inline void transpose8(const Vec256 *x, Byte *out) {
Vec256 a[8], b[8];
#pragma GCC unroll 4
for (U32 i = 0; i < 4; ++i) {
a[2 * i] = _mm256_unpacklo_epi8(x[2 * i], x[2 * i + 1]);
a[2 * i + 1] = _mm256_unpackhi_epi8(x[2 * i], x[2 * i + 1]);
}
#pragma GCC unroll 2
for (U32 i = 0; i < 2; ++i) {
#pragma GCC unroll 2
for (U32 j = 0; j < 2; ++j) {
b[4 * i + 2 * j] = _mm256_unpacklo_epi16(a[4 * i + j], a[4 * i + j + 2]);
b[4 * i + 2 * j + 1] = _mm256_unpackhi_epi16(a[4 * i + j], a[4 * i + j + 2]);
}
}
#pragma GCC unroll 4
for (U32 i = 0; i < 4; ++i) {
Vec256 lo = _mm256_unpacklo_epi32(b[i], b[i + 4]);
Vec256 hi = _mm256_unpackhi_epi32(b[i], b[i + 4]);
_mm256_store_si256((Vec256 *)(out + 32 * i), _mm256_permute2x128_si256(lo, hi, 0x20));
_mm256_store_si256((Vec256 *)(out + 128 + 32 * i), _mm256_permute2x128_si256(lo, hi, 0x31));
}
}
[[gnu::always_inline]] inline void sort16_columns(const Byte *src, Byte *out, Byte *out2) {
Vec256 x[16];
repeat<16>([&](auto j) __attribute__((always_inline)) {
x[j] = _mm256_load_si256((const Vec256 *)(src + 256 * j));
});
#define CMP16(A,B) { Vec256 lo = _mm256_min_epu8(x[A], x[B]); \
x[B] = _mm256_max_epu8(x[A], x[B]); x[A] = lo; }
CMP16(0, 13)
CMP16(1, 12)
CMP16(2, 15)
CMP16(3, 14)
CMP16(4, 8)
CMP16(5, 6)
CMP16(7, 11)
CMP16(9, 10)
CMP16(0, 5)
CMP16(1, 7)
CMP16(2, 9)
CMP16(3, 4)
CMP16(6, 13)
CMP16(8, 14)
CMP16(10, 15)
CMP16(11, 12)
CMP16(0, 1)
CMP16(2, 3)
CMP16(4, 5)
CMP16(6, 8)
CMP16(7, 9)
CMP16(10, 11)
CMP16(12, 13)
CMP16(14, 15)
CMP16(0, 2)
CMP16(1, 3)
CMP16(4, 10)
CMP16(5, 11)
CMP16(6, 7)
CMP16(8, 9)
CMP16(12, 14)
CMP16(13, 15)
CMP16(1, 2)
CMP16(3, 12)
CMP16(4, 6)
CMP16(5, 7)
CMP16(8, 10)
CMP16(9, 11)
CMP16(13, 14)
CMP16(1, 4)
CMP16(2, 6)
CMP16(5, 8)
CMP16(7, 10)
CMP16(9, 13)
CMP16(11, 14)
CMP16(2, 4)
CMP16(3, 6)
CMP16(9, 12)
CMP16(11, 13)
CMP16(3, 5)
CMP16(6, 8)
CMP16(7, 9)
CMP16(10, 12)
CMP16(3, 4)
CMP16(5, 6)
CMP16(7, 8)
CMP16(9, 10)
CMP16(11, 12)
CMP16(6, 7)
CMP16(8, 9)
#undef CMP16
transpose8(x, out);
transpose8(x + 8, out2);
}
[[gnu::always_inline]] inline U32 restore_counts(U32 *count, U32 *exceptional) {
U32 *end = count + 256; U32 over;
const Vec256 limit16 = _mm256_set1_epi32(16), limit32 = _mm256_set1_epi32(32);
__asm__ volatile("vpxor %%ymm4, %%ymm4, %%ymm4\n\t"
".p2align 4\n\t1:\n\t"
".irp reg,0,1,2,3; vmovdqa 32*\\reg(%[count]), %%ymm\\reg; .endr\n\t"
".irp reg,0,1,2,3; vpsrld $8, %%ymm\\reg, %%ymm\\reg; .endr\n\t"
".irp reg,0,1,2,3; vmovdqa %%ymm\\reg, 32*\\reg(%[count]); .endr\n\t"
"vpmaxud %%ymm1, %%ymm0, %%ymm0; vpmaxud %%ymm3, %%ymm2, %%ymm2\n\t"
"vpmaxud %%ymm2, %%ymm0, %%ymm0\n\t"
"vpcmpgtd %[limit16], %%ymm0, %%ymm1\n\t"
"vmovmskps %%ymm1, %%eax; movl %%eax, (%[flags])\n\t"
"vpmaxud %%ymm0, %%ymm4, %%ymm4\n\t"
"addq $128, %[count]; addq $4, %[flags]\n\t"
"cmpq %[end], %[count]; jb 1b\n\t"
"vpcmpgtd %[limit32], %%ymm4, %%ymm0; vmovmskps %%ymm0, %k[over]\n\t"
: [count] "+&r"(count), [flags] "+&r"(exceptional), [over] "=&r"(over)
: [end] "r"(end), [limit16] "x"(limit16), [limit32] "x"(limit32)
: "rax", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "cc", "memory");
return over;
}
[[gnu::always_inline]] inline U32 *emit32_buckets16(const Byte *first, const Byte *second, const U32 *count,
U32 *dst, Vec256 prefix) {
const Vec256 step = _mm256_set1_epi32(256); const U32 *end = count + 32;
__asm__ volatile(".p2align 4\n\t1:\n\t"
".irp lane,0,1,2,3,4,5,6,7\n\t"
"movl 4*\\lane(%[count]), %%eax\n\t"
"vpmovzxbd 8*\\lane(%[first]), %%ymm1; vpmovzxbd 8*\\lane(%[second]), %%ymm2\n\t"
"vpor %[prefix], %%ymm1, %%ymm1; vpor %[prefix], %%ymm2, %%ymm2\n\t"
"vmovdqu %%ymm1, (%[dst]); vmovdqu %%ymm2, 32(%[dst])\n\t"
"leaq (%[dst],%%rax,4), %[dst]\n\t"
"vpaddd %[step], %[prefix], %[prefix]\n\t"
".endr\n\t"
"addq $64, %[first]; addq $64, %[second]; addq $32, %[count]\n\t"
"cmpq %[end], %[count]; jb 1b\n\t"
: [dst] "+&r"(dst), [prefix] "+&x"(prefix), [count] "+&r"(count), [first] "+&r"(first),
[second] "+&r"(second)
: [end] "r"(end), [step] "x"(step)
: "rax", "ymm1", "ymm2", "cc", "memory");
return dst;
}
inline void sort_leaf(const U16 *src, U32 n, U32 *dst, U32 base, U32 *end) {
if (n < 64) {
U32 copy[64];
for (U32 i = 0; i < n; ++i) copy[i] = base | load<U16>(src + i);
std::sort(copy, copy + n);
std::memcpy(dst, copy, n * 4); return;
}
if (n > kLeafLimit) return counting_sort(src, n, dst, base);
alignas(64) static Byte columns[kRadix * kLeafLimit];
alignas(64) U32 count[256];
std::memset(columns, 255, 16 * kRadix);
#ifndef AMPT
#define AMPT 1
#endif
{for (U32 i = 0; i < 256; ++i) count[i] = i;
#pragma GCC unroll 16
for (U32 i = 0; i < n; ++i) {
U32 v = load<U16>(src + i), h = v >> 8, slot = count[h];
count[h] = slot + 256; columns[slot] = Byte(v);
}}
for (int amp_t = 1; amp_t < AMPT; ++amp_t) {
for (U32 i = 0; i < 256; ++i) count[i] = i;
#pragma GCC unroll 16
for (U32 i = 0; i < n; ++i) {
U32 v = load<U16>(src + i), h = v >> 8, slot = count[h];
count[h] = slot + 256; columns[slot] = Byte(v);
}
}
U32 exceptional[8];
if (restore_counts(count, exceptional)) return counting_sort(src, n, dst, base);
const Vec256 idx = _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7);
for (U32 block = 0; block < 256; block += 32) {
alignas(32) Byte first[256], second[256];
#pragma GCC unroll 16
for(unsigned pf=0;pf<16;++pf)asm volatile("prefetchw %0"::"m"(*(const char*)((uintptr_t)dst+2048+64*pf)));
sort16_columns(columns + block, first, second);
Vec256 prefix = _mm256_set1_epi32(base | (block << 8));
if (exceptional[block >> 5] == 0 && end - dst >= 512) {
dst = emit32_buckets16(first, second, count + block, dst, prefix);
continue;
}
for (U32 i = 0; i < 32; ++i) {
U32 h = block + i, c = count[h];
if (__builtin_expect(c <= 16, 1)) {
Vec256 w0 = _mm256_or_si256(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const Vec128 *)(first + 8 * i))), prefix);
Vec256 w1 = _mm256_or_si256(_mm256_cvtepu8_epi32(_mm_loadl_epi64((const Vec128 *)(second + 8 * i))), prefix);
if (__builtin_expect(end - dst >= 16, 1)) {
_mm256_storeu_si256((Vec256 *)dst, w0);
_mm256_storeu_si256((Vec256 *)(dst + 8), w1);
} else {
_mm256_maskstore_epi32((int *)dst, _mm256_cmpgt_epi32(_mm256_set1_epi32((int)(c < 8u ? c : 8u)), idx), w0);
_mm256_maskstore_epi32((int *)(dst + 8), _mm256_cmpgt_epi32(_mm256_set1_epi32((int)(c > 8u ? c - 8u : 0u)), idx), w1);
}
} else {
Byte copy[32];
for (U32 j = 0; j < c; ++j) copy[j] = columns[256 * j + h];
std::sort(copy, copy + c);
for (U32 j = 0; j < c; ++j) dst[j] = base | (h << 8) | copy[j];
}
dst += c;
prefix = _mm256_add_epi32(prefix, _mm256_set1_epi32(256));
}
}
}
}
static void small(unsigned char *__restrict s,unsigned *__restrict d,unsigned n,unsigned prefix,unsigned *array_end){
if(n<65536){small_original(s,d,n,prefix,array_end);return;}
unsigned mu=(n+255)/256;unsigned need=mu+(unsigned)(6.5*__builtin_sqrt((double)mu)+0.5)+8;unsigned odd=(need+31)/32;if(!(odd&1))odd++;unsigned cap=32u*odd;
if(cap>4096){small_original(s,d,n,prefix,array_end);return;}
unsigned mp[256];
#ifndef AMPM
#define AMPM 1
#endif
for(unsigned k=0;k<256;k++)mp[k]=k*cap;
{const unsigned char *r=s;const unsigned char *rend=s+3u*n;
#pragma GCC unroll 16
for(;r<rend;r+=3){
unsigned k=r[2];uint16_t value;memcpy(&value,r,2);
middle[mp[k]++]=value;
}}
for(int amp_m=1;amp_m<AMPM;amp_m++){
for(unsigned k=0;k<256;k++)mp[k]=k*cap;
{const unsigned char *r=s;const unsigned char *rend=s+3u*n;
#pragma GCC unroll 16
for(;r<rend;r+=3){
unsigned k=r[2];uint16_t value;memcpy(&value,r,2);
middle[mp[k]++]=value;
}}
}
for(unsigned k=0;k<256;k++)if(mp[k]-k*cap>cap){small_original(s,d,n,prefix,array_end);return;}
for(unsigned k=0;k<256;k++){
unsigned count=mp[k]-k*cap,pre=prefix|(k<<16);
const uint16_t *src=middle+k*cap;
byte_sort::sort_leaf(src,count,d,pre,array_end);
d+=count;
}
}
template<bool Reserved> static bool partition(unsigned *__restrict a,unsigned n,unsigned capacity=0){
for(unsigned k=0;k<K;k++)byte_cursor[k]=k*320;
for(unsigned page=0;page<n;page+=1024){
if(page+5120<n)asm volatile("orl $0, %0" : "+m"(a[page+4112]));
unsigned stop=std::min(n,page+1024);
#pragma GCC unroll 16
for(unsigned i=page;i<stop;i++){
unsigned x=a[i],k=x>>24,o=byte_cursor[k]+3;
memcpy(staging_bytes+o-3,&x,4);
if(__builtin_expect((o&63)==0,0)){
if(Reserved && pos[k]+64-starts[k]>capacity){_mm_sfence();return false;}
unsigned char *d=workspace+3*pos[k];
const unsigned char *src=staging_bytes+o-192;
#pragma GCC unroll 6
for(unsigned j=0;j<192;j+=32)
_mm256_stream_si256((__m256i*)(d+j),_mm256_load_si256((const __m256i*)(src+j)));
pos[k]+=64;o-=192;
}
byte_cursor[k]=o;
}
}
for(unsigned k=0;k<K;k++)pos[k]+=(byte_cursor[k]-k*320)/3;
_mm_sfence();
if(Reserved)for(unsigned k=0;k<K;k++)if(pos[k]-starts[k]>capacity)return false;
return true;
}
void sort(unsigned *a,int n){
{ static int done=0; if(!done){ done=1; syscall(28,(void*)workspace,sizeof(workspace),23);
  syscall(28,(void*)middle,sizeof(middle),14); } }
if(n<256){std::sort(a,a+n);return;}
unsigned capacity=((unsigned(n)+K-1)/K*17/16+BS-1)&~(BS-1);
for(unsigned k=0;k<K;k++)starts[k]=pos[k]=k*capacity;
if(!partition<true>(a,n,capacity)){
memset(hist,0,sizeof(hist));
for(int i=0;i<n;i++)++hist[a[i]>>(32-TB)];
unsigned z=0;
for(unsigned k=0;k<K;k++){starts[k]=pos[k]=z;z=(z+hist[k]+BS-1)&~(BS-1);}
partition<false>(a,n);
}
#ifndef AMPP
#define AMPP 4
#endif
/* [AMP] partition is deterministic and reads nothing but `a`; resetting pos[] and re-running it
   leaves workspace/pos[] bit-identical to a single run => output unchanged (Accepted). */
for(int amp_p=1;amp_p<AMPP;amp_p++){ for(unsigned k=0;k<K;k++)pos[k]=starts[k]; partition<false>(a,n); }
unsigned z=0;
for(unsigned k=0;k<K;k++){
unsigned size=pos[k]-starts[k];
memcpy(workspace+3*(pos[k]&~(BS-1)),staging_bytes+k*320,3*(pos[k]&(BS-1)));
small(workspace+3*starts[k],a+z,size,k<<24,a+n);
z+=size;
}
}

评测结果