// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2,提交 #102895 <https://duck.ac/submission/102895>
// 用途:**本文件正文的基底**(30 位打包键 c0..c5、按 (c0,c1) MSD 散射 + 桶内两趟 10 位 LSD、
// 融合 walk(`clz(kprev^k)` 经 `D5B` 得 LCP、等键 run 用 `suf_less` 精修)、
// `T4a/T3b/T3x` 两级分块写入器、SA-IS 回退)。本发只改「桶内 LSD 与 walk 的融合」一处(见思路)。
// [2] duck.ac 用户 saffah_cc_v41_260924,提交 #101791 / #100017 / #99905
// <https://duck.ac/submission/101791> <https://duck.ac/submission/100017>
// <https://duck.ac/submission/99905>
// 用途:上述引擎骨架(打包键、桶内 LSD、`D5B` 长度表、`prefix_excl1024` 的 AVX2 写法)
// 与 `suf_less` refine 的原始来源(经 [1] 转抄)。
// [3] 本账号 saffah_cc_v41_agg1,提交 #102957 <https://duck.ac/submission/102957>
// (文件 problems/1006e6/work/k6_1.cpp)用途:**本发相对它做单变量改动** —— 它是我们上一代
// AC 最好件(9.851890 ms),本发的对照基线,正文与它逐字节相同(除下述一处)。
// [4] 本账号 saffah_cc_v41_agg1,提交 #102784 / #103944
// <https://duck.ac/submission/102784> <https://duck.ac/submission/103944>
// 用途:[3] 的直接上游,以及"key 加宽/形状闸门"那一支的实验记录(本题未采用)。
// [5] /home/yjp/duck.ac/BRIEF.md §2.18.605(先画像再下刀)与 problems/1006e6/notes.md。
// 各提交的公开源码中未附许可证声明;此处已按账号 / 原始提交 URL / 所用内容逐条列出。
// ======================
// ===== 思路 =====
// 【本发改了什么(相对 #102957,单变量)】**把 walk 融进最后一趟桶内 LSD**。
// 原结构:桶内 LSD 两趟(pass1 `RB[lo..hi)` → 寄存器级 scratch `SCR`;pass2 `SCR` → **写回 `RB`**),
// 全部桶做完之后 walk 再**从 DRAM 顺序读 `RB`**(8 MB)并发射。
// 新结构:pass2 改写进**第二个 scratch `SCR2`**(仍是桶内 ~977 条 = 7.8 KB,L1 常驻),
// 然后**立刻用 walk 处理该桶**(`W = SCR2 - lo`,让 `W[x]` 仍是全局排序位置 x 的记录)。
// ⇒ 省掉 **pass2 对 RB 的 8 MB 写 + 8 MB RFO** 与 **walk 对 RB 的 8 MB 顺序读**,共约 24 MB 的 DRAM 流量。
// 正确性依据:等键 run **不可能跨桶**(桶由键的高 10 位确定,等键 ⇒ 同桶)⇒ 每个桶的 walk 边界
// (`x .. xe = START[bk+1]`) 是精确的;`kprev` 仍然跨桶携带,所以桶边界处的 LCP 与原来一致;
// m==1 的桶直接用 `W = RB`(那条记录本来就在位);m==0 跳过。
// run 超 RUNMAX 时仍然 `g_fallback=1; o=ostart`(丢弃已发射的输出、整段回退 SA-IS)——回退路径
// 与 #102957 逐字节相同。
// 【本发是**探针式提交**】本族"记账类"改动**本地会低估判题机**(BRIEF §2.18.4),而本机 L2 是判题机
// 的 4 倍,8 MB 的数组在本机可能被 L3 吸收 ⇒ 本地只量到 **−0.94 %**(n=1e6,min-of-5),
// 判题机应更大。逐点明细会直接给出这笔账的兑现率。
// 【正确性闸门】`-DBENCH`:n=2e5 与 n=1e6 的 `olen`/`ohash` 与基线 **逐字节相同**
// (n=1e6: olen=8888894、ohash=16407443965629840480);另用 stdin 驱动对 8 种**结构化/回退**
// 输入(a^k、(ab)^k、(abc)^k、Thue-Morse、Fibonacci、长 run、二元伪随机、26 字母伪随机)
// 逐字节对比,**全部相同**(其中 7 种触发 SA-IS 回退路径)。
// ================
#pragma GCC optimize("O3","unroll-loops","rename-registers")
#pragma GCC target("arch=skylake,avx2,bmi,bmi2,popcnt,lzcnt")
#define JUDGE 1
// 1006e6 v3 (lane e_variant_lane): SA + LCP for a uniformly-random lowercase string.
// stage 1 build u64 records [key30(c0..c5, 5 bits each) | idx20] (BMI2 pext, no gathers)
// stage 2 3 LSD radix passes, 10 bits each, 1024 buckets (L1 counters)
// stage 3 one walk: emit SA digits, derive LCP from packed-key xor (clz), refine equal-key runs
// stage 4 emit the LCP line (byte array -> digits)
// Fallback to SA-IS when a 6-char-key run exceeds RUNMAX (non-random data).
#include <cstdint>
#include <cstring>
#include <immintrin.h>
#include <immintrin.h>
#define MAXN 1000005
#define PADB 64
static unsigned char str[MAXN + PADB];
__attribute__((aligned(64))) static uint64_t RA[MAXN], RB[MAXN];
static uint32_t CNT[1024];
static uint32_t CNT2[1024];
static uint32_t START[1025];
#define SCRSZ 65536
static uint64_t SCR[SCRSZ];
static uint64_t SCR2[SCRSZ]; // y2 FUSE: pass2 target -- stays in L1 (a bucket is ~977 records)
static unsigned char LCPB[MAXN + 8];
static uint32_t D4[10000];
static unsigned char L10[64];
/* LANE l1006e6b, copied from #100017: `(clz(x)-2)/5` is a multiply-high chain ON the
critical path of the four LCPB byte stores in the walk; clz of a non-zero 30-bit key
is 2..31, so the quotient is a 32-entry table (and clz(0) = 32 maps to 6, which is what
the old `(__builtin_clz(0)-2)/5` produced too). */
static const unsigned char D5B[64] = {0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6};
static unsigned char CLA[64];
static uint64_t PAIR[256];
static const unsigned PW[12] = {0, 1, 10, 100, 1000, 10000, 100000, 1000000,
10000000, 100000000, 1000000000u, 0};
/* AVX2 exclusive prefix sum over the 1024 radix counters. Directly copied from the
rival's #99905, whose comment reads: the scalar form was a serial `s += c` chain of
1024 dependent adds plus 1024 single-word stores, i.e. ~1024 cycles per call, and it
is called twice per non-empty bucket (~676 of them at n=1e6). `acc - original`
recovers the exclusive prefix, so the only shuffles left are the two vpermd broadcasts
of lane 3 and lane 7. */
static inline void prefix_excl864(uint32_t *C) {
__m256i carry = _mm256_setzero_si256();
const __m256i msk = _mm256_setr_epi32(0, 0, 0, 0, -1, -1, -1, -1);
const __m256i id3 = _mm256_set1_epi32(3), id7 = _mm256_set1_epi32(7);
for (int d = 0; d < 864; d += 8) {
__m256i o = _mm256_loadu_si256((const __m256i *)(C + d));
__m256i v = o;
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 4));
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 8));
v = _mm256_add_epi32(v, _mm256_and_si256(_mm256_permutevar8x32_epi32(v, id3), msk));
__m256i acc = _mm256_add_epi32(v, carry);
_mm256_storeu_si256((__m256i *)(C + d), _mm256_sub_epi32(acc, o));
carry = _mm256_permutevar8x32_epi32(acc, id7);
}
}
static inline uint64_t pext64(uint64_t v, uint64_t m) {
uint64_t r; __asm__("pext %2,%1,%0" : "=r"(r) : "r"(v), "r"(m)); return r;
}
#define M5X8 0x1F1F1F1F1F1F1F1FULL
#define REC_KEY(r) ((uint32_t)((r) >> 20))
#define REC_IDX(r) ((uint32_t)((r) & 0xFFFFFu))
static inline uint32_t key_at(int i) {
uint64_t v; memcpy(&v, str + i, 8);
return (uint32_t)(pext64(__builtin_bswap64(v), M5X8) >> 10);
}
/* [x6 cbd] the MSD bucket of suffix i is exactly bits 20..29 of key_at(i), i.e.
c1 | c0<<5 with cj = str[i+j] & 31 -- obtainable from a single two-byte load with no
bswap and no pext. Used by the COUNT pass only; the scatter still needs the full
30-bit key for the record it writes. */
static inline uint32_t bk_dig(int i) {
uint16_t v; memcpy(&v, str + i, 2);
return ((uint32_t)(v >> 8) & 31u) | (((uint32_t)v & 31u) << 5);
}
static inline int len10c(uint32_t v) { return CLA[31 - __builtin_clz(v | 1)]; }
static inline int len10(uint32_t v) {
int b = 31 - __builtin_clz(v | 1);
int l = L10[b];
l -= (v < PW[l]);
return l;
}
// pack v (<= 9999999) as digits followed by a space in an 8-byte word (left-aligned, low byte first)
/* LANE pragma_lane: TWO-LEVEL CHUNK WRITER. The old form built the number from a 4+4 digit split
and then LENGTH-ARITHMETIC'd it into place: 2 magic divisions, 6 compares, a cmov, and TWO variable
shifts (`>>= (8-l)*8` and `' ' << (8*l)`). This form looks the digits up in two 8 KB L1-resident
tables of PRE-SHIFTED, pre-spaced 8-byte words (`T4a[A]` = A's digits + space + its length in byte
7; `T3b[B]` = B as 3 zero-padded digits + space) and places the low part with ONE shift by
8*len(A). Byte-identical: the same digits in the same bytes, and the bytes ABOVE the space are
always overwritten by the next store (the last number of the line still goes through
`putu_exact`). MEASURED in a faithful mikro: -19.3 % instructions at unchanged IPC (3.29 -> 3.24),
i.e. -18.1 % cycles, on a loop that IS throughput-bound. */
static uint64_t T3x[1024], T4a[1024]; static uint32_t T3b[1024];
static unsigned char LENA[1024];
static void build_t4(void) {
for (int v = 0; v < 1024; v++) {
char t[8] = {0,0,0,0,0,0,0,0}; int l;
if (v < 10) { t[0] = (char)('0' + v); t[1] = ' '; l = 2; }
else if (v < 100) { t[0] = (char)('0' + v / 10); t[1] = (char)('0' + v % 10); t[2] = ' '; l = 3; }
else { t[0] = (char)('0' + v / 100); t[1] = (char)('0' + (v / 10) % 10); t[2] = (char)('0' + v % 10); t[3] = ' '; l = 4; }
uint64_t u = 0; memcpy(&u, t, 8);
T3x[v] = u | ((uint64_t)(l - 1) << 56); /* for the A == 0 case */
char z[8] = {(char)('0' + (v / 100) % 10), (char)('0' + (v / 10) % 10), (char)('0' + v % 10), ' ', 0,0,0,0};
uint64_t uz = 0; memcpy(&uz, z, 8); T3b[v] = (uint32_t)uz; /* 3 zero-padded digits + space */
char w[8] = {0,0,0,0,0,0,0,0}; int la;
if (v < 10) { w[0] = (char)('0' + v); w[1] = ' '; la = 1; }
else if (v < 100) { w[0] = (char)('0' + v / 10); w[1] = (char)('0' + v % 10); w[2] = ' '; la = 2; }
else if (v < 1000) { w[0] = (char)('0' + v / 100); w[1] = (char)('0' + (v / 10) % 10); w[2] = (char)('0' + v % 10); w[3] = ' '; la = 3; }
else { w[0] = '1'; w[1] = '0'; w[2] = '0'; w[3] = '0'; w[4] = ' '; la = 4; }
uint64_t uw = 0; memcpy(&uw, w, 8); T4a[v] = uw | ((uint64_t)la << 56); LENA[v] = (unsigned char)la;
}
}
static inline uint64_t pack_num(uint32_t v, int &adv) {
uint32_t A = v / 1000, B = v - A * 1000;
uint64_t ta = T4a[A], tb = (uint64_t)T3b[B];
/* LANE l1006e6b (copied from #100017): T3x[B] was loaded unconditionally and
selected with a cmov although it is only needed when A == 0 (v < 1000);
A's length now sits in T4a's byte 7 instead of the LENA[] array. */
if (!A) { uint64_t tn = T3x[B]; adv = (int)((tn >> 56) + 1); return tn; }
int la = (int)(ta >> 56);
uint64_t u = (ta & 0x00FFFFFFFFFFFFFFULL) | (tb << (8 * la));
adv = la + 4;
return u;
}
static inline char *pack_emit(char *o, uint32_t v) {
int adv;
uint64_t u = pack_num(v, adv);
memcpy(o, &u, 8);
return o + adv;
}
// write v (<= 9999999) followed by a space; 8-byte store, advances by len+1
// byte-exact writer (no overrun) for the final element of each line
static inline char *putu_exact(char *p, unsigned v) {
char tmp[8];
int i = 8;
while (v >= 100) { unsigned q = v / 100, r = v - q * 100; tmp[--i] = '0' + r % 10; tmp[--i] = '0' + r / 10; v = q; }
if (v >= 10) { tmp[--i] = (char)('0' + v % 10); tmp[--i] = (char)('0' + v / 10); }
else tmp[--i] = (char)('0' + v);
memcpy(p, tmp + i, 8 - i);
return p + (8 - i);
}
static inline void build_D4_if_needed() {
static int ready = 0;
if (ready) return;
for (int i = 0; i < 10000; i++) {
char t[4];
t[0] = (char)('0' + i / 1000); t[1] = (char)('0' + (i / 100) % 10);
t[2] = (char)('0' + (i / 10) % 10); t[3] = (char)('0' + i % 10);
memcpy(&D4[i], t, 4);
}
ready = 1;
}
static inline char *wr(char *p, uint32_t v) {
build_D4_if_needed();
uint32_t hi = v / 10000, lo = v % 10000;
uint64_t u = (uint64_t)D4[hi] | ((uint64_t)D4[lo] << 32);
int l = len10(v);
u >>= (8 - l) * 8;
u |= (uint64_t)' ' << (8 * l);
memcpy(p, &u, 8);
return p + l + 1;
}
// small-value writer for the LCP line (v < 1000 typical)
static uint64_t T3[1024];
static void build_T3() {
for (int v = 0; v < 1024; v++) {
char t[8]; int l;
if (v < 10) { t[0] = (char)('0' + v); t[1] = ' '; l = 2; }
else if (v < 100) { t[0] = (char)('0' + v / 10); t[1] = (char)('0' + v % 10); t[2] = ' '; l = 3; }
else { t[0] = (char)('0' + v / 100); t[1] = (char)('0' + (v / 10) % 10); t[2] = (char)('0' + v % 10); t[3] = ' '; l = 4; }
uint64_t u = 0; memcpy(&u, t, 8);
T3[v] = u; // left-aligned: low byte = first char
}
}
static inline char *wrs(char *p, uint32_t v) {
if (v >= 1024) return wr(p, v); // rare: long LCP (non-random data)
uint64_t u = T3[v];
memcpy(p, &u, 8);
return p + ((v < 10) ? 2 : (v < 100) ? 3 : 4);
}
// ---------------- SA-IS fallback (unchanged behaviour) ----------------
#define TPOOL 24
static int g_bkt[MAXN + 2];
static int g_cntc[MAXN + 2];
static unsigned char t_pool[TPOOL][(MAXN + 7) / 8];
static inline int tget(const unsigned char *t, int i) { return (t[i >> 3] >> (i & 7)) & 1; }
static inline void tset(unsigned char *t, int i, int b) {
if (b) t[i >> 3] |= (unsigned char)(1 << (i & 7));
else t[i >> 3] &= (unsigned char)~(1 << (i & 7));
}
static inline int isLMS(const unsigned char *t, int i) { return i > 0 && tget(t, i) && !tget(t, i - 1); }
template <typename T> static void getCounts(const T *s, int *C, int n, int K) {
for (int i = 0; i <= K; i++) C[i] = 0;
for (int i = 0; i < n; i++) C[(int)s[i]]++;
}
static void getBuckets(const int *C, int *B, int K, bool end) {
int sum = 0;
if (end) { for (int i = 0; i <= K; i++) { sum += C[i]; B[i] = sum; } }
else { for (int i = 0; i <= K; i++) { B[i] = sum; sum += C[i]; } }
}
template <typename T> static void induceSAl(const unsigned char *t, int *SA, const T *s, int *bkt, int n, int K) {
getBuckets(g_cntc, bkt, K, false);
for (int i = 0; i < n; i++) { int j = SA[i] - 1; if (j >= 0 && !tget(t, j)) SA[bkt[(int)s[j]]++] = j; }
}
template <typename T> static void induceSAs(const unsigned char *t, int *SA, const T *s, int *bkt, int n, int K) {
getBuckets(g_cntc, bkt, K, true);
for (int i = n - 1; i >= 0; i--) { int j = SA[i] - 1; if (j >= 0 && tget(t, j)) SA[--bkt[(int)s[j]]] = j; }
}
template <typename T> static void sais(const T *s, int *SA, int n, int K, int depth) {
unsigned char *t = t_pool[depth];
int *bkt = g_bkt;
int i, j;
tset(t, n - 2, 0); tset(t, n - 1, 1);
for (i = n - 3; i >= 0; i--) tset(t, i, (s[i] < s[i + 1] || (s[i] == s[i + 1] && tget(t, i + 1))));
getCounts(s, g_cntc, n, K);
getBuckets(g_cntc, bkt, K, true);
for (i = 0; i < n; i++) SA[i] = -1;
for (i = 1; i < n; i++) if (isLMS(t, i)) SA[--bkt[(int)s[i]]] = i;
induceSAl(t, SA, s, bkt, n, K);
induceSAs(t, SA, s, bkt, n, K);
int n1 = 0;
for (i = 0; i < n; i++) if (isLMS(t, SA[i])) SA[n1++] = SA[i];
for (i = n1; i < n; i++) SA[i] = -1;
int name = 0, prev = -1;
for (i = 0; i < n1; i++) {
int pos = SA[i];
bool diff = false;
int dmax = n - (pos > prev ? pos : prev);
for (int d = 0; d < dmax; d++)
if (prev == -1 || s[pos + d] != s[prev + d] || tget(t, pos + d) != tget(t, prev + d)) { diff = true; break; }
else if (d > 0 && (isLMS(t, pos + d) || isLMS(t, prev + d))) break;
if (diff) { name++; prev = pos; }
pos = (pos % 2 == 0) ? pos / 2 : (pos - 1) / 2;
SA[n1 + pos] = name - 1;
}
for (i = n - 1, j = n - 1; i >= n1; i--) if (SA[i] >= 0) SA[j--] = SA[i];
int *SA1 = SA, *s1 = SA + n - n1;
if (name < n1) sais<int>(s1, SA1, n1, name - 1, depth + 1);
else for (i = 0; i < n1; i++) SA1[s1[i]] = i;
getCounts(s, g_cntc, n, K);
getBuckets(g_cntc, bkt, K, true);
for (i = 1, j = 0; i < n; i++) if (isLMS(t, i)) s1[j++] = i;
for (i = 0; i < n1; i++) SA1[i] = s1[SA1[i]];
for (i = n1; i < n; i++) SA[i] = -1;
for (i = n1 - 1; i >= 0; i--) { j = SA[i]; SA[i] = -1; SA[--bkt[(int)s[j]]] = j; }
induceSAl(t, SA, s, bkt, n, K);
induceSAs(t, SA, s, bkt, n, K);
}
static int g_sa[MAXN];
// ---------------- helpers: full suffix compare / lcp ----------------
static inline bool suf_less(uint32_t i, uint32_t j) {
const unsigned char *p = str + i, *q = str + j;
for (;;) {
uint64_t a; memcpy(&a, p, 8); a = __builtin_bswap64(a);
uint64_t b; memcpy(&b, q, 8); b = __builtin_bswap64(b);
if (a != b) return a < b;
p += 8; q += 8;
}
}
static inline int lcp_len(uint32_t i, uint32_t j) {
const unsigned char *p = str + i, *q = str + j;
int h = 0;
for (;;) {
uint64_t a; memcpy(&a, p, 8); a = __builtin_bswap64(a);
uint64_t b; memcpy(&b, q, 8); b = __builtin_bswap64(b);
if (a != b) return h + ((int)(__builtin_clzll(a ^ b) >> 3));
h += 8; p += 8; q += 8;
}
}
#define RUNMAX 48
static int g_fallback = 0;
static char *g_sa_out_start;
static long g_osal;
#ifdef PHASE
static unsigned long long P[16];
#define MARK(k) P[k] = __builtin_ia32_rdtsc()
#else
#define MARK(k)
#endif
static int solve_string(const unsigned char *q, long sn, char **opo) {
int n = 0;
while (n < sn && (q[n] == ' ' || q[n] == '\n' || q[n] == '\r' || q[n] == '\t')) { n++; }
{ /* [x6 parse] was: byte-at-a-time store with four comparisons per byte. Locate the
token end with 32-byte SIMD compares, then copy the whole token in one memcpy.
Same token: the loop stops at the first byte in {0, ' ', '\n', '\r', '\t'}. */
long st = n, p = n;
int done = 0;
const __m256i sp8 = _mm256_set1_epi8(' '), nl8 = _mm256_set1_epi8('\n');
const __m256i cr8 = _mm256_set1_epi8('\r'), tb8 = _mm256_set1_epi8('\t');
const __m256i zz8 = _mm256_setzero_si256();
while (p + 32 <= sn) {
__m256i v; memcpy(&v, q + p, 32);
__m256i hit = _mm256_or_si256(
_mm256_or_si256(_mm256_cmpeq_epi8(v, sp8), _mm256_cmpeq_epi8(v, nl8)),
_mm256_or_si256(_mm256_cmpeq_epi8(v, cr8),
_mm256_or_si256(_mm256_cmpeq_epi8(v, tb8), _mm256_cmpeq_epi8(v, zz8))));
unsigned msk = (unsigned)_mm256_movemask_epi8(hit);
if (msk) { p += (long)__builtin_ctz(msk); done = 1; break; }
p += 32;
}
if (!done)
while (p < sn && q[p] != 0 && q[p] != ' ' && q[p] != '\n' && q[p] != '\r' && q[p] != '\t') p++;
n = (int)(p - st);
if (n > 0) memcpy(str, q + st, (size_t)n);
}
if (n <= 0) return 0;
memset(str + n, 0, PADB);
// L10[b] = decimal digits of the largest value with bit-length b+1 (= digits(2^(b+1)-1))
for (int b = 0; b < 32; b++) {
unsigned long long hi = (b == 31) ? 0xFFFFFFFFull : ((1ull << (b + 1)) - 1);
int l = 1; unsigned long long pp = 1;
while (pp * 10 <= hi) { pp *= 10; l++; }
L10[b] = (unsigned char)l;
}
for (int d = 0; d < 256; d++) {
uint64_t u = (d < 10) ? ((uint64_t)('0' + d) | ((uint64_t)' ' << 8)) : 0;
PAIR[d] = u;
}
build_T3();
build_t4();
MARK(0);
// ---- build records + MSD count (fused), then MSD scatter by chars c0,c1 ----
memset(CNT, 0, sizeof(CNT));
for (int i = 0; i < n; i++) CNT[bk_dig(i)]++;
{
uint32_t sum = 0;
for (int d = 0; d < 1024; d++) { uint32_t c = CNT[d]; START[d] = sum; CNT[d] = sum; sum += c; }
START[1024] = sum;
}
MARK(1);
for (int i = 0; i < n; i++) {
uint32_t k = key_at(i);
RB[CNT[k >> 20]++] = ((uint64_t)k << 20) | (uint64_t)(uint32_t)i;
}
MARK(2);
{ // oversized bucket => non-random data: bail out to SA-IS
int bad = 0;
for (int d = 0; d < 1024; d++) if (START[d + 1] - START[d] > SCRSZ) { bad = 1; break; }
if (bad) {
g_fallback = 1;
sais<unsigned char>(str, g_sa, n + 1, 256, 0);
for (int x = 0; x < n; x++) g_sa[x] = g_sa[x + 1];
char *o = *opo;
for (int x = 0; x < n - 1; x++) o = wr(o, (uint32_t)g_sa[x] + 1);
o = putu_exact(o, (uint32_t)g_sa[n - 1] + 1);
*o++ = '\n';
if (n > 1) {
for (int x = 1; x < n - 1; x++) o = wrs(o, (uint32_t)lcp_len((uint32_t)g_sa[x - 1], (uint32_t)g_sa[x]));
o = putu_exact(o, (uint32_t)lcp_len((uint32_t)g_sa[n - 2], (uint32_t)g_sa[n - 1]));
}
*o++ = '\n';
*opo = o;
return 1;
}
}
MARK(3);
// ---- per-bucket LSD by (c4,c5) then (c2,c3), FUSED with the walk ----
// The last LSD pass used to write its result back into RB (8 MB of stores + 8 MB of RFO) and
// the walk then read those 8 MB back from DRAM. Both legs are pure overhead: the bucket's
// sorted records are only ever needed in sorted order, and inside the bucket they already sit
// in an L1-resident scratch. So pass2 now targets a second scratch and the walk runs on it
// immediately, per bucket. Runs of equal keys never straddle a bucket boundary (the bucket is
// determined by the top 10 key bits), so the walk's per-bucket bounds are exact.
char *o = *opo;
{
const char *ostart = o;
uint32_t kprev = 0;
int x = 0;
LCPB[0] = 0;
for (int bk = 0; bk < 864; bk++) {
uint32_t lo = START[bk], hi = START[bk + 1];
uint32_t m = hi - lo;
if (m == 0) continue;
uint64_t *W;
if (m == 1) {
W = RB; // single record: already in place, walk it from RB
} else {
uint64_t *cur = RB + lo, *nx = SCR;
{
int sh = 20;
memset(CNT, 0, sizeof(uint32_t) * 864);
for (uint32_t i = 0; i < m; i++) CNT[(uint32_t)(cur[i] >> sh) & 1023]++;
prefix_excl864(CNT);
memset(CNT2, 0, sizeof(uint32_t) * 864);
for (uint32_t i = 0; i < m; i++) {
uint64_t record = cur[i];
nx[CNT[(uint32_t)(record >> sh) & 1023]++] = record;
CNT2[(uint32_t)(record >> 30) & 1023]++;
}
}
prefix_excl864(CNT2);
W = SCR2 - lo; // W[x] is the record at global sorted position x
for (uint32_t i = 0; i < m; i++) {
uint64_t record = nx[i];
W[lo + CNT2[(uint32_t)(record >> 30) & 1023]++] = record;
}
}
const int xe = (int)hi;
/* LANE pragma_lane: THE FAST PATH USED TO EXIT ON THE FIRST REPEATED 6-GRAM AND NEVER RESUME.
`if (k1==k0 || ...) break;` fires as soon as two ADJACENT sorted suffixes share a 30-bit
key, i.e. the first repeated 6-gram -- and a random 1e6 string contains ~1600 of them
(n^2/2/26^6), the first at x=732. So the 4-at-a-time path processed 732 of 1e6 elements
and the per-element run loop did the other 999 268. Restructure: the fast loop is now the
OUTER loop and the run loop handles EXACTLY ONE equal-key run before control returns to it.
Byte-exact: same key, same `kprev` carry, same `g_fallback` bail (the `break` inside the
plain block exits the outer `while (x < n)`), same RUNMAX semantics. */
while (x < xe) {
while (x + 5 <= xe) { // fast path: 4 at a time; requires no run touching x..x+4
uint64_t r0 = W[x], r1 = W[x + 1], r2 = W[x + 2], r3 = W[x + 3];
uint32_t k0 = (uint32_t)(r0 >> 20), k1 = (uint32_t)(r1 >> 20);
uint32_t k2 = (uint32_t)(r2 >> 20), k3 = (uint32_t)(r3 >> 20);
uint32_t k4 = (uint32_t)(W[x + 4] >> 20);
if (k1 == k0 || k2 == k1 || k3 == k2 || k4 == k3) break; // run inside or just past the batch
if (x && k0 == kprev) break; // run straddling the previous batch boundary
if (x) LCPB[x] = (unsigned char)(D5B[__builtin_clz(kprev ^ k0)]);
LCPB[x + 1] = (unsigned char)(D5B[__builtin_clz(k0 ^ k1)]);
LCPB[x + 2] = (unsigned char)(D5B[__builtin_clz(k1 ^ k2)]);
LCPB[x + 3] = (unsigned char)(D5B[__builtin_clz(k2 ^ k3)]);
kprev = k3;
o = pack_emit(o, (uint32_t)(r0 & 0xFFFFFu) + 1);
o = pack_emit(o, (uint32_t)(r1 & 0xFFFFFu) + 1);
o = pack_emit(o, (uint32_t)(r2 & 0xFFFFFu) + 1);
o = pack_emit(o, (uint32_t)(r3 & 0xFFFFFu) + 1);
x += 4;
}
if (g_fallback) break;
{ // ONE equal-key run, then back to the fast path
uint32_t k = (uint32_t)(W[x] >> 20);
int j = x + 1;
while (j < xe && (uint32_t)(W[j] >> 20) == k) j++;
if (j - x > RUNMAX) { g_fallback = 1; o = (char *)ostart; break; }
if (j - x > 1) {
for (int y = x + 1; y < j; y++) {
uint64_t r = W[y];
int z = y - 1;
while (z >= x && suf_less((uint32_t)(r & 0xFFFFFu), (uint32_t)(W[z] & 0xFFFFFu))) { W[z + 1] = W[z]; z--; }
W[z + 1] = r;
}
}
for (int y = x; y < j; y++) {
uint32_t id = (uint32_t)(W[y] & 0xFFFFFu);
uint32_t ky = (uint32_t)(W[y] >> 20);
if (y != 0) LCPB[y] = (ky != kprev) ? (unsigned char)((__builtin_clz(kprev ^ ky) - 2) / 5)
: (unsigned char)lcp_len((uint32_t)(W[y - 1] & 0xFFFFFu), id);
kprev = ky;
if (y == n - 1) o = putu_exact(o, id + 1);
else {
int width;
uint64_t packed = pack_num(id + 1, width);
memcpy(o, &packed, 8); o += width;
}
}
x = j;
} // (E)
} // while (x < xe)
if (g_fallback) break; // bail out of the whole bucket loop
} // bucket loop
if (!g_fallback) {
MARK(5);
*o++ = '\n';
if (n > 1) {
int y = 1;
/* [x6 lcp] eight LCP bytes -> sixteen output bytes per iteration. PAIR[d] for
d < 10 is exactly the u16 ('0'+d) | (' '<<8), i.e. byte-identical to
`_mm_unpacklo_epi8(_mm_add_epi8(d, '0'), ' ')`, and PAIR[d] == 0 for d >= 10,
which is what makes the scalar loop bail out. Any byte >= 9 (unsigned test,
so >= 128 included) therefore falls back to the original scalar loop. */
const __m128i nine = _mm_set1_epi8(9);
const __m256i digit_space = _mm256_set1_epi16(0x2030);
while (y + 16 <= n - 1) {
__m128i v = _mm_loadu_si128((const __m128i *)(LCPB + y));
if (_mm_movemask_epi8(_mm_cmpeq_epi8(_mm_min_epu8(v, nine), nine))) break;
__m256i words = _mm256_cvtepu8_epi16(v);
_mm256_storeu_si256((__m256i *)o, _mm256_add_epi16(words, digit_space));
o += 32; y += 16;
}
for (; y < n - 1; y++) o = wrs(o, LCPB[y]);
o = putu_exact(o, LCPB[n - 1]);
}
*o++ = '\n';
}
}
if (g_fallback) {
sais<unsigned char>(str, g_sa, n + 1, 256, 0);
for (int x = 0; x < n; x++) g_sa[x] = g_sa[x + 1];
char *o = *opo;
for (int x = 0; x < n - 1; x++) o = wr(o, (uint32_t)g_sa[x] + 1);
o = putu_exact(o, (uint32_t)g_sa[n - 1] + 1);
*o++ = '\n';
if (n > 1) {
for (int x = 1; x < n - 1; x++) o = wrs(o, (uint32_t)lcp_len((uint32_t)g_sa[x - 1], (uint32_t)g_sa[x]));
o = putu_exact(o, (uint32_t)lcp_len((uint32_t)g_sa[n - 2], (uint32_t)g_sa[n - 1]));
}
*o++ = '\n';
*opo = o;
return 1;
}
*opo = o;
return 1;
}
#ifdef BENCH
#include <cstdio>
#include <cstdlib>
#define OUTMAX (12u << 20)
static char OUTB[OUTMAX];
int main(int argc, char **argv) {
int n = (argc > 1) ? atoi(argv[1]) : 1000000;
unsigned long long seed = 123456789ULL;
for (int i = 0; i < n; i++) {
seed ^= seed << 13; seed ^= seed >> 7; seed ^= seed << 17;
str[i] = (unsigned char)('a' + (seed % 26));
}
str[n] = 0;
char *o = OUTB;
unsigned long long t0 = __builtin_ia32_rdtsc();
solve_string(str, n, &o);
unsigned long long t1 = __builtin_ia32_rdtsc();
unsigned long long oh = 1469598103934665603ULL;
for (long i = 0; i < (long)(o - OUTB); i++) oh = (oh ^ (unsigned char)OUTB[i]) * 1099511628211ULL;
// same checksum recipe as problems/1006e6/sol_key.cpp (0-based sa_, lcp_ 1..n-1)
unsigned long long sum = 0;
if (!g_fallback) {
for (int i = 0; i < n; i++) sum = sum * 131 + (unsigned)REC_IDX(RB[i]);
for (int i = 1; i < n; i++) sum = sum * 131 + (unsigned)LCPB[i];
} else {
for (int i = 0; i < n; i++) sum = sum * 131 + (unsigned)g_sa[i];
for (int i = 1; i < n; i++) sum = sum * 131 + (unsigned)sum; // fallback: not comparable
}
printf("n=%d total=%llu fb=%d olen=%ld ohash=%llu chk=%llu\n", n, t1 - t0, g_fallback, (long)(o - OUTB), oh, sum);
#ifdef PHASE
{ unsigned long long prev = t0; const char *nm[6]={"parse","build","radix","scanruns","emit+LCPB","lcpline"};
for (int k = 0; k < 6; k++) { printf("%s=%llu ", nm[k], (P[k]-prev)/1000); prev = P[k]; }
printf("perbucket=%llu walk=%llu afterwalk=%llu\n",(P[3]-P[2])/1000,(P[4]-P[3])/1000,(t1-P[5])/1000); }
#endif
if (argc > 2 && argv[2][0] == 'd') { fwrite(OUTB, 1, (size_t)(o - OUTB), stdout); }
return 0;
}
#elif defined(JUDGE)
struct DI {
unsigned long abi;
const char *s; unsigned long sn;
char *o; unsigned long ol; unsigned long os;
char *e; unsigned long el; unsigned long es;
const char *IB; unsigned long IBl;
char *OB; unsigned long OBl;
unsigned long tsc;
} __attribute__((packed));
extern "C" void __libc_start_main(void *m, int argc, char **argv) {
(void)m;
unsigned long *p = (unsigned long *)(argv + argc + 1);
while (*p) p++;
p++;
DI *d = 0;
for (int i = 0; i < 32 && p[0]; i++, p += 2)
if (p[0] == 0x6b637564UL) { d = (DI *)p[1]; break; }
if (d) {
char *o = d->o;
if (solve_string((const unsigned char *)d->s, (long)d->sn, &o)) d->os = (unsigned long)(o - d->o);
}
__asm__ volatile("syscall" ::"a"(60), "D"(0) : "rcx", "r11", "memory");
for (;;);
}
int main() { return 0; }
#endif