// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #103006 <https://duck.ac/submission/103006>
// （t=5806.714 us，本题当前榜首）用途：**整份复刻为基座** —— 本文件正文逐字节等于该提交。
// ★ 全量 diff 已核（本账号本轮实测）：该提交 = **本账号 #102869
// <https://duck.ac/submission/102869> 的正文 + 两处改动**（16 行差异、2 个实质 hunk）：
// (a) `TGT`/`FT` 的 target 从 `"avx2,fma"` 改成 **`"avx2,fma,tune=skylake"`**；
// (b) `IM` 偏移从 `RE_ + ARRSZ + 256` 改成 **`+ 112`**（他自述出自他自己的 #102991）。
// [2] 本账号（saffah_cc_v41_agg1）提交 #102869 <https://duck.ac/submission/102869>
// （t=5885.598 us，`work/s21_sd16.cpp`）= 本文件去掉上面两处改动后的正文；再往上是
// #102866 <https://duck.ac/submission/102866>（STRIDE k=40）与 #102860
// <https://duck.ac/submission/102860>（预取 4）。
// [3] /home/yjp/duck.ac/BRIEF.md §2.18.499①（`mine > T` 时先整份复刻）、§2.18.609
// （同族/对手在替我们做移植）、§2.18.347（复刻前先确认没有"其实是我们先做的"东西）。
// [5] duck.ac 用户 saffah_codex_6s_agg2, 提交 #102813 <https://duck.ac/submission/102813>
// 用途：直接复制了它相对本账号 #99814 的全部三处几何常量（其 diff 仅此三行，且逐函数 md5 全同）：
// `#define CBW2 8→16`、`#define STRIDE (NC+8)→(NC+24)`、`IM = RE_+ARRSZ+128→+256`。
// ======================
// ===== 思路 =====
// 局势：`T = 5806.714`（对手 #103006）、**目标全文 `0.99*T + 1us = 5749.647`**、
// 我方 mine = 5885.598（#102869）⇒ **还需快 135.951 us = 2.31%**。
// ⚠ `exact.py` 当前显示 `[1.005*T]`（对手的提交比我方新），显示缺口只有 48.850 us；
// 但**一旦我方刷新最好件口径立刻翻严** ⇒ **按 `0.99*T` 设计**（`§2.18.588`）。
//
// ★ 本发**只做复刻、不做任何改动**：`rivalcopy.py 1004` 报 #103006 与本账号 #102869 的
// fwd/rev 重合 **99.61%**，全量 diff 确认只有两处：
// (a) `tune=skylake`（编入 target 属性，改调度）；（b）**`IM` 偏移 256 → 112**。
// ★ (b) 是一条我方 notes 未扫过的轴：旧 notes 只测过 `IM` +128 / +256（+384 越界），
// 并把它记为"上限就是 +256" ⇒ **"上限"其实是"我扫过的最大值"，不是真正的边界**。
//
// 【正确性】正文与站点上 #103006 的副本逐字节相同（本脚本的 md5 断言）；
// 本地逐字节对拍作为提交前置条件。
// ===== 【本次提交（代理 zj_）】单变量：文件级 optimize 列表 3 项 -> 20 项 =====
// 【为什么这是本题唯一没被扫过的"纯 codegen"轴】本题的 pragma 只有 `O3,unroll-loops,
// rename-registers` 三项；而 `1002i` 的榜首族经过对手 **16 发只改 pragma** 的逐项扫描后收敛到
// 下面这 20 项（notes 第八/九轮，全部有判题机读数，且逐项是**加上去**的）。
// 注：本题 03:0x 那一轮用 **probe** 试过"强化 pragma"并报"无效果"，但那是**旧引擎 + probe**
// （本族 probe 与正式提交会反号，见 RULES 1.6e），**没有判题机的正式读数**。
// 【正确性】只改编译器开关 ⇒ 算术与地址不变；本地 drv + Python 大整数对拍照跑。
// 【预期（可证伪）】对方在他的基座上把 pragma 栈量到 -43 µs；我方在 1002i 上只兑现 -5 µs
// （§2.18.484① 同族刀收益重叠）。本题基底是 3 项 ⇒ 未兑现空间更大，但**符号未知**。
// 判据：< 5.813008 ms。
// ===== 【本次提交（代理 zj_）】单变量：给 `cw()` 加 `TGT`（使它可被内联） =====
// 【发现】`cw()`（twiddle 查表）是本文件里**唯一没有 `TGT`/`FT` 属性的函数**，而它被
// `diagonal_row`（每行 8 次）与 `diagonal18_row`（每行 8 次）在**最内层 j 循环里**调用。
// 本地 `g++ -O2 -S` 实测：它**没有被内联**，成了一个独立符号 `_ZL2cwjRdS_i.constprop.0`，
// 全文件 8 个调用点（`grep -c 'call.*cw'` = 8）⇒ 运行期约 **12 300 次函数调用**
// （1024 行 ×2 + 512 行 ×2，每次 8 个 j），且它是以**基线 ISA（SSE2）**编译的。
// 按每次调用 ~20~25 cycles 估 ⇒ **~0.12M cycles ≈ 34 µs**，正好是本題的缺口量级。
// 【改动】只给它的定义加 `TGT`（= `__attribute__((target("avx2,fma,tune=skylake")))`）⇒
// 调用方与它同 target ⇒ gcc 会内联（函数体小）。**算术一字不改**。
// 【正确性风险与检查】加了 target 之后理论上可能开启 FMA 收缩而改变舍入 —— 但本文件里
// `cw` 的表达式是 `a*b - c*d`，而判题端用 `-std=c++17`（ISO 模式，`-ffp-contract` 默认 off）
// ⇒ 收缩不生效。**以本地逐字节对拍为准**（判据：与基线 FNV `h=14084918943729329413` 相同）。
// 【预期（可证伪）】mine 5.783880；判据 < 5.783880 ms。缺口 **34.23 µs**。
// ===== 【本次提交（代理 zj_）】单变量：把 `pointwise_vec` 强制内联 =====
// 【同一类】`#104222`（`cw` 可内联化）= **−12.97 µs**。本发处理**最后一个**还在热路径上
// out-of-line 的函数：本地 `g++ -O2 -S` 清点，全文件只剩 `row_dif` / `pointwise_vec` /
// `kb_scan`(回退路径，不跑) / `parse_digits2`(2 次) 四个非库调用；`row_dif` 的强制内联
// 已试过（`#104223` = **+20.3 µs**，2 KB 函数体复制到 4 个点 ⇒ 取指代价压过调用开销）。
// `pointwise_vec` **比 row_dif 小得多**（每次一行对的一趟复数乘）⇒ 本发量它。
// 【改动】加 `inline __attribute__((always_inline))`。**算术一字不改**。
// 【正确性】本地 `drv` + Python 大整数对拍：与基线逐字节 MATCH（`h=14084918943729329413`）。
// 【预期（可证伪）】mine 5.770911；判据 < 5.770911 ms。缺口 **21.264 µs**。
// ===== 【本次提交（代理 zj_）】单变量：`TGT`/`FT` 的 target 追加 `bmi,bmi2,popcnt,lzcnt` =====
// 【为什么】本題**没有文件级 `#pragma GCC target`**，ISA 只来自 `TGT`/`FT`，而这两个宏只有
// `avx2,fma,tune=skylake`。对照：**1002i 的榜首件自带 `#pragma GCC target("avx2,bmi,bmi2,popcnt,lzcnt,tune=haswell")`**
// —— 同一族的另一题认为 `bmi/bmi2/popcnt/lzcnt` 值得开。本題 `parse_digits2`/`kb_scan` 里有
// `movemask + 位扫描` ⇒ 值得一次判题机定价。⚠ 它们只影响**标量**位操作（热循环是 intrinsics）。
// 【正确性】只改可用指令集（算术一字不改）⇒ 输出不变；本地对拍照跑（判据 `h=14084918943729329413`）。
// 【预期（可证伪）】mine 5.770385；判据 < 5.770385 ms。缺口 **20.738 µs**。
// ===== 【本次提交（代理 zj_）】单变量：把 `fused_row_loop` 强制内联回 `run_job` =====
// 【发现】本地 `g++ -O2 -S` 清点调用表：`cw` 内联化（`#104222` = **−12.97 µs**）与 `pointwise_vec`
// 内联化之后，**`fused_row_loop` 反而变成了 out-of-line 的 `_ZL14fused_row_loopid.constprop.0`**
// （`-finline-functions-called-once` 没吃下它，因为它把 `reformat18_row`/`row_dit`/`diagonal18_row`
// 三个内核都内联进去了）。它是**单调用点**（只从 `run_job` 调一次），且占全程序 **19.6%**。
// ⇒ 把它合回 `run_job` 可以让编译器**跨这个边界**统一调度（regalloc / 指令排布 / 常量提升）。
// 【诚实标注】同族的 `row_dif` 强制内联是 **+20.3 µs**（2 KB × 4 个调用点 ⇒ 取指代价），
// 但本函数**只有一个调用点**、且合并后 `run_job` 本来就已经包含了 `fft_forward`/`col_dif`/`col512_dit`
// ⇒ 预期不亏、小幅为正。判据就是判题机读数。
// 【正确性】只改内联 ⇒ 算术一字不改；本地对拍照跑（判据 `h=14084918943729329413`）。
// 【预期（可证伪）】mine 5.759391；判据 < 5.759391 ms。缺口 **9.744 µs**。
// ================
/* Compile-time switch defaults. These lines are part of the rank-1 source
#103006; our replica generator truncated the downloaded file at the first
`#pragma GCC optimize` line, which sits AFTER them, so `#if FUSE16>0` and
`#if NTPARSE>0` below quietly evaluated `0>0`. Restored verbatim (in the
original the first `#ifndef` wins; the later duplicate 0-valued defaults are
dead there and are omitted here, which is preprocessor-identical). */
#ifndef FUSE16
#define FUSE16 1
#endif
#ifndef NTPARSE
#define NTPARSE 1
#endif
#ifndef NOMEMSET
#define NOMEMSET 1
#endif
#ifndef PFV
#define PFV 4
#endif
#ifndef PFW
#define PFW 2
#endif
#ifndef ROWPF
#define ROWPF 0
#endif
#ifndef NTCARRY
#define NTCARRY 0
#endif
#pragma GCC optimize("O3","unroll-loops","rename-registers","live-range-shrinkage","ira-loop-pressure","modulo-sched","web","peel-loops","unswitch-loops","split-paths","gcse-after-reload","tree-vectorize","predictive-commoning","schedule-insns2","no-stack-protector","omit-frame-pointer","reorder-blocks-and-partition","sched-pressure","sched-spec-load","modulo-sched-allow-regmoves")
// 1004 -- multiply two 1,000,000-digit decimals. base 10^4 limbs, packed z=a+ib.
#include <cstdint>
#include <cstring>
#include <cmath>
#include <immintrin.h>
#ifndef CBW
#define CBW 16
#endif
#ifndef ZSKIP
#define ZSKIP 1
#endif
#ifndef CBW2
#define CBW2 16
#endif
#ifndef PF
#define PF 0
#endif
#define LOGN 19
#define NN (1<<LOGN)
#define NR (1<<10)
#define NC (1<<9)
#define STRIDE (NC+16)
#define ARRSZ ((size_t)NR*STRIDE)
typedef uint32_t u32; typedef uint64_t u64; typedef int64_t i64;
typedef double f64;
struct DI {
unsigned long abi;
const char *s; unsigned long sn;
char *o; unsigned long ol; unsigned long os;
char *e; unsigned long el; unsigned long es;
const char *IB; unsigned long IBl;
char *OB; unsigned long OBl;
unsigned long tsc;
} __attribute__((packed));
#define TGT __attribute__((target("avx2,fma,tune=skylake,bmi,bmi2,popcnt,lzcnt")))
#define FT __attribute__((target("avx2,fma,tune=skylake,bmi,bmi2,popcnt,lzcnt")))
// RE and IM must NOT be congruent mod 4096: the butterfly touches 4 rows x 8 lines of
// RE and the same of IM, and with a page-congruent IM those 64 lines fall into the same
// 8 L1 sets (8 ways -- the entire L1). Measured 29.4 -> 19.4 cyc/butterfly when the two
// arrays are 1 KB apart mod 4 KB instead of 0.
alignas(4096) static f64 RE_[2*ARRSZ + 256];
static f64 *const RE = RE_;
static f64 *const IM = RE_ + ARRSZ + 112; // +1024 B mod 4096
/* reformat18 twiddles: CJ0[j]=j/2 and CJ1[j]=256+j/2 for even j, so the destination is
plain sequential and only the DQ lookup is permuted. Precompute per kq so the vector
kernel loads four in a row. */
/* ================= compile-time tables =================
Every twiddle, bit-reversal and digit table below is a pure function of the
transform size, and the FFT sizes are compile-time constants, so the whole
~190 KB of tables is materialised by a constexpr constructor and lands in
.rodata: build_r4() + build_all() + build_tab4() together cost ~316 k cycles
judge-side (1.3% of the run, measured by rdtsc phase map), i.e. ~9500 libm
sincos calls plus the stores. Angles are evaluated in long double inside the
compiler (64-bit significand, additionally reduced Cody-Waite style), so the
materialised values are at least as accurate as the runtime libm ones (the
runtime form rounds the *angle* to double first, which alone is a ~1e-15
relative error on the twiddle - the compile-time form reduces that to ~1e-19).
Verified: the 2*10^6-digit product is byte-identical to the runtime-table
build for random inputs. */
#ifndef NN2
#define NN2 ((1<<18))
#endif
namespace ct {
constexpr long double PIL = 3.14159265358979323846264338327950288L;
constexpr long double TAU = 2.0L * PIL;
constexpr long double HPI1 = 1.5707963267948966L; /* pi/2, double part */
constexpr long double HPI2 = 6.12323399573676603587e-17L; /* pi/2 - HPI1 */
constexpr long double ffloor(long double v) {
long long t = (long long)v;
return (t > v) ? (long double)(t - 1) : (long double)t;
}
constexpr long double sin_series(long double x) { /* |x| <= pi/4 */
long double x2 = x * x, t = x, s = x;
for (int k = 1; k < 13; ++k) { t *= -x2 / (long double)((2*k) * (2*k + 1)); s += t; }
return s;
}
constexpr long double cos_series(long double x) { /* |x| <= pi/4 */
long double x2 = x * x, t = 1.0L, s = 1.0L;
for (int k = 1; k < 13; ++k) { t *= -x2 / (long double)((2*k - 1) * (2*k)); s += t; }
return s;
}
constexpr void sincos_l(long double x, long double &sn, long double &cs) {
long double k = ffloor(x / PIL * 2.0L + 0.5L);
long double r = (x - k * HPI1) - k * HPI2;
int q = (int)(k - 4.0L * ffloor(k * 0.25L));
long double s = sin_series(r), c = cos_series(r);
if (q == 0) { sn = s; cs = c; }
else if (q == 1) { sn = c; cs = -s; }
else if (q == 2) { sn = -s; cs = -c; }
else { sn = -c; cs = s; }
}
constexpr long double SN(long double x) { long double s = 0, c = 0; sincos_l(x, s, c); return s; }
constexpr long double CS(long double x) { long double s = 0, c = 0; sincos_l(x, s, c); return c; }
}
struct Tables {
f64 TWR[2][NC], TWI[2][NC];
f64 RDr[NC], RDi[NC], RIr[NC], RIi[NC];
f64 CDr[NR], CDi[NR], CIr[NR], CIi[NR];
f64 DQr[NR], DQi[NR], DTr[NC], DTi[NC];
u32 BREV[NR], BREVC[NC], BREV9[512], CJ0[NC], CJ1[NC];
f64 DQ18[512], DQ18i[512], DT18[NC], DT18i[NC];
u32 RVIX[NR], CVIX0[NC], CVIX1[NC];
f64 T1r[400], T1i[400], T2r[400], T2i[400], T3r[400], T3i[400];
f64 RT1r[176],RT1i[176],RT2r[176],RT2i[176],RT3r[176],RT3i[176];
f64 RDT1r[168],RDT1i[168],RDT2r[168],RDT2i[168],RDT3r[168],RDT3i[168];
f64 CDT1r[85],CDT1i[85],CDT2r[85],CDT2i[85],CDT3r[85],CDT3i[85];
u32 TAB4[10000];
constexpr Tables()
: TWR(), TWI(), RDr(), RDi(), RIr(), RIi(), CDr(), CDi(), CIr(), CIi(),
DQr(), DQi(), DTr(), DTi(), BREV(), BREVC(), BREV9(), CJ0(), CJ1(),
DQ18(), DQ18i(), DT18(), DT18i(), RVIX(), CVIX0(), CVIX1(),
T1r(), T1i(), T2r(), T2i(), T3r(), T3i(),
RT1r(),RT1i(),RT2r(),RT2i(),RT3r(),RT3i(),
RDT1r(),RDT1i(),RDT2r(),RDT2i(),RDT3r(),RDT3i(),
CDT1r(),CDT1i(),CDT2r(),CDT2i(),CDT3r(),CDT3i(), TAB4() {
using ct::PIL; using ct::TAU; using ct::SN; using ct::CS;
/* ---- build_r4 ---- */
int k = 0;
for (int h = NR>>2; h >= 1; h >>= 2) {
long double w = -TAU / (long double)(4*h);
for (int j = 0; j < h; j++) {
T1r[k+j] = (f64)CS(w*j); T1i[k+j] = (f64)SN(w*j);
T2r[k+j] = (f64)CS(2*w*j); T2i[k+j] = (f64)SN(2*w*j);
T3r[k+j] = (f64)CS(3*w*j); T3i[k+j] = (f64)SN(3*w*j);
}
k += h;
}
k = 0;
for (int h = NC>>2; h >= 8; h >>= 2) {
long double w = -TAU / (long double)(4*h);
for (int j = 0; j < h; j++) {
RT1r[k+j] = (f64)CS(w*j); RT1i[k+j] = (f64)SN(w*j);
RT2r[k+j] = (f64)CS(2*w*j); RT2i[k+j] = (f64)SN(2*w*j);
RT3r[k+j] = (f64)CS(3*w*j); RT3i[k+j] = (f64)SN(3*w*j);
}
k += h;
}
{ int m = 0;
for (int h = 8; h < NC; h <<= 2) {
long double aa = PIL/(long double)h;
for (int j = 0; j < h; j++) {
long double t = aa*(long double)j, u = t*0.5L;
RDT1r[m+j]=(f64)CS(t); RDT1i[m+j]=(f64)SN(t);
RDT2r[m+j]=(f64)CS(u); RDT2i[m+j]=(f64)SN(u);
RDT3r[m+j]=(f64)(-SN(u)); RDT3i[m+j]=(f64)CS(u);
}
m += h; } }
{ int m = 0;
for (int h = 1; h*4 <= 512; h <<= 2) {
long double aa = PIL/(long double)h;
for (int j = 0; j < h; j++) {
long double t = aa*(long double)j, u = t*0.5L;
CDT1r[m+j]=(f64)CS(t); CDT1i[m+j]=(f64)SN(t);
CDT2r[m+j]=(f64)CS(u); CDT2i[m+j]=(f64)SN(u);
CDT3r[m+j]=(f64)(-SN(u)); CDT3i[m+j]=(f64)CS(u);
}
m += h; } }
/* ---- build_all ---- */
for (int q = 0; q < NR; q++) { long double a = -TAU*(long double)(q*NC)/(long double)NN;
DQr[q]=(f64)CS(a); DQi[q]=(f64)SN(a); }
for (int t = 0; t < NC; t++) { long double a = -TAU*(long double)t/(long double)NN;
DTr[t]=(f64)CS(a); DTi[t]=(f64)SN(a); }
for (int i = 0; i < NR; i++) { u32 r = 0; for (int b = 0; b < 10; b++) if (i & (1u<<b)) r |= 1u << (9-b); BREV[i] = r; }
for (int i = 0; i < 512; i++) { u32 r = 0; for (int b = 0; b < 9; b++) if (i & (1u<<b)) r |= 1u << (8-b); BREV9[i] = r; }
for (int q = 0; q < 512; q++) { long double a = -TAU*(long double)(q*512)/(long double)NN2;
DQ18[q]=(f64)CS(a); DQ18i[q]=(f64)SN(a); }
for (int t = 0; t < NC; t++) { long double a = -TAU*(long double)t/(long double)NN2;
DT18[t]=(f64)CS(a); DT18i[t]=(f64)SN(a); }
for (int j = 0; j < NC; j++) { u32 r = 0; for (int b = 0; b < 9; b++) if (j & (1u<<b)) r |= 1u << (8-b); BREVC[j] = r; }
for (int i = 0; i < NR; i++) { u32 k1 = BREV[i]; u32 k1n = (NR-k1)%NR; RVIX[i] = BREV[k1n]; }
for (int j = 0; j < NC; j++) { u32 k2 = BREVC[j];
CVIX0[j] = BREVC[(NC-k2)%NC]; CVIX1[j] = BREVC[(NC-k2+NC-1)%NC]; }
for (int j = 0; j < NC; j += 2) { u32 k2 = BREVC[j];
for (int par = 0; par < 2; par++) { u32 k2p = par + 2*k2; u32 cc = BREVC[k2p & 511];
if (par) CJ1[j] = cc; else CJ0[j] = cc; } }
{ int n = NC, kk = 0;
for (int h = n>>1; h >= 1; h >>= 1) { for (int j = 0; j < h; j++) { long double a = -PIL*(long double)j/(long double)h;
RDr[kk+j]=(f64)CS(a); RDi[kk+j]=(f64)SN(a);} kk += h; }
kk = n-1;
for (int h = 1; h < n; h <<= 1) { kk -= h;
for (int j = 0; j < h; j++) { long double a = PIL*(long double)j/(long double)h;
RIr[kk+j]=(f64)CS(a); RIi[kk+j]=(f64)SN(a);} } }
for (int kq = 0; kq < 2; kq++)
for (int c = 0; c < NC; c++) { u32 eh = kq + 2u*BREVC[(2u*c) & 511u]; TWR[kq][c] = DQr[eh]; TWI[kq][c] = DQi[eh]; }
{ int n = NR, kk = 0;
for (int h = n>>1; h >= 1; h >>= 1) { for (int j = 0; j < h; j++) { long double a = -PIL*(long double)j/(long double)h;
CDr[kk+j]=(f64)CS(a); CDi[kk+j]=(f64)SN(a);} kk += h; }
kk = n-1;
for (int h = 1; h < n; h <<= 1) { kk -= h;
for (int j = 0; j < h; j++) { long double a = PIL*(long double)j/(long double)h;
CIr[kk+j]=(f64)CS(a); CIi[kk+j]=(f64)SN(a);} } }
/* ---- build_tab4 ---- */
for (u32 d = 0; d < 10000; d++)
TAB4[d] = (u32)('0'+d/1000) | ((u32)('0'+(d/100)%10) << 8)
| ((u32)('0'+(d/10)%10) << 16) | ((u32)('0'+d%10) << 24);
}
};
static constexpr Tables TB{};
#define TWR TB.TWR
#define TWI TB.TWI
#define RDr TB.RDr
#define RDi TB.RDi
#define RIr TB.RIr
#define RIi TB.RIi
#define CDr TB.CDr
#define CDi TB.CDi
#define CIr TB.CIr
#define CIi TB.CIi
#define DQr TB.DQr
#define DQi TB.DQi
#define DTr TB.DTr
#define DTi TB.DTi
#define BREV TB.BREV
#define BREVC TB.BREVC
#define BREV9 TB.BREV9
#define CJ0 TB.CJ0
#define CJ1 TB.CJ1
#define DQ18 TB.DQ18
#define DQ18i TB.DQ18i
#define DT18 TB.DT18
#define DT18i TB.DT18i
#define RVIX TB.RVIX
#define CVIX0 TB.CVIX0
#define CVIX1 TB.CVIX1
#define T1r TB.T1r
#define T1i TB.T1i
#define T2r TB.T2r
#define T2i TB.T2i
#define T3r TB.T3r
#define T3i TB.T3i
#define RT1r TB.RT1r
#define RT1i TB.RT1i
#define RT2r TB.RT2r
#define RT2i TB.RT2i
#define RT3r TB.RT3r
#define RT3i TB.RT3i
#define RDT1r TB.RDT1r
#define RDT1i TB.RDT1i
#define RDT2r TB.RDT2r
#define RDT2i TB.RDT2i
#define RDT3r TB.RDT3r
#define RDT3i TB.RDT3i
#define CDT1r TB.CDT1r
#define CDT1i TB.CDT1i
#define CDT2r TB.CDT2r
#define CDT2i TB.CDT2i
#define CDT3r TB.CDT3r
#define CDT3i TB.CDT3i
#define TAB4 TB.TAB4
FT static void build_r4(void);
FT static void build_all(void) { }
FT static inline void low3_dif(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0,-0.7071067811865475244,-1.0, -0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d sr=_mm256_add_pd(r0,r1), si=_mm256_add_pd(i0,i1);
__m256d dr=_mm256_sub_pd(r0,r1), di=_mm256_sub_pd(i0,i1);
__m256d t2=_mm256_mul_pd(di,wi), drr=_mm256_fmsub_pd(dr,wr,t2);
__m256d t3=_mm256_mul_pd(dr,wi), dii=_mm256_fmadd_pd(di,wr,t3);
r0=sr;i0=si;r1=drr;i1=dii; }
#define H2D(V,W) do{ \
__m256d tv=_mm256_permute2f128_pd(V,V,1), tw=_mm256_permute2f128_pd(W,W,1); \
__m256d sV=_mm256_add_pd(V,tv), dV=_mm256_sub_pd(V,tv); \
__m256d sW=_mm256_add_pd(W,tw), dW=_mm256_sub_pd(W,tw); \
__m256d M=_mm256_permute2f128_pd(dV,dW,0x20); \
__m256d P1=_mm256_permute4x64_pd(M,0xC0); \
__m256d P2=_mm256_permute4x64_pd(M,_MM_SHUFFLE(1,2,0,0)); \
__m256d sg=_mm256_set_pd(-1.0,1.0,1.0,1.0); \
V=_mm256_blend_pd(sV,P1,0xC); \
W=_mm256_blend_pd(sW,_mm256_mul_pd(P2,sg),0xC); \
}while(0)
H2D(r0,i0); H2D(r1,i1);
#undef H2D
#define H1D(V,W) do{ \
__m256d tv=_mm256_permute_pd(V,5), tw=_mm256_permute_pd(W,5); \
V=_mm256_blend_pd(_mm256_add_pd(V,tv),_mm256_permute4x64_pd(_mm256_sub_pd(V,tv),0x80),0xA); \
W=_mm256_blend_pd(_mm256_add_pd(W,tw),_mm256_permute4x64_pd(_mm256_sub_pd(W,tw),0x80),0xA); \
}while(0)
H1D(r0,i0); H1D(r1,i1);
#undef H1D
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
FT static inline void low3_dit(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
#define H1I(V,W) do{ \
__m256d tv=_mm256_permute_pd(V,5), tw=_mm256_permute_pd(W,5); \
V=_mm256_blend_pd(_mm256_add_pd(V,tv),_mm256_permute4x64_pd(_mm256_sub_pd(V,tv),0x80),0xA); \
W=_mm256_blend_pd(_mm256_add_pd(W,tw),_mm256_permute4x64_pd(_mm256_sub_pd(W,tw),0x80),0xA); \
}while(0)
H1I(r0,i0); H1I(r1,i1);
#undef H1I
#define H2I(V,W) do{ \
__m256d tv=_mm256_permute2f128_pd(V,V,1), tw=_mm256_permute2f128_pd(W,W,1); \
__m256d Ur=_mm256_blend_pd(tv,_mm256_sub_pd(_mm256_setzero_pd(),tw),0x2); \
__m256d Ui=_mm256_blend_pd(tw,tv,0x2); \
__m256d sV=_mm256_add_pd(V,Ur), dV=_mm256_sub_pd(V,Ur); \
__m256d sW=_mm256_add_pd(W,Ui), dW=_mm256_sub_pd(W,Ui); \
V=_mm256_blend_pd(sV,_mm256_permute4x64_pd(dV,0x40),0xC); \
W=_mm256_blend_pd(sW,_mm256_permute4x64_pd(dW,0x40),0xC); \
}while(0)
H2I(r0,i0); H2I(r1,i1);
#undef H2I
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0, 0.7071067811865475244, 1.0, 0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d t2=_mm256_mul_pd(i1,wi), prr=_mm256_fmsub_pd(r1,wr,t2);
__m256d t3=_mm256_mul_pd(r1,wi), pii=_mm256_fmadd_pd(i1,wr,t3);
__m256d s=_mm256_add_pd(r0,prr), d=_mm256_sub_pd(r0,prr);
__m256d s2=_mm256_add_pd(i0,pii), d2=_mm256_sub_pd(i0,pii);
r0=s;i0=s2;r1=d;i1=d2; }
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
// Radix-4 DIF butterfly with outputs written as soon as each is complete. The
// monolithic form (compute all of y0,y1,y2,y3 then store) needs ~24 live YMM values
// against 16 architectural registers, so gcc spilled and the loop ran at IPC ~0.9.
#define R4BODY(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
#define R4BODY0(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_setzero_pd(),x2i=x2r,x3r=x2r,x3i=x2r; /* rows 512..1023 are exactly +0.0 */ \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
// Radix-4 DIT butterfly, fusing the radix-2 DIT stages of half-sizes h and 2h.
// Rows are (s+j, s+j+h, s+j+2h, s+j+3h); W1 = exp(i*pi*j/h) is the stage-h twiddle
// (used twice, on x1 and x3); W2 = exp(i*pi*j/(2h)) is the stage-2h twiddle for the
// (A,C) pair and W3 = i*W2 the one for the (B,D) pair, whose within-group offset is
// j+h, i.e. exp(i*pi*j/(2h))*exp(i*pi/2).
#define R4DIT(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d a1r=_mm256_fmsub_pd(x1r,W1R,_mm256_mul_pd(x1i,W1I)); \
__m256d a1i=_mm256_fmadd_pd(x1r,W1I,_mm256_mul_pd(x1i,W1R)); \
__m256d a3r=_mm256_fmsub_pd(x3r,W1R,_mm256_mul_pd(x3i,W1I)); \
__m256d a3i=_mm256_fmadd_pd(x3r,W1I,_mm256_mul_pd(x3i,W1R)); \
__m256d Ar=_mm256_add_pd(x0r,a1r), Ai=_mm256_add_pd(x0i,a1i); \
__m256d Br=_mm256_sub_pd(x0r,a1r), Bi=_mm256_sub_pd(x0i,a1i); \
__m256d Cr=_mm256_add_pd(x2r,a3r), Ci=_mm256_add_pd(x2i,a3i); \
__m256d Dr=_mm256_sub_pd(x2r,a3r), Di=_mm256_sub_pd(x2i,a3i); \
__m256d c2r=_mm256_fmsub_pd(Cr,W2R,_mm256_mul_pd(Ci,W2I)); \
__m256d c2i=_mm256_fmadd_pd(Cr,W2I,_mm256_mul_pd(Ci,W2R)); \
__m256d d2r=_mm256_fmsub_pd(Dr,W3R,_mm256_mul_pd(Di,W3I)); \
__m256d d2i=_mm256_fmadd_pd(Dr,W3I,_mm256_mul_pd(Di,W3R)); \
_mm256_store_pd(R0+c,_mm256_add_pd(Ar,c2r)); \
_mm256_store_pd(I0+c,_mm256_add_pd(Ai,c2i)); \
_mm256_store_pd(R2+c,_mm256_sub_pd(Ar,c2r)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(Ai,c2i)); \
_mm256_store_pd(R1+c,_mm256_add_pd(Br,d2r)); \
_mm256_store_pd(I1+c,_mm256_add_pd(Bi,d2i)); \
_mm256_store_pd(R3+c,_mm256_sub_pd(Br,d2r)); \
_mm256_store_pd(I3+c,_mm256_sub_pd(Bi,d2i)); }while(0)
#define R4BODYT(R0,R1,R2,R3,I0,I1,I2,I3) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
_mm256_store_pd(R1+c,_mm256_sub_pd(ar,cr)); \
_mm256_store_pd(I1+c,_mm256_sub_pd(ai,ci)); \
_mm256_store_pd(R2+c,_mm256_add_pd(br,di)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(bi,dr)); \
_mm256_store_pd(R3+c,_mm256_sub_pd(br,di)); \
_mm256_store_pd(I3+c,_mm256_add_pd(bi,dr)); }while(0)
#ifndef NST
#define NST 5
#endif
FT static void row_dif(f64 *pr, f64 *pi) {
const int n=NC; int k=0;
// Radix-4 DIF: fuses the radix-2 stages (256,128),(64,32),(16,8) into three passes;
// low3_dif still covers (4,2,1), so the total is the same 9 radix-2 equivalents and the
// output ordering is bit-for-bit what the radix-2 kernel produced.
for (int h=n>>2;h>=8;k+=h,h>>=2)
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RT1r+k+j), w1i=_mm256_load_pd(RT1i+k+j);
__m256d w2r=_mm256_load_pd(RT2r+k+j), w2i=_mm256_load_pd(RT2i+k+j);
__m256d w3r=_mm256_load_pd(RT3r+k+j), w3i=_mm256_load_pd(RT3i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4BODY(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int s=0;s<n;s+=16) { low3_dif(pr+s,pi+s); low3_dif(pr+s+8,pi+s+8); }
}
// Row-pass radix-4 twiddles. Fusing the radix-2 DIF stages (256,128),(64,32),(16,8) of the
// 512-point row transform needs W_t = exp(-2*pi*i*t*j/(2*h2)) with h2 the LARGER radix-2
// half-size, i.e. h_code=h2/2 and w=-2*pi/(4*h_code) -- the same convention build_r4 uses
// for the column pass, so the tables cannot be shared.
FT static void build_r4(void) { }
FT static void col_dif(int j0){ const int n=NR;
// Four-step column pass. Twiddle table offsets: the radix-4 passes fuse the radix-2
// half-sizes (512,256),(128,64),(32,16),(8,4),(2,1); build_r4 lays them out at
// cumulative offsets 0,256,320,336,340.
// The tail passes are TILED so that each pass works on the smallest window its
// butterflies need: h=64 over a 256-row window, h=16 over 64 rows, h=4 over 16 rows
// and h=1 over 4 rows -- those last two are L1-resident, so their re-reads never
// reach L2/L3.
{
const int h=n>>2; const int K=0;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K+j]),w1i=_mm256_set1_pd(T1i[K+j]);
__m256d w2r=_mm256_set1_pd(T2r[K+j]),w2i=_mm256_set1_pd(T2i[K+j]);
__m256d w3r=_mm256_set1_pd(T3r[K+j]),w3i=_mm256_set1_pd(T3i[K+j]);
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
f64 *r0=RE+(size_t)j*STRIDE+j0,*i0=IM+(size_t)j*STRIDE+j0;
/* FIXED row prefetch (taken from the rival's current #99680, see REFERENCES):
the previous form advanced a char* by st1 DOUBLES (byte distance 8x too
small) so it named lines far from anything the pass reads, and `o<CBW`
(16) covered only the first of the two 64 B lines of a row segment.
Corrected: byte stride, both 64 B lines, lookahead 4 rows. */
if (j+4<h){
const char *p0=(const char*)(RE+(size_t)(j+4)*STRIDE+j0);
const char *q0=(const char*)(IM+(size_t)(j+4)*STRIDE+j0);
const size_t stp=(size_t)h*STRIDE*8;
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); } }
for(int s=0;s<n;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
#if ZSKIP
for(int c=0;c<CBW;c+=4) R4BODY0(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#else
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#endif
r0+=st4; i0+=st4; }
}
}
const int GW=n>>2, K64=n>>2, K16=K64+(n>>4), K4=K16+(n>>6);
for (int g=0; g<n; g+=GW) {
{ const int h=GW>>2; // 64
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K64+j]),w1i=_mm256_set1_pd(T1i[K64+j]);
__m256d w2r=_mm256_set1_pd(T2r[K64+j]),w2i=_mm256_set1_pd(T2i[K64+j]);
__m256d w3r=_mm256_set1_pd(T3r[K64+j]),w3i=_mm256_set1_pd(T3i[K64+j]);
f64 *r0=RE+(size_t)(g+j)*STRIDE+j0,*i0=IM+(size_t)(g+j)*STRIDE+j0;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
for(int s=g;s<g+GW;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
r0+=st4; i0+=st4; }
}
}
for (int g2=g; g2<g+GW; g2+=64) {
{ const int h=16;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K16+j]),w1i=_mm256_set1_pd(T1i[K16+j]);
__m256d w2r=_mm256_set1_pd(T2r[K16+j]),w2i=_mm256_set1_pd(T2i[K16+j]);
__m256d w3r=_mm256_set1_pd(T3r[K16+j]),w3i=_mm256_set1_pd(T3i[K16+j]);
f64 *r0=RE+(size_t)(g2+j)*STRIDE+j0,*i0=IM+(size_t)(g2+j)*STRIDE+j0;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g2+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g2+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g3=g2; g3<g2+64; g3+=16) {
{ const int h=4;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K4+j]),w1i=_mm256_set1_pd(T1i[K4+j]);
__m256d w2r=_mm256_set1_pd(T2r[K4+j]),w2i=_mm256_set1_pd(T2i[K4+j]);
__m256d w3r=_mm256_set1_pd(T3r[K4+j]),w3i=_mm256_set1_pd(T3i[K4+j]);
f64 *r0=RE+(size_t)(g3+j)*STRIDE+j0,*i0=IM+(size_t)(g3+j)*STRIDE+j0;
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g4=g3; g4<g3+16; g4+=4) {
f64 *r0=RE+(size_t)g4*STRIDE+j0,*i0=IM+(size_t)g4*STRIDE+j0;
f64 *r1=r0+STRIDE,*i1=i0+STRIDE, *r2=r1+STRIDE,*i2=i1+STRIDE, *r3=r2+STRIDE,*i3=i2+STRIDE;
for(int c=0;c<CBW;c+=4) R4BODYT(r0,r1,r2,r3,i0,i1,i2,i3);
}
}
}
}
}
FT static void row_dit(f64 *pr, f64 *pi) {
const int n=NC;
for (int s=0;s<n;s+=16) { low3_dit(pr+s,pi+s); low3_dit(pr+s+8,pi+s+8); }
// Fused radix-4 DIT over the remaining stages (8,16),(32,64),(128,256).
int k=0;
for (int h=8; h*4<=n; h<<=2) {
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RDT1r+k+j), w1i=_mm256_load_pd(RDT1i+k+j);
__m256d w2r=_mm256_load_pd(RDT2r+k+j), w2i=_mm256_load_pd(RDT2i+k+j);
__m256d w3r=_mm256_load_pd(RDT3r+k+j), w3i=_mm256_load_pd(RDT3i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4DIT(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
}
k+=h;
}
}
FT static void col_dit(int j0) {
const int n=NR; int k=n-1;
for (int h=1;h<n;h<<=1) { k-=h;
for (int s=0;s<n;s+=2*h)
for (int j=0;j<h;j++) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
__m256d wr=_mm256_set1_pd(CIr[k+j]), wi=_mm256_set1_pd(CIi[k+j]);
for (int c=0;c<128;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
TGT static inline void cw(u32 e, f64 &r, f64 &i, int sign) {
f64 rr = DQr[e>>9]*DTr[e&(NC-1)] - DQi[e>>9]*DTi[e&(NC-1)];
f64 ii = DQr[e>>9]*DTi[e&(NC-1)] + DQi[e>>9]*DTr[e&(NC-1)];
if (sign<0) ii = -ii;
r=rr; i=ii;
}
FT static inline void diagonal_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
/* four independent twiddle chains: the base advances by r16=r4^4 per 16 elements,
and chains 1..3 start at base*r4, base*r4^2, base*r4^3. A single chain was
latency-bound (one dependent complex multiply per 4 elements). */
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i; /* r4^2 */
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i; /* r4^4 */
for (int j=0;j<NC;j+=64) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN-1));
f64 base_r, base_i; cw(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
#if PFDIAG>0
if (j+64<NC) { _mm_prefetch((const char*)(pr+j+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+j+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+96),_MM_HINT_T0); }
#endif
for (int t=0;t<64;t+=16) {
f64 br_,bi_; int off;
#define DCHAIN(BR,BI,OFF) do{ \
__m256d br=_mm256_set1_pd(BR), bi=_mm256_set1_pd(BI); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
DCHAIN(base_r,base_i,0); DCHAIN(b1r,b1i,4);
DCHAIN(b2r,b2i,8); DCHAIN(b3r,b3i,12);
#undef DCHAIN
(void)br_;(void)bi_;(void)off;
f64 n0r=base_r*r16r-base_i*r16i, n0i=base_r*r16i+base_i*r16r;
f64 n1r=b1r*r16r-b1i*r16i, n1i=b1r*r16i+b1i*r16r;
f64 n2r=b2r*r16r-b2i*r16i, n2i=b2r*r16i+b2i*r16r;
f64 n3r=b3r*r16r-b3i*r16i, n3i=b3r*r16i+b3i*r16r;
base_r=n0r;base_i=n0i; b1r=n1r;b1i=n1i; b2r=n2r;b2i=n2i; b3r=n3r;b3i=n3i;
}
}
}
}
FT static void diagonal(int sign, f64 sc) {
for (int i=0;i<NR;i++) diagonal_row(i, sign, sc);
}
FT static inline void cw18(u32 e, f64 &r, f64 &i, int sign) {
r = DQ18[e>>9]*DT18[e&(NC-1)] - DQ18i[e>>9]*DT18i[e&(NC-1)];
i = DQ18[e>>9]*DT18i[e&(NC-1)] + DQ18i[e>>9]*DT18[e&(NC-1)];
if (sign<0) i = -i;
}
FT static inline void diagonal18_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV9[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw18(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i;
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i;
for (int j=0;j<NC;j+=64) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN2-1));
f64 base_r, base_i; cw18(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
for (int t=0;t<64;t+=16) {
#define ECHAIN(BR,BI,OFF) do{ \
__m256d br=_mm256_set1_pd(BR), bi=_mm256_set1_pd(BI); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
ECHAIN(base_r,base_i,0); ECHAIN(b1r,b1i,4);
ECHAIN(b2r,b2i,8); ECHAIN(b3r,b3i,12);
#undef ECHAIN
f64 n0r=base_r*r16r-base_i*r16i, n0i=base_r*r16i+base_i*r16r;
f64 n1r=b1r*r16r-b1i*r16i, n1i=b1r*r16i+b1i*r16r;
f64 n2r=b2r*r16r-b2i*r16i, n2i=b2r*r16i+b2i*r16r;
f64 n3r=b3r*r16r-b3i*r16i, n3i=b3r*r16i+b3i*r16r;
base_r=n0r;base_i=n0i; b1r=n1r;b1i=n1i; b2r=n2r;b2i=n2i; b3r=n3r;b3i=n3i;
}
}
}
}
FT static void diagonal18(int sign, f64 sc) {
for (int i=0;i<512;i++) diagonal18_row(i, sign, sc);
}
// Radix-4 DIT along the strided column axis: fuse the stage pairs (1,2),(4,8),(16,32),
// (64,128), leaving only h=256 as a radix-2 pass -- 5 passes instead of 9 and ~35% fewer
// instructions. The CDT tables (built for row_dit) have exactly the (h, j) layout needed.
// Unlike row_dit, the vector here spans 4 COLUMNS (broadcast twiddles), because a column
// block is many columns wide and only the position index j carries a twiddle.
FT static void col512_dit(int j0) {
const int n=512;
int m=0;
for (int h=1; h*4<=n; h<<=2) {
for (int j=0;j<h;j++) {
__m256d w1r=_mm256_set1_pd(CDT1r[m+j]), w1i=_mm256_set1_pd(CDT1i[m+j]);
__m256d w2r=_mm256_set1_pd(CDT2r[m+j]), w2i=_mm256_set1_pd(CDT2i[m+j]);
__m256d w3r=_mm256_set1_pd(CDT3r[m+j]), w3i=_mm256_set1_pd(CDT3i[m+j]);
for (int s=0;s<n;s+=4*h) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
f64 *cr=br+(size_t)h*STRIDE, *ci=bi+(size_t)h*STRIDE;
f64 *dr=cr+(size_t)h*STRIDE, *di=ci+(size_t)h*STRIDE;
for (int c=0;c<CBW2;c+=4)
R4DIT(ar,br,cr,dr,ai,bi,ci,di,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
m+=h;
}
{ const int h=256; // last stage, twiddles live at RIr[0..255]
for (int j=0;j<h;j++) {
__m256d wr=_mm256_set1_pd(RIr[j]), wi=_mm256_set1_pd(RIi[j]);
f64 *ar=RE+(size_t)j*STRIDE+j0, *ai=IM+(size_t)j*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
for (int c=0;c<CBW2;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
alignas(64) static f64 SVr[2][NC], SVi[2][NC];
#ifndef R18VEC
#define R18VEC 1
#endif
#ifndef PFR18
#define PFR18 1
#endif
#ifndef PFDIAG
#define PFDIAG 1
#endif
#ifndef PFPAIR
#define PFPAIR 0
#endif
/* reformat18 twiddles: CJ0[j]=j/2, CJ1[j]=256+j/2 for even j, so the destination is
plain sequential and only the DQ lookup is permuted. Precompute it per kq so the
vector kernel can load four in a row. */
#if R18VEC>0
/* Vectorised reformat18. ccol is exactly j/2 (CJ0) or 256+j/2 (CJ1), so the writes are
sequential; the only irregular part is the DQ lookup, which TWR/TWI linearise. The
FMA contraction order is the one gcc-9 emits for the scalar form (verified from asm),
so this is bit-identical to it. */
FT static inline void reformat18_row(int i) {
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
f64 *wr = RE+(size_t)r*STRIDE + (par?256:0), *wi = IM+(size_t)r*STRIDE + (par?256:0);
u32 k1 = BREV[i]; u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
const f64 *TR = TWR[kq], *TI = TWI[kq];
const __m256d sg=_mm256_set1_pd(-0.0);
const __m256d Vr=_mm256_set1_pd(dr), Vi=_mm256_set1_pd(di);
for (int c=0;c<(NC>>1);c+=4) {
#if PFR18>0
if (c+32 < (NC>>1)) { _mm_prefetch((const char*)(pr+2*c+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+2*c+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+96),_MM_HINT_T0); }
#endif
__m256d A=_mm256_load_pd(pr+2*c), B=_mm256_load_pd(pr+2*c+4);
// even/odd deinterleave with 2 permutes instead of 3 (identical values):
// P=perm4x64(A,0xD8)=[a0,a2,a1,a3]; the 128-halves of (P,Q) give even and odd.
__m256d P=_mm256_permute4x64_pd(A,0xD8), Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1r=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2r=_mm256_permute2f128_pd(P,Q,0x31);
A=_mm256_load_pd(pi+2*c); B=_mm256_load_pd(pi+2*c+4);
P=_mm256_permute4x64_pd(A,0xD8); Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1i=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2i=_mm256_permute2f128_pd(P,Q,0x31);
__m256d qr0=_mm256_load_pd(TR+c), qi0=_mm256_load_pd(TI+c);
__m256d er=_mm256_add_pd(z1r,z2r), dr2=_mm256_sub_pd(z1r,z2r);
__m256d ei=_mm256_add_pd(z1i,z2i), di2=_mm256_sub_pd(z1i,z2i);
__m256d t0=_mm256_mul_pd(qr0,Vr), tr=_mm256_fnmadd_pd(qi0,Vi,t0);
__m256d t1=_mm256_mul_pd(qr0,Vi), ti=_mm256_fnmadd_pd(qi0,Vr,_mm256_xor_pd(t1,sg));
__m256d t2=_mm256_mul_pd(tr,dr2), orr=_mm256_fnmadd_pd(ti,di2,t2);
__m256d t3=_mm256_mul_pd(tr,di2), oii=_mm256_fmadd_pd(ti,dr2,t3);
_mm256_store_pd(wr+c,_mm256_sub_pd(er,oii));
_mm256_store_pd(wi+c,_mm256_add_pd(ei,orr));
}
}
#else
FT static inline void reformat18_row_scalar(int i);
FT static inline void reformat18_row(int i) { reformat18_row_scalar(i); }
#endif
FT static inline void reformat18_row_scalar(int i) {
{
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
const u32 *CJ = par ? CJ1 : CJ0;
f64 *wr = RE+(size_t)r*STRIDE, *wi = IM+(size_t)r*STRIDE;
u32 k1 = BREV[i];
u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
/* rows 0,1 were clobbered by the pointwise before their own pass, so they come from
the copies taken at entry; select the base pointers once instead of testing i<2
in the inner loop (the two branches have identical layout). */
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
for (int j=0;j<NC;j+=2) {
f64 z1r,z1i,z2r,z2i;
{ z1r=pr[j]; z1i=pi[j]; z2r=pr[j+1]; z2i=pi[j+1]; }
u32 eh = kq + 2u*BREVC[j];
f64 qr0=DQr[eh], qi0=DQi[eh];
f64 tr = qr0*dr - qi0*di;
f64 ti = -(qr0*di + qi0*dr);
f64 er=z1r+z2r, ei=z1i+z2i;
f64 dr=z1r-z2r, di=z1i-z2i;
f64 orr=tr*dr-ti*di, oii=tr*di+ti*dr;
int ccol = (int)CJ[j];
wr[ccol] = er - oii;
wi[ccol] = ei + orr;
}
}
}
FT static void reformat18(void) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) reformat18_row(i);
}
FT static int parse_digits2(const char *s, long len, f64 *dst, int stride) {
// v43 -> c2: the destination pointer is now advanced incrementally (wp += 8, plus
// a row jump when the column wraps). The old form recomputed
// (base>>9)*stride + (base&511) from scratch every 32 digits, which cost ~10 of
// the loop's 22 instructions (shift, imul by stride, and, two lea) for what is
// really a constant +8 step. Identical stores, identical order.
int m = (int)((len + 3) / 4);
long lo = len & 3;
if (lo) { u32 v=0; for(long j=0;j<lo;j++) v=v*10+(u32)(s[j]-'0'); long kk=m-1;
dst[(size_t)(kk>>9)*stride + (kk & 511)] = (f64)v; }
const __m256i c10 = _mm256_setr_epi8(10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1,
10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1);
const __m256i c100 = _mm256_setr_epi16(100,1,100,1,100,1,100,1,100,1,100,1,100,1,100,1);
const __m256i adj = _mm256_set1_epi32(53328);
const __m256i rv = _mm256_setr_epi32(7,6,5,4,3,2,1,0);
long i = len;
f64 *wp = dst;
int col = 0;
// two chunks (64 digits) per iteration: halves the loop overhead and lets the two
// (load, maddubs, madd, permute, cvt, cvt, store, store) chains overlap.
while (i - 64 >= lo) {
__m256i d0 = _mm256_loadu_si256((const __m256i*)(s+i-32));
__m256i d1 = _mm256_loadu_si256((const __m256i*)(s+i-64));
i -= 64;
// the scan runs BACKWARDS, which the hardware prefetchers largely ignore; pull
// the next few lines in by hand
{ long pa = i-320; if (pa < 0) pa = 0; _mm_prefetch((const char*)(s+pa), _MM_HINT_T0);
long pb = i-448; if (pb < 0) pb = 0; _mm_prefetch((const char*)(s+pb), _MM_HINT_T0); }
__m256i v0 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d0,c10),c100), adj);
__m256i v1 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d1,c10),c100), adj);
v0 = _mm256_permutevar8x32_epi32(v0, rv);
v1 = _mm256_permutevar8x32_epi32(v1, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v0));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v0,1));
__m256d q2=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v1));
__m256d q3=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v1,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_stream_pd(wp,q2); _mm256_stream_pd(wp+4,q3);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_store_pd(wp,q2); _mm256_store_pd(wp+4,q3);
#endif
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
}
while (i - 32 >= lo) {
i -= 32;
__m256i v = _mm256_sub_epi32(_mm256_madd_epi16(
_mm256_maddubs_epi16(_mm256_loadu_si256((const __m256i*)(s+i)),c10),c100), adj);
v = _mm256_permutevar8x32_epi32(v, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
#endif
wp += 8; col += 8;
if (col == 512) { col = 0; wp += stride - 512; }
}
#if NTPARSE>0
_mm_sfence();
#endif
int k = (int)((i - lo) / 4);
while (i - 4 >= lo) {
i -= 4;
long kk = --k;
dst[(size_t)(kk>>9)*stride + (kk & 511)] =
(f64)((((s[i]-'0')*10 + (s[i+1]-'0'))*10 + (s[i+2]-'0'))*10 + (s[i+3]-'0'));
}
return m;
}
/* Vectorised boundary scan. The scalar original ran at ~1.09 cycles per input
byte over the whole 2 MB input: measured judge-side at 2.18 M cycles = 8.5 % of
the run, and it was absent from every earlier phase map. AVX2 here: 32 bytes
per iteration, exact same semantics (including "no digit at all" -> 0 return).
ALIGNED loads only -- gcc-9 splits every unaligned 256-bit access into
vmovdqu xmm + vinserti128 (verified in the .s), which is what the note in
JUDGE_MEM warns about. A <=31-byte scalar prologue gets us to alignment. */
TGT static inline long kb_scan(const char *s, long i, long sn, int want_digit) {
const __m256i lo = _mm256_set1_epi8('0'-1), hi = _mm256_set1_epi8('9'+1);
if (i < sn) {
while ((((uintptr_t)(s + i)) & 31) != 0) {
if (i >= sn) break;
if (((s[i] >= '0' && s[i] <= '9') ? 1 : 0) == want_digit) return i;
i++;
}
long n = sn - 32;
if (want_digit) {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (m) return i + __builtin_ctz(m);
}
while (i < sn && (s[i] < '0' || s[i] > '9')) i++;
} else {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (~m) return i + __builtin_ctz(~m);
}
while (i < sn && s[i] >= '0' && s[i] <= '9') i++;
}
}
return i;
}
/* FAST PATH for the split. The statement pins BOTH factors at exactly 10^6
digits, so the two digit runs are known to be [0,10^6) and [10^6+1, 2*10^6+1)
and the 2 MB digit/non-digit sweep of split_lines (measured judge-side at
~230 k cycles warm, ~480 k cycles with the input cold in DRAM - i.e. ~1..2 % of
the whole run) can be replaced by six byte tests. The sweep reads the input
once and parse_digits2 reads it again; this path reads it once. Any layout
that does not match (different digit counts, CRLF, leading blank, trailing
garbage) fails a test and falls through to the general sweep, so the fast path
can only ever be a pure win in time, never a correctness risk. */
TGT static inline int split_lines_fast(const char *s, long sn, const char **b0, long *n0,
const char **b1, long *n1) {
if (sn < 2000001) return 0;
unsigned c0 = (unsigned)(s[0] - '0'), c1 = (unsigned)(s[999999] - '0');
unsigned c2 = (unsigned)(s[1000001] - '0'),c3 = (unsigned)(s[1999999] - '0');
if (c0 > 9 || c1 > 9 || c2 > 9 || c3 > 9) return 0;
if (s[1000000] != '\n') return 0;
if (sn > 2000001 && s[2000001] != '\n') return 0;
*b0 = s; *n0 = 1000000; *b1 = s + 1000001; *n1 = 1000000;
return 1;
}
TGT static inline int split_lines(const char *s, long sn, const char **b0, long *n0,
const char **b1, long *n1) {
long i = kb_scan(s, 0, sn, 1);
long j = kb_scan(s, i, sn, 0);
if (j == i) return 0;
*b0 = s + i; *n0 = j - i;
long p = kb_scan(s, j, sn, 1);
long k = kb_scan(s, p, sn, 0);
if (k == p) return 0;
*b1 = s + p; *n1 = k - p;
return 1;
}
#define COEF(k) RE[(size_t)((k)>>9)*STRIDE + ((k) & 511)]
FT static void build_tab4(void) { }
FT static long carry_format_m(int tot, char *o, char *end) {
// ---- v43 -> c1 changes (all provably exact integer transformations) ----
// (1) double->u64 via cvtsd2si (round-to-nearest, default MXCSR). The old
// (u64)(x+0.5) is floor(x+0.5); the two agree whenever |x-n| < 0.5 for the
// true integer n, which is exactly the condition under which the old form
// was correct, so this is bit-identical. It also drops the addsd and the
// branchy >=2^63 fixup gcc emits for a double->u64 conversion.
// (2) hi/lo are taken from t4 = total/1e4 rather than from dg = total mod 1e8.
// total/1e4 does not depend on ex, so BOTH magic divisions now hang off the
// loop-carried chain instead of one behind the other. Digit group is
// unchanged: (total mod 1e8) = ((total/1e4) mod 1e4)*1e4 + total mod 1e4.
// (3) two 32-bit ASCII stores instead of salq/orq + one 64-bit store.
// (4) the tail's 429497 (ceil) magic had a fixup that only handled the
// under-estimate; the over-estimate case (dg % 1e4 >= 9938) indexed TAB4
// with a wrapped u32 -> out-of-bounds read. The tail is now plain / and %.
const int nm = (tot + 1) >> 1;
int nfull = ((tot - 1) + 1) >> 1;
if (nfull > nm) nfull = nm;
if (nfull < 0) nfull = 0;
u64 carry = 0;
char *p = end;
const f64 *pr = RE, *pi = IM;
int m = 0;
while (m < nfull) {
int cnt = nfull - m; if (cnt > NC) cnt = NC;
const f64 *qr = pr, *qi = pi;
for (int k = 0; k < cnt; k++) {
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qr + k)))
+ (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qi + k))) * 10000ull;
u64 total = C + carry;
u64 t4 = total / 10000ull; // independent of the carry chain
u64 ex = total / 100000000ull; // the chain's next carry
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
}
pr += STRIDE; pi += STRIDE; m += cnt;
}
pr = RE + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
pi = IM + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
for (; m < nm; m++) { // tail: no IM term
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd(pr));
u64 total = C + carry;
u64 t4 = total / 10000ull;
u64 ex = total / 100000000ull;
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
pr++; pi++;
}
while (carry) { u32 d = (u32)(carry % 100000000ull); carry /= 100000000ull;
u32 hi = d / 10000u, lo = d % 10000u;
p -= 8; *(u32*)p = TAB4[hi]; *(u32*)(p+4) = TAB4[lo]; }
while (p < end && *p == '0') p++;
if (p == end) *(--p) = '0';
(void)o;
return (long)(end - p);
}
FT static inline void pointwise_row(int i);
static unsigned char FWD_DONE[NR];
/* Fused forward: the diagonal twiddle and the pointwise product are both per-row
operations, so they run while their rows are still L1/L2-resident straight out of
row_dif instead of as two extra full sweeps of the 8.5 MB array. Rows are visited in
pairs (i, RVIX[i]) -- RVIX is an involution -- so both operands of every pointwise
are transformed before it runs, and the pair is always entered at its smaller index,
which keeps the pointwise application order (and hence the rounding) identical to the
old separate loop. */
#if FUSE16>0
FT static void fft_forward2(int sign, f64 sc);
#endif
FT static void fft_forward(void) {
#ifdef REVCD
for (int j0=NC-CBW;j0>=0;j0-=CBW) col_dif(j0);
#else
for (int j0=0;j0<NC;j0+=CBW) col_dif(j0);
#endif
memset(FWD_DONE, 0, sizeof FWD_DONE);
for (int i=0;i<NR;i++) {
if (FWD_DONE[i]) continue;
int m = (int)RVIX[i];
#if ROWPF>0
{ int n1=i+1; if(n1<NR && !FWD_DONE[n1]){ const char*q0=(const char*)(RE+(size_t)n1*STRIDE); const char*w0=(const char*)(IM+(size_t)n1*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(w0+o,_MM_HINT_T0); } }
int n2=(int)RVIX[i<NR-1?i+1:i]; if(n2!=i+1 && !FWD_DONE[n2]){ const char*q1=(const char*)(RE+(size_t)n2*STRIDE); const char*w1=(const char*)(IM+(size_t)n2*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q1+o,_MM_HINT_T0); _mm_prefetch(w1+o,_MM_HINT_T0); } } }
#endif
diagonal_row(i, +1, 1.0);
row_dif(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
if (m != i) {
diagonal_row(m, +1, 1.0);
row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE);
FWD_DONE[m] = 1;
}
FWD_DONE[i] = 1;
pointwise_row(i);
if (m != i) pointwise_row(m);
}
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
}
FT static void fft_inverse(void) {
for (int i=0;i<NR;i++) row_dit(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
diagonal(-1, 1.0/(f64)NN);
for (int j0=0;j0<NC;j0+=128) col_dit(j0);
}
/* Fused inverse. reformat18 finalises destination row r at source index i=2r+1 and
nothing later in the sweep touches row r again, so row_dit and the diagonal18 twiddle
run on row r immediately, while it is still cache-resident, instead of as two extra
sweeps of the 4 MB half-array. */
#if FUSE16>0
/* Fused: reformat18 for source rows 2r,2r+1 finalises destination row r and nothing
later in the sweep touches row r again, so row_dit + diagonal18 run on it immediately
while it is still L1/L2-resident -- and because the fold writes row r only after both
of its source rows were pointwise-complete (and every later pair reads only rows >= i),
this is exactly as safe as the separate sweep. */
FT static inline __attribute__((always_inline)) void fused_row_loop(int sign, f64 sc) {
for (int i=0;i<NR;i++) {
if (!FWD_DONE[i]) {
int m = (int)RVIX[i];
diagonal_row(i, +1, 1.0);
row_dif(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
if (m != i) { diagonal_row(m, +1, 1.0); row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE); }
FWD_DONE[i] = 1; FWD_DONE[m] = 1;
pointwise_row(i);
if (m != i) pointwise_row(m);
}
if (i & 1) {
int r = i>>1;
if (i == 1) { f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
reformat18_row(i-1);
reformat18_row(i);
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc);
}
}
}
#endif
FT static void fft_inverse18(int sign, f64 sc) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) {
reformat18_row(i);
if (i & 1) { int r = i>>1;
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc); }
}
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
}
#ifndef PWVEC
#define PWVEC 1
#endif
/* Vectorised pointwise. With U = row i at the 8 positions [8a,8a+8) and
S/V/T = row i / row ri at the partner block [8b,8b+8), b=63-a, the pair map is
exactly p -> 511-p, so block a partners block 63-a reversed. Every position's
input is read before it is written (even lanes of row i and odd lanes of row ri
are written by this call, the other half by the partner call, which runs after),
so the read-modify-write merge of the untouched lanes below is safe. Each output
is computed by the same expression from the same two inputs as the scalar form,
hence bit-identical (verified against it). */
#if PWVEC>0
#define BLD(X,Y,M) _mm256_blend_pd(X,Y,M)
#define PWC(z1r,z1i,z2r,z2i,cr,ci) do{ \
__m256d a1r=_mm256_add_pd(z1r,z2r), a1i=_mm256_sub_pd(z1i,z2i); \
__m256d b1r=_mm256_add_pd(z1i,z2i), b1i=_mm256_sub_pd(z2r,z1r); \
/* gcc-9 contracts the scalar form as cr = fma(a1r,b1r, rnd(a1i*b1i)) and
ci = fma(a1r,b1i, rnd(a1i*b1r)) -- i.e. the OTHER product is the exact one in
each. Mirroring that exactly is what makes these two bit-identical. */ \
__m256d u=_mm256_mul_pd(a1i,b1i); cr=_mm256_fmsub_pd(a1r,b1r,u); \
__m256d v=_mm256_mul_pd(a1i,b1r); ci=_mm256_fmadd_pd(a1r,b1i,v); }while(0)
FT static inline __attribute__((always_inline)) void pointwise_vec(int i) {
/* Row i's FIRST half (positions 0..255) pairs with row ri's SECOND half:
position x pairs with 511-x. Each pair is computed exactly once here and
both of its outputs are written with full-width 256-bit stores, so there is
no read-modify-write lane merge and the partner call pointwise_row(ri)
writes only the complementary halves (row ri's first half, row i's second
half) -- the two calls touch disjoint memory and their reads are never
affected by each other's stores. Per-element arithmetic is the scalar
form's, so this is bit-identical. */
unsigned ri = RVIX[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
const __m256d sg=_mm256_set1_pd(-0.0);
for(int A=0;A<256;A+=8){
int B=504-A; /* 511-(A+l) == B+7-l */
__m256d Ur0=_mm256_load_pd(pr+A), Ur1=_mm256_load_pd(pr+A+4);
__m256d Ui0=_mm256_load_pd(pi+A), Ui1=_mm256_load_pd(pi+A+4);
__m256d Vr0=_mm256_load_pd(qr+B), Vr1=_mm256_load_pd(qr+B+4);
__m256d Vi0=_mm256_load_pd(qi+B), Vi1=_mm256_load_pd(qi+B+4);
__m256d z2r0=_mm256_permute4x64_pd(Vr1,0x1B),z2r1=_mm256_permute4x64_pd(Vr0,0x1B);
__m256d z2i0=_mm256_permute4x64_pd(Vi1,0x1B),z2i1=_mm256_permute4x64_pd(Vi0,0x1B);
__m256d cr0,ci0,cr1,ci1;
PWC(Ur0,Ui0,z2r0,z2i0,cr0,ci0);
PWC(Ur1,Ui1,z2r1,z2i1,cr1,ci1);
_mm256_store_pd(pr+A, cr0); _mm256_store_pd(pr+A+4,cr1);
_mm256_store_pd(pi+A, ci0); _mm256_store_pd(pi+A+4,ci1);
_mm256_store_pd(qr+B, _mm256_permute4x64_pd(cr1,0x1B));
_mm256_store_pd(qr+B+4, _mm256_permute4x64_pd(cr0,0x1B));
_mm256_store_pd(qi+B, _mm256_xor_pd(_mm256_permute4x64_pd(ci1,0x1B),sg));
_mm256_store_pd(qi+B+4, _mm256_xor_pd(_mm256_permute4x64_pd(ci0,0x1B),sg));
}
}
#endif
FT static inline void pointwise_row_scalar(int i);
FT static inline void pointwise_row(int i) {
#if PWVEC>0
/* vec() writes row i and row ri in the same pass; if they are the same row the
two read-modify-write merges clobber one another. RVIX[i]==i only for i=0,1. */
unsigned ri = RVIX[i];
if (BREV[i] && ri != (unsigned)i) { pointwise_vec(i); return; }
#endif
pointwise_row_scalar(i);
}
FT static inline void pointwise_row_scalar(int i) {
{
unsigned ri = RVIX[i];
const u32 *CX = BREV[i] ? CVIX1 : CVIX0;
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
for (int j=0;j<NC;j+=2) {
int jr = (int)CX[j];
f64 z1r=pr[j], z1i=pi[j], z2r=qr[jr], z2i=qi[jr];
f64 a1r=(z1r+z2r), a1i=(z1i-z2i);
f64 b1r=(z1i+z2i), b1i=-(z1r-z2r);
f64 cr=a1r*b1r - a1i*b1i, ci=a1r*b1i + a1i*b1r;
pr[j]=cr; pi[j]=ci;
qr[jr]=cr; qi[jr]=-ci;
}
}
}
FT static void pointwise(void) {
for (int i=0;i<NR;i++) pointwise_row(i);
/* self-conjugate bin: raw at this index, and the loop now produces 4*A*B, so scale by 4 */
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
}
FT static void run_job(DI *d) {
const char *s; long sn;
s = d->s; sn = (long)d->sn;
const char *b0,*b1; long n0,n1;
if (!split_lines_fast(s, sn, &b0, &n0, &b1, &n1) &&
!split_lines(s, sn, &b0, &n0, &b1, &n1)) return;
build_all();
const int m0 = (int)((n0 + 3) / 4), m1 = (int)((n1 + 3) / 4);
/* The pre-fault pass that used to live here was PURE WASTE (measured judge-side,
2026-09-25): a READ fault on an anonymous page maps the shared zero page
read-only, and the first WRITE to it then takes a SECOND fault (COW). So the
prefault paid one 846-cycle fault per page and the kernels still paid the COW
fault afterwards. Deleting it: prefault phase 1 022 960 -> 1 612 cycles,
col_dif unchanged (+23 k), judge COLD total 23 659 776 -> 22 687 674, output
byte-identical. Let the kernels take the single write fault instead. */
parse_digits2(b0, n0, RE, STRIDE);
parse_digits2(b1, n1, IM, STRIDE);
int tot = m0 + m1 - 1;
#if FUSE16>0
fft_forward();
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
fused_row_loop(-1, 1.0/(4.0*(f64)NN));
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
#else
fft_forward();
fft_inverse18(-1, 1.0/(4.0*(f64)NN));
#endif
/* Write the digit groups backwards from a FIXED offset (the answer is at most
2*10^6 digits, so d->o[0..2000000] is exactly the region the final memmove
needs) instead of from the far end of the judge's buffer: with a generously
sized d->ol the old form dirtied a second, disjoint 2 MB region. */
/* Cap the backward write at the region the final memmove needs (the product of two
10^6-digit factors has at most 2*10^6 digits), instead of the far end of the
judge's output buffer: with a generously sized d->ol the old form dirtied a
second, disjoint 2 MB / 489-page region (RSS 12460 -> 10460 KB, mem_kb 12416 ->
10460). Clamped so a short d->ol degrades to the old bound instead of running
off the front. */
long olim = (long)d->ol - 1; if (olim > 2000000) olim = 2000000;
char *cap = (char *)d->o + olim;
build_tab4();
long L = carry_format_m(tot, cap, cap);
/* with cap at the fixed offset the source of the shift is already d->o whenever
the product has the full 2*10^6 digits, so skip the 2 MB self-copy */
if (cap - L != (char *)d->o) memmove(d->o, cap - L, (size_t)L);
((char *)d->o)[L] = '\n';
d->os = L + 1;
}
#ifndef NO_LOCAL_TEST
extern "C" void __libc_start_main(void *m, int argc, char **argv) {
(void)m;
unsigned long *p = (unsigned long *)(argv + argc + 1);
while (*p) p++;
p++;
DI *d = 0;
for (int i = 0; i < 32 && p[0]; i++, p += 2)
if (p[0] == 0x6b637564UL) { d = (DI *)p[1]; break; }
if (d) run_job(d);
__asm__ volatile("syscall" ::"a"(60), "D"(0) : "rcx", "r11", "memory");
for (;;);
}
int main(){ return 0; }
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 5.715 ms | 10 MB + 192 KB | Accepted | Score: 100 | 显示更多 |