/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121157: taken verbatim as this file's body --
the complete current best piece including its independently-balanced decimal parser (the carry
chosen from the raw limb, telescoping correction) and all of its own reference notices to
saffah_cc_v41_agg1 (#120960) and the inherited HintFFT/TSKY (MIT) sources, retained below.
- saffah_cc_v41_agg1 (this account), https://duck.ac/submission/120960: the forward-mid fusion /
staggered prefetch / bigint engine that #121157 carries.
Idea [E4W-ZQ]: the forward FFT top stage (fwd4's slice loop) loads the four quarters of the padded
array, but for this size regime the two upper quarters are identically zero: the operand data
(used <= 2*child doubles) occupies only the first two quarters, quarters 2 and 3 being untouched
zero padding. Skip those two loads and the two add/sub pairs that combine them whenever
used <= child<<1, writing the algebraically identical result:
t0=x0+x2 -> x0, t1=x0-x2 -> x0, t2=x1+x3 -> x1, t3=x1-x3 -> x1, so
y0=x0+x1, y1=x0+negI(x1), y2=x0-x1, y3=x0+posI(x1).
All four outputs, the twiddles, the store positions and every later stage are unchanged; the
small-N / huge-N path keeps the original four-load form.
Experiment purpose: drop 2*16.7 MB of pointless zero-reads from the forward top stage (memory
traffic is the binding resource on the judge) with every output byte identical.
*/
/* References:
- saffah_codex_6a_agg3, https://duck.ac/submission/121144: our independently
balanced decimal parser and its cited FFT engine; all notices retained.
Idea: Keep the preceding balancing masks in vector registers. Byte alignment
constructs all four incoming carries together, eliminating per-iteration
scalar extract/insert and the separate low-lane addition.
Experiment purpose: Reduce parser instructions with exactly the same limbs.
*/
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/120960: forward-mid fusion,
staggered prefetch and full bigint implementation; original notices below.
- TSKY/WithSky, HintFFT: inherited MIT FFT code; full license retained below.
Idea: Choose each base100000 balancing carry from the raw decimal limb,
independently of the previous carry. Let g[i]=(raw[i]>=50000), and store
raw[i]+g[i-1]-100000*g[i]. The carry terms telescope exactly, so the integer
is unchanged; balanced limbs lie in[-50000,50000]. This eliminates the
exceptional49999 carry-propagation branches and all inter-vector carry
dependencies in the decision flags. Scalar remainder remains exact.
Experiment purpose: Reduce input conversion instructions and dependency chains;
validate complete products, including repeated threshold-boundary limbs.
*/
// ===== REFERENCES =====
// 作者: saffah_cc_v41_agg1
// 原始提交: https://duck.ac/submission/120625 （65.945398 ms，我方现役最好件）
// **本文件正文 = 其正文逐字节 + 三处 [E4Z-*] 改动（见下"本次思路"），其余未改一字** ✓
// 该件的来源链（本队历史提交，原样保留在本文件下方各段头注中，未改一字）:
// #111416（[E4Y-TL2] 前向顶层↔递归首级 分片融合）· #111350（[E7Y2-PFS]）· #110398（[E7Y2-PREFAULT]）
// · 对手 #110391 · #110191（[E7X-R4] radix-5 基座 → 纯 2 幂次 radix-4 基座）· #109940 · #109800
// · 对手 #109582 · #109394 · pdoom #89173 · TSKY(WithSky) HintFFT（MIT，许可证全文见下方头注）
// 另引用同族账号 saffah_codex_6a_agg3 的件（其三个 token 被正文逐字节复刻，见下）:
// https://duck.ac/submission/120673 （65.499619 ms）· https://duck.ac/submission/120613
// 其自述思想原文: "Retain the smaller FFT leaves and compact loop bodies. Add explicit
// prefetches to the paired child pointwise-product loop"（LOG_SHORT 11→9 /
// #pragma optimize("O3,no-unroll-loops") / 点乘配对环四路预取）
// 本题地址: https://duck.ac/problem/1004e7
// 使用: 正文复刻 + 本席三处 [E4Z-*] 改动 + 复刻上述三 token；许可: 与本队历史提交一致
// ======================
// ===== 思路 =====
// 【进场口径】`python3 tools/exact.py 1004e7`：mine 65.945398（#120625）· T 65.499619（#120673）
// ⇒ **严支** 0.99T+1µs = 64.845623 ⇒ 需再快 **1.099775 ms（1.668%）** ✓
// （宽支 1.005T=65.828117 不可达：任何新提交都会把 T_sub<mine_sub 变真 ⇒ 锁严支）
// 【本发三处改动（单变量：只动调度/预取，不动任何浮点表达式、表、下标）】
// ① [E4Z-M23] **前向 `difRecMid` 的前两层（stride cs/4 与 cs/16）融合成一趟"垂直(列)通道"**：
// 列 j0 内先对 4 个 (k2,象限) 位置各做一次 dif4（层1，twiddle `MT(cs/2)[j>>3]`），
// 再用其 4 个输出做层2（twiddle `MT(cs/8)[j0>>3]`）并写回——**中间结果只驻寄存器**，
// 原先"层1 全数组写一遍 → 层2 全数组读回来"的一整趟 DRAM 往返消失 ✓
// 残下递归 `difRecMid(块大小 cs/16)` 与叶子逐字节不动；小 N/超大 N 走原路径 ✓
// ② [E4Z-PF1] inv4 收尾合并环 4 路**错位**预取（+256/+264/+272/+280；同余 4096 的 4 条流错开一行）
// ③ [E4Z-PF2] ①的列通道 4 路错位预取（+192/+200/+208）
// 另复刻同族账号三 token（LOG_SHORT=9 / no-unroll-loops / 点乘配对环四路预取，见 REFERENCES）✓
// 【判题机探针（tools/probe.py 立即体验；程序内自造 2×10^7 位输入 + rdtscp 相位表 + chk 校验和；
// 臂序轮转；每次 --reps 3）】全部臂 chk = d9765fc410434a23（同一输出 ⇒ 逐位等价）✓
// | 件 | fwd_a+fwd_b | dot | inv | tot | times_ns (min) |
// | base = #120625 正文 | 104.82M | 18.35M | 49.88M | 237.86M | 66.24 ms |
// | v1 = base + ① | 103.29M | 18.03M | 49.76M | 236.3M | 65.84 ms |
// | v5 = v1 + 同族三token| 102.87M | 16.04M | 49.75M | 233.51M | 65.13 ms |
// | 本发 v6 = v5 + ②③ | **100.33M** | 16.02M | **47.99M** | **229.16M** | **63.96 ms** |
// ⇒ 判题机探针两轮独立复测：base→v6 = **−2.28 ms（−3.7%）**，两轮读数一致（±0.02 ms）✓
// ★ 本机（Skylake-X，L3 38.5 MB ⇒ 数组全驻 L3）对本发 ① 的读数 ≈ 0 甚至为负，判题机为 −0.38 ms；
// 对逆侧融合（负项）本机为 −6%、判题机为 +13.5% ⇒ **本题"本机不可外推"再次兑现** ✓
// 【判负项（勿重走，均已过逐位闸门但判题机为负）】
// · [E4Z-ABC] 顶层+首级+次首级 三级融合同一趟（64 流 × 2 KB 块）：judge `fwd` **+9.4%** ✗
// · [E4Z-INV2] 逆侧前两级合并的垂直通道（16 流 × 8 KB 片 / 单趟全扫两种形态）：judge `inv` **+13.5%** ✗
// （16 条流 stride 均为 4096 倍数 ⇒ 同一 L1/L2 组 8/16 路冲突；AVX2 下 16 个 C64X4 同时活跃 ⇒ 溢出）
// 【闸门（判题同款 ref/gcc9/g9.sh -O2 -std=c++17 -static -U_FORTIFY_SOURCE）】
// · 生产形 work/e7a_in.txt（2×10^7 位）⇒ md5 **0ec796190f0afec508c63b9327368b5e = 基线预言机** ✓
// · 8 形状（100×101 / 1750×1750 / 2000×2000 / 200k / 525k / 1M / 2.1M×1.9M / 4.2M 位）⇒
// 与 #120625 正文输出**逐字节相同** ✓（其中 cs=2^15/2^16/2^17/2^18 四档均走到本发新通道 ✓）
// · 200k 位与 1M 位 ⇒ 对 Python 大整数**逐位相同** ✓
// 【预期】65.945398 − 2.28 ≈ **63.67 ms** ⇒ 过严支线 64.845623 ✓（余量 ≈ 1.1 ms）
// ======================
// ===== REFERENCES =====
// 作者: saffah_cc_v41_agg1
// 原始提交: https://duck.ac/submission/111416 <https://duck.ac/submission/111416>
// （69626.386 µs，我方现役最好件）：**本文件正文 = 其正文逐字节 + 一处 [E4Y-TL2] 改动** ✓
// 该件的来源链（本队历史提交，原样保留在本文件下方各段头注中，未改一字）:
// #111350（[E7Y2-PFS] 递推层预取错位）· #110398（[E7Y2-PREFAULT] 页预触）· 对手 #110391
// · #110191（[E7X-R4] radix-5 基座 → 纯 2 幂次 radix-4 基座）· #109940 · #109800
// · 对手 #109582 · #109394 · pdoom #89173 · TSKY(WithSky) HintFFT（MIT，许可证全文见下方头注）
// 本题地址: https://duck.ac/problem/1004e7
// 使用: 正文复刻 + 一处 [E4Y-TL2]（前向 FFT 顶层与递归首级的分片融合）改动；许可: 与本队历史提交一致
// ======================
// ===== 思路 =====
// 【进场口径】`python3 tools/exact.py 1004e7`：mine 69.626386（#111416）· T 69.255130（#120488）
// ⇒ 宽支差 23.98 µs 但**不可达**（发射即锁严支）⇒ **严支** 0.99T+1µs = 68.563579 ⇒ 需再快 **1.062807 ms（1.526%）** ✓
// 【本发单变量】[E4Y-TL2] **把 `fwd4` 的 radix-4 顶层与其后 `difRecLong` 递归的第一级，融合成一趟分片(blocked)循环**
// · 病灶（判题机相位表，本席探针实测）：前向两趟 `fwd_a+fwd_b = 118.2M 拍 / 全集 251.2M 拍 = 47%`，是整条流水线最大项；
// 顶层的 4 条流在 `stride = child`（8 MB）上，紧跟的递归首级 4 条流在 `stride = child/4`（2 MB）上 ——
// **顶层把整个数组写一遍、递归首级立刻把它读回来**，两次都是 DRAM 往返（每操作数 67 MB，两操作数 134 MB）✓
// · 处置：以 **SL 双精度/流 为片**，片内先做顶层（4 子带 × 4 流，twiddle 用 `roots[(s*cs+j0+j)/8]`），
// 再做 4 个子带的递归首级（4 流，twiddle 用 `FFTAVX::sqrt_table_1[log2(child/2)]` 的原下标 `(j0+j)/8`）；
// 片的中间结果（16 流 × SL × 8 B）只驻留 L2/L3，**不再回 DRAM** ⇒ 省掉一整趟读 + 一整趟写 ✓
// 顶层原先对子带 2/3 用的 NT store 也一并消失（其余子带的下游读者就在同一片内）✓
// · **算术逐位不动**：同一表达式序列（顶层 = 原 6 条加减 + 3 次复乘；首级 = 引擎 `dif4` + 3 次复乘）、
// 同一 twiddle 表对象、同一 twiddle 下标语义（均为"块内绝对下标"，与原来逐位相同）✓
// · 小 N（`child/2 ≤ MID_LEN` 或 `≥ 2^22`，即 `difRecLong` 不进入 4 子带层：小输入直接 `difRecMid`，
// 超大输入走 `fusedLong`）走**原路径逐字节不动** ✓
// 【判题机实测（本席探针：程序内自造 2×10^7 位输入 + rdtscp 相位表；臂序轮转；两轮共 44 次调用）】
// | 件 | fwd_a+fwd_b (min) | dot | inv | 全程 tot (min) | vs base |
// | base（#111416 正文） | 118.12M | 18.03M | 49.78M | 250.51M | — |
// | 本发 [E4Y-TL2] **SL=1024** | **104.75M** | 17.87M | 49.76M | **237.29M** | **−13.22M 拍（−5.28%）≈ −3.68 ms** |
// | （片长对照臂）SL=256 / 512 / 2048 / 4096 | 104.97M / 104.72M / 106.07M / 108.43M | — | — | 237.87M / 237.60M / 238.65M / 240.99M | −5.05% / −5.15% / −4.99% / −3.94% |
// ★ **内置正控**：未改动的 `inv`/`dot` 相位两臂读数一致（49.78M → 49.76M、18.03M → 17.87M，差 < 1%）⇒ 效应属于被改的 fwd ✓
// ★ 同源空对照：同一 base 源重复提交 4 次，tot 极差 0.26%（250.51M vs 251.18M）⇒ 上述 5% 级差异远超噪声 ✓
// 【闸门（判题同款 `ref/gcc9/g9.sh -O2 -std=c++17 -static -U_FORTIFY_SOURCE`）】
// · 生产形 `work/e7a_in.txt`（2×10^7 位）⇒ md5 **`0ec796190f0afec508c63b9327368b5e` = 基线预言机** ✓
// · 另 2 个全新 2×10^7 位随机输入 ⇒ 与 #111416 正文输出**逐字节相同**（20,000,001 B）✓
// · N=2^21（3M 位）与 N=2^22（6M 位）⇒ 与正文**逐字节相同**，且 3M 位**对 Python 大整数逐位相同** ✓
// · 4 个小形状（2000×2000 / 1750×1750 / 100×101 / 199932×199932 位）⇒ 与正文逐字节相同 + 对 Python ✓
// 【预期】判题机 ≈ 69.63 ms − 3.77 ms ≈ **65.9 ms** ⇒ 过严支线 68.563579 ✓（余量 ≈ 2.7 ms）
// ======================
// ===== REFERENCES =====
// 作者: saffah_cc_v41_agg1
// 原始提交: https://duck.ac/submission/111350 <https://duck.ac/submission/111350>
// （70783.555 µs，我方现役最好件）：**本文件正文 = 其正文逐字节 + 一处 [E7Y2-PFS] 改动** ✓
// 该件的来源: https://duck.ac/submission/110398（其正文 = 对手 https://duck.ac/submission/110391
// 的逐字节复刻）· 本题 https://duck.ac/problem/1004e7
// 上游继承链（原样保留在本文件下方各段头注中，未改一字）: #110191 · #109940 · #109800 · #109582
// · #109394 · pdoom #89173 · TSKY(WithSky) HintFFT（MIT，许可证全文见下方头注）
// 使用: 正文复刻 + 一处 [E7Y2-PFS]（radix-4 递推层预取错位）改动；许可: 与本队历史提交一致
// ======================
// ===== 思路 =====
// 【进场口径】mine 70.783555（#111350）· T 70.928420（#110391）· 严支 0.99T+1µs = 70.220136
// ⇒ 需再快 **0.563419 ms** ✓
// 【判题机相位表（本发先做的那一步）】自造数探针（`tools/probe.py`，程序内生成 2×1e7 位）实测（拍，3.59 GHz）：
// parse 15.9M · tabs 0.04M · **fwd_a top 13.6M + rec 46.7M** · **fwd_b top 14.1M + rec 46.9M**
// · dot 18.2M · **inv rec 37.1M + comb 13.6M** · fmt 24.7M ⇒ 合计 230.9M 拍
// ⇒ ★ **递推（difRecLong/iditRecLong）占整条流水线的 52~57%**，而本机相位表把首触页错误当成了主角
// （本机 fwd top 44M ↔ 判题机 13.6M）⇒ **本机相位表在本题不可用，必须以判题机相位表为准** ✓
// 【本发单变量】[E7Y2-PFS] **radix-4 递推层的 4 路预取错位**：
// `difLayer`/`iditLayer` 的 4 条流在 `stride = float_len/4` 上，而 **stride 恒为 4096 字节的倍数**
// ⇒ 4 路（载 4 + 存 4 = 8 流）**全部落在同一个 L1/L2 set**（该判据由引擎自己的
// `[A12X-PFS]` 注释在 `fusedLong` 上写明，只是从未应用到递推层）⇒ 把 4 条预取的距离从
// 一律 `+320` 改成 **`+192/+256/+320/+384`**（各差 64 B = 1 条 cache line ⇒ 落进 4 个不同 set）✓
// 【判题机实测（探针装置：同一程序内 前向 4 子块 / 逆向 4 子块，交错 2 轮，各程序 2 次）】
// | 件 | fwd1 | fwd2 | inv1 | inv2 |
// | 原版(320×4) | 45.71M | 48.35M | 36.94M | 36.81M |
// | 仅前向错位 | **44.33M(−2.99%)** | **46.78M(−3.28%)** | 36.89M(**不变**) | 36.84M(**不变**) |
// | 前向+逆向错位 | 44.43M | 46.72M | **35.79M(−3.2%)** | **35.74M(−3.2%)** |
// ⇒ ★ **内置正控**：只改前向时**逆向臂读数纹丝不动**（36.94→36.89），改前向臂同步 −3.0%
// ⇒ 排除机器漂移，效应属于被改的那段代码 ✓
// 【预期】递推段 130.8M 拍 ×≈3% ≈ **−4.0M 拍 = 整条流水线 −1.7% ⇒ 判题机约 −1.2 ms**
// （判据 0.563419 ms ⇒ 2× 余量）✓
// 【闸门（判题同款 gcc9）】编译 rc=0 ✓ · 生产形 `work/e7a_in.txt` ⇒ md5
// `0ec796190f0afec508c63b9327368b5e` 与 #111350 逐字节相同 ✓ · 4 形状（4000/3501/201/399864 B）逐字节相同 ✓
// （`prefetch` 不参与架构状态，等价性由构造保证，仍按全闸门复核）✓
// ======================
// ===== REFERENCES =====
// 作者: saffah_cc_v41_agg1
// 原始提交: https://duck.ac/submission/110398 <https://duck.ac/submission/110398>
// （70806.501 µs，我方现役最好件）：**本文件正文 = 其正文逐字节 + 一处 [E7Y2-PREFAULT] 改动** ✓
// 同一代码的对手提交: https://duck.ac/submission/110391
// （70828.420 µs，saffah_codex_6s_agg2 的当代最好件；#110398 是其正文的逐字节复刻）
// 本题地址: https://duck.ac/problem/1004e7
// 上游继承链（原样保留在本文件下方各段头注中，未改一字）:
// https://duck.ac/submission/110191（本账号 r4 换引擎：radix-5 -> 纯 2 幂次 radix-4）
// · https://duck.ac/submission/109940 · https://duck.ac/submission/109800
// · https://duck.ac/submission/109582（对手件）· https://duck.ac/submission/109394
// · pdoom #89173 · TSKY(WithSky) HintFFT（MIT 许可，许可证全文见下方头注）
// 使用: 正文复刻 + 一处 [E7Y2-PREFAULT]（页预触）改动；许可: 与本队历史提交一致（站点公开提交）
// ======================
// ===== 思路 =====
// 【进场口径】`python3 tools/exact.py 1004e7`：mine = 70.806501（#110398）· T = 70.928420（#110391）
// 两件号仅差 7 ⇒ 我方提交更新 ⇒ **严支** 0.99T+1µs = **70.220136** ⇒ 需再快 **0.586365 ms（0.828%）** ✓
// 【本发单变量】[E7Y2-PREFAULT] **把"首次触碰"的页错误从 FFT 热环里挪到紧凑顺序循环**：
// ① 诊断（本席自建 rdtscp 相位表 + 故意错的重复臂）：对 `fwd4` 的顶阶同一段代码、同一批地址跑三遍 ⇒
// 第 1 遍 top = 43.6M~49.5M 拍，第 2/3 遍 = **7.3M / 6.9M 拍**（差 6 倍）⇒ 差额 = **页错误**
// （两个 32 MB `_mm_malloc` 数组里，解析器只写了各自低半；FFT 顶阶是**第一段写高半**的代码，
// 在 4 流 gather/scatter 环里吃 ~16 MB/操作数的"写时分配零页"错误）✓
// ② 处置：FFT 前对**从未被触碰的区间**（`A[na..N)` / `B[nb..N)`）每 4 KiB 页写一个 `0.0`
// ⇒ 页错误在紧凑顺序循环内完成（实测该循环 = 40.4M 拍；FFT 侧省下 ~74M 拍）✓
// ③ 等价性：写入值 `0.0` 与该处**既有语义**（新鲜匿名映射的零页，原码注释已声明依赖它）逐位相同
// ⇒ **输出逐字节不变**（含"不写"的那 511/512 个 double——内核分配时整页清零）✓
// 【闸门（本席自建，判题同款：`ref/gcc9/g++-9 -O2 -std=c++17 -static -U_FORTIFY_SOURCE`）】
// ① 编译 rc=0（零 error/零 warning）✓
// ② 生产形 `work/e7a_in.txt`（2×1e7 位）⇒ 输出 **20000001 字节 / md5
// `0ec796190f0afec508c63b9327368b5e` = 与 #110398 逐字节相同** ✓
// ③ 形状差分（见下）：小/中形状与 #110398 输出逐字节相同 ✓
// 【实测（本机，3 轮交错 A/B，同窗口）】总拍数 base 270.6/268.8/271.8M ↔ 本发 **247.0/235.9/231.4M**
// ⇒ **−12.0%**（本机只作粗筛，量级 ≫ 5% 门槛；定价以判题机为准）✓
// ======================
// References:
// [1] saffah_cc_v41_agg1, https://duck.ac/submission/110191 : directly
// reused the accepted FFT engine, parser, formatter and public stdio driver.
// [2] saffah_codex_6s_agg2, https://duck.ac/submission/109582 : parts of
// the parser and big integer driver reached [1] through our earlier code.
// [3] TSKY (WithSky), HintFFT, https://github.com/With-Sky/HintFFT : inherited
// FFT kernels are MIT licensed; full copyright/permission notice follows.
// The duck.ac submissions provide no further license notice. Authors,
// submission addresses, and reused components are identified above.
// Approach:
// Load four consecutive complex base twiddles with two AVX loads, transpose
// real and imaginary lanes, and multiply by the segment constant in SIMD.
// This replaces four scalar complex multiplies and per-lane branch checks.
// Purpose:
// Experimental official evaluation of vectorized twiddle generation.
/*
MIT License -- HintFFT by TSKY (WithSky)
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #109940 <https://duck.ac/submission/109940>
// （77868.066 µs，我方现役最好件 = work/e7x_G1_sub.cpp）：**本文件正文 = 其正文逐字节 + 一处
// [E7X-R4] 引擎级改动**（见下「思路」）✓
// [2] #109940 的正文来源链（本账号 #109800 · 对手 saffah_codex_6s_agg2 #109582 · 本账号 #109394 ·
// pdoom #89173 · TSKY HintFFT(MIT) 等）**原样保留在本文件下半部分的头注中**，未改一字 ✓
// 许可证合规：duck.ac 题面明示「你提交的代码将会被公开，所有人都可见」⇒ 站上提交按站点规则公开；
// 本发逐条标注原作者账号与原提交地址，未宣称任何部分为原创 ✓
// ======================
// ===== 思路 =====
// 【本发单变量】[E7X-R4] **换引擎类：radix-5 实 FFT 基座 -> 纯 2 幂次 radix-4 基座**
// （N = 2 幂次，超配 4.8%；现役 radix-5 = 5·2^k，超配 31%）⇒ 变换长度约少 20%，
// FFT 相位（前向×2 + 点乘 + 逆）随之省下；判题机独立原型上界实测 20.98%（≈11.7 ms，15× 缺口）✓
// · 前向：自写 radix-4 顶层（4 子带）＋引擎 `difRecLong` 递归；**打包融进顶层**（`load(p+j,true_type)`）✓
// · 点乘：自写向量化 unpack/repack（公式与已过闸的标量版逐字一致），twiddle 用**定长 1024 段表**
// （ts[t]=e^{-2πi·bitrev10(t)/2048}，与 Lc 无关）+ 每段一个常数（1 次复乘/lane）；
// 子带 0 直接用引擎自带 `real_dot_binrev4`（逐元素已对参考正控）✓
// · 逆：自写 radix-4 顶层 ＋ 引擎 `iditRecLong` 递归；**缩放+取整+交错融进顶层**（`*scale -> transToI64`）✓
// 【等价性/正控】生产形 md5 = 0ec796190f0afec508c63b9327368b5e（= 基线 oracle ✓）；
// 6 个形状（85x115 / 2000 / 4000 / 50000 / 200k / 20 MB 位）逐个对 **Python 大整数**参考逐字节一致 ✓；
// ★ 200k 位输入上现役件反而错（399860 位 vs 参考 399863 位），本发正确 ✓
// 【与前序改动的关系】[E7X-G1]（末级逆 radix-5 第 5 子带"写了从不被读"的两路特化）所消除的冗余
// 在 radix-4 顶层里天然不存在（顶层只写 4 个子带）⇒ 本发已包含该收益 ✓
// 【实测（本机筛选；律一：本机不作量级判据，生产件一律判题机实测定价）】
// user CPU（20 MB 输入，10 次求和）：本发 0.95 s vs 现役 #109940 1.00 s ⇒ **−5.0%** ✓
// ======================
// ===== 思路（本席）=====
// 【进场口径】mine 63.708706 (#120960) · T 62.829802 (#121157, 同族 saffah_codex_6a_agg3)
// ⇒ 宽支 thr 63.144951 数值上只差 0.563755 ms 但**不可达**（我方任何新最快提交 ⇒ 翻严支）
// ⇒ 真靶 = 严支 0.99·T+1µs = 62.202504 ms（需 −1.506202 ms）
// 【§2.18.349 抄件分解】#121157 正文与我方 #120960 正文逐行 diff 只有 2 个 hunk：
// ① 自撰头注（无作用）；② `parse_balanced_base100000_dual` 的进位逻辑重写（g[i]=(raw[i]>=50000)，
// 存 raw[i]+g[i-1]-100000*g[i]，进位项完全望远镜 ⇒ 整数值不变；去掉 49999 边界传播分支与向量内
// 进位依赖）。判题机相位表实测 parse 21.02M → 17.83M 拍（−3.19M ≈ −0.879 ms，与板面
// 63.708706 → 62.829802 精确吻合）⇒ 抄件上限 = 对手件本身 = 62.829802 ⇒ 必须另补一刀 ✓
// 【本发单刀 [E4W-ZQ]】前向 FFT 顶层四路载入中**上两象限恒为零**（used=2 000 001 ≤ 2·child=2 097 152，
// Q2/Q3 全落在零填充区）⇒ 跳过这两路载入与 t0/t1/t2/t3 组合，直接写代数等价结果；四个输出、
// twiddle、写出位置及其后全部各级逐字节不动。判题机探针实测：fa+fb 100.33M → 99.42M 拍 ·
// tot 225.89M → 224.99M（−0.90M 拍 = −0.25 ms，−0.40%）✓
// 【闸门】判题同款 gcc9 -O2 -static：生产形 md5 `0ec796190f0afec508c63b9327368b5e` = 基线预言机 ✓
// · 10 形状（100×101 / 1750 / 2000 / 200k / 525k / 1M / 2.1M×1.9M / 4.2M / 3M / 4e6+1 位）
// 与基座输出**逐字节相同** ✓ · 25 例对 Python 大整数**逐位相同**（含 8 个 pos=1 构造）✓
// ======================
#include <complex>
#include <iostream>
#include <type_traits>
#include <cstdint>
#include <climits>
#include <cstring>
#include <immintrin.h>
#pragma GCC target("avx2,fma,tune=skylake")
#pragma GCC optimize("O3,no-unroll-loops","rename-registers")
namespace hint
{
template <typename T, size_t ALIGN = 64>
class AlignMem{
public:
using Ptr = T *;
using ConstPtr = const T *;
~AlignMem(){
if (ptr){
_mm_free(ptr);}};
AlignMem() : ptr(nullptr), len(0) {}
AlignMem(size_t n) : ptr(reinterpret_cast<Ptr>(_mm_malloc(n * sizeof(T), ALIGN))), len(n) {}
AlignMem(const AlignMem &) = delete;
AlignMem &operator=(const AlignMem &) = delete;
T &operator[](size_t i){
return ptr[i];}
const T &operator[](size_t i) const{
return ptr[i];}
Ptr begin(){
return ptr;}
Ptr end(){
return ptr + len;}
ConstPtr begin() const{
return ptr;}
ConstPtr end() const{
return ptr + len;}
size_t size() const{
return len;}
private:
T *ptr;
size_t len;};
template <typename YMM>
inline void transpose64_2X4(YMM &row0, YMM &row1){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
row0 = YMM(_mm256_permute2f128_pd(t0, t1, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t0, t1, 0x31));}
template <typename YMM>
inline void transpose64_4X2(YMM &row0, YMM &row1){
auto t0 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x20);
auto t1 = _mm256_permute2f128_pd(__m256d(row0), __m256d(row1), 0x31);
row0 = YMM(_mm256_unpacklo_pd(t0, t1));
row1 = YMM(_mm256_unpackhi_pd(t0, t1));}
template <typename YMM>
inline void transpose64_4X4(YMM &row0, YMM &row1, YMM &row2, YMM &row3){
auto t0 = _mm256_unpacklo_pd(__m256d(row0), __m256d(row1));
auto t1 = _mm256_unpackhi_pd(__m256d(row0), __m256d(row1));
auto t2 = _mm256_unpacklo_pd(__m256d(row2), __m256d(row3));
auto t3 = _mm256_unpackhi_pd(__m256d(row2), __m256d(row3));
row0 = YMM(_mm256_permute2f128_pd(t0, t2, 0x20));
row1 = YMM(_mm256_permute2f128_pd(t1, t3, 0x20));
row2 = YMM(_mm256_permute2f128_pd(t0, t2, 0x31));
row3 = YMM(_mm256_permute2f128_pd(t1, t3, 0x31));}
class Float64X4{
public:
using F64 = double;
using F64X4 = Float64X4;
Float64X4() : data(_mm256_setzero_pd()) {}
Float64X4(__m256d in_data) : data(in_data) {}
Float64X4(F64 in_data) : data(_mm256_set1_pd(in_data)) {}
Float64X4(const F64 *in_data) : data(_mm256_load_pd(in_data)) {}
F64X4 operator+(const F64X4 &other) const{
return _mm256_add_pd(data, other.data);}
F64X4 operator-(const F64X4 &other) const{
return _mm256_sub_pd(data, other.data);}
F64X4 operator*(const F64X4 &other) const{
return _mm256_mul_pd(data, other.data);}
F64X4 operator/(const F64X4 &other) const{
return _mm256_div_pd(data, other.data);}
static F64X4 fmadd(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmadd_pd(a.data, b.data, c.data);}
static F64X4 fmsub(const F64X4 &a, const F64X4 &b, const F64X4 &c){
return _mm256_fmsub_pd(a.data, b.data, c.data);}
template <int N>
F64X4 permute4x64() const{
return _mm256_permute4x64_pd(data, N);}
F64X4 reverse() const{
return permute4x64<0b00011011>();}
void load(const F64 *p){
data = _mm256_load_pd(p);}
void load1(const F64 *p){
data = _mm256_broadcast_sd(p);}
void store(F64 *p) const{
_mm256_store_pd(p, data);}
void store_nt(F64 *p) const{ /* p8c */
_mm256_stream_pd(p, data);}
operator __m256d() const{
return data;}
__m256i toI64X4() const{
constexpr uint64_t mask = (uint64_t(1) << 52) - 1;
constexpr uint64_t offset = (uint64_t(1) << 10) - 1;
const __m256i f64bits = _mm256_castpd_si256(data);
__m256i tail = _mm256_and_si256(f64bits, _mm256_set1_epi64x(mask));
tail = _mm256_or_si256(tail, _mm256_set1_epi64x(mask + 1));
__m256i exp = _mm256_srli_epi64(f64bits, 52);
exp = _mm256_sub_epi64(_mm256_set1_epi64x(offset + 52), exp);
return _mm256_srlv_epi64(tail, exp);}
private:
__m256d data;};
struct Complex64X4{
using C64X4 = Complex64X4;
using F64X4 = Float64X4;
using F64 = double;
Complex64X4() {}
Complex64X4(F64X4 real, F64X4 imag) : real(real), imag(imag) {}
Complex64X4(const F64 *p) : real(p), imag(p + 4) {}
Complex64X4(const F64 *p_real, const F64 *p_imag) : real(p_real), imag(p_imag) {}
C64X4 operator+(const C64X4 &other) const{
return C64X4(real + other.real, imag + other.imag);}
C64X4 operator-(const C64X4 &other) const{
return C64X4(real - other.real, imag - other.imag);}
C64X4 operator*(const F64X4 &other) const{
return C64X4(real * other, imag * other);}
C64X4 mul(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmsub(real, other.real, ii);
const F64X4 i = F64X4::fmadd(imag, other.real, ri);
return C64X4(r, i);}
C64X4 mulConj(const C64X4 &other) const{
const F64X4 ii = imag * other.imag;
const F64X4 ri = real * other.imag;
const F64X4 r = F64X4::fmadd(real, other.real, ii);
const F64X4 i = F64X4::fmsub(imag, other.real, ri);
return C64X4(r, i);}
C64X4 reverse() const{
return C64X4(real.reverse(), imag.reverse());}
void set1(F64 real_in, F64 imag_in){
real = F64X4(real_in);
imag = F64X4(imag_in);}
template <typename T>
void load(const T *p, std::false_type){
this->load(p);}
template <typename T>
void load(const T *p, std::true_type){
this->load(p);
*this = this->toRRIIPermu();}
template <typename T>
void load(const T *p){
real.load(reinterpret_cast<const F64 *>(p));
imag.load(reinterpret_cast<const F64 *>(p) + 4);}
void load1(const F64 *real_p, const F64 *imag_p){
real.load1(real_p);
imag.load1(imag_p);}
template <typename T>
void store(T *p, std::false_type) const{
this->store(p);}
template <typename T>
void store(T *p, std::true_type) const{
this->toRIRIPermu().store(p);}
template <typename T>
void store(T *p) const{
real.store(reinterpret_cast<F64 *>(p));
imag.store(reinterpret_cast<F64 *>(p) + 4);}
template <typename T>
void store_nt(T *p) const{ /* p8c */
real.store_nt(reinterpret_cast<F64 *>(p));
imag.store_nt(reinterpret_cast<F64 *>(p) + 4);}
C64X4 toRIRIPermu() const{
C64X4 res = *this;
transpose64_2X4(res.real, res.imag);
return res;}
C64X4 toRRIIPermu() const{
C64X4 res = *this;
transpose64_4X2(res.real, res.imag);
return res;}
C64X4 transToI64(std::false_type) const{
return *this;}
C64X4 transToI64(std::true_type) const{
const __m256d magic=_mm256_set1_pd(6755399441055744.0);
const __m256i bits=_mm256_castpd_si256(magic);
auto real_i64=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(real,magic)),bits);
auto imag_i64=_mm256_sub_epi64(_mm256_castpd_si256(_mm256_add_pd(imag,magic)),bits);
return C64X4(__m256d(real_i64),__m256d(imag_i64));}
F64X4 real, imag;};
using Float32 = float;
using Float64 = double;
constexpr Float64 HINT_PI = 3.141592653589793238462643;
constexpr Float64 HINT_2PI = HINT_PI * 2;
constexpr Float64 COS_PI_8 = 0.707106781186547524400844;
template <typename T>
constexpr T int_floor2(T n){
constexpr int bits = sizeof(n) * 8;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return (n >> 1) + 1;}
template <typename T>
constexpr T int_ceil2(T n){
constexpr int bits = sizeof(n) * 8;
n--;
for (int i = 1; i < bits; i *= 2){
n |= (n >> i);}
return n + 1;}
template <typename IntTy>
constexpr bool is_2pow(IntTy n){
return n != 0 && (n & (n - 1)) == 0;}
template <typename T>
constexpr int hint_log2(T n){
constexpr int bits = sizeof(n) * 8;
int l = -1, r = bits;
while ((l + 1) != r){
int mid = (l + r) / 2;
if ((T(1) << mid) > n){
r = mid;}
else{
l = mid;}}
return l;}
constexpr uint32_t crc32(const char *str){
uint32_t crc = 0xFFFFFFFF;
while (*str != '\0'){
crc ^= *str;
for (int i = 0; i < 8; ++i){
crc = (crc >> 1) ^ (0 - (crc & 1)) & 0xEDB88320;}
str++;}
return ~crc;}
namespace transform{
template <typename T>
inline void transform2(T &sum, T &diff){
T temp0 = sum, temp1 = diff;
sum = temp0 + temp1;
diff = temp0 - temp1;}
namespace fft{
using F64 = Float64;
using C64 = std::complex<F64>;
using F64X4 = Float64X4;
using C64X4 = Complex64X4;
template <typename Float, size_t OMEGA_LEN>
class TableFix{
alignas(64) Float table[OMEGA_LEN * 2];
public:
TableFix(size_t theta_divider, size_t factor, size_t stride){
const Float theta = -HINT_2PI * factor / theta_divider;
for (size_t begin = 0, index = 0; begin < OMEGA_LEN * 2; begin += stride * 2)
{
for (size_t j = 0; j < stride; j++, index++)
{
table[begin + j] = std::cos(theta * index);
table[begin + j + stride] = std::sin(theta * index);}}}
constexpr const Float &operator[](size_t index) const{
return table[index];}};
void initOmegaX4(F64 *arr, size_t fft_len, int table_len, int factor){
table_len /= 4;
const F64 theta = -HINT_2PI * factor / fft_len;
auto arrx4 = reinterpret_cast<C64X4 *>(arr);
arr[0] = 1, arr[4] = 0;
arr[1] = std::cos(theta), arr[5] = std::sin(theta);
arr[2] = std::cos(theta * 2), arr[6] = std::sin(theta * 2);
arr[3] = std::cos(theta * 3), arr[7] = std::sin(theta * 3);
for (size_t begin = 1; begin < table_len; begin *= 2){
size_t nth = begin * 4;
C64X4 unit;
unit.set1(std::cos(theta * nth), std::sin(theta * nth));
for (size_t i = 0; i < begin; i++)
{
arrx4[begin + i] = arrx4[i].mul(unit);}}}
template <typename Float, int LOG_BEGIN, int LOG_END, int DIV>
class TableFixMulti{
static_assert(LOG_END >= LOG_BEGIN);
static_assert(is_2pow(DIV));
static constexpr size_t TABLE_CPX_LEN = (size_t(1) << (LOG_END + 1)) / DIV;
alignas(64) Float table[TABLE_CPX_LEN * 2];
public:
TableFixMulti(size_t factor, size_t stride = 4){
initBottomUp(factor, stride);}
void initBottomUp(size_t factor, size_t stride){
static_assert(std::is_same<Float, Float64>::value);
size_t len = size_t(1) << LOG_BEGIN, cpx_len = len / DIV;
auto it = getBeginLog(LOG_BEGIN);
initOmegaX4(it, len, cpx_len, factor);
for (int log_len = LOG_BEGIN + 1; log_len <= LOG_END; log_len++)
{
len = size_t(1) << log_len, cpx_len = len / DIV;
Float theta = -HINT_2PI * factor / len;
auto it = getBeginLog(log_len), it_last = getBeginLog(log_len - 1);
C64X4 unit(std::cos(theta), std::sin(theta));
for (auto end = it + cpx_len * 2; it < end; it += 16, it_last += 8)
{
C64X4 omega0, omega1;
omega0.load(it_last);
omega1 = omega0.mul(unit);
transpose64_2X4(omega0.real, omega1.real);
transpose64_2X4(omega0.imag, omega1.imag);
omega0.store(it), omega1.store(it + 8);}}}
constexpr const Float *getBeginLog(int log_rank) const{
return getBegin(size_t(1) << log_rank);}
constexpr Float *getBeginLog(int log_rank){
return getBegin(size_t(1) << log_rank);}
constexpr const Float *getBegin(size_t rank) const{
return &table[rank * 2 / DIV];}
constexpr Float *getBegin(size_t rank){
return &table[rank * 2 / DIV];}};
template <int CACHE_LOG_LEN>
class FFTSqrtTableC64X4{
public:
using F64 = double;
using C64 = std::complex<double>;
using C64X4 = hint::Complex64X4;
static constexpr size_t CACHE_LEN = size_t(1) << CACHE_LOG_LEN;
static constexpr size_t MASK = CACHE_LEN - 1;
static constexpr size_t C4_COUNT = sizeof(C64X4) / sizeof(C64);
~FFTSqrtTableC64X4(){
if (high)
{
delete[] high;}}
FFTSqrtTableC64X4() {}
FFTSqrtTableC64X4(size_t fft_len, int len_div, int factor){
init(fft_len, len_div, factor);}
void init(size_t fft_len, int len_div, int factor){
size_t table_len = fft_len / len_div;
size_t low_len = CACHE_LEN * C4_COUNT, high_len = table_len / low_len;
if (high != nullptr)
{
delete[] high;}
high = new C64[high_len];
auto p = reinterpret_cast<F64 *>(&low[0]);
initOmegaX4(p, fft_len, low_len, factor);
const F64 theta = -HINT_2PI * factor / fft_len;
high[0] = C64(1, 0);
for (size_t begin = 1; begin < high_len; begin *= 2)
{
C64 unit = std::polar<F64>(1.0, theta * begin * low_len);
for (size_t i = 0; i < begin; i++)
{
high[i + begin] = high[i] * unit;}}}
C64X4 operator[](size_t i) const{
C64X4 hi;
auto p = reinterpret_cast<const F64 *>(&high[i >> CACHE_LOG_LEN]);
hi.load1(p, p + 1);
return low[i & MASK].mul(hi);}
private:
alignas(64) C64X4 low[CACHE_LEN];
C64 *high = nullptr;};
template <int DIV, int LOG_BEGIN, int LOG_MAX, int CACHE_LOG_LEN>
class FFTTableSqrt{
using TableLong = FFTSqrtTableC64X4<CACHE_LOG_LEN>;
static constexpr size_t SHORT_LEN = size_t(1) << LOG_BEGIN;
static constexpr size_t TABLE_LEN = LOG_MAX - LOG_BEGIN + 1;
public:
FFTTableSqrt(int factor){
for (int i = 0; i < TABLE_LEN; i++)
{
size_t fft_len = SHORT_LEN << i;
table[i].init(fft_len, DIV, factor);}}
const TableLong &operator[](int log_len) const{
log_len -= LOG_BEGIN;
return table[log_len];}
TableLong &operator[](int log_len){
log_len -= LOG_BEGIN;
return table[log_len];}
private:
TableLong table[TABLE_LEN];};
struct FFT{
template <typename Float>
static void trans2MulI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 + i1;
r0 = r0 - i1;
i1 = i0 - temp;
i0 = i0 + temp;}
template <typename Float>
static void trans2MulNegI(Float &r0, Float &i0, Float &r1, Float &i1){
auto temp = r1;
r1 = r0 - i1;
r0 = r0 + i1;
i1 = i0 + temp;
i0 = i0 - temp;}
template <typename Float>
static void dif4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
difSplit(r0, i0, r1, i1, r2, i2, r3, i3);
transform2(r0, r1);
transform2(i0, i1);}
template <typename Float>
static void idit4(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r1);
transform2(i0, i1);
iditSplit(r0, i0, r1, i1, r2, i2, r3, i3);}
template <typename Float>
static void difSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r0, r2);
transform2(i0, i2);
transform2(r1, r3);
transform2(i1, i3);
trans2MulNegI(r2, i2, r3, i3);}
template <typename Float>
static void iditSplit(Float &r0, Float &i0, Float &r1, Float &i1, Float &r2, Float &i2, Float &r3, Float &i3){
transform2(r2, r3);
transform2(i2, i3);
transform2(r0, r2);
transform2(i0, i2);
trans2MulI(r1, i1, r3, i3);}};
struct FFTAVX : public FFT{
static constexpr int LOG_SHORT = 9, LOG_MID = 17, LOG_MAX = 22, LOG_CACHE = 8;
static constexpr size_t SHORT_LEN = size_t(1) << LOG_SHORT, MID_LEN = size_t(1) << LOG_MID, MAX_LEN = size_t(1) << LOG_MAX;
using TableFix4 = const TableFix<Float64, 4>;
using TableFix8 = const TableFix<Float64, 8>;
using TableMulti1 = const TableFixMulti<Float64, LOG_SHORT + 1, LOG_MID, 4>;
using TableMulti2 = const TableFixMulti<Float64, 6, LOG_SHORT + 1, 4>;
using TableMulti3 = const TableFixMulti<Float64, 6, LOG_SHORT, 4>;
using TableSqrt = const FFTTableSqrt<4, LOG_MID + 1, LOG_MAX, LOG_CACHE>;
static TableFix4 table_8, table_16_1, table_16_3;
static TableFix8 table_32_1, table_32_3;
static TableMulti2 multi_table_2;
static TableMulti3 multi_table_3;
static TableMulti1 multi_table_1;
static TableSqrt sqrt_table_1;
static constexpr const Float64 *it8 = &table_8[0], *it16_1 = &table_16_1[0], *it16_3 = &table_16_3[0], *it32_1 = &table_32_1[0], *it32_3 = &table_32_3[0];
static void dif4x4(F64X4 &r0, F64X4 &i0, F64X4 &r1, F64X4 &i1, F64X4 &r2, F64X4 &i2, F64X4 &r3, F64X4 &i3){
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);
dif4(r0, i0, r1, i1, r2, i2, r3, i3);
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);}
static void idit4x4(F64X4 &r0, F64X4 &i0, F64X4 &r1, F64X4 &i1, F64X4 &r2, F64X4 &i2, F64X4 &r3, F64X4 &i3){
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);
idit4(r0, i0, r1, i1, r2, i2, r3, i3);
transpose64_4X4(r0, r1, r2, r3);
transpose64_4X4(i0, i1, i2, i3);}
static void dif8x2(C64X4 &c0, C64X4 &c1, C64X4 &c2, C64X4 &c3){
C64X4 omega(it8);
transform2(c0, c1);
transform2(c2, c3);
c1 = c1.mul(omega), c3 = c3.mul(omega);
dif4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);}
static void idit8x2(C64X4 &c0, C64X4 &c1, C64X4 &c2, C64X4 &c3){
C64X4 omega(it8);
idit4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mulConj(omega), c3 = c3.mulConj(omega);
transform2(c0, c1);
transform2(c2, c3);}
static void dif16(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[1], c2 = p[2], c3 = p[3];
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(it8)), c2 = c2.mul(C64X4(it16_1)), c3 = c3.mul(C64X4(it16_3));
dif4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[0] = c0, p[1] = c1, p[2] = c2, p[3] = c3;}
static void idit16(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[1], c2 = p[2], c3 = p[3], omega;
idit4x4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mulConj(C64X4(it8)), c2 = c2.mulConj(C64X4(it16_1)), c3 = c3.mulConj(C64X4(it16_3));
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[0] = c0, p[1] = c1, p[2] = c2, p[3] = c3;}
static void dif32(Float64 in_out[]){
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[0], c1 = p[2], c2 = p[4], c3 = p[6];
difSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c2 = c2.mul(C64X4(it32_1)), c3 = c3.mul(C64X4(it32_3));
p[0] = c0, p[2] = c1, p[4] = c2, p[6] = c3;
c0 = p[1], c1 = p[3], c2 = p[5], c3 = p[7];
difSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c2 = c2.mul(C64X4(it32_1 + 8)), c3 = c3.mul(C64X4(it32_3 + 8));
p[1] = c0, p[3] = c1, c0 = p[4], c1 = p[6];
dif8x2(c0, c2, c1, c3);
p[4] = c0, p[5] = c2, p[6] = c1, p[7] = c3;
dif16(in_out);}
static void idit32(Float64 in_out[]){
idit16(in_out);
auto p = reinterpret_cast<C64X4 *>(in_out);
C64X4 c0 = p[4], c1 = p[5], c2 = p[6], c3 = p[7];
idit8x2(c0, c1, c2, c3);
p[5] = c1, p[7] = c3, c1 = p[0], c3 = p[2];
c0 = c0.mulConj(C64X4(it32_1)), c2 = c2.mulConj(C64X4(it32_3));
iditSplit(c1.real, c1.imag, c3.real, c3.imag, c0.real, c0.imag, c2.real, c2.imag);
p[0] = c1, p[2] = c3, p[4] = c0, p[6] = c2;
c0 = p[1], c1 = p[3], c2 = p[5], c3 = p[7];
c2 = c2.mulConj(C64X4(it32_1 + 8)), c3 = c3.mulConj(C64X4(it32_3 + 8));
iditSplit(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
p[1] = c0, p[3] = c1, p[5] = c2, p[7] = c3;}
template <typename F32, typename F16>
static void fftTiny(Float64 in_out[], size_t float_len, F32 &&func32, F16 &&func16){
if (hint_log2(float_len / 2) % 2 == 0)
{
for (auto end = in_out + float_len; in_out < end; in_out += 32)
{
func16(in_out);}}
else
{
for (auto end = in_out + float_len; in_out < end; in_out += 64)
{
func32(in_out);}}}
static void difIter(Float64 in_out[], size_t float_len){
size_t fft_len = float_len / 2;
C64X4 c0, c1, c2, c3;
size_t stride = fft_len / 2;
auto it0 = in_out, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (size_t rank = fft_len; rank >= 64; rank /= 4)
{
stride = rank / 2;
for (auto begin = in_out, end = in_out + float_len; begin < end; begin += rank * 2)
{
auto table1 = multi_table_2.getBegin(rank * 2), table2 = multi_table_2.getBegin(rank), table3 = multi_table_3.getBegin(rank);
it0 = begin, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (; it0 < begin + stride; it0 += 8, it1 += 8, it2 += 8, it3 += 8, table1 += 8, table2 += 8, table3 += 8)
{
c0 = it0, c1 = it1, c2 = it2, c3 = it3;
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c1 = c1.mul(C64X4(table2)), c2 = c2.mul(C64X4(table1)), c3 = c3.mul(C64X4(table3));
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3);
}}}
fftTiny(in_out, float_len, dif32, dif16);}
static void iditIter(Float64 in_out[], size_t float_len){
size_t fft_len = float_len / 2;
size_t rank = hint_log2(fft_len) % 2 == 0 ? 64 : 128;
fftTiny(in_out, float_len, idit32, idit16);
for (; rank <= fft_len; rank *= 4)
{
const size_t stride = rank / 2;
for (auto begin = in_out, end = in_out + float_len; begin < end; begin += rank * 2)
{
auto table1 = multi_table_2.getBegin(rank * 2), table2 = multi_table_2.getBegin(rank), table3 = multi_table_3.getBegin(rank);
auto it0 = begin, it1 = it0 + stride, it2 = it1 + stride, it3 = it2 + stride;
for (; it0 < begin + stride; it0 += 8, it1 += 8, it2 += 8, it3 += 8, table1 += 8, table2 += 8, table3 += 8)
{
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3;
c1 = c1.mulConj(C64X4(table2)), c2 = c2.mulConj(C64X4(table1)), c3 = c3.mulConj(C64X4(table3));
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3);
}}}}
#define difLayer(dif_func, in_out, stride, table) \
do \
{ \
auto it0 = in_out, it1 = in_out + stride, it2 = it1 + stride, it3 = it2 + stride; \
size_t indx = 0; \
for (auto end = it1; FROM_RIRI_PERM && it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
C64X4 c0, c1, c2, c3, omega1, omega2; \
c0.load(it0, FromRIRI{}), c1.load(it1, FromRIRI{}), c2 = c0,c3 = c1; \
transform2(c0,c1); \
trans2MulNegI(c2.real,c2.imag,c3.real,c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
for (auto end = it1; (!FROM_RIRI_PERM) && it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
_mm_prefetch((const char*)(it0+192),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+256),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+320),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+384),_MM_HINT_T0); \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
omega1 = table[indx], c2 = c2.mul(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2); \
c3 = c3.mul(omega2.mul(omega1)); \
c0.store(it0), c1.store(it1), c2.store(it2), c3.store(it3); \
} \
dif_func(in_out, stride); \
dif_func(in_out + stride, stride); \
dif_func(in_out + stride * 2, stride); \
dif_func(in_out + stride * 3, stride); \
} while (0)
#define iditLayer(idit_func, in_out, stride, table) \
do \
{ \
idit_func(in_out, stride); \
idit_func(in_out + stride, stride); \
idit_func(in_out + stride * 2, stride); \
idit_func(in_out + stride * 3, stride); \
auto it0 = in_out, it1 = in_out + stride, it2 = it1 + stride, it3 = it2 + stride; \
size_t indx = 0; \
_Pragma("GCC unroll 4") \
for (auto end = it1; it0 < end; it0 += 8, it1 += 8, it2 += 8, it3 += 8, indx++) \
{ \
_mm_prefetch((const char*)(it0+192),_MM_HINT_T0); \
_mm_prefetch((const char*)(it1+256),_MM_HINT_T0); \
_mm_prefetch((const char*)(it2+320),_MM_HINT_T0); \
_mm_prefetch((const char*)(it3+384),_MM_HINT_T0); \
C64X4 c0 = it0, c1 = it1, c2 = it2, c3 = it3, omega1, omega2; \
omega1 = table[indx], c2 = c2.mulConj(omega1); \
omega2 = omega1.mul(omega1), c1 = c1.mulConj(omega2); \
c3 = c3.mulConj(omega2.mul(omega1)); \
idit4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag); \
c0 = c0.transToI64(ToI64{}), c1 = c1.transToI64(ToI64{}), c2 = c2.transToI64(ToI64{}), c3 = c3.transToI64(ToI64{}); \
c0.store(it0, ToRIRI{}), c1.store(it1, ToRIRI{}), c2.store(it2, ToRIRI{}), c3.store(it3, ToRIRI{}); \
} \
} while (0)
template <bool FROM_RIRI_PERM = false>
static void difRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
difIter(in_out, float_len);
return;}
using FromRIRI = std::integral_constant<bool, FROM_RIRI_PERM>;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
const size_t stride = float_len / 4;
difLayer(difRecMid, in_out, stride, table1);}
template <bool TO_RIRI_PERM = false, bool TO_INT64 = false>
static void iditRecMid(Float64 in_out[], size_t float_len){
const size_t fft_len = float_len / 2;
if (fft_len <= SHORT_LEN)
{
iditIter(in_out, float_len);
return;}
using ToRIRI = std::integral_constant<bool, TO_RIRI_PERM>;
using ToI64 = std::integral_constant<bool, TO_INT64>;
const size_t stride = float_len / 4;
auto table1 = reinterpret_cast<const C64X4 *>(multi_table_1.getBegin(fft_len));
iditLayer(iditRecMid, in_out, stride, table1);}
// Fuse two long radix-4 stages while a 16-vector tile is in L1.
// Preserves the original butterfly formulas and bit-reversed ordering.
static void fusedLong(Float64* a,size_t len,bool inverse){
size_t stride=len/16;
const auto &t0=sqrt_table_1[hint_log2(len/2)];
const auto &t1=sqrt_table_1[hint_log2(len/8)];
if(inverse)for(int k=0;k<16;k++)iditRecLong<false,false>(a+k*stride,stride);
for(size_t j=0;j<stride;j+=8){
// a12x_ [A12X-PF]: the 16 panel rows are `stride` doubles apart (stride = len/16);
// each j-step reads 16 x 32 B from 16 different pages -> the HW prefetcher cannot
// track 16 streams, so the loads are latency-bound. Prefetch the whole tile 128
// doubles (= 16 j-steps) ahead; row stride and j are loop-invariant so the addresses
// are one lea each. Values are untouched -> output must be bit-identical.
/* a12x_ [A12X-PFS] set-conflict fix. `stride = len/16` doubles is ALWAYS a
multiple of 4096 bytes (len = 5*2^m => stride = 5*2^(m-4) doubles = 5*2^(m-1) bytes),
so all 16 rows of the tile are congruent mod 4096 => every row's line lands in the
SAME L1 (64-set) and L2 set. Prefetching 16 lines into one 8-way set means they
evict each other before use. Staggering the per-row distance by 8 doubles (one
cache line = set index +1) puts the 16 prefetches in 16 DIFFERENT sets, keeping the
same lead (64..184 doubles ahead). */
if(j>=64 && j+416<stride){
for(int k=0;k<16;k++) _mm_prefetch((const char*)(a+(size_t)k*stride+j-64+32*k),_MM_HINT_T0);
}
C64X4 v[16];
auto butterfly=[&](C64X4 &c0,C64X4 &c1,C64X4 &c2,C64X4 &c3,const C64X4 &w){
C64X4 w2=w.mul(w),w3=w2.mul(w);
if(inverse){
c1=c1.mulConj(w2);c2=c2.mulConj(w);c3=c3.mulConj(w3);
idit4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
}else{
dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
c1=c1.mul(w2);c2=c2.mul(w);c3=c3.mul(w3);
}
};
if(inverse){
#pragma GCC unroll 4
for(int q=0;q<4;q++){
v[4*q]=C64X4(a+(4*q)*stride+j);
v[4*q+1]=C64X4(a+(4*q+1)*stride+j);
v[4*q+2]=C64X4(a+(4*q+2)*stride+j);
v[4*q+3]=C64X4(a+(4*q+3)*stride+j);
butterfly(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3],t1[j/8]);
}
#pragma GCC unroll 4
for(int k=0;k<4;k++){
butterfly(v[k],v[k+4],v[k+8],v[k+12],t0[(k*stride+j)/8]);
v[k].store(a+k*stride+j);
v[k+4].store(a+(k+4)*stride+j);
v[k+8].store(a+(k+8)*stride+j);
v[k+12].store(a+(k+12)*stride+j);
}
}else{
#pragma GCC unroll 4
for(int k=0;k<4;k++){
v[k]=C64X4(a+k*stride+j);
v[k+4]=C64X4(a+(k+4)*stride+j);
v[k+8]=C64X4(a+(k+8)*stride+j);
v[k+12]=C64X4(a+(k+12)*stride+j);
butterfly(v[k],v[k+4],v[k+8],v[k+12],t0[(k*stride+j)/8]);
}
#pragma GCC unroll 4
for(int q=0;q<4;q++){
butterfly(v[4*q],v[4*q+1],v[4*q+2],v[4*q+3],t1[j/8]);
v[4*q].store(a+(4*q)*stride+j);
v[4*q+1].store(a+(4*q+1)*stride+j);
v[4*q+2].store(a+(4*q+2)*stride+j);
v[4*q+3].store(a+(4*q+3)*stride+j);
}
}
}
if(!inverse)for(int k=0;k<16;k++)difRecLong<false>(a+k*stride,stride);
}
template <bool FROM_RIRI_PERM = false>
static void difRecLong(Float64 in_out[], size_t float_len){
if constexpr(!FROM_RIRI_PERM) if(float_len/2>=(1ULL<<22)){fusedLong(in_out,float_len,false);return;}
const size_t fft_len = float_len / 2;
if (fft_len <= MID_LEN)
{
difRecMid<FROM_RIRI_PERM>(in_out, float_len);
return;}
using FromRIRI = std::integral_constant<bool, FROM_RIRI_PERM>;
const auto &table1 = sqrt_table_1[hint_log2(fft_len)];
const size_t stride = float_len / 4;
difLayer(difRecLong, in_out, stride, table1);}
template <bool TO_RIRI_PERM = false, bool TO_INT64 = false>
static void iditRecLong(Float64 in_out[], size_t float_len){
if constexpr(!TO_RIRI_PERM && !TO_INT64) if(float_len/2>=(1ULL<<22)){fusedLong(in_out,float_len,true);return;}
const size_t fft_len = float_len / 2;
if (fft_len <= MID_LEN)
{
iditRecMid<TO_RIRI_PERM, TO_INT64>(in_out, float_len);
return;}
using ToRIRI = std::integral_constant<bool, TO_RIRI_PERM>;
using ToI64 = std::integral_constant<bool, TO_INT64>;
const size_t stride = float_len / 4;
const auto &table1 = sqrt_table_1[hint_log2(fft_len)];
iditLayer(iditRecLong, in_out, stride, table1);}};
#undef difLayer
#undef iditLayer
constexpr int FFTAVX::LOG_SHORT, FFTAVX::LOG_MID, FFTAVX::LOG_MAX, FFTAVX::LOG_CACHE;
constexpr size_t FFTAVX::SHORT_LEN, FFTAVX::MID_LEN, FFTAVX::MAX_LEN;
FFTAVX::TableFix4 FFTAVX::table_8(8, 1, 4), FFTAVX::table_16_1(16, 1, 4), FFTAVX::table_16_3(16, 3, 4);
FFTAVX::TableFix8 FFTAVX::table_32_1(32, 1, 4), FFTAVX::table_32_3(32, 3, 4);
FFTAVX::TableMulti2 FFTAVX::multi_table_2(2);
FFTAVX::TableMulti3 FFTAVX::multi_table_3(3);
FFTAVX::TableMulti1 FFTAVX::multi_table_1(1);
FFTAVX::TableSqrt FFTAVX::sqrt_table_1(1);
constexpr uint32_t bitrev32(uint32_t n){
constexpr uint32_t mask55 = 0x55555555;
constexpr uint32_t mask33 = 0x33333333;
constexpr uint32_t mask0f = 0x0f0f0f0f;
constexpr uint32_t maskff = 0x00ff00ff;
n = ((n & mask55) << 1) | ((n >> 1) & mask55);
n = ((n & mask33) << 2) | ((n >> 2) & mask33);
n = ((n & mask0f) << 4) | ((n >> 4) & mask0f);
n = ((n & maskff) << 8) | ((n >> 8) & maskff);
return (n << 16) | (n >> 16);}
constexpr uint32_t bitrev(uint32_t n, int len){
return bitrev32(n) >> (32 - len);}
template <int MAX_LOG_LEN, int DIV>
class BinRevTableC64X4HP{
public:
static constexpr int LOG_BLOCK = 2, BLOCK = 1 << LOG_BLOCK;
static constexpr size_t MAX_LEN = size_t(1) << MAX_LOG_LEN;
struct Unit{
C64 units[MAX_LOG_LEN]{};
F64 block[BLOCK * 2]{};
Unit()
{
constexpr F64 factor = F64(1) / DIV;
for (int i = 0; i < MAX_LOG_LEN; i++)
{
units[i] = getOmega(size_t(1) << (i + 1), 1, factor);}
block[0] = 1, block[BLOCK] = 0;
for (int i = 1; i < BLOCK; i++)
{
C64 omega = getOmega(BLOCK, bitrev(i, LOG_BLOCK), factor);
block[i] = omega.real(), block[i + BLOCK] = omega.imag();}}};
BinRevTableC64X4HP() : index(0), pop(0){
std::memcpy(table, unit_table.block, sizeof(unit_table.block));}
void reset(size_t i = 0){
if (i == 0)
{
pop = 0, index = i;
return;}
pop = 1, index = i / BLOCK;
int zero = __builtin_ctzll(index);
auto fp = reinterpret_cast<const F64 *>(&unit_table.units[zero + 2]);
table[1].load1(fp, fp + 1);
table[1] = table[1].mul(table[0]);}
C64X4 iterate(){
C64X4 res = table[pop], unit4;
index++;
int zero = __builtin_ctzll(index);
auto fp = reinterpret_cast<const F64 *>(&unit_table.units[zero + 2]);
unit4.load1(fp, fp + 1);
pop -= zero;
table[pop + 1] = table[pop].mul(unit4);
pop++;
return res;}
static C64 getOmega(size_t n, size_t index, F64 factor = 1){
F64 theta = -HINT_2PI * index / n;
return std::polar<F64>(1, theta * factor);}
private:
alignas(64) static const Unit unit_table;
alignas(64) C64X4 table[MAX_LOG_LEN];
size_t index;
int pop;
int log_max_iter, log_fft_len;};
template <int MAX_LOG_LEN, int DIV>
alignas(64) const typename BinRevTableC64X4HP<MAX_LOG_LEN, DIV>::Unit BinRevTableC64X4HP<MAX_LOG_LEN, DIV>::unit_table;
template <size_t RI_DIFF = 1, typename FloatTy>
inline void dot_rfft(FloatTy *inout0, FloatTy *inout1, const FloatTy *in0, const FloatTy *in1,
const std::complex<FloatTy> &omega, const FloatTy inv = 1){
using Complex = std::complex<FloatTy>;
auto addConj = [](Complex c0, Complex c1){ return Complex(c0.real() + c1.real(), c0.imag() - c1.imag()); };
Complex x0(inout0[0], inout0[RI_DIFF]), x1(inout1[0], inout1[RI_DIFF]),
y0(in0[0], in0[RI_DIFF]), y1(in1[0], in1[RI_DIFF]);
auto t0 = x0 * y0, t1 = x1 * y1, xy0 = addConj(x0, x1), xy1 = addConj(y0, y1);
auto t2 = xy0 * xy1;
y1 = addConj(t0, t1);
x1 = (y1 + y1 - t2) * omega * omega;
const auto inv2 = inv + inv;
x0 = (t2 - x1) * inv, x1 = Complex(t0.real() - t1.real(), t0.imag() + t1.imag()) * inv2;
Complex out0 = x0 + x1, out1(x0.real() - x1.real(), x1.imag() - x0.imag());
inout0[0] = out0.real(), inout0[RI_DIFF] = out0.imag();
inout1[0] = out1.real(), inout1[RI_DIFF] = out1.imag();}
inline void dot_rfftX4(F64 *inout0, F64 *inout1, const F64 *in0, const F64 *in1, const C64X4 &omega, const F64X4 &inv){
auto addConj = [](const C64X4 &x0, const C64X4 &x1){
return C64X4(x0.real + x1.real, x0.imag - x1.imag);};
C64X4 x0 = inout0, x1 = inout1, y0 = in0, y1 = in1;
x1 = x1.reverse();
y1 = y1.reverse();
C64X4 t0 = x0.mul(y0), t1 = x1.mul(y1);
C64X4 xy0 = addConj(x0, x1), xy1 = addConj(y0, y1);
C64X4 t2 = xy0.mul(xy1);
y1 = addConj(t0, t1);
x1 = (y1 + y1 - t2).mul(omega);
const F64X4 inv2 = inv + inv;
x0 = (t2 - x1) * inv, x1 = C64X4(t0.real - t1.real, t0.imag + t1.imag) * inv2;
C64X4 out0 = x0 + x1, out1(x0.real - x1.real, x1.imag - x0.imag);
out0.store(inout0), out1.reverse().store(inout1);}
inline void real_dot_binrev4(Float64 in_out[], Float64 in[], size_t float_len,Float64 scale=1){
Float64 inv = 2.0 * scale / float_len;{
auto r0 = in_out[0], i0 = in_out[4], r1 = in[0], i1 = in[4];
transform2(r0, i0);
transform2(r1, i1);
r0 *= r1, i0 *= i1;
transform2(r0, i0);
in_out[0] = r0 * 0.5 * inv, in_out[4] = i0 * 0.5 * inv;}
auto temp = C64(in_out[1], in_out[5]) * C64(in[1], in[5]) * inv;
in_out[1] = temp.real(), in_out[5] = temp.imag();
inv /= 4;
dot_rfft<4>(&in_out[2], &in_out[3], &in[2], &in[3], C64(COS_PI_8, -COS_PI_8), inv);
constexpr Float64 COS_16_1 = 0.92387953251128675612818318939;
constexpr Float64 SIN_16_1 = 0.38268343236508977172845998403;
dot_rfft<4>(&in_out[8], &in_out[11], &in[8], &in[11], C64(COS_16_1, -SIN_16_1), inv);
dot_rfft<4>(&in_out[9], &in_out[10], &in[9], &in[10], C64(-SIN_16_1, -COS_16_1), inv);
const Float64X4 inv4 = F64X4(0.5 * scale / float_len);
BinRevTableC64X4HP<28, 1> table;
for (size_t begin = 16; begin < float_len; begin *= 2){
table.reset(begin / 2);
auto it0 = in_out + begin, it1 = it0 + begin - 8, it2 = in + begin, it3 = it2 + begin - 8;
for (; it0 < it1; it0 += 8, it1 -= 8, it2 += 8, it3 -= 8)
{
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(it2+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it3-256),_MM_HINT_T0);
dot_rfftX4(it0, it1, it2, it3, table.iterate(), inv4);}}}
template <bool TO_INT = false>
inline void real_conv_avx(F64 *in_out1, F64 *in2, size_t float_len){
FFTAVX::difRecLong<true>(in_out1, float_len);
FFTAVX::difRecLong<true>(in2, float_len);
real_dot_binrev4(in_out1, in2, float_len);
FFTAVX::iditRecLong<true, TO_INT>(in_out1, float_len);}}}}
#include <cstdio>
#include <string>
#include <memory>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/auxv.h>
#include <cstddef>
// Public standard-stream ABI 0.04 (version 40), exact 64-bit field layout:
// https://github.com/JudgeDuck/JudgeDuck-OS/blob/d4df797bad6dc9b66c00312e97c05ca33adc3abd/inc/abi.hpp
// The public stdin and stdout fields implement zero-copy standard streams.
namespace duck_public_stdio {
constexpr unsigned long AT_DUCK = 0x6b637564;
struct DuckInfo_t {
uint64_t abi_version;
const char *stdin_ptr;
uint64_t stdin_size;
char *stdout_ptr;
uint64_t stdout_limit;
uint64_t stdout_size;
} __attribute__((packed));
static_assert(sizeof(void*)==8 && sizeof(DuckInfo_t)==48,"64-bit Duck ABI required");
static_assert(offsetof(DuckInfo_t,stdin_ptr)==8 && offsetof(DuckInfo_t,stdin_size)==16,
"Public stdin fields must match ABI 0.04");
static_assert(offsetof(DuckInfo_t,stdout_ptr)==24 && offsetof(DuckInfo_t,stdout_size)==40,
"Public stdout fields must match ABI 0.04");
}
#include <cstddef>
#include <cstdint>
#include <cstring>
namespace bigint_io_opt {
inline void parse_balanced_base100000_dual(const char* s0,std::size_t len0,double* out0,
const char* s1,std::size_t len1,double* out1) {
size_t k0=0,k1=0; int carry0=0,carry1=0;
__m128i prev0=_mm_setzero_si128(),prev1=_mm_setzero_si128();
const __m128i zeros=_mm_set1_epi8('0'),w1=_mm_set1_epi16(0x010a),w2=_mm_set1_epi32(0x00010064);
const __m128i sh=_mm_setr_epi8(11,12,13,14,6,7,8,9,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i tail=_mm_setr_epi8(15,-128,-128,-128,10,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i threshold=_mm_set1_epi32(49999),minusbase=_mm_set1_epi32(-100000);
while(len0>=26 && len1>=26) {
__m128i a0=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s0+len0-16)),zeros);
__m128i b0=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s0+len0-26)),zeros);
__m128i a1=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s1+len1-16)),zeros);
__m128i b1=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s1+len1-26)),zeros);
__m128i x0=_mm_unpacklo_epi64(_mm_shuffle_epi8(a0,sh),_mm_shuffle_epi8(b0,sh));
__m128i y0=_mm_unpacklo_epi64(_mm_shuffle_epi8(a0,tail),_mm_shuffle_epi8(b0,tail));
__m128i x1=_mm_unpacklo_epi64(_mm_shuffle_epi8(a1,sh),_mm_shuffle_epi8(b1,sh));
__m128i y1=_mm_unpacklo_epi64(_mm_shuffle_epi8(a1,tail),_mm_shuffle_epi8(b1,tail));
x0=_mm_madd_epi16(_mm_maddubs_epi16(x0,w1),w2);
x1=_mm_madd_epi16(_mm_maddubs_epi16(x1,w1),w2);
x0=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x0,3),_mm_slli_epi32(x0,1)),y0);
x1=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x1,3),_mm_slli_epi32(x1,1)),y1);
__m128i g0=_mm_cmpgt_epi32(x0,threshold);
__m128i g1=_mm_cmpgt_epi32(x1,threshold);
x0=_mm_add_epi32(_mm_sub_epi32(x0,_mm_alignr_epi8(g0,prev0,12)),_mm_and_si128(g0,minusbase));
x1=_mm_add_epi32(_mm_sub_epi32(x1,_mm_alignr_epi8(g1,prev1,12)),_mm_and_si128(g1,minusbase));
prev0=g0;prev1=g1;
_mm256_stream_pd(out0+k0,_mm256_cvtepi32_pd(x0));k0+=4;len0-=20; // [b12q NT-PARSE]
_mm256_stream_pd(out1+k1,_mm256_cvtepi32_pd(x1));k1+=4;len1-=20; // [b12q NT-PARSE]
}
_mm_sfence();
carry0=_mm_extract_epi32(prev0,3)&1;carry1=_mm_extract_epi32(prev1,3)&1;
while(len0>=5) { int v=0;for(int j=5;j>0;--j)v=v*10+s0[len0-j]-'0';v+=carry0; carry0=v>=50000;out0[k0++]=v-carry0*100000;len0-=5; }
{ int v=0;for(std::size_t j=0;j<len0;++j)v=v*10+s0[j]-'0';out0[k0]=v+carry0; }
while(len1>=5) { int v=0;for(int j=5;j>0;--j)v=v*10+s1[len1-j]-'0';v+=carry1; carry1=v>=50000;out1[k1++]=v-carry1*100000;len1-=5; }
{ int v=0;for(std::size_t j=0;j<len1;++j)v=v*10+s1[j]-'0';out1[k1]=v+carry1; }
}
inline void parse_balanced_base100000(const char* s,std::size_t len,double* out) {
size_t k=0;int carry=0;
const __m128i zeros=_mm_set1_epi8('0'),w1=_mm_set1_epi16(0x010a),w2=_mm_set1_epi32(0x00010064);
const __m128i sh=_mm_setr_epi8(11,12,13,14,6,7,8,9,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i tail=_mm_setr_epi8(15,-128,-128,-128,10,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128,-128);
const __m128i threshold=_mm_set1_epi32(49999),minusbase=_mm_set1_epi32(-100000);
while(len>=26) {
__m128i a=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s+len-16)),zeros);
__m128i b=_mm_sub_epi8(_mm_loadu_si128((const __m128i*)(s+len-26)),zeros);
__m128i x=_mm_unpacklo_epi64(_mm_shuffle_epi8(a,sh),_mm_shuffle_epi8(b,sh));
__m128i y=_mm_unpacklo_epi64(_mm_shuffle_epi8(a,tail),_mm_shuffle_epi8(b,tail));
x=_mm_madd_epi16(_mm_maddubs_epi16(x,w1),w2);
x=_mm_add_epi32(_mm_add_epi32(_mm_slli_epi32(x,3),_mm_slli_epi32(x,1)),y);
x=_mm_add_epi32(x,_mm_cvtsi32_si128(carry));
__m128i g=_mm_cmpgt_epi32(x,threshold),p=_mm_cmpeq_epi32(x,threshold);
if(!_mm_testz_si128(p,p)){
g=_mm_or_si128(g,_mm_and_si128(p,_mm_slli_si128(g,4)));
p=_mm_and_si128(p,_mm_slli_si128(p,4));
g=_mm_or_si128(g,_mm_and_si128(p,_mm_slli_si128(g,8)));
}
x=_mm_add_epi32(_mm_sub_epi32(x,_mm_slli_si128(g,4)),_mm_and_si128(g,minusbase));
carry=_mm_extract_epi32(g,3)&1;
_mm256_store_pd(out+k,_mm256_cvtepi32_pd(x));k+=4;len-=20;
}
while(len>=5) {
int v=0;for(int j=5;j>0;--j)v=v*10+s[len-j]-'0';v+=carry;
carry=v>=50000;out[k++]=v-carry*100000;len-=5;
}
int v=0;for(size_t j=0;j<len;++j)v=v*10+s[j]-'0';out[k]=v+carry;
}
} // namespace bigint_io_opt
#include <cerrno>
static bool write_all_stdout(const char* data,size_t remaining) {
while(remaining){
const ssize_t sent=write(STDOUT_FILENO,data,remaining);
if(sent>0){data+=sent;remaining-=(size_t)sent;}
else if(sent<0 && errno==EINTR)continue;
else return false;
}
return true;
}
static inline int64_t floor_base(int64_t v) {
return int64_t((uint64_t(v)+4000000000000000000ULL)/100000ULL)-40000000000000LL;
}
static inline int64_t boundary_carry(const int64_t* c,size_t end,int64_t bound) {
size_t width=6;
for(;;) {
size_t start=end>width?end-width:0;
int64_t lo=start?-bound:0,hi=start?bound:0;
for(size_t j=start;j<end;++j) {lo=floor_base(lo+c[j]);hi=floor_base(hi+c[j]);}
if(lo==hi)return lo;width*=2;
}
}
static size_t format_parallel(const int64_t* c,size_t conv,char* out,size_t digits,uint64_t) {
uint16_t pair[100];uint32_t triple[1000];
for(unsigned v=0;v<100;++v)pair[v]=('0'+v/10)|(('0'+v%10)<<8);
for(unsigned v=0;v<1000;++v){triple[v]=((uint32_t)('0'+v/100)<<8)|((uint32_t)('0'+v/10%10)<<16)|((uint32_t)('0'+v%10)<<24);}
enum { NC = 2 };
size_t m=conv-1,step=m/NC;
int64_t carry[NC]={0,boundary_carry(c,step,1LL<<50)};
char* p[NC]={out+digits,out+digits-5*step};
for(size_t i=0;i<step;++i) {
#pragma GCC unroll 8
for(int r=0;r<NC;++r) {
int64_t x=carry[r]+c[r*step+i],q=floor_base(x);
unsigned v=x-q*100000;carry[r]=q;
char* d=p[r]-5;p[r]=d;
*(uint32_t*)(d+1)=triple[v%1000];
memcpy(d,pair+v/1000,2);
}
}
size_t pos=digits-5*step*NC;
for(size_t i=(size_t)NC*step;i<m;++i) {
int64_t x=carry[NC-1]+c[i],q=floor_base(x);
unsigned v=x-q*100000;carry[NC-1]=q;pos-=5;
*(uint32_t*)(out+pos+1)=triple[v%1000];
memcpy(out+pos,pair+v/1000,2);
}
int64_t top=carry[NC-1]+c[m];
while(top>0) {out[--pos]='0'+top%10;top/=10;}
if(pos==digits)out[--pos]='0';
while(pos+1<digits&&out[pos]=='0')++pos;
return pos;
}
namespace hint { namespace transform { namespace fft {
static void real_conv5(F64* a,F64* b,size_t total,size_t used_a,size_t used_b) {
const size_t child=total/5,L=child/2;
FFTSqrtTableC64X4<6> roots(total/2,5,1);
const F64X4 c1(0.30901699437494742410229341718282),c2(-0.80901699437494742410229341718282),
s1(0.95105651629515357211643933337938),s2(0.58778525229247312916870595463907),zero(0.0);
auto negI=[&](const C64X4& z){return C64X4(z.imag,zero-z.real);};
auto posI=[&](const C64X4& z){return C64X4(zero-z.imag,z.real);};
// The second active input child ends at each operand's limb count.
// Beyond that boundary, the third radix-5 input is identically zero.
auto forward_stage=[&](auto has2_tag,F64* x,size_t j0,size_t jend){
constexpr bool HAS2=decltype(has2_tag)::value;
for(size_t j=j0;j<jend;j+=8) {
_mm_prefetch((const char*)(x+j+512),_MM_HINT_T0);
_mm_prefetch((const char*)(x+child+j+512),_MM_HINT_T0);
if constexpr(HAS2)_mm_prefetch((const char*)(x+2*child+j+512),_MM_HINT_T0);
C64X4 x0,x1,x2;
x0.load(x+j,std::true_type{});
x1.load(x+child+j,std::true_type{});
if constexpr(HAS2)x2.load(x+2*child+j,std::true_type{});
C64X4 t1,t2,u1,u2;
if constexpr(HAS2) {
t1=x0+x1*c1+x2*c2;
t2=x0+x1*c2+x2*c1;
u1=negI(x1*s1+x2*s2);
u2=negI(x1*s2-x2*s1);
} else {
t1=x0+x1*c1;
t2=x0+x1*c2;
u1=negI(x1*s1);
u2=negI(x1*s2);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
if constexpr(HAS2)(x0+x1+x2).store(x+j);
else (x0+x1).store(x+j);
(t1+u1).mul(w).store(x+child+j);
(t2+u2).mul(w2).store_nt(x+2*child+j);
(t2-u2).mul(w3).store_nt(x+3*child+j);
(t1-u1).mul(w4).store_nt(x+4*child+j);
}
};
auto forward=[&](F64* x,size_t used){
size_t active2=used>2*child?used-2*child:0;
size_t split=std::min(child,(active2+7)&~size_t(7));
forward_stage(std::true_type{},x,0,split);
forward_stage(std::false_type{},x,split,child);
_mm_sfence();
for(int r=0;r<5;++r)FFTAVX::difRecLong<false>(x+r*child,child);
};
forward(a,used_a);forward(b,used_b);
real_dot_binrev4(a,b,child,0.2);
const F64X4 inv(0.5/total);
for(int r=1;r<=2;++r) {
double theta=-HINT_2PI*r/(5*L);C64X4 offset(std::cos(theta),std::sin(theta));
BinRevTableC64X4HP<28,1> table;
for(size_t j=0;j<child;j+=8) {
if(__builtin_expect(j + 192 < child,1)) {
_mm_prefetch((const char*)(a+r*child+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(a+(6-r)*child-8-j-192),_MM_HINT_T0);
_mm_prefetch((const char*)(b+r*child+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(b+(6-r)*child-8-j-192),_MM_HINT_T0);
}
dot_rfftX4(a+r*child+j,a+(6-r)*child-8-j,b+r*child+j,b+(6-r)*child-8-j,table.iterate().mul(offset),inv);
}
}
for(int r=0;r<5;++r)FFTAVX::iditRecLong<false,false>(a+r*child,child);
/* [E7X-G1] 表示级：末级逆 radix-5 的第 5 个子带（a+4*child ..）只有落在最终卷积区间
[0, conv = na+nb-1) 内的部分才会被 `format_parallel` 读到；jlim = max(0, conv-4*child)
（向下取 8 的倍数）之外的 **5th store 与它的预取**是"写了从不被读"的纯浪费。
两路特化（`§2.19.1292` 的规程：把原路径整段照抄、条件在循环外判定）⇒ 语义逐位不变 ✓ */
const size_t conv0=used_a+used_b-1; /* [E7X-G1] 需要的输出系数个数（= format_parallel 读的上界） */
size_t jlim=0;
if(conv0>4*child){ jlim=(conv0-4*child<child)?(conv0-4*child):child; jlim&=~size_t(7); }
for(size_t j=0;j<jlim;j+=8) {
if(__builtin_expect(j + 496 < child,1)) {
_mm_prefetch((const char*)(a+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(a+child+j+316),_MM_HINT_T0);
_mm_prefetch((const char*)(a+2*child+j+376),_MM_HINT_T0);
_mm_prefetch((const char*)(a+3*child+j+436),_MM_HINT_T0);
_mm_prefetch((const char*)(a+4*child+j+496),_MM_HINT_T0);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
C64X4 x0=a+j,x1=C64X4(a+child+j).mulConj(w),x2=C64X4(a+2*child+j).mulConj(w2),
x3=C64X4(a+3*child+j).mulConj(w3),x4=C64X4(a+4*child+j).mulConj(w4);
C64X4 sum1=x1+x4,sum2=x2+x3,diff1=x1-x4,diff2=x2-x3;
C64X4 t1=x0+sum1*c1+sum2*c2,t2=x0+sum1*c2+sum2*c1,u1=posI(diff1*s1+diff2*s2),u2=posI(diff1*s2-diff2*s1);
(x0+sum1+sum2).transToI64(std::true_type{}).store(a+j,std::true_type{});
(t1+u1).transToI64(std::true_type{}).store(a+child+j,std::true_type{});
(t2+u2).transToI64(std::true_type{}).store(a+2*child+j,std::true_type{});
(t2-u2).transToI64(std::true_type{}).store(a+3*child+j,std::true_type{});
(t1-u1).transToI64(std::true_type{}).store(a+4*child+j,std::true_type{});
}
for(size_t j=jlim;j<child;j+=8) {
if(__builtin_expect(j + 496 < child,1)) {
_mm_prefetch((const char*)(a+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(a+child+j+316),_MM_HINT_T0);
_mm_prefetch((const char*)(a+2*child+j+376),_MM_HINT_T0);
_mm_prefetch((const char*)(a+3*child+j+436),_MM_HINT_T0);
}
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w),w4=w2.mul(w2);
C64X4 x0=a+j,x1=C64X4(a+child+j).mulConj(w),x2=C64X4(a+2*child+j).mulConj(w2),
x3=C64X4(a+3*child+j).mulConj(w3),x4=C64X4(a+4*child+j).mulConj(w4);
C64X4 sum1=x1+x4,sum2=x2+x3,diff1=x1-x4,diff2=x2-x3;
C64X4 t1=x0+sum1*c1+sum2*c2,t2=x0+sum1*c2+sum2*c1,u1=posI(diff1*s1+diff2*s2),u2=posI(diff1*s2-diff2*s1);
(x0+sum1+sum2).transToI64(std::true_type{}).store(a+j,std::true_type{});
(t1+u1).transToI64(std::true_type{}).store(a+child+j,std::true_type{});
(t2+u2).transToI64(std::true_type{}).store(a+2*child+j,std::true_type{});
(t2-u2).transToI64(std::true_type{}).store(a+3*child+j,std::true_type{});
}
}
}}}
#include <vector>
#include <complex>
#include <cmath>
namespace r4 {
using namespace hint::transform::fft;
static constexpr size_t RSYNC = 1024; /* twiddle segment length (lanes) */
// [E7X-R4] segmented twiddle table: ts[t] = omega_{2Lc}^{bitrev10(t)*2^{L-10}} = e^{-2*pi*i*bitrev10(t)/2048}
// (independent of Lc -- build once; segment constant folds in the rest).
static inline const std::vector<std::complex<double>>& ts_table(){
static const std::vector<std::complex<double>> t = [](){
std::vector<std::complex<double>> v(1024);
for(size_t i=0;i<1024;i++){ size_t r=0; for(size_t j=0;j<10;j++) r=(r<<1)|((i>>j)&1);
double th=-3.14159265358979323846*2.0*(double)r/2048.0;
v[i]=std::complex<double>(std::cos(th),std::sin(th)); }
return v; }();
return t; }
static inline size_t brev(size_t v, size_t w){ size_t r=0; for(size_t i=0;i<w;i++) r=(r<<1)|((v>>i)&1); return r; }
using F64=double; using C64X4=hint::Complex64X4; using F64X4=hint::Float64X4;
static inline C64X4 negI(const C64X4& z){ return C64X4(z.imag, F64X4(0.0)-z.real); }
static inline C64X4 posI(const C64X4& z){ return C64X4(F64X4(0.0)-z.imag, z.real); }
// [E4Z-23] Vertical (column) pass fusing difRecMid's layer-1 (stride sa=cs/4, table MT(cs/2))
// with its layer-2 (stride sb=cs/16, table MT(cs/8)); the rest of the recursion
// (legacy difRecMid on the sub-blocks of size sb) is byte-for-byte the original.
// The two layers and the top-stage->layer-1 fusion share the same trick: the intermediate
// never leaves L1/L2 by construction -- here it never even leaves registers (A/B arrays).
// Same dif4 + same three complex multiplies, same tables, same indices => bit-identical.
static void mid23(F64* a, size_t cs){
const size_t sa = cs>>2, sb = cs>>4;
const C64X4* TA = reinterpret_cast<const C64X4*>(FFTAVX::multi_table_1.getBegin(cs>>1));
const C64X4* TB = reinterpret_cast<const C64X4*>(FFTAVX::multi_table_1.getBegin(cs>>3));
for(size_t j0=0;j0<sb;j0+=8){
C64X4 A[4][4]; /* A[k2][t] = layer-1 output t at position k2*sb+j0 */
for(int k2=0;k2<4;k2++){
const size_t j = (size_t)k2*(sb>>3) + (j0>>3);
F64* pj = a + ((size_t)k2*sb + j0);
_mm_prefetch((const char*)(pj+sa+192+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+2*sa+200+8*k2),_MM_HINT_T0);
_mm_prefetch((const char*)(pj+3*sa+208+8*k2),_MM_HINT_T0);
C64X4 c0(pj), c1(pj+sa), c2(pj+2*sa), c3(pj+3*sa);
FFTAVX::dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
C64X4 w = TA[j], w2 = w.mul(w);
c2 = c2.mul(w); c1 = c1.mul(w2); c3 = c3.mul(w2.mul(w));
A[k2][0]=c0; A[k2][1]=c1; A[k2][2]=c2; A[k2][3]=c3;
}
for(int k1=0;k1<4;k1++){
C64X4 c0=A[0][k1], c1=A[1][k1], c2=A[2][k1], c3=A[3][k1];
FFTAVX::dif4(c0.real,c0.imag,c1.real,c1.imag,c2.real,c2.imag,c3.real,c3.imag);
C64X4 w = TB[j0>>3], w2 = w.mul(w);
c2 = c2.mul(w); c1 = c1.mul(w2); c3 = c3.mul(w2.mul(w));
F64* q = a + (size_t)k1*sa + j0;
c0.store(q); c1.store(q+sb); c2.store(q+2*sb); c3.store(q+3*sb);
}
}
for(int k1=0;k1<4;k1++) for(int k2=0;k2<4;k2++)
FFTAVX::difRecMid<false>(a + (size_t)k1*sa + (size_t)k2*sb, sb);
}
static void fwd4(F64* p, size_t child, FFTSqrtTableC64X4<6>& roots, size_t used){
// [E4Y-TL2] top stage + first recursion level fused, slice-blocked (see header).
const size_t cs = child >> 2; /* recursion level-1 stride */
const bool zq = (used <= (child<<1)); /* [E4W-ZQ] quarters 2/3 carry no operand data */
const size_t flen = child >> 1;
if(flen <= FFTAVX::MID_LEN || flen >= (size_t(1)<<22)){
// difRecLong(child) does not enter its radix-4 layer here (small inputs go straight to
// difRecMid, huge ones to fusedLong) -> keep the original two passes verbatim.
for(size_t j=0;j<child;j+=8){
C64X4 x0; x0.load(p+j, std::true_type{});
C64X4 x1; x1.load(p+child+j, std::true_type{});
C64X4 x2; x2.load(p+2*child+j, std::true_type{});
C64X4 x3; x3.load(p+3*child+j, std::true_type{});
C64X4 t0=x0+x2,t1=x0-x2,t2=x1+x3,t3=x1-x3;
C64X4 y0=t0+t2,y2=t0-t2,y1=t1+negI(t3),y3=t1+posI(t3);
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w);
y0.store(p+j); y1.mul(w).store(p+child+j); y2.mul(w2).store_nt(p+2*child+j); y3.mul(w3).store_nt(p+3*child+j);
}
_mm_sfence();
for(int r=0;r<4;++r) FFTAVX::difRecLong<false>(p+r*child,child);
return;
}
const auto &t2 = FFTAVX::sqrt_table_1[hint::hint_log2(flen)]; /* == difRecLong(child) table */
constexpr size_t SL = 1024; /* doubles per stream per slice */
for(size_t j0=0;j0<cs;j0+=SL){
size_t Tj = cs-j0; if(Tj>SL) Tj=SL;
for(int s=0;s<4;++s){
F64 *b0=p+(size_t)s*cs+j0, *b1=b0+child, *b2=b1+child, *b3=b2+child;
for(size_t j=0;j<Tj;j+=8){
// [E7X-R4] FUSED pack (unchanged): plain interleaved limbs -> split lanes in-register
C64X4 x0; x0.load(b0+j, std::true_type{});
C64X4 x1; x1.load(b1+j, std::true_type{});
C64X4 y0,y1,y2,y3;
if(!zq){
C64X4 x2; x2.load(b2+j, std::true_type{});
C64X4 x3; x3.load(b3+j, std::true_type{});
C64X4 r0=x0+x2,r1=x0-x2,r2=x1+x3,r3=x1-x3;
y0=r0+r2;y2=r0-r2;y1=r1+negI(r3);y3=r1+posI(r3);
} else {
y0=x0+x1; y1=x0+negI(x1); y2=x0-x1; y3=x0+posI(x1);
}
C64X4 w=roots[(s*cs+j0+j)/8],w2=w.mul(w),w3=w2.mul(w);
y0.store(b0+j); y1.mul(w).store(b1+j); y2.mul(w2).store(b2+j); y3.mul(w3).store(b3+j);
}
}
for(int r=0;r<4;++r){
F64* q=p+(size_t)r*child+j0;
for(size_t j=0;j<Tj;j+=8){
_mm_prefetch((const char*)(q+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+cs+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+2*cs+j+192),_MM_HINT_T0);
_mm_prefetch((const char*)(q+3*cs+j+192),_MM_HINT_T0);
C64X4 c0(q+j), c1(q+cs+j), c2(q+2*cs+j), c3(q+3*cs+j), omega1, omega2;
FFTAVX::dif4(c0.real, c0.imag, c1.real, c1.imag, c2.real, c2.imag, c3.real, c3.imag);
omega1 = t2[(j0+j)/8], c2 = c2.mul(omega1);
omega2 = omega1.mul(omega1), c1 = c1.mul(omega2);
c3 = c3.mul(omega2.mul(omega1));
c0.store(q+j), c1.store(q+cs+j), c2.store(q+2*cs+j), c3.store(q+3*cs+j);
}
}
}
_mm_sfence();
for(int r=0;r<4;++r) for(int s=0;s<4;++s){
F64* q = p+(size_t)r*child+(size_t)s*cs;
if(cs>=(size_t(1)<<15) && cs<=(size_t(1)<<18)) mid23(q, cs);
else FFTAVX::difRecMid<false>(q, cs);
}
}
template<bool TO_INT=false>
static void inv4(F64* p, size_t child, FFTSqrtTableC64X4<6>& roots, double scale=1.0){
for(int r=0;r<4;++r) FFTAVX::iditRecLong<false,false>(p+r*child,child);
for(size_t j=0;j<child;j+=8){
_mm_prefetch((const char*)(p+j+256),_MM_HINT_T0);
_mm_prefetch((const char*)(p+child+j+264),_MM_HINT_T0);
_mm_prefetch((const char*)(p+2*child+j+272),_MM_HINT_T0);
_mm_prefetch((const char*)(p+3*child+j+280),_MM_HINT_T0);
C64X4 w=roots[j/8],w2=w.mul(w),w3=w2.mul(w);
C64X4 x0=p+j,x1=C64X4(p+child+j).mulConj(w),x2=C64X4(p+2*child+j).mulConj(w2),x3=C64X4(p+3*child+j).mulConj(w3);
C64X4 t0=x0+x2,t1=x0-x2,t2=x1+x3,t3=x1-x3;
if constexpr(TO_INT){
const F64X4 sc(scale);
(((t0+t2))* sc).transToI64(std::true_type{}).store(p+j, std::true_type{});
(((t1+posI(t3)))* sc).transToI64(std::true_type{}).store(p+child+j, std::true_type{});
(((t0-t2))* sc).transToI64(std::true_type{}).store(p+2*child+j, std::true_type{});
(((t1-posI(t3)))* sc).transToI64(std::true_type{}).store(p+3*child+j, std::true_type{});
} else {
(t0+t2).store(p+j); (t1+posI(t3)).store(p+child+j); (t0-t2).store(p+2*child+j); (t1-posI(t3)).store(p+3*child+j);
}
}
}
static void dot_vec(F64* a, F64* b, size_t child, size_t total){
// [E7X-R4 route乙] vectorized radix-4 dot, FULLY SELF-CONTROLLED algebra (no black box).
// child 0 : engine's generic block routine, rescaled to this file's unscaled convention
// (fit was 2/child with scale=1 => pass scale=Lc=R/8 to get unscaled).
// children 1<->3 : flat full-mirror pair loop, MY algebra (same formulas as the verified
// scalar dot_only) + MY incremental bit-reverse twiddle generator.
// child 2 : same, self-paired full mirror.
// Twiddles: w(m) = omega_R^{4*bitrev_L(m)+r} = omega_{2Lc}^{bitrev_L(m)} * omega_R^r
// (generator validated lane-by-lane in work/e7x_series.cpp, 16/16).
// Scaling: all three blocks emit the UNSCALED product convention of dot_only
// => unpack_out uses the same 2/total as the scalar path.
const size_t R=total, Lc=child/2;
const double PI=3.14159265358979323846;
real_dot_binrev4(a, b, child, (double)Lc);
size_t L=0; while(((size_t)1<<L) < Lc) ++L;
std::vector<std::complex<double>> unit(L); // unit[j] = omega_{2Lc}^{2^j}
for(size_t j=0;j<L;j++){ double th=-2*PI*(double)((size_t)1<<j)/(double)(2.0*(double)Lc);
unit[j]=std::complex<double>(std::cos(th),std::sin(th)); }
const F64X4 half(0.5);
const C64X4 negHalfI(F64X4(0.0), F64X4(-0.5)); // 1/(2i)
auto conjC=[&](const C64X4& v){ return C64X4(v.real, F64X4(0.0)-v.imag); };
auto mulI =[&](const C64X4& v){ return C64X4(F64X4(0.0)-v.imag, v.real); };
auto run_child=[&](int r, bool selfpair){
std::complex<double> cur(1,0); size_t m=0;
double th=-2*PI*(double)r/(double)R; std::complex<double> off(std::cos(th),std::sin(th));
const std::complex<double>* TS=&ts_table()[0];
size_t seg=0; std::complex<double> segc=off;
auto next4=[&](){
if((m & (RSYNC-1))==0){
std::complex<double> sc=std::polar(1.0, -3.14159265358979323846*2.0*(double)brev(seg, L-10)/(double)(2*Lc));
segc=sc*off; ++seg;
}
const double* tp=reinterpret_cast<const double*>(TS+(m & (RSYNC-1)));
__m256d q0=_mm256_loadu_pd(tp), q1=_mm256_loadu_pd(tp+4);
__m256d tr=_mm256_permute4x64_pd(_mm256_unpacklo_pd(q0,q1),0xd8);
__m256d ti=_mm256_permute4x64_pd(_mm256_unpackhi_pd(q0,q1),0xd8);
const __m256d sr=_mm256_set1_pd(segc.real()), si=_mm256_set1_pd(segc.imag());
m+=4;
return C64X4(F64X4(_mm256_sub_pd(_mm256_mul_pd(tr,sr),_mm256_mul_pd(ti,si))),
F64X4(_mm256_add_pd(_mm256_mul_pd(tr,si),_mm256_mul_pd(ti,sr))));
};
// Paired children (1<->3): the Lc pairs (1,m)<->(3,Lc-1-m) are all DISTINCT for m in [0,Lc),
// so the sweep must cover the FULL child (child/8 blocks). A self-paired child needs only
// the first half, because one step there covers both members.
const size_t steps = selfpair ? (child/2)/8 : child/8;
// BUG FIX: the partner block must start at the END of the partner child (4-r),
// i.e. doubles (5-r)*child-8 -- the old (4-r)*child-8 landed in the wrong child
// and scribbled over it (that is what wrecked the self-paired child 2).
F64* it1=(selfpair? a+3*child-8 : a+(5-r)*child-8);
F64* jt1=(selfpair? b+3*child-8 : b+(5-r)*child-8);
F64* it0=a+r*child; F64* jt0=b+r*child;
for(size_t st=0; st<steps; st++, it0+=8, it1-=8, jt0+=8, jt1-=8){
_mm_prefetch((const char*)(it0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(it1-256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt0+256),_MM_HINT_T0);
_mm_prefetch((const char*)(jt1-256),_MM_HINT_T0);
C64X4 w = next4();
C64X4 x0 = it0, x1 = C64X4(it1).reverse();
C64X4 y0 = jt0, y1 = C64X4(jt1).reverse();
C64X4 cx1=conjC(x1), cy1=conjC(y1);
C64X4 Ea=(x0+cx1)*half, Oa=(x0-cx1).mul(negHalfI);
C64X4 Eb=(y0+cy1)*half, Ob=(y0-cy1).mul(negHalfI);
C64X4 wo=w.mul(Oa), wb=w.mul(Ob);
C64X4 Xak=Ea+wo, Xakc=Ea-wo, Xbk=Eb+wb, Xbkc=Eb-wb;
C64X4 Pk=Xak.mul(Xbk), Pkc=Xakc.mul(Xbkc);
C64X4 Ep=(Pk+Pkc)*half;
C64X4 Op=(Pk-Pkc).mulConj(w)*half; // (Pk-Pkc)/(2w), |w|=1
C64X4 Zp=Ep+mulI(Op);
C64X4 Zq=conjC(Ep-mulI(Op)); // partner position value
Zp.store(it0); Zq.reverse().store(it1); } };
run_child(1,false);
run_child(2,true);
}
static void pack_inplace(F64* p, size_t total, size_t used){
// plain real limbs (padding already zero) -> packed complex z[n]=x[2n]+i*x[2n+1],
// written back in the internal split-lane form, in place, via the engine's own permutation.
(void)used;
for(size_t j=0;j<total;j+=8){
C64X4 v; v.load(p+j, std::true_type{});
v.store(p+j); }
}
static void unpack_int64(F64* p, size_t total, double scale){
const F64X4 sc(scale);
int64_t* o=reinterpret_cast<int64_t*>(p);
for(size_t j=0;j<total;j+=8){
C64X4 v(p+j);
v = v * sc;
v.transToI64(std::true_type{}).store(o+j, std::true_type{}); }
}
static void real_conv4(F64* a, F64* b, size_t total, size_t used_a, size_t used_b){
const size_t child=total/4;
FFTSqrtTableC64X4<6> roots(total/2,4,1);
fwd4(a,child,roots,used_a); fwd4(b,child,roots,used_b); /* [E4W-ZQ] */
dot_vec(a,b,child,total);
inv4<true>(a,child,roots,2.0/(double)total); /* scale+round+interleave fused */
}
} // namespace r4
int main(){
struct stat st{};
std::unique_ptr<char[]> storage;
std::string fallback;
const char* input=nullptr;
size_t input_size=0;
auto* duck=reinterpret_cast<duck_public_stdio::DuckInfo_t*>(getauxval(duck_public_stdio::AT_DUCK));
if(duck && duck->abi_version==40 && duck->stdin_ptr){
input=duck->stdin_ptr;input_size=duck->stdin_size;
}else if(fstat(STDIN_FILENO,&st)==0 && st.st_size>0){
size_t capacity=(size_t)st.st_size;
storage.reset(new char[capacity]);
while(input_size<capacity){
size_t got=fread(storage.get()+input_size,1,capacity-input_size,stdin);
if(!got)break;
input_size+=got;
}
input=storage.get();
}else{
char block[1<<16];size_t got;
while((got=fread(block,1,sizeof block,stdin)))fallback.append(block,got);
input=fallback.data();input_size=fallback.size();
}
const char* newline=(input_size>10000000 && input[10000000]=='\n') ? input+10000000 : (const char*)memchr(input,'\n',input_size);
if(!newline)return 2;
size_t cut=(size_t)(newline-input);
size_t l1=cut;while(l1&&input[l1-1]<'0')--l1;
size_t off=cut+1;while(off<input_size&&input[off]<'0')++off;
size_t l2=input_size-off;while(l2&&input[off+l2-1]<'0')--l2;
size_t na=l1/5+1,nb=l2/5+1,conv=na+nb-1,N=std::max<size_t>(16384,hint::int_ceil2(na+nb+4)); /* [E7X-R4] pure 2-power base */
hint::AlignMem<double>A(N),B(N+256); /* +128 doubles: (A-B) % 4096 = 3072 != 0 */
// Both inputs occupy the first half; the radix-5 layer supplies the zero padding.
// A's unwritten padding is likewise zero from its fresh large Linux anonymous allocation.
// The fresh large _mm_malloc allocation is backed by zero-filled Linux
// anonymous pages in the tested environment. Leave the unused second
// half of B untouched until the FFT reads it as zeros.
bigint_io_opt::parse_balanced_base100000_dual(input,l1,A.begin(),input+off,l2,B.begin()+256);
storage.reset();std::string().swap(fallback);
// [E7Y2-PREFAULT] Explicit page pre-touch of the never-touched half of both operands.
// The two 32 MB _mm_malloc'd arrays are zero-filled Linux anonymous mappings; the parser
// writes only the low half of each (na / nb limbs). The FFT top stage is therefore the
// FIRST code to write the high half, i.e. it takes ~16 MB of write (COW/zero) faults per
// operand *in the middle of a 4-stream gather/scatter loop*, where each fault flushes the
// pipeline and resets the HW prefetchers. rdtscp phase table (local, production shape):
// first fwd4 top-stage pass 43.6M-49.5M ticks vs 7.3M/6.9M ticks for the 2nd/3rd pass over
// the very same addresses and the very same code (only the page state differs).
// One store per 4 KiB page (512 doubles) makes every page present in a tight sequential
// loop; the values written are 0.0, which is exactly what the code already requires the
// padding to be (fresh anonymous mapping), so every output byte is unchanged.
{
double* pa=A.begin(); double* pb=B.begin()+256;
const size_t PG=512; /* doubles per 4 KiB page */
for(size_t i=na;i<N;i+=PG) pa[i]=0.0;
for(size_t i=nb;i<N;i+=PG) pb[i]=0.0;
}
r4::real_conv4(A.begin(),B.begin()+256,N,na,nb);
int64_t* coeff=(int64_t*)A.begin();
const size_t digits=l1+l2;
bool direct=duck && duck->abi_version==40 && duck->stdout_ptr && duck->stdout_limit>=digits+1;
char* out=direct?duck->stdout_ptr:(char*)B.begin();
size_t pos=format_parallel(coeff,conv,out,digits,std::min(na,nb)*999ULL+1);
size_t len=digits-pos;
if(direct) {
if(pos)memmove(out,out+pos,len);
out[len]='\n';duck->stdout_size=len+1;
} else {
if(!write_all_stdout(out+pos,len) || !write_all_stdout("\n",1))return 3;
}
return 0;
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 62.535 ms | 84 MB + 352 KB | Accepted | Score: 100 | 显示更多 |