「快速数列变换这个词怎么造出来的」的一些阅读感悟——2025.7.3 鲜花
「快速数列变换这个词怎么造出来的」的一些阅读感悟
冠世一战
虽千万人逆我之 我仍执着
侠客行遍九州 恩仇奔波
刀光剑影里是谁的明眸粲粲如星火
纵使我双掌倚天 苍龙伏获
奈何降不住这人心如魔
所谓天下第一 怎囿我 天高海阔
飞雪连天 扬沙满地华山巅 谁论剑
盖世绝学 宏图业 究竟能有多执念
狭路相逢 离离分分 游丝难织就双飞燕
天之道 损有余 所求仍不得圆
斗胆教日月射落
白针将剑招拆破
谁逐鹿 谁流离 谁问鼎 谁失所
此战 何必笑我
虽千万人逆我之 我仍执着
侠客行遍九州 恩仇奔波
刀光剑影里是谁的明眸粲粲如星火
纵使我双掌倚天 苍龙伏获
奈何降不住这人心如魔
所谓天下第一 怎囿我 天高海阔
左右逢源 大悟大彻黑白颠 谁周旋
书剑快意 仁无敌越是冠绝 越决绝
痴缠眷侣 生死相许 渺万里层云 千山雪
侠之道 为国为民 敢向天下先
内力将乾坤腾挪
轻功点水上烟波
阴与阳 柔与刚 强与弱 对与错
此战 我战与我
虽千万人逆我之 我仍执着
侠客行遍九州 恩仇奔波
刀光剑影里是谁的明眸粲粲如星火
纵使我双掌倚天 苍龙伏获
奈何降不住这人心如魔
所谓天下第一 怎囿我 天高海阔
所幸敝屣荣华且把酒
风雨偏洗岳阳楼
冠世侠者 先天下之忧
此战 再出手
虽千万人逆我之 我仍执着
寒铓再出江湖 付身家国
踏破贺兰山缺 君自有还我山河胆魄
纵使我双掌倚天 苍龙伏获
摒此生 尚不知 是功是过
若非冠世一战 后人何处听传说
我不会代数，很难的，我看不太懂，简单说一下我懂的部分，如果有不懂的可以看原文，写的很细的。
首先会一下朴素的 FWT，可以直接看 位运算卷积(FWT) & 集合幂级数，这里简单说一下：
就是考虑高维前缀和：
for(int i = 0; i < n; ++i)
for(int j = 0; j < (1 << n); ++j)
if(!(j >> i & 1)) A[j | 1<< i] += A[j];
我们扩展到 \(a \circ b\)，其中 \(\circ\) 是位运算。
可以发现过程中都是线性操作，我们用一个矩阵刻画内层操作，具体的，我们构造一个矩阵
可以证明，只要满足限制：\(c_{i, j}c_{i, k} = c_{i, j \circ k}\)，我们用这个矩阵做 FWT 即可满足限制。
容易发现这个约束是充分的，但其确实不是很必要，因为已经有不能用这个做的题了：
题意：求
其中 \(mex_3\) 表示转成三进制后每位做 \(mex\)。
如果用上面的方法，我们应该能立即给出一个 \(3^nn\) 的做法，但事实上是不太可能的。
这里我们以 \(h_k = \sum_{i \oplus j = k} f_ig_j\) 为例，给出一个分治乘法的做法：
这里我本来想写「这本来有个绝妙的例子，但笔者写完这是电脑死机没存，所以不再写了」的 QwQ
设其长为 \(2^n\)，我们假设我们已经会了 \(n - 1\) 的做法 \(W_{n - 1}(f, g)\)，尝试递归调用。
边界是显然的，对于 \(n\) 足够小暴力即可。复杂度跟每层调用的 \(W_{n - 1}\) 次数有关，设其为 \(r\)，则复杂度是 \(T(n) = rT(n - 1) + \mathcal{O}(2^n)\)
我们将 \(f, g, h\) 按照下标最后一位拆成 \(f_0, f_1, g_0, g_1, h_0, h_1\)，我们用形式幂级数给出一个例子：设 \(f = \{3, x, 5x^2, 7x^3\}\)，则 \(f_0 = \{3, 5x^2\}\)，\(f_1 = \{x, 7x^3\}\)。
我们显然有：
如果暴力做 \(r = 4\)，此时复杂度 \(T(n) = \mathcal{O}(4 ^ n)\) 和暴力无异。
但是我们发现：
于是我们只用调用 \(2\) 次，复杂度是 \(T(n) = \mathcal{O}(2^nn)\)。
我们用矩阵刻画一下刚才的变换，并将进制扩展到 \(v\) 进制：
考虑其限制，容易发现是：
熟悉张量的同学应该已经发现了，这个是三维张量 \(\{[i \circ j = k]\}\) 的 CP 分解，我们可以手玩甚至人类智慧甚至爆搜甚至贺一个 CP 分解板子并在实数范围里调参。（爆搜是基于 [集训队互测 2024] 欧伊昔 告诉我们在操作矩阵是 \(3 \times 3\) 时 \(r \le 4\) 且用 \(-1, 0, 1\) 就能搜出解，剪个枝就能跑）。
我们用 Tritwise Mex 来感受一下如何人类智慧：
我们将操作表打出来（\((i, j)\) 的值表示 \(mex_3(i, j)\)）：
发现 \(h_0 = W(f_1 + f_2, g_1 + g_2)\)，\(h_2 = W(f_1, g_0) + W(f_0, g_1)\)，\(h_1 = W(f_1 + f_2 + f_3, g_1 + g_2 + g_3) - h_0 - h_2\)，于是 \(r = 4\) 解决。
给出实现：
Tritwise Mex
/*
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra && time ./%< && size %<
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra -fsanitize=address,undefined -g && time ./%< && size %<
echo && cat out.out && echo
*/
#include <bits/stdc++.h>
using namespace std;
using llt = long long;
using llf = long double;
using ull = unsigned long long;
#define endl '\n'
#ifdef LOCAL
FILE *InFile = freopen("in.in", "r", stdin), *OutFile = freopen("out.out", "w", stdout);
#endif
#define int llt
const int LG = 12, N = pow(3, LG) + 3;
int _a[N], _b[N], _c[N], n;
template <const int V, const int R> class FWT{
private:
int pw[max(V, R) + 1][LG + 1];
public:
using Tag = vector<vector<int>>;
Tag ta, tb, tc;
FWT(){
for(auto i : {V, R}){
pw[i][0] = 1;
for(int j = 1; j <= LG; ++j) pw[i][j] = pw[i][j - 1] * i;
}
}
void Tim(int *f, int *g, int *h, int n){
if(!n) return h[0] = f[0] * g[0], void();
int l = pw[V][n], nl = pw[V][n - 1], a[nl] = {}, b[nl] = {}, c[nl] = {};
memset(h, 0, sizeof(*h) * l);
for(int k = 0; k < R; ++k){
for(int i = 0; i < nl; ++i){
a[i] = b[i] = 0;
for(int j = 0; j < V; ++j)
a[i] = a[i] + f[i * V + j] * ta[k][j],
b[i] = b[i] + g[i * V + j] * tb[k][j];
}
Tim(a, b, c, n - 1);
for(int i = 0; i < nl; ++i) for(int j = 0; j < V; ++j)
h[i * V + j] = h[i * V + j] + c[i] * tc[k][j];
}
}
};
FWT<3, 4> Fwt;
int ca[N], cb[N];
signed main(){
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n; int l = pow(3, n);
for(int i = 0; i < l; ++i) cin >> _a[i];
for(int i = 0; i < l; ++i) cin >> _b[i];
Fwt.ta = {{1, 1, 1}, {0, 1, 1}, {1, 0, 0}, {0, 1, 0}};
Fwt.tb = {{1, 1, 1}, {0, 1, 1}, {0, 1, 0}, {1, 0, 0}};
Fwt.tc = {{0, 1, 0}, {1, -1, 0}, {0, -1, 1}, {0, -1, 1}};
Fwt.Tim(_a, _b, _c, n);
for(int i = 0; i < l; ++i) cout << _c[i] << ' ';
cout << endl;
}
我们称这个东西几乎覆盖了所有快速数列变换的场景，所以我们让它支持子集卷积。
考虑子集卷积：
好消息是暴力做就是 \(3^n\) 的，坏消息是这个东西不太能直接用上面的办法做了。
我们发现其长得像一个真正的卷积，具体的，我们整个形式元 \(x\)，则有：
也就是 \(h_0 + h_1x + W(f_0 + f_1, g_0 + g_1)x^2\)。
我们考虑 \(h_0\) 可以直接做，只要将 \(h_1\) 设置成最低次项系数就可以完美保留，于是就简单了：
这里 \(h_1\) 是一个多项式，我们最后的答案是提取其最低项也就是 \(\mathrm{popcount}(n)\) 项系数。
一个小优化：我们发现其每次加法都是对应位的，我们初始时插入 \(\mathrm{popcount}(n)\) 项系数即可不用 \(\times x\)，也就省掉了一次平移。
若想卡到 \(\mathcal{O}(2^n)\) 空间其实可以直接拉插，原文有更简单的做法，我没仔细看。
给出实现：
子集卷积
/*
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra && time ./%< && size %<
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra -fsanitize=address,undefined -g && time ./%< && size %<
echo && cat out.out && echo
*/
#include <bits/stdc++.h>
using namespace std;
using llt = long long;
using llf = long double;
using ull = unsigned long long;
#define endl '\n'
#ifdef LOCAL
FILE *InFile = freopen("in.in", "r", stdin), *OutFile = freopen("out.out", "w", stdout);
#endif
const int LG = 20, N = 1 << LG | 3, M = LG + 3, MOD = 1e9 + 9;
int _a[N][M], _b[N][M], _c[N][M], n;
vector<vector<int>> ta, tb, tc;
void Tim(int (*f)[M], int (*g)[M], int (*h)[M], int n){
if(!n){
memset(h[0], 0, sizeof h[0]);
for(int i = 0; i <= ::n; ++i) for(int j = 0; i + j <= ::n; ++j)
h[0][i + j] = (h[0][i + j] + 1ll * f[0][i] * g[0][j]) % MOD;
return ;
}
int l = 1 << n, nl = l >> 1, a[nl][M], b[nl][M], c[nl][M];
memset(h, 0, sizeof(*h) * l);
for(int w = 0; w < 2; ++w){
for(int i = 0; i < nl; ++i){
memset(a[i], 0, sizeof a[i]), memset(b[i], 0, sizeof b[i]);
for(int j = 0; j < 2; ++j) for(int k = 0; k <= ::n; ++k)
a[i][k] = (a[i][k] + 1ll * f[i << 1 | j][k] * ta[w][j]) % MOD,
b[i][k] = (b[i][k] + 1ll * g[i << 1 | j][k] * tb[w][j]) % MOD;
}
Tim(a, b, c, n - 1);
for(int i = 0; i < nl; ++i) for(int j = 0; j < 2; ++j) for(int k = 0; k <= ::n; ++k)
h[i << 1 | j][k] = (h[i << 1 | j][k] + 1ll * c[i][k] * tc[w][j]) % MOD;
}
}
int main(){
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n; int l = 1 << n;
for(int i = 0; i < l; ++i) cin >> _a[i][__builtin_popcount(i)];
for(int i = 0; i < l; ++i) cin >> _b[i][__builtin_popcount(i)];
ta = tb = {{1, 1}, {1, 0}}, tc = {{0, 1}, {1, -1}};
Tim(_a, _b, _c, n);
for(int i = 0; i < l; ++i) cout << (_c[i][__builtin_popcount(i)] + MOD) % MOD << ' ';
}
但是说实话，可能是我实现的问题，这个玩意的常数实在是不敢恭维。
进一步！
我们发现其限制有着天然的优美，即其形式高度一致，但 \(a, b\) 依赖输入，\(c\) 依赖输出。
这也就意味着我们可以轻松实现差卷积，即：
具体的，我们只需要交换其 CP 分解的第 \(1, 3\) 维即可，也就是交换 \(a, c\)。
但是我们有时想要的是一种「变换-逆变换」的形式，我们也可以用 FWT 在同样的时间复杂度内实现。
先考虑 \(r = v\)，我们发现我们的矩阵简直就是位矩阵，所以我们直接做即可，唯一要注意的是做 IFWT 时要对 \(c\) 转置。
考虑 \(r > v\) 时，我们容易想到将 \(v\) 补齐到 \(r\)，但是这样复杂度会多 \(n\)。
我们考虑在做 FWT 时将 \(v\) 进制转化成 \(r\) 进制，相当于时一个扩张，在 IFWT 时将 \(r\) 进制转回来，相当于一个缩减。
容易将复杂度证到和原来一样。
给出实现：
【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
/*
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra && time ./%< && size %<
ulimit -s 128000 && clear && g++ % -o %< -O2 -std=c++14 -DLOCAL -Wall -Wextra -fsanitize=address,undefined -g && time ./%< && size %<
echo && cat out.out && echo
*/
#include <bits/stdc++.h>
using namespace std;
using llt = long long;
using llf = long double;
using ull = unsigned long long;
#define endl '\n'
#ifdef LOCAL
FILE *InFile = freopen("in.in", "r", stdin), *OutFile = freopen("out.out", "w", stdout);
#endif
const int LG = 17, N = 1 << LG | 3, MOD = 998244353, I2 = (MOD + 1) >> 1;
int _a[N], _b[N], _c[N], n, _u;
template <const int V, const int R> class FWT{
private:
int pw[max(V, R) + 1][LG + 1], g[N];
public:
using Tag = vector<vector<int>>;
Tag ta, tb, tc;
FWT(){
for(auto i : {V, R}){
pw[i][0] = 1;
for(int j = 1; j <= LG; ++j) pw[i][j] = pw[i][j - 1] * i;
}
}
void Fwt(int *f, int n, const Tag &a, int v = V, int r = R){
int l = pw[v][n];
for(int i = 0; i < n; ++i){
int nl = l / v * r, p = pw[r][i]; memset(g, 0, sizeof(*g) * nl);
for(int j = 0; j < l; ++j) if(j / p % v == 0){
int t[r] = {0};
for(int k = 0; k < r; ++k) for(int x = 0; x < v; ++x)
t[k] = (t[k] + 1ll * f[j + x * p] * a[k][x]) % MOD;
int nj = j / p / v * p * v + j % p;
for(int k = 0; k < r; ++k) g[nj + k * p] = t[k];
}
memcpy(f, g, sizeof(*g) * nl);
}
}
void Ifw(int *f, int n, const Tag &a){
Tag x(V, vector<int>(R));
for(int i = 0; i < V; ++i) for(int j = 0; j < R; ++j) x[i][j] = a[j][i];
Fwt(f, n, x, R, V);
}
};
FWT<2, 2> Fwt;
int ca[N], cb[N];
int main(){
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n; int l = 1 << n;
for(int i = 0; i < l; ++i) cin >> _a[i];
for(int i = 0; i < l; ++i) cin >> _b[i];
memcpy(ca, _a, sizeof ca), memcpy(cb, _b, sizeof cb);
Fwt.Fwt(ca, n, {{1, 1}, {1, 0}}), Fwt.Fwt(cb, n, {{1, 1}, {1, 0}});
for(int i = 0; i < l; ++i) _c[i] = 1ll * ca[i] * cb[i] % MOD;
Fwt.Ifw(_c, n, {{0, 1}, {1, -1}});
for(int i = 0; i < l; ++i) cout << (_c[i] + MOD) % MOD << ' ';
cout << endl;
memcpy(ca, _a, sizeof ca), memcpy(cb, _b, sizeof cb);
Fwt.Fwt(ca, n, {{1, 1}, {0, 1}}), Fwt.Fwt(cb, n, {{1, 1}, {0, 1}});
for(int i = 0; i < l; ++i) _c[i] = 1ll * ca[i] * cb[i] % MOD;
Fwt.Ifw(_c, n, {{1, 0}, {-1, 1}});
for(int i = 0; i < l; ++i) cout << (_c[i] + MOD) % MOD << ' ';
cout << endl;
memcpy(ca, _a, sizeof ca), memcpy(cb, _b, sizeof cb);
Fwt.Fwt(ca, n, {{1, 1}, {1, -1}}), Fwt.Fwt(cb, n, {{1, 1}, {1, -1}});
for(int i = 0; i < l; ++i) _c[i] = 1ll * ca[i] * cb[i] % MOD;
Fwt.Ifw(_c, n, {{I2, I2}, {I2, -I2}});
for(int i = 0; i < l; ++i) cout << (_c[i] + MOD) % MOD << ' ';
cout << endl;
}
较为可惜的是，我现在还不知道一种能够在一般情况下扩展 CF1119H Triple 的方法，目前只能在操作满足交换结合时将 \(v\) 进制转化成 \(r\) 进制，并用原来的 \(c_{i, j}c_{i, k} = c_{i, j \circ k}\) 限制做，做法就和正常的一样了，复杂度依然要乘 \(n\)，但是逆变换好像可以一样的缩减，所以逆变换的复杂度好像是对的。
P
本文来自博客园，作者：xrlong，转载请注明原文链接：https://www.cnblogs.com/xrlong/p/18962183
版权声明：本作品采用 「署名-非商业性使用-相同方式共享 4.0 国际」许可协议（CC BY-NC-SA 4.0） 进行许可。
浙公网安备 33010602011771号