代码

// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #105885
//     <https://duck.ac/submission/105885>
//     用途：**本文件正文 = #105885**（`work/a13e_sub2.cpp`，逐字节相同，已核 md5）。
//     1.1 <https://duck.ac/submission/106354> saffah_codex_6s_agg2 —— 当代对手 T。
//     1.2 <https://duck.ac/problem/1002e8> 题目页（n=m=1e8，系数 <10）。
// ======================
// ===== 思路 =====
// **单变量：把 twiddle 的 6 条并发取数流交织成 1 条。**
//
// 病灶（真 gcc 9.3 `-O2` 汇编普查 + 同进程双臂定价）：
//   `row_dif_avx` / `row_dit_avx` 的融合 radix-4 内环每 32 点读 **6 条 32 B twiddle 向量**
//   （t1, ts1, t2, ts2, t3, ts3），它们位于 2 个数组的 6 个相距很远的区段
//   ⇒ 该环**同时有 10 条输入流在飞**（4 条数据流 = 蝶形本征 + 6 条 twiddle）。
//   内环 101 条指令、端口/前端下界 **27.67 拍**（绑定项 = 所有向量指令共享 {p0,p1,p5}：83/3），
//   实测 **≈55 拍**；串行链 ≈30 拍 ⇒ **缺口 ≈25 拍在访存**，其中 twiddle 侧占大头。
//   twiddle 读流量占该环读流量的 **~60%**。
//
// 改法（**一处、语义不变**）：为每个 len 级预建一张**交织表**，把 `(t1,ts1,t2,ts2,t3,ts3)`
//   在同一个 j 上的 6 条向量**相邻摆放**（每 j 前进 192 B）。
//   * 内环 **6 条流 → 1 条顺序流**；**字节数不变、取值不变、使用次序不变**；
//   * 蝶形/下标/输出/主循环 **一字未动**；只有 twiddle 的取址表达式改了；
//   * 表很小（全部 len 级合计 ≤ 256 KB/表），**在 `eng_init` 里建一次**，不进热路径；
//   * `twint_for()` 只对 **C 维**两张表返回交织缓冲；**R 维调用与所有非融合路径原样走旧码**
//     （模板参数 `TI` 编译期折叠 ⇒ R 维与旧码**逐字节相同**）。
//
// 闸门（本地，生产 n = 1e8）：
//   * **逐位相同**：`CHK=7586622636533660264` 与基座一致 ✓
//   * **逐元素对拍**（`2.19.107`：checksum 对方向/置换盲）：把两臂的 `c[0..2n]` 各 800 MB
//     落盘后 `cmp` ⇒ **完全一致** ✓
//   * 建表正确性另有**内建自检**（逐 (len,j,k) 比对交织缓冲与直接取值，>8M 次 0 bad）✓
//
// 定价（同进程交替 A/B，本机 n=5e7；**方向以判题机为准**）：
//   * 形态上界臂（把 6 条流压到 1 条 cache line）**−11.6%**；
//   * 忠实臂（1 条顺序流、同字节、400 MB 足迹）**−13.4%**（4 次同号）；
//   * **真刀**（本文件）：3 次交替配对 **−8.5% / −1.6% / −4.1%（均值 ≈ −4.7%）**。
// ================
#define LGR_  12
#define LGC_  15
#define LGH_  26
#define LGR2_ 13
#define LGC2_ 13
#define EXACT 1
#define XPAD 768
#define NRW_ (1 << LGR_)
#define NCW_ (1 << LGC_)
#define XST_ (NCW_ + XPAD)
#include <cstdint>
#include <cstring>
#ifndef NO_AVX2
#include <immintrin.h>
#define TIK(j, k) ((size_t)((j) >> 3) * 48 + (size_t)(k) * 8)
#endif
typedef uint16_t u16;
typedef uint32_t u32;
#define NR_  (1u << LGR_)
#define NC_  (1u << LGC_)
#define N1_  (1u << (LGR_ + LGC_))
#define NR2_ (1u << LGR2_)
#define NC2_ (1u << LGC2_)
#define N2_  (1u << (LGR2_ + LGC2_))
static u32 t1fR[NR_], t1iR[NR_], t1fC[NC_], t1iC[NC_];
static u16 b1R[NR_], b1C[NC_];
static u32 w1R[2 * NR_], w1C[2 * NC_];
static u32 t2fR[NR2_], t2iR[NR2_], t2fC[NC2_], t2iC[NC2_];
static u16 b2R[NR2_], b2C[NC2_];
static u32 w2R[2 * NR2_], w2C[2 * NC2_];
static u32 t1FI[2 * NC_], t1II[2 * NC_];
static u32 t2FI[2 * NC2_], t2II[2 * NC2_];
static inline const u32 *twint_for(const u32 *tw) {
  if (tw == t1fC) return t1FI;
  if (tw == t1iC) return t1II;
  if (tw == t2fC) return t2FI;
  if (tw == t2iC) return t2II;
  return 0;                       
}
typedef uint64_t u64;
typedef int64_t  i64;
typedef int32_t  i32;
typedef unsigned __int128 u128;
#pragma GCC push_options
#pragma GCC target("avx2")
static inline __attribute__((always_inline)) __m256i ldu256(const void *p) {
  __m256i r; __asm__("vmovdqu %1, %0" : "=x"(r) : "m"(*(const __m256i *)p)); return r;
}
static inline __attribute__((always_inline)) void stu256(void *p, __m256i v) {
  __asm__ volatile("vmovdqu %1, %0" : "=m"(*(__m256i *)p) : "x"(v));
}
static inline __attribute__((always_inline)) void stnt256(void *p, __m256i v) {
  __asm__ volatile("vmovntdq %1, %0" :: "m"(*(__m256i *)p), "x"(v));
}
static inline __attribute__((always_inline)) __m256i vaddm256(__m256i a, __m256i b, __m256i Pv) {
  __m256i s = _mm256_add_epi32(a, b);
  return _mm256_min_epu32(s, _mm256_sub_epi32(s, Pv));
}
static inline __attribute__((always_inline)) __m256i vsubm256(__m256i a, __m256i b, __m256i Pv) {
  __m256i d = _mm256_sub_epi32(a, b);
  return _mm256_min_epu32(d, _mm256_add_epi32(d, Pv));
}
#pragma GCC pop_options
#pragma GCC optimize("O3","rename-registers","live-range-shrinkage","ira-loop-pressure","no-stack-protector")
#define TS 64
static const u32 P = 2013265921u;
static const u32 G = 31u;
static u64 IM;
static u32 TBUF[TS * TS];
static u32 TW[TS * TS];
static inline u32 mm(u32 a, u32 b) {
  u64 z = (u64)a * b;
  u64 q = (u64)(((u128)z * IM) >> 64);
  u64 r = z - q * (u64)P;
  return (u32)(r >= P ? r - P : r);
}
static inline u32 shoup_c(u32 w) { return (u32)(((u64)w << 32) / (u64)P); }
static u32 pw(u32 a, u32 e) {
  u32 r = 1;
  while (e) { if (e & 1) r = (u32)((u64)r * a % P); a = (u32)((u64)a * a % P); e >>= 1; }
  return r;
}
static void build_tab(u32 *tw, int L, u32 w) {
  int off = 0;
  for (int len = 2; len <= L; len <<= 1) {
    u32 wl = pw(w, (u32)(L / len));
    u32 c = 1;
    int h = len >> 1;
    for (int j = 0; j < h; j++) { tw[off + j] = c; c = mm(c, wl); }
    off += h;
  }
}
static void build_br(u16 *br, int L) {
  int bits = 0; while ((1 << bits) < L) bits++;
  for (int i = 0; i < L; i++) {
    int r = 0;
    for (int b = 0; b < bits; b++) if (i & (1 << b)) r |= 1 << (bits - 1 - b);
    br[i] = (u16)r;
  }
}
#ifdef NO_AVX2
#define HAVE_AVX2 0
#else
#define HAVE_AVX2 1
#pragma GCC push_options
#pragma GCC target("avx2")
static inline __m256i vmulhi8(__m256i a, __m256i b) {
  __m256i e = _mm256_mul_epu32(a, b);
  __m256i o = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
  e = _mm256_srli_epi64(e, 32);
  return _mm256_blend_epi32(e, o, 0xAA);
}
static inline __m256i shoup8(__m256i v, __m256i w, __m256i ws) {
  const __m256i Pv = _mm256_set1_epi32((int)P);
  __m256i t = _mm256_mullo_epi32(v, w);
  __m256i q = vmulhi8(v, ws);
  __m256i r = _mm256_sub_epi32(t, _mm256_mullo_epi32(q, Pv));
  return _mm256_min_epu32(r, _mm256_sub_epi32(r, Pv));
}
#define SP(w, ws) _mm256_and_si256(w,M32), _mm256_srli_epi64(w,32), _mm256_and_si256(ws,M32), _mm256_srli_epi64(ws,32)
template<int TI> static void row_dit_avx(u32 *x, int L, const u32 *tw, const u32 *tws, const u32 *IF) {
  __m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
  int off = 0;
  if (L >= 8) {
    const __m256i V4  = _mm256_setr_epi32((int)tw[1], (int)tw[2], 1, 1, (int)tw[1], (int)tw[2], 1, 1);
    const __m256i V4s = _mm256_setr_epi32((int)tws[1], (int)tws[2], 2, 2, (int)tws[1], (int)tws[2], 2, 2);
    const __m256i V8  = _mm256_setr_epi32((int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6], 1, 1, 1, 1);
    const __m256i V8s = _mm256_setr_epi32((int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6], 2, 2, 2, 2);
    const __m256i IDX8  = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
    const __m256i IDXQ4 = _mm256_setr_epi32(0, 1, 0, 1, 4, 5, 4, 5);
    const __m256i IDXQ8 = _mm256_setr_epi32(0, 1, 2, 3, 0, 1, 2, 3);
    for (int i = 0; i < L; i += 8) {
      __m256i v = ldu256((const void *)((x + i)));
      {  
        __m256i y = _mm256_shuffle_epi32(v, 0xB1);        
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_sub_epi32(y, v);               
        d = _mm256_add_epi32(d, _mm256_and_si256(_mm256_cmpgt_epi32(Z, d), Pv));
        v = _mm256_blend_epi32(sa, d, 0xAA);
      }
      {  
        __m256i y = _mm256_shuffle_epi32(v, 0x4E);        
        y = shoup8(y, V4, V4s);
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_min_epu32(_mm256_sub_epi32(v, y), _mm256_add_epi32(_mm256_sub_epi32(v, y), Pv));
        d = _mm256_permutevar8x32_epi32(d, IDXQ4);        
        v = _mm256_blend_epi32(sa, d, 0xCC);
      }
      {  
        __m256i y = _mm256_permutevar8x32_epi32(v, IDX8); 
        y = shoup8(y, V8, V8s);
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_min_epu32(_mm256_sub_epi32(v, y), _mm256_add_epi32(_mm256_sub_epi32(v, y), Pv));
        d = _mm256_permutevar8x32_epi32(d, IDXQ8);        
        v = _mm256_blend_epi32(sa, d, 0xF0);
      }
      stu256((void *)((x + i)), v);
    }
    off = 7;
    int len = 16;
    size_t ioff = 0;
    for (; (len << 1) <= L; len <<= 2) {
      int h = len >> 1;                                   
      const u32 *t1 = tw + off,  *ts1 = tws + off;        
      const u32 *t2 = t1 + h,    *ts2 = ts1 + h;          
      const u32 *t3 = t2 + h,    *ts3 = ts2 + h;          
      off += 3 * h;
      for (int i = 0; i < L; i += (len << 1)) {
        u32 *p = x + i;
        for (int j = 0; j < h; j += 8) {
          __m256i A = ldu256((const void *)(p + j));
          __m256i B = ldu256((const void *)(p + j + h));
          __m256i C = ldu256((const void *)(p + j + 2 * h));
          __m256i D = ldu256((const void *)(p + j + 3 * h));
          __m256i w1 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 0)) : (t1 + j))), ws1 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 1)) : (ts1 + j)));
          __m256i Bm = shoup8(B, w1, ws1), Dm = shoup8(D, w1, ws1);
          __m256i Ap = vaddm256(A, Bm, Pv), Bp = vsubm256(A, Bm, Pv);
          __m256i Cp = vaddm256(C, Dm, Pv), Dp = vsubm256(C, Dm, Pv);
          __m256i w2 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 2)) : (t2 + j))), ws2 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 3)) : (ts2 + j)));
          __m256i w3 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 4)) : (t3 + j))), ws3 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 5)) : (ts3 + j)));
          __m256i Cm = shoup8(Cp, w2, ws2), Dn = shoup8(Dp, w3, ws3);
          stu256((void *)(p + j),         vaddm256(Ap, Cm, Pv));
          stu256((void *)(p + j + 2 * h), vsubm256(Ap, Cm, Pv));
          stu256((void *)(p + j + h),     vaddm256(Bp, Dn, Pv));
          stu256((void *)(p + j + 3 * h), vsubm256(Bp, Dn, Pv));
        }
      }
      ioff += TIK(h, 0);
    }
    if (len <= L) {                                       
      int h = len >> 1;
      const u32 *t = tw + off, *ts = tws + off; off += h;
      for (int i = 0; i < L; i += len) {
        u32 *p = x + i, *q = x + i + h;
        for (int j = 0; j < h; j += 8) {
          __m256i u = ldu256((const void *)(p + j));
          __m256i v = shoup8(ldu256((const void *)(q + j)),
                             ldu256((const void *)(t + j)), ldu256((const void *)(ts + j)));
          stu256((void *)(p + j), vaddm256(u, v, Pv));
          stu256((void *)(q + j), vsubm256(u, v, Pv));
        }
      }
    }
  } else {
    for (int len = 2; len <= L; len <<= 1) {
      int h = len >> 1;
      const u32 *t = tw + off; off += h;
      for (int i = 0; i < L; i += len) {
        u32 *p = x + i, *q = x + i + h;
        for (int j = 0; j < h; j++) {
          u32 u = p[j], v = mm(q[j], t[j]);
          u32 sa = u + v; if (sa >= P) sa -= P;
          u32 d = u + P - v; if (d >= P) d -= P;
          p[j] = sa; q[j] = d;
        }
      }
    }
  }
}
template<int TI> static void row_dif_avx(u32 *x, int L, const u32 *tw, const u32 *tws, const u32 *IF) {
  __m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
  if (L >= 8) {
    int len = L;
    size_t ioff = 0;
    for (; len >= 32; len >>= 2) {
      int h = len >> 1, h2 = len >> 2;
      const u32 *t1 = tw + h - 1,  *ts1 = tws + h - 1;      
      const u32 *t2 = t1 + h2,    *ts2 = ts1 + h2;          
      const u32 *t3 = tw + h2 - 1, *ts3 = tws + h2 - 1;     
      for (int i = 0; i < L; i += len) {
        u32 *p = x + i;
        for (int j = 0; j < h2; j += 8) {
          __m256i v0 = ldu256((const void *)((p + j)));
          __m256i v1 = ldu256((const void *)((p + j + h2)));
          __m256i v2 = ldu256((const void *)((p + j + h)));
          __m256i v3 = ldu256((const void *)((p + j + h + h2)));
          __m256i s02 = _mm256_min_epu32(_mm256_add_epi32(v0, v2), _mm256_sub_epi32(_mm256_add_epi32(v0, v2), Pv));
          __m256i d02 = _mm256_sub_epi32(v0, v2);
          __m256i s13 = _mm256_min_epu32(_mm256_add_epi32(v1, v3), _mm256_sub_epi32(_mm256_add_epi32(v1, v3), Pv));
          __m256i d13 = _mm256_sub_epi32(v1, v3);
          d02 = _mm256_min_epu32(d02, _mm256_add_epi32(d02, Pv));
          d13 = _mm256_min_epu32(d13, _mm256_add_epi32(d13, Pv));
          d02 = shoup8(d02, ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 0)) : (t1 + j))), ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 1)) : (ts1 + j))));
          d13 = shoup8(d13, ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 2)) : (t2 + j))), ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 3)) : (ts2 + j))));
          __m256i w3 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 4)) : (t3 + j)));
          __m256i ws3 = ldu256((const void *)((TI) ? (IF + ioff + TIK(j, 5)) : (ts3 + j)));
          __m256i ra = _mm256_min_epu32(_mm256_add_epi32(s02, s13), _mm256_sub_epi32(_mm256_add_epi32(s02, s13), Pv));
          __m256i rb = _mm256_min_epu32(_mm256_sub_epi32(s02, s13), _mm256_add_epi32(_mm256_sub_epi32(s02, s13), Pv));
          rb = shoup8(rb, w3, ws3);
          __m256i rc = _mm256_min_epu32(_mm256_add_epi32(d02, d13), _mm256_sub_epi32(_mm256_add_epi32(d02, d13), Pv));
          __m256i rd = _mm256_min_epu32(_mm256_sub_epi32(d02, d13), _mm256_add_epi32(_mm256_sub_epi32(d02, d13), Pv));
          rd = shoup8(rd, w3, ws3);
          stu256((void *)((p + j)), ra);
          stu256((void *)((p + j + h2)), rb);
          stu256((void *)((p + j + h)), rc);
          stu256((void *)((p + j + h + h2)), rd);
        }
      }
      ioff += TIK(h2, 0);
    }
    if (len >= 16) {                              
      int h = len >> 1;
      const u32 *t = tw + h - 1, *ts = tws + h - 1;
      for (int i = 0; i < L; i += len) {
        u32 *p = x + i, *q = x + i + h;
        for (int j = 0; j < h; j += 8) {
          __m256i u = ldu256((const void *)((p + j)));
          __m256i v = ldu256((const void *)((q + j)));
          __m256i w = ldu256((const void *)((t + j)));
          __m256i ws = ldu256((const void *)((ts + j)));
          __m256i sa = _mm256_min_epu32(_mm256_add_epi32(u, v), _mm256_sub_epi32(_mm256_add_epi32(u, v), Pv));
          __m256i d = _mm256_min_epu32(_mm256_sub_epi32(u, v), _mm256_add_epi32(_mm256_sub_epi32(u, v), Pv));
          d = shoup8(d, w, ws);
          stu256((void *)((p + j)), sa);
          stu256((void *)((q + j)), d);
        }
      }
    }
    const __m256i V8  = _mm256_setr_epi32(1, 1, 1, 1, (int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6]);
    const __m256i V8s = _mm256_setr_epi32(2, 2, 2, 2, (int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6]);
    const __m256i V4  = _mm256_setr_epi32(1, 1, (int)tw[1], (int)tw[2], 1, 1, (int)tw[1], (int)tw[2]);
    const __m256i V4s = _mm256_setr_epi32(2, 2, (int)tws[1], (int)tws[2], 2, 2, (int)tws[1], (int)tws[2]);
    const __m256i IDX8 = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
    for (int i = 0; i < L; i += 8) {
      __m256i v = ldu256((const void *)((x + i)));
      {  
        __m256i y = _mm256_permutevar8x32_epi32(v, IDX8);
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_sub_epi32(y, v);            
        d = _mm256_add_epi32(d, _mm256_and_si256(_mm256_cmpgt_epi32(Z, d), Pv));
        d = shoup8(d, V8, V8s);
        v = _mm256_blend_epi32(sa, d, 0xF0);
      }
      {  
        __m256i y = _mm256_shuffle_epi32(v, 0x4E);
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_sub_epi32(y, v);            
        d = _mm256_min_epu32(d, _mm256_add_epi32(d, Pv));
        d = shoup8(d, V4, V4s);
        v = _mm256_blend_epi32(sa, d, 0xCC);
      }
      {  
        __m256i y = _mm256_shuffle_epi32(v, 0xB1);
        __m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
        __m256i d = _mm256_min_epu32(_mm256_sub_epi32(y, v), _mm256_add_epi32(_mm256_sub_epi32(y, v), Pv));
        v = _mm256_blend_epi32(sa, d, 0xAA);
      }
      stu256((void *)((x + i)), v);
    }
  } else {
    for (int len = L; len >= 2; len >>= 1) {
      int h = len >> 1;
      const u32 *t = tw + h - 1;
      for (int i = 0; i < L; i += len) {
        u32 *p = x + i, *q = x + i + h;
        for (int j = 0; j < h; j++) {
          u32 u = p[j], v = q[j];
          u32 sa = u + v; if (sa >= P) sa -= P;
          u32 d = u + P - v; if (d >= P) d -= P;
          p[j] = sa; q[j] = mm(d, t[j]);
        }
      }
    }
  }
}
static inline __m256i shoup_of(__m256i w) {
  const __m256i ML = _mm256_set1_epi32(572662301);
  __m256i lo = _mm256_mul_epu32(w, ML);
  __m256i hi = _mm256_mul_epu32(_mm256_srli_epi64(w, 32), ML);
  __m256i m = _mm256_blend_epi32(_mm256_srli_epi64(lo, 32), _mm256_slli_epi64(_mm256_srli_epi64(hi, 32), 32), 0xAA);
  return _mm256_add_epi32(_mm256_add_epi32(w, w), m);
}
static inline void t8x8(__m256i *r) {
  __m256i t0 = _mm256_unpacklo_epi32(r[0], r[1]), t1 = _mm256_unpackhi_epi32(r[0], r[1]);
  __m256i t2 = _mm256_unpacklo_epi32(r[2], r[3]), t3 = _mm256_unpackhi_epi32(r[2], r[3]);
  __m256i t4 = _mm256_unpacklo_epi32(r[4], r[5]), t5 = _mm256_unpackhi_epi32(r[4], r[5]);
  __m256i t6 = _mm256_unpacklo_epi32(r[6], r[7]), t7 = _mm256_unpackhi_epi32(r[6], r[7]);
  __m256i u0 = _mm256_unpacklo_epi64(t0, t2), u1 = _mm256_unpackhi_epi64(t0, t2);
  __m256i u2 = _mm256_unpacklo_epi64(t1, t3), u3 = _mm256_unpackhi_epi64(t1, t3);
  __m256i u4 = _mm256_unpacklo_epi64(t4, t6), u5 = _mm256_unpackhi_epi64(t4, t6);
  __m256i u6 = _mm256_unpacklo_epi64(t5, t7), u7 = _mm256_unpackhi_epi64(t5, t7);
  r[0] = _mm256_permute2x128_si256(u0, u4, 0x20); r[1] = _mm256_permute2x128_si256(u1, u5, 0x20);
  r[2] = _mm256_permute2x128_si256(u2, u6, 0x20); r[3] = _mm256_permute2x128_si256(u3, u7, 0x20);
  r[4] = _mm256_permute2x128_si256(u0, u4, 0x31); r[5] = _mm256_permute2x128_si256(u1, u5, 0x31);
  r[6] = _mm256_permute2x128_si256(u2, u6, 0x31); r[7] = _mm256_permute2x128_si256(u3, u7, 0x31);
}
static inline u32 mont1(u32 a, u32 b);
static inline __m256i mont8(__m256i a, __m256i b);
static void pointwise_mul_avx(u32 *A, const u32 *B, u32 n, u32 sc) {
  (void)sc;
  u32 i = 0;
  for (; i + 8 <= n; i += 8) {
    __m256i a = ldu256((const void *)(A + i));
    __m256i b = ldu256((const void *)(B + i));
    stu256((void *)(A + i), mont8(a, b));
  }
  for (; i < n; i++) A[i] = mont1(A[i], B[i]);
}
static u32 NPI = 0;                              
static const u32 R32 = (u32)((1ull << 32) % P);
static inline u32 mont1(u32 a, u32 b) {
  u64 z = (u64)a * (u64)b;
  u32 m = (u32)((u64)(u32)z * (u64)NPI);
  u64 t = (z + (u64)m * (u64)P) >> 32;
  return (u32)(t >= P ? t - P : t);
}
static inline __m256i mont8(__m256i a, __m256i b) {
  __m256i P64 = _mm256_set1_epi64x((long long)P), NP = _mm256_set1_epi64x((long long)NPI);
  __m256i Pv = _mm256_set1_epi32((int)P);
  __m256i ahi = _mm256_srli_epi64(a, 32), bhi = _mm256_srli_epi64(b, 32);
  __m256i tlo = _mm256_mul_epu32(a, b), thi = _mm256_mul_epu32(ahi, bhi);
  __m256i e2 = _mm256_mul_epu32(tlo, NP), o2 = _mm256_mul_epu32(thi, NP);
  __m256i e4 = _mm256_mul_epu32(e2, P64), o4 = _mm256_mul_epu32(o2, P64);
  __m256i e5 = _mm256_add_epi64(tlo, e4), o5 = _mm256_add_epi64(thi, o4);
  __m256i ue = _mm256_srli_epi64(e5, 32), uo = _mm256_srli_epi64(o5, 32);
  __m256i u = _mm256_blend_epi32(ue, _mm256_slli_epi64(uo, 32), 0xAA);
  return _mm256_min_epu32(u, _mm256_sub_epi32(u, Pv));
}
static void scale_row_mont(u32 *dst, const u32 *src, u32 kM) {
  const __m256i kv = _mm256_set1_epi32((int)kM);   
  for (int v = 0; v < TS; v += 8)
    stu256((void *)((dst + v)), mont8(ldu256((const void *)((src + v))), kv));
}
static void twiddle_rows_mont(u32 *TW, const u32 *RVM) {
  for (int u = 1; u < TS; u++) {
    const u32 *prv = TW + (u - 1) * TS;
    u32 *cur = TW + u * TS;
    for (int v = 0; v < TS; v += 8)
      stu256((void *)((cur + v)), mont8(ldu256((const void *)((prv + v))), ldu256((const void *)((RVM + v)))));
  }
}
static void scale_row_avx(u32 *dst, const u32 *src, u32 k) {
  const __m256i kv = _mm256_set1_epi32((int)k), ks = _mm256_set1_epi32((int)shoup_c(k));
  for (int v = 0; v < TS; v += 8)
    stu256((void *)((dst + v)), shoup8(ldu256((const void *)((src + v))), kv, ks));
}
static void twiddle_rows_avx(u32 *TW, const u32 *RV) {
  __m256i RVS[TS / 8];
  for (int v = 0; v < TS; v += 8) RVS[v >> 3] = shoup_of(ldu256((const void *)((RV + v))));
  for (int u = 1; u < TS; u++) {
    const u32 *prv = TW + (u - 1) * TS;
    u32 *cur = TW + u * TS;
    for (int v = 0; v < TS; v += 8) {
      __m256i r = ldu256((const void *)((RV + v)));
      stu256((void *)((cur + v)), shoup8(ldu256((const void *)((prv + v))), r, RVS[v >> 3]));
    }
  }
}
static void corner_tile_avx(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut) {
  for (int u0 = 0; u0 < TS; u0 += 64) {
    u32 cb[TS * 64] __attribute__((aligned(64)));
    for (int u = 0; u < 64; u += 8) {
      for (int v0 = 0; v0 < TS; v0 += 8) {
        __m256i r[8];
        for (int k = 0; k < 8; k++) {
          const u32 *s = in + (u64)(A0 + u0 + u + k) * strIn + B0 + v0;
          const u32 *tw = TW + (u0 + u + k) * TS + v0;
          __m256i x = ldu256((const void *)(s));
          __m256i w = ldu256((const void *)(tw));      
          r[k] = mont8(x, w);
        }
        t8x8(r);
        for (int j = 0; j < 8; j++)
          stu256((void *)(cb + (u64)(v0 + j) * 64 + u), r[j]);
      }
    }
    for (int v = 0; v < TS; v++) {
      u32 *d = out + (u64)(B0 + v) * strOut + A0 + u0;
      for (int q = 0; q < 64; q += 8)
        stnt256((void *)(d + q), ldu256((const void *)(cb + (u64)v * 64 + q)));
    }
  }
}
static void load_tile_avx(const u32 *src, u32 *S, int r0, int c0, int R, int C, int ctr, int allzero) {
  for (int rb0 = 0; rb0 < TS; rb0 += 64) {
    u32 buf[TS * 64] __attribute__((aligned(64)));
    for (int rb = 0; rb < 64; rb += 8) {
      for (int cb0 = 0; cb0 < TS; cb0 += 8) {
        __m256i r[8];
        if (allzero) {
          for (int k = 0; k < 8; k++) r[k] = _mm256_setzero_si256();
        } else {
          for (int k = 0; k < 8; k++) {
            const u32 *s = src + (u64)(r0 + rb0 + rb + k) * C + c0 + cb0;
            __m256i x = ldu256((const void *)(s));
            if (ctr) {
              __m256i y = _mm256_sub_epi32(x, _mm256_set1_epi32(4));
              __m256i z = _mm256_add_epi32(x, _mm256_set1_epi32((int)(P - 4)));
              x = _mm256_min_epu32(y, z);
            }
            r[k] = x;
          }
        }
        t8x8(r);
        for (int j = 0; j < 8; j++)
          stu256((void *)(buf + (u64)(cb0 + j) * 64 + rb), r[j]);
      }
    }
    for (int m = 0; m < TS; m++) {
      u32 *d = S + (u64)(c0 + m) * R + r0 + rb0;
      for (int q = 0; q < 64; q += 8)
        stnt256((void *)(d + q), ldu256((const void *)(buf + (u64)m * 64 + q)));
    }
  }
}
static void untr_tile_avx(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut) {
  for (int u0 = 0; u0 < TS; u0 += 8) {
    for (int v0 = 0; v0 < TS; v0 += 8) {
      __m256i r[8];
      for (int k = 0; k < 8; k++)
        r[k] = ldu256((const void *)((in + (u64)(A0 + u0 + k) * strIn + B0 + v0)));
      t8x8(r);
      for (int k = 0; k < 8; k++)
        stu256((void *)((out + (u64)(B0 + v0 + k) * strOut + A0 + u0)), r[k]);
    }
  }
}
static void untr_tile_nt(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut) {
  for (int u0 = 0; u0 < TS; u0 += 64) {
    u32 buf[TS * 64] __attribute__((aligned(64)));
    for (int u = 0; u < 64; u += 8) {
      for (int v = 0; v < TS; v += 8) {
        __m256i r[8];
        for (int k = 0; k < 8; k++)
          r[k] = ldu256((const void *)(in + (u64)(A0 + u0 + u + k) * strIn + B0 + v));
        t8x8(r);
        for (int j = 0; j < 8; j++)
          stu256((void *)(buf + (u64)(v + j) * 64 + u), r[j]);
      }
    }
    for (int v = 0; v < TS; v++) {
      u32 *d = out + (u64)(B0 + v) * strOut + A0 + u0;
      for (int q = 0; q < 64; q += 8)
        stnt256((void *)(d + q), ldu256((const void *)(buf + (u64)v * 64 + q)));
    }
  }
}
static void build_tab2(u32 *tw, u32 *tws, int L, u32 w) {
  int off = 0;
  for (int len = 2; len <= L; len <<= 1) {
    u32 wl = pw(w, (u32)(L / len));
    u32 c = 1;
    int h = len >> 1;
    for (int j = 0; j < h; j++) { tw[off + j] = c; tws[off + j] = shoup_c(c); c = mm(c, wl); }
    off += h;
  }
}
#pragma GCC pop_options
#endif
static void row_dit(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
  { const u32 *IF = twint_for(tw);
    if (IF) row_dit_avx<1>(x, L, tw, tws, IF); else row_dit_avx<0>(x, L, tw, tws, 0); }
  return;
#else
  (void)tws;
#endif
  int off = 0;
  for (int len = 2; len <= L; len <<= 1) {
    int h = len >> 1;
    const u32 *t = tw + off; off += h;
    for (int i = 0; i < L; i += len) {
      u32 *p = x + i, *q = x + i + h;
      for (int j = 0; j < h; j++) {
        u32 u = p[j], v = mm(q[j], t[j]);
        u32 s = u + v; if (s >= P) s -= P;
        u32 d = u + P - v; if (d >= P) d -= P;
        p[j] = s; q[j] = d;
      }
    }
  }
}
static void row_dif(u32 *x, int L, const u32 *tw, const u32 *tws);
static void pointwise_dit(u32 *A, u32 *B, u32 n, u32 sc,
                          const u32 *tfC, const u32 *wsCf,
                          const u32 *tiC, const u32 *wsCi) {
  (void)n;
  for (int k0 = 0; k0 < (int)NRW_; k0++) {
    u32 *Aa = A + (u64)k0 * (u64)XST_;
    u32 *Bb = B + (u64)k0 * (u64)XST_;
    row_dif(Aa, NCW_, tfC, wsCf);   
    row_dif(Bb, NCW_, tfC, wsCf);
#if HAVE_AVX2
    pointwise_mul_avx(Aa, Bb, NCW_, sc);
#else
    for (u32 i = 0; i < NCW_; i++) Aa[i] = mm(Aa[i], Bb[i]);
#endif
    row_dit(Aa, NCW_, tiC, wsCi);
  }
}
static void row_dif(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
  { const u32 *IF = twint_for(tw);
    if (IF) row_dif_avx<1>(x, L, tw, tws, IF); else row_dif_avx<0>(x, L, tw, tws, 0); }
  return;
#else
  (void)tws;
#endif
  for (int len = L; len >= 2; len >>= 1) {
    int h = len >> 1;
    const u32 *t = tw + h - 1;
    for (int i = 0; i < L; i += len) {
      u32 *p = x + i, *q = x + i + h;
      for (int j = 0; j < h; j++) {
        u32 u = p[j], v = q[j];
        u32 s = u + v; if (s >= P) s -= P;
        u32 d = u + P - v; if (d >= P) d -= P;
        p[j] = s; q[j] = mm(d, t[j]);
      }
    }
  }
}
static inline void bitrev_pass(u32 *x, int L, const u16 *br) {
  for (int i = 0; i < L; i++) {
    int k = br[i];
    if (i < k) { u32 t = x[i]; x[i] = x[k]; x[k] = t; }
  }
}
struct Eng {
  int R, C, N;
  u32 w;
  u32 *tfR; u32 *tiR; u32 *tfC; u32 *tiC;
  u32 *wsR; u32 *wsC;
  u16 *brR; u16 *brC;
};
static void load_T(Eng &e, const u32 *src, u32 nsrc, u32 *S, int ctr) {
  const int R = e.R, C = e.C;
#if HAVE_AVX2
  if (nsrc <= (u32)R * (u32)C) {
    for (int r0 = 0; r0 < R; r0 += TS) {
      u64 tilerow = (u64)r0 * C;
      for (int c0 = 0; c0 < C; c0 += TS) {
        u64 firstidx = tilerow + c0;
        u64 lastidx  = (u64)(r0 + TS - 1) * C + c0 + TS - 1;
        if (lastidx < nsrc)        { load_tile_avx(src, S, r0, c0, R, C, ctr, 0); continue; }
        if (firstidx >= nsrc)      { load_tile_avx(src, S, r0, c0, R, C, ctr, 1); continue; }
        for (int rb = 0; rb < TS; rb++) {
          u64 base = (u64)(r0 + rb) * C + c0;
          for (int cb = 0; cb < TS; cb++) {
            u64 id = base + cb;
            u32 v = 0;
            if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
            TBUF[cb * TS + rb] = v;
          }
        }
        for (int cb = 0; cb < TS; cb++) {
          u32 *d = S + (u64)(c0 + cb) * R + r0;
          const u32 *s = TBUF + cb * TS;
          for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
        }
      }
    }
    return;
  }
#endif
  for (int r0 = 0; r0 < R; r0 += TS) {
    for (int c0 = 0; c0 < C; c0 += TS) {
      u64 lastidx = (u64)(r0 + TS - 1) * C + (c0 + TS - 1);
      if (lastidx < nsrc) {
        for (int rb = 0; rb < TS; rb++) {
          const u32 *s = src + (u64)(r0 + rb) * C + c0;
          for (int cb = 0; cb < TS; cb++) {
            u32 v = s[cb];
            if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u));
            TBUF[cb * TS + rb] = v;
          }
        }
      } else {
        for (int rb = 0; rb < TS; rb++) {
          u64 base = (u64)(r0 + rb) * C + c0;
          for (int cb = 0; cb < TS; cb++) {
            u64 id = base + cb;
            u32 v = 0;
            if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
            TBUF[cb * TS + rb] = v;
          }
        }
      }
      for (int cb = 0; cb < TS; cb++) {
        u32 *d = S + (u64)(c0 + cb) * R + r0;
        const u32 *s = TBUF + cb * TS;
        for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
      }
    }
  }
}
static u32 REV6[TS];
static u32 g_inv_scale = 1u;
static u32 POWT[1024], POWT64[1024];
static void build_rev6(void) {
  for (int i = 0; i < TS; i++) {
    int r = 0;
    for (int b = 0; b < 6; b++) if (i & (1 << b)) r |= 1 << (5 - b);
    REV6[i] = (u32)r;
  }
}
static inline u32 revbits(u32 x, int lg) {
  u32 r = 0;
  for (int b = 0; b < lg; b++) if (x & (1u << b)) r |= 1u << (lg - 1 - b);
  return r;
}
static void corner_generic(Eng &e, const u32 *in, u32 *out, int nA, int nB, int strIn, int strOut, u32 root, int perm, int fuse_row) {
  u32 wv[TS];
  { u32 c = 1; for (int v = 0; v < TS; v++) { wv[v] = c; c = mm(c, root); } }
  const u32 w64 = pw(root, TS);
  int lgB = 0, lgA = 0;
  while ((1 << lgB) < nB) lgB++;
  while ((1 << lgA) < nA) lgA++;
  const int S1 = nB / TS, S2 = nA / TS;          
  const int HB = lgB - 6, HA = lgA - 6;          
  u32 RV2[TS], RV2b[TS], PW2[TS], cu[TS], cu64[TS], base[TS];
  if (perm == 1) {
    u32 g = pw(root, (u32)S1);                   
    { u32 gM = mm(g, R32), c = R32; for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; c = mont1(c, gM); } }
    { u32 rootM = mm(root, R32), p = R32; for (int t = 0; t < S1; t++) { POWT[t] = p; p = mont1(p, rootM); } }
  } else if (perm == 2) {
    u32 g = pw(root, (u32)S2), g64 = pw(g, 64);
    const u32 rootM2 = mm(root, R32), gM = mm(g, R32), g64M = mm(g64, R32);
    { u32 c = R32, c64 = R32;
      for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; RV2b[REV6[t]] = c64; c = mont1(c, gM); c64 = mont1(c64, g64M); } }
    { u32 p = R32; for (int t = 0; t < S2; t++) { POWT[t] = p; p = mont1(p, rootM2); } }
    { u32 q = pw(root, 64), qM = mm(q, R32), p = R32; for (int t = 0; t < S2; t++) { POWT64[t] = p; p = mont1(p, qM); } }
  }
  u32 wr = 1;                                   
  const u32 INVSCALE = g_inv_scale;             
  for (int A0 = 0; A0 < nA; A0 += TS) {
    if (fuse_row) {
      u32 *Sw = (u32 *)in + (u64)A0 * (u64)strIn;
      for (int u = 0; u < TS; u++) row_dif(Sw + (u64)u * (u64)strIn, nB, e.tfR, e.wsR);
    }
    u32 wr64 = pw(wr, TS);                      
    if (perm == 1) {
      u32 g = pw(wr, (u32)S1), gM2 = mm(g, R32), c = R32;
      for (int t = 0; t < TS; t++) { PW2[REV6[t]] = c; c = mont1(c, gM2); }
    } else if (perm == 2) {
      int hi = A0 / TS;
      u32 rH = POWT[revbits((u32)hi, HA)], r64H = POWT64[revbits((u32)hi, HA)];
      const u32 INVSCALE_M = mm(INVSCALE, R32);
      for (int u = 0; u < TS; u++) { cu[u] = mont1(RV2[u], rH); cu64[u] = mont1(RV2b[u], r64H); base[u] = INVSCALE_M; }
    } else {
      u32 c = 1; for (int v = 0; v < TS; v++) { PW2[v] = c; c = mm(c, wr); }
    }
    u32 Vloc = 1;                               
    u32 Vr = 1;                                 
    for (int B0 = 0; B0 < nB; B0 += TS) {
      if (perm == 1) {
        u32 H = revbits((u32)(B0 / TS), HB);
        u32 wrH = pw(wr, H), rootH = POWT[H];
        u32 RV[TS];
#if HAVE_AVX2
        scale_row_mont(RV, RV2, rootH);                                 
        scale_row_mont(TW, PW2, mm(wrH, R32));                          
        twiddle_rows_mont(TW, RV);
#else
        for (int v = 0; v < TS; v++) RV[v] = mont1(RV2[v], rootH);
        for (int v = 0; v < TS; v++) TW[v] = mont1(PW2[v], mm(wrH, R32));
        for (int u = 1; u < TS; u++) {
          const u32 *prv = TW + (u - 1) * TS;
          u32 *cur = TW + u * TS;
          for (int v = 0; v < TS; v++) cur[v] = mont1(prv[v], RV[v]);
        }
#endif
      } else if (perm == 2) {
        if (B0 == 0) {
          for (int u = 0; u < TS; u++) {
            u32 *cur = TW + u * TS;
            cur[0] = base[u];
            for (int v = 1; v < TS; v++) cur[v] = mont1(cur[v - 1], cu[u]);
          }
        } else {
          for (int u = 0; u < TS; u++) {
            u32 *cur = TW + u * TS;
#if HAVE_AVX2
            scale_row_mont(cur, cur, cu64[u]);
#else
            u32 d = cu64[u];
            for (int v = 0; v < TS; v++) cur[v] = mont1(cur[v], d);
#endif
          }
        }
        for (int u = 0; u < TS; u++) base[u] = mont1(base[u], cu64[u]);
      } else {
        u32 RV[TS];
        for (int v = 0; v < TS; v++) RV[v] = mm(Vloc, wv[v]);
        for (int v = 0; v < TS; v++) TW[v] = mm(Vr, PW2[v]);
        for (int u = 1; u < TS; u++) {
          const u32 *prv = TW + (u - 1) * TS;
          u32 *cur = TW + u * TS;
          for (int v = 0; v < TS; v++) cur[v] = mm(prv[v], RV[v]);
        }
      }
#if HAVE_AVX2
      corner_tile_avx(in, out, A0, B0, strIn, strOut);
#else
      for (int u = 0; u < TS; u++) {
        const u32 *s = in + (u64)(A0 + u) * strIn + B0;
        const u32 *tw = TW + u * TS;
        for (int v = 0; v < TS; v++) TBUF[u * TS + v] = mont1(s[v], tw[v]);
      }
      for (int v = 0; v < TS; v++) {
        u32 *d = out + (u64)(B0 + v) * strOut + A0;
        for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
      }
#endif
      Vr = mm(Vr, wr64);
      Vloc = mm(Vloc, w64);
    }
    wr = mm(wr, w64);
  }
}
static inline void corner_fwd(Eng &e, const u32 *S, u32 *X, u32 root, int fuse_row) {
  corner_generic(e, S, X, e.C, e.R, e.R, XST_, root, 1, fuse_row);   
}
static inline void corner_inv(Eng &e, const u32 *X, u32 *S, u32 root) {
  corner_generic(e, X, S, e.R, e.C, XST_, e.R, root, 2, 0);   
}
static void untranspose(Eng &e, const u32 *S, u32 *dst) {
  const int R = e.R, C = e.C;
#if HAVE_AVX2
  for (int c0 = 0; c0 < C; c0 += TS) {        
    { 
      u32 *Sw = (u32 *)S + (u64)c0 * (u64)R;
      for (int u = 0; u < TS; u++) row_dit(Sw + (u64)u * (u64)R, R, e.tiR, e.wsR + e.R); }
    for (int r0 = 0; r0 < R; r0 += TS)        
      untr_tile_nt(S, dst, c0, r0, R, XST_);   
  }
  return;
#endif
  for (int c0 = 0; c0 < C; c0 += TS) {          
    for (int r0 = 0; r0 < R; r0 += TS) {        
      for (int u = 0; u < TS; u++) {
        const u32 *s = S + (u64)(c0 + u) * R + r0;
        u32 *d = TBUF + u * TS;
        for (int v = 0; v < TS; v++) d[v] = s[v];
      }
      for (int v = 0; v < TS; v++) {
        u32 *d = dst + (u64)(r0 + v) * (u64)XST_ + c0;   
        for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
      }
    }
  }
}
static void untr_window(Eng &e, const u32 *S, u32 *dst, long lo, long hi) {
  const int R = e.R, C = e.C;
  for (int c0 = 0; c0 < C; c0 += TS) {
    { 
      u32 *Sw = (u32 *)S + (u64)c0 * (u64)R;
      for (int u = 0; u < TS; u++) row_dit(Sw + (u64)u * (u64)R, R, e.tiR, e.wsR + e.R); }
    for (int r0 = 0; r0 < R; r0 += TS) {
      long tlo = (long)r0 * C + c0;
      long thi = (long)(r0 + TS - 1) * C + c0 + TS - 1;
      if (thi < lo || tlo > hi) continue;
      if (tlo >= lo && thi <= hi) { untr_tile_avx(S, dst, c0, r0, R, C); continue; }
      for (int u = 0; u < TS; u++) {
        long b = (long)r0 * C + c0 + u;
        for (int v = 0; v < TS; v++) {
          long idx = b + (long)v * C;
          if (idx >= lo && idx <= hi) dst[idx] = S[(u64)(c0 + u) * R + r0 + v];
        }
      }
    }
  }
}
static void fwd(Eng &e, const u32 *src, u32 nsrc, u32 *S, u32 *X, int ctr,
                bool defer_final_row = false) {
  load_T(e, src, nsrc, S, ctr);
  corner_fwd(e, S, X, e.w, 1);
  if (!defer_final_row) for (int k0 = 0; k0 < e.R; k0++) {
    u32 *row = X + (u64)k0 * (u64)XST_;
    row_dif(row, e.C, e.tfC, e.wsC);
  }
}
static void inv(Eng &e, u32 *X, u32 *S) {
  corner_inv(e, X, S, pw(e.w, P - 2u));
}
static u32 Aarr[NRW_ * XST_], Barr[NRW_ * XST_ + 1352];   
static Eng E1, E2;
static int inited = 0;
static void eng_init(Eng &e, int lgR, int lgC, u32 w,
                     u32 *tfR, u32 *tiR, u32 *tfC, u32 *tiC, u16 *brR, u16 *brC,
                     u32 *wsR, u32 *wsC) {
  e.R = 1 << lgR; e.C = 1 << lgC; e.N = e.R * e.C; e.w = w;
  e.tfR = tfR; e.tiR = tiR; e.tfC = tfC; e.tiC = tiC; e.brR = brR; e.brC = brC;
  e.wsR = wsR; e.wsC = wsC;
  u32 wR = pw(w, (u32)e.C), wC = pw(w, (u32)e.R);
#if HAVE_AVX2
  build_tab2(tfR, e.wsR, e.R, wR); build_tab2(tiR, e.wsR + e.R, e.R, pw(wR, P - 2u));
  build_tab2(tfC, e.wsC, e.C, wC); build_tab2(tiC, e.wsC + e.C, e.C, pw(wC, P - 2u));
#else
  build_tab(tfR, e.R, wR); build_tab(tiR, e.R, pw(wR, P - 2u));
  build_tab(tfC, e.C, wC); build_tab(tiC, e.C, pw(wC, P - 2u));
#endif
  { u32 *fo = (u32 *)twint_for(tfC);
    size_t off = 0;
    for (int len = e.C; len >= 32; len >>= 2) {
      int h = len >> 1, h2 = len >> 2;
      const u32 *t1 = tfC + h - 1,   *ts1 = e.wsC + h - 1;
      const u32 *t2 = t1 + h2,       *ts2 = ts1 + h2;
      const u32 *t3 = tfC + h2 - 1,  *ts3 = e.wsC + h2 - 1;
      for (int j = 0; j < h2; j += 8) {
        u32 *d = fo + off + TIK(j, 0);
        memcpy(d +  0, t1  + j, 32); memcpy(d +  8, ts1 + j, 32);
        memcpy(d + 16, t2  + j, 32); memcpy(d + 24, ts2 + j, 32);
        memcpy(d + 32, t3  + j, 32); memcpy(d + 40, ts3 + j, 32);
      }
      off += TIK(h2, 0);
    }
    u32 *io = (u32 *)twint_for(tiC);
    size_t o2 = 0; int toff = 7;
    for (int len = 16; (len << 1) <= e.C; len <<= 2) {
      int h = len >> 1;
      const u32 *t1 = tiC + toff,      *ts1 = e.wsC + e.C + toff;
      const u32 *t2 = t1 + h,          *ts2 = ts1 + h;
      const u32 *t3 = t2 + h,          *ts3 = ts2 + h;
      toff += 3 * h;
      for (int j = 0; j < h; j += 8) {
        u32 *d = io + o2 + TIK(j, 0);
        memcpy(d +  0, t1  + j, 32); memcpy(d +  8, ts1 + j, 32);
        memcpy(d + 16, t2  + j, 32); memcpy(d + 24, ts2 + j, 32);
        memcpy(d + 32, t3  + j, 32); memcpy(d + 40, ts3 + j, 32);
      }
      o2 += TIK(h, 0);
    }
  }
  build_br(brR, e.R); build_br(brC, e.C);
  build_rev6();
}
#pragma GCC push_options
#pragma GCC target("avx2")
static inline __attribute__((always_inline)) __m256i scan8add(__m256i v) {
  __m256i t = _mm256_slli_si256(v, 4);
  v = _mm256_add_epi32(v, t);
  t = _mm256_slli_si256(v, 8);
  v = _mm256_add_epi32(v, t);
  t = _mm256_permute2x128_si256(v, v, 0x08);
  t = _mm256_shuffle_epi32(t, 0xFF);
  return _mm256_add_epi32(v, t);
}
#pragma GCC pop_options
#pragma GCC push_options
#pragma GCC target("avx2")
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
  const int c32 = (int)((((uintptr_t)(const void *)c) & 31u) == 0u);
  IM = ~0ull / (u64)P;
  if (!inited) {
    inited = 1;
    { u64 inv = 1; for (int i = 0; i < 5; i++) inv = inv * (2 - (u64)P * inv); NPI = (u32)(-(u64)inv); }
    u32 w1 = pw(G, (u32)((P - 1) / N1_));
    u32 w2 = pw(G, (u32)((P - 1) / N2_));
    eng_init(E1, LGR_, LGC_, w1, t1fR, t1iR, t1fC, t1iC, b1R, b1C, w1R, w1C);
    eng_init(E2, LGR2_, LGC2_, w2, t2fR, t2iR, t2fC, t2iC, b2R, b2C, w2R, w2C);
  }
  u32 *S = (u32 *)c;                 
  u32 *A = Aarr, *B = Barr;
  const int ctr = EXACT;
  const u32 sc = pw((u32)(N1_ % P), P - 2u);     
#if HAVE_AVX2
  g_inv_scale = mm(sc, R32);  
#else
  g_inv_scale = sc;
#endif
  const long Sa = (long)N1_ - (long)m;
  const long Sb = (long)N1_ - (long)n;
  const long M  = (long)n + (long)m - (long)N1_;
  if (M >= 0) {
    fwd(E1, (const u32 *)a + Sa, (u32)(M + 1), S, A, ctr, true);
    fwd(E1, (const u32 *)b + Sb, (u32)(M + 1), S, B, ctr, true);
    pointwise_dit(A, B, N1_, sc, E1.tfC, E1.wsC, E1.tiC, E1.wsC + E1.C);
    inv(E1, A, S);
    untr_window(E1, S, c + (N1_ - M), M, 2 * M);
  }
  fwd(E1, (const u32 *)a, (u32)n + 1, S, A, ctr, true);
  fwd(E1, (const u32 *)b, (u32)m + 1, S, B, ctr, true);
  pointwise_dit(A, B, N1_, sc, E1.tfC, E1.wsC, E1.tiC, E1.wsC + E1.C);
  inv(E1, A, S);
  untranspose(E1, S, A);
  long nmax = (long)n + (long)m;
  i64 WA = 0, WB = 0;
  if (n == m && M >= 0 && (long)M < (long)n && (long)n < (long)N1_ && (long)N1_ < nmax) {
    const long N = (long)n;
    const __m256i Pvm = _mm256_set1_epi32((int)P);
    const __m256i Phm = _mm256_set1_epi32((int)(P >> 1));
    const __m256i Zv = _mm256_setzero_si256();
    __m256i WAv = Zv, WBv = Zv;
    const long rb[5] = {0, M + 1, N + 1, (long)N1_, nmax + 1};
    for (int j1 = 0; j1 < (int)NRW_; j1++) {
      const long rbase = (long)j1 * (long)NCW_;
      const long rend  = rbase + (long)NCW_;
      u32 *const Arow = A + (u64)j1 * (u64)XST_;
      for (int rg = 0; rg < 3; rg++) {
        long s = rb[rg];       if (s  < rbase) s  = rbase;
        long k1 = rb[rg + 1];  if (k1 > rend)  k1 = rend;
        if (s >= k1) continue;
        WAv = _mm256_set1_epi32((int)WA);
        WBv = _mm256_set1_epi32((int)WB);
        long vend = k1 - (rg == 0 ? 16 : 8);
        if (vend < s) vend = s;
        long k = s;
        if (vend > k) {
          const int incPos = (rg < 2);
          const int ntA = c32 & (((k & 7) == 0) ? 1 : 0);   
          for (; k + 8 <= vend; k += 8) {
            __m256i dav, dbv;
            if (incPos) { dav = ldu256(a + k); dbv = ldu256(b + k); }
            else { dav = ldu256(a + k - N - 1); dbv = ldu256(b + k - N - 1);
                   dav = _mm256_sub_epi32(Zv, dav); dbv = _mm256_sub_epi32(Zv, dbv); }
            __m256i wav = _mm256_add_epi32(WAv, scan8add(dav));
            __m256i wbv = _mm256_add_epi32(WBv, scan8add(dbv));
            WAv = _mm256_permutevar8x32_epi32(wav, _mm256_set1_epi32(7));
            WBv = _mm256_permutevar8x32_epi32(wbv, _mm256_set1_epi32(7));
            __m256i tv;
            if (incPos) {
              tv = _mm256_setr_epi32((int)(k + 1), (int)(k + 2), (int)(k + 3), (int)(k + 4),
                                     (int)(k + 5), (int)(k + 6), (int)(k + 7), (int)(k + 8));
            } else {
              tv = _mm256_setr_epi32((int)(2 * N - k + 1), (int)(2 * N - k), (int)(2 * N - k - 1), (int)(2 * N - k - 2),
                                     (int)(2 * N - k - 3), (int)(2 * N - k - 4), (int)(2 * N - k - 5), (int)(2 * N - k - 6));
            }
            __m256i tab;
            if (rg == 0) {
              tab = _mm256_sub_epi32(ldu256(Arow + (k - rbase)), ldu256(c + N1_ + k));
              tab = _mm256_min_epu32(tab, _mm256_add_epi32(tab, Pvm));
            } else {
              tab = ldu256(Arow + (k - rbase));
            }
            __m256i gt = _mm256_cmpgt_epi32(tab, Phm);
            __m256i ts = _mm256_sub_epi32(tab, _mm256_and_si256(gt, Pvm));
            __m256i s4 = _mm256_slli_epi32(_mm256_add_epi32(wav, wbv), 2);
            __m256i v16 = _mm256_slli_epi32(tv, 4);
            if (ntA) stnt256(c + k, _mm256_sub_epi32(_mm256_add_epi32(ts, s4), v16));
            else     stu256(c + k, _mm256_sub_epi32(_mm256_add_epi32(ts, s4), v16));
            WA = _mm256_extract_epi32(wav, 7);
            WB = _mm256_extract_epi32(wbv, 7);
          }
        }
        for (; k < k1; k++) {
          if (rg < 2) { WA += a[k]; WB += b[k]; } else { WA -= a[k - N - 1]; WB -= b[k - N - 1]; }
          long t = (k <= N) ? (k + 1) : (2 * N - k + 1);
          u32 tab;
          if (rg == 0) { u32 x = Arow[k - rbase], y = c[N1_ + k]; tab = (x >= y) ? (x - y) : (x + (P - y)); }
          else tab = Arow[k - rbase];
          i32 ts = (i32)tab;
          if (ts > (i32)(P >> 1)) ts -= (i32)P;
          c[k] = (u32)((i64)ts + 4 * (WA + WB) - 16 * t);
        }
      }
    }
    {
      WAv = _mm256_set1_epi32((int)WA);
      WBv = _mm256_set1_epi32((int)WB);
      long k1 = nmax + 1;
      long vend = k1 - 16;
      long k = (long)N1_;
      if (vend < k) vend = k;
      const int ntB = c32;   
      for (; k + 8 <= vend; k += 8) {
        __m256i dav = ldu256(a + k - N - 1), dbv = ldu256(b + k - N - 1);
        dav = _mm256_sub_epi32(Zv, dav); dbv = _mm256_sub_epi32(Zv, dbv);
        __m256i wav = _mm256_add_epi32(WAv, scan8add(dav));
        __m256i wbv = _mm256_add_epi32(WBv, scan8add(dbv));
        WAv = _mm256_permutevar8x32_epi32(wav, _mm256_set1_epi32(7));
        WBv = _mm256_permutevar8x32_epi32(wbv, _mm256_set1_epi32(7));
        __m256i tv = _mm256_setr_epi32((int)(2 * N - k + 1), (int)(2 * N - k), (int)(2 * N - k - 1), (int)(2 * N - k - 2),
                                       (int)(2 * N - k - 3), (int)(2 * N - k - 4), (int)(2 * N - k - 5), (int)(2 * N - k - 6));
        __m256i tab = ldu256(c + k);
        __m256i gt = _mm256_cmpgt_epi32(tab, Phm);
        __m256i ts = _mm256_sub_epi32(tab, _mm256_and_si256(gt, Pvm));
        __m256i s4 = _mm256_slli_epi32(_mm256_add_epi32(wav, wbv), 2);
        __m256i v16 = _mm256_slli_epi32(tv, 4);
        if (ntB) stnt256(c + k, _mm256_sub_epi32(_mm256_add_epi32(ts, s4), v16));
        else     stu256(c + k, _mm256_sub_epi32(_mm256_add_epi32(ts, s4), v16));
        WA = _mm256_extract_epi32(wav, 7);
        WB = _mm256_extract_epi32(wbv, 7);
      }
      for (; k < k1; k++) {
        WA -= a[k - N - 1]; WB -= b[k - N - 1];
        long t = 2 * N - k + 1;
        u32 tab = c[k];
        i32 ts = (i32)tab;
        if (ts > (i32)(P >> 1)) ts -= (i32)P;
        c[k] = (u32)((i64)ts + 4 * (WA + WB) - 16 * t);
      }
    }
    if (c32) _mm_sfence();
    (void)E2;
    return;
  }
  u32 *Arow = A;
  long rbase = 0;
  for (long k = 0; k <= nmax; k++) {
    if (k >= rbase + (long)NCW_ && rbase + (long)NCW_ < (long)N1_) {
      rbase += (long)NCW_; Arow += (u64)XST_;
    }
    if (k <= n) WA += a[k];
    if (k <= m) WB += b[k];
    if (k > (long)m) WA -= a[k - 1 - m];
    if (k > (long)n) WB -= b[k - 1 - n];
    long lo = (k - (long)m > 0) ? (k - (long)m) : 0;
    long hi = (k < (long)n) ? k : (long)n;
    long t = hi - lo + 1;
    u32 tab;
    if (k <= M) {
      u32 x = Arow[k - rbase], y = c[N1_ + k];   
      tab = (x >= y) ? (x - y) : (x + (P - y));
    } else if (k < (long)N1_) {
      tab = Arow[k - rbase];
    } else {
      tab = c[k];                                
    }
#if EXACT
    i32 ts = (i32)tab;
    if (ts > (i32)(P >> 1)) ts -= (i32)P;
    i64 val = (i64)ts + 4 * (WA + WB) - 16 * t;
    c[k] = (u32)val;
#else
    c[k] = tab;
#endif
  }
  if (c32) _mm_sfence();
  (void)E2;
}
#pragma GCC pop_options

评测结果