快速沃尔什变换解决的卷积问题
快速沃尔什变换（FWT）是解决这样一类卷积问题：
其中，(odot) 是位运算的一种。举个例子，给定数列 (a,b)，求：
FWT 的思想
看到 FWT 的名字，我们可以联想到之前学过的 FFT（很可惜，我没有写过 FFT 的笔记，所以没有链接），先看看 FFT 的原理：
- 把 (a,b) 变换为 (A,B)，(O(nlog n))；
- 通过 (C_i=A_iB_i) 计算，(O(n));
- 把 (C) 变换回 (c)，(O(nlog n))。
综上，时间复杂度是 (O(nlog n)) 的。
在 FFT 中，我们构造了 (A,B) 为 (a,b) 的点值表示法，这么做满足 (C_i=A_iB_i) 且容易变换。
其实 FWT 的思想也是一样的，主要也是需要构造 (A,B)，使得其满足 (C_i=A_iB_i) 且可以快速变换。下面我们举 (cup)（按位或）、(cap)（按位与）和 (oplus)（按位异或）为例。
因为数列长度是 (2) 的幂会更好处理，所以下文认为数列长度为 (2^n)。
按位或
我们可以构造 (A_i=sum_{icup j=i} a_i)。看看为什么需要这么构造。
首先，它满足 (C_i=A_iB_i)：
其次，它可以快速变换。举顺变换的例子。类比 FFT 的步骤，我们采用分治的方法来处理它。假设目前考虑到第 (i) 位，其中 (A_0) 和 (A_1) 是 (i-1) 位分治的结果：
其中，(A_0) 是数列 (A) 的左半部分，(A_1) 是 (A) 的右半部分。(text{merge}) 函数就是把两个数列像拼接字符串一样拼接起来。(+) 则是将两个数列对应相加。
这么做为什么是正确的呢？容易发现，(A_0) 恰好是当前处理到的二进制位为 (0) 的子数列，(A_1) 则是当前处理到的二进制位为 (1) 的子数列。若当前位为 (0)，则只能取二进制位为 (0) 的子数列 (A_0) 才能使得 (icup j=i)。而若当前位为 (1)，则两种序列都能取。
考虑逆变换，则是将加上的 (A_0) 减回去：
下面我们给出代码实现。容易发现顺变换和逆变换可以合并为一个函数，顺变换时 (text{type}=1)，逆变换时 (text{type}=-1)。
void Or(ll *a, ll type) { // 迭代实现，常数更小 for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j + k] += a[i + j] * type) %= P; } } } } 按位与
同理构造 (A_i=sum_{icap j=i} a_i)。(C_i=A_iB_i) 的正确性不证了。
容易发现，(A_0) 恰好是当前处理到的二进制位为 (0) 的子数列，(A_1) 则是当前处理到的二进制位为 (1) 的子数列。若当前位为 (1)，则只能取二进制位为 (1) 的子数列 (A_0) 才能使得 (icap j=i)。而若当前位为 (0)，则两种序列都能取。
下面我们给出代码实现。顺变换时 (text{type}=1)，逆变换时 (text{type}=-1)。
void And(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j] += a[i + j + k] * type) %= P; } } } } 按位异或
发现异或有点难搞，但这怎么会难倒沃尔什大佬呢？我们引入一个新的运算符 (circ)。定义 (xcirc y=text{popcnt}(xcap y)bmod 2)，其中 (text{popcnt}) 表示二进制下 (1) 的个数，并重申一下 (cap) 表示按位与。
不用慌，我们也不需要你真正实现一个 (text{popcnt})，它仅仅只是作为一个理解的辅助罢了。
我们发现它满足 ((xcirc y)oplus (xcirc z)=xcirc(yoplus z))。（重申一下 (oplus) 表示按位异或）
感性证明：发现这个新的运算符 (circ) 其实就是 (x) 与 (y) 相同位数的奇偶性。若 ((xcirc y)oplus (xcirc z)=0)，则 (x) 与 (y)、(x) 与 (z) 相同位数个数奇偶性相同，所以 (yoplus z) 和 (x) 相同位数个数奇偶性也是相同的 ；若 ((xcirc y)oplus (xcirc z)=1)，则 (x) 与 (y)、(x) 与 (z) 相同位数个数奇偶性不同，所以 (yoplus z) 和 (x) 相同位数个数奇偶性也是不同的。
设 (A_i=sum_{icirc j=0}a_j-sum_{icirc j=1}a_j)。我们来证一下 (C_i=A_iB_i) 的正确性：
来看看怎么快速计算 (A,B) 的值，依旧是分治：
对于 (i) 在当前位为 (0) 的子数列 (A_0)，进行 (circ) 运算时发现它和 (0) 计算或和 (1) 计算结果都不会变（因为 (0cap 0=0,0cap1=0)），所以 (A_i=sum_{icirc j=0}a_j-sum_{icirc j=1}a_j) 中的 (sum_{icirc j=1}a_j=0)。
对于 (i) 在当前位为 (1) 的子数列 (A_1)，进行 (circ) 运算时发现它和 (0) 计算结果是 (0)，和 (1) 计算结果是 (1)（因为 (1cap 0=0,1cap1=1)）。
综上，有：
也就是：
逆变换易得：
给出代码，顺变换时 (text{type}=1)，逆变换时 (text{type}=frac{1}{2})。
void Xor(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j] += a[i + j + k]) %= P; (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P; (a[i + j] *= type) %= P; (a[i + j + k] *= type) %= P; } } } } 现在大家能去切前两道模板例题，并挑战一下后面的几道题目了。
从另一个角度看待 FWT
我们设 (c(i,j)) 是 (a_j) 对 (A_i) 的贡献系数。我们可以重新描述 FWT 变换的过程：
因为有：
所以我们可以通过简单的证明得到：(c(i,j)c(i,k)=c(i,jodot k))。其中 (odot) 是任意一种位运算。
同时，(c) 函数还有一个重要的性质，它可以按位处理。
举个例子，我们变换的时候：
这么做是比较劣的，我们将其拆分：
考虑前面的式子和后面的式子 (i,j) 的区别，发现只有最高位不同。
所以我们将 (i,j) 去除最高位的值为 (i',j')，并记 (i_0) 为 (i) 的最高位。有：
如果 (i_0=0)，则有：
(i_0=1) 则有：
也就是说，我们只需要：
四个数就可以完成变换了。我们称这个矩阵为位矩阵。
如果我们要进行逆变换，则需要上面的位矩阵的逆矩阵。
若逆矩阵为 (c^{-1})，可以通过类似操作得到原数：
逆矩阵不一定存在，比如如果有一排 (0) 或者一列 (0) 那么这个矩阵就没有逆，我们在构造时需要格外小心。
按位或
我们可以构造：
这样满足 (c(i,j)c(i,k)=c(i,jcup k))。我们发现，这和我们前面推出的 (A=text{merge}(A_0, A_0+A_1)) 一模一样！同理，下面也是一个满足这个条件的矩阵，但我们一般使用上面这个：
虽然下面这个矩阵也满足 (c(i,j)c(i,k)=c(i,jcup k))，但这个矩阵存在一排 (0)，不存在逆，所以不合法：
如果我们要进行逆变换，则需要对矩阵求逆，以最上面这个矩阵为例，得：
然后按照顺变换的方法，把逆变换矩阵代入即可。
按位与
我们可以构造：
这样满足 (c(i,j)c(i,k)=c(i,jcap k))。
逆矩阵：
按位异或
我们可以构造：
这样满足 (c(i,j)c(i,k)=c(i,joplus k))。
逆矩阵：
FWT 的性质
FWT 是线性变换。
若 (FWT(X)) 是 (X) 的 FWT 变换，则有：
以及：
这样就可以实现快速卷积，参考第四道例题。
K 维 FWT
max 运算
我们重新看看我们的 (cup) 运算，发现他实际上就是二进制下的取 (max)。我们将其拓展到 (K) 进制，有：
若 (j=k)，那么上式又是：
也就是说，每一行的 (1) 必定只能在 (0) 的前面，如果在后面则不合法了。手玩一下可以发现一组合法构造：
求逆可得：
min 运算
我们重新看看我们的 (cap) 运算，发现他实际上就是二进制下的取 (min)。我们将其拓展到 (K) 进制，有：
若 (j=k)，那么上式又是：
也就是说，每一行的 (1) 必定只能在 (0) 的后面，如果在后面则不合法了。手玩一下可以发现一组合法构造：
求逆可得：
前两者用得比较少，用得比较多的是：
不进位加法
我们重新看看我们的 (oplus) 运算，发现他实际上就是二进制下的不进位加法。我们将其拓展到 (K) 进制，有：
我们构造 (c(i,j)=omega_{K}^j)，就可以满足要求了：
但是每一行都一样矩阵也没有逆，所以我们可以构造 (c(i,j)=omega_{K}^{(i-1)j}) 即可。
有下面这个矩阵：
这不就是我们熟悉的范德蒙德矩阵吗？
现在我们也知道矩阵的逆了：
如果我们题目给出的模数是存在单位根的，我们就可以简单实现，可以参考第六道例题。
但是单位根在模意义下可能不存在，所以我们考虑扩域，就是人为地定义一个 (x)，满足 (x^K=1)，然后直接把 (x) 代入计算，这样每个数都是一个关于 (x) 的 (k-1) 次多项式。我们只需要在 (pmod {x^K-1}) 下计算即可。那么矩阵可以这么表示：
但是这么做可能会存在零因子，也就是一个数有多种表示方法，我们无法确定一个数的真实值。
我们考虑不 (pmod {x^K-1}) 了，我们 (bmod) 分圆多项式 (Phi_{K}(x))，他满足 (x) 的阶为 (k)，且在 (Q) 上不可约。所以我们定义上面的计算是在 (pmod {Phi_{K}(x)}) 下进行即可。
另一方面，如何求分圆多项式，这一点可以在因式分解这道题的题解区里了解。这里给出分圆多项式的表：
还有一个问题是，(bmod Phi_{K}(x)) 常数大（因为 (Phi) 本身就是一个多项式）。但是因为 (Phi_{K}(x)mid x^k-1)，我们只需要在计算时 (bmod x^k -1)，最后再 (bmod Phi_{K}(x)) 即可。
具体实现参考第七道例题。
例题
「洛谷 P4717」 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
求 (cup)、(cap)、(oplus) 的三种卷积。
(nle17)
这题也就是模板题了，下文直接给出代码：
#include <bits/stdc++.h> using namespace std; #define ll long long #define P 998244353 const ll N = 1 << 18; ll n; ll A[N], B[N]; ll a[N], b[N]; void init() { for(ll i = 0; i < n; i ++) a[i] = A[i], b[i] = B[i]; } void Or(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j + k] += a[i + j] * type) %= P; } } } } void And(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j] += a[i + j + k] * type) %= P; } } } } void Xor(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j] += a[i + j + k]) %= P; (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P; (a[i + j] *= type) %= P; (a[i + j + k] *= type) %= P; } } } } void calc() { for(ll i = 0; i < n; i ++) (a[i] *= b[i]) %= P; } void print() { for(ll i = 0; i < n; i ++) printf("%lld ", (a[i] % P + P) % P); printf("n"); } int main() { scanf("%lld", &n); n = 1 << n; for(ll i = 0; i < n; i ++) scanf("%lld", &A[i]); for(ll i = 0; i < n; i ++) scanf("%lld", &B[i]); init(); Or(a, 1); Or(b, 1); calc(); Or(a, P - 1); print(); init(); And(a, 1); And(b, 1); calc(); And(a, P - 1); print(); init(); Xor(a, 1); Xor(b, 1); calc(); Xor(a, 499122177); print(); } 「洛谷 P6097」 【模板】子集卷积
求：
[c_k=sum_{substack{{i cap j=0}\{icup j=k}}} a_i b_j ](nle20)
首先，下半部分是我们喜闻乐见的 FWT 常见形式，而上半部分我们可以看成是 (i) 与 (j) 不交。有：
所以我们可以构造：
可以枚举 (text{popcnt}) 的值，分开考虑。
那么求 (C) 的时候有 (C_{i,j}=sum_{j=0}^n A_{i,k}B_{i,j-k})。
然后就可以做了。
#include <bits/stdc++.h> using namespace std; #define popcnt(x) __builtin_popcountll(x) #define ll long long const ll M = 20, N = 1 << M, P = 1e9 + 9; ll n, m; ll a[M + 1][N], b[M + 1][N], c[M + 1][N]; void Or(ll *a, ll type) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j + k] += a[i + j] * type) %= P; } } } } int main() { scanf("%lld" ,&m); n = 1 << m; for(ll i = 0; i < n; i ++) { scanf("%lld", &a[popcnt(i)][i]); } for(ll i = 0; i < n; i ++) { scanf("%lld", &b[popcnt(i)][i]); } for(ll i = 0; i <= m; i ++) { Or(a[i], 1); Or(b[i], 1); } for(ll i = 0; i <= m; i ++) { for(ll j = 0; j <= i; j ++) { for(ll k = 0; k < n; k ++) { (c[i][k] += a[j][k] * b[i - j][k]) %= P; } } } for(ll i = 0; i <= m; i ++) { Or(c[i], -1); } for(ll i = 0; i < n; i ++) { printf("%lld ", (c[popcnt(i)][i] % P + P) % P); } } 「牛客 881D」Parity of Tuples
给定 (ntimes m) 的矩阵 (a)，定义 (text{cnt}(x)) 为矩阵中有多少行对于 (x) 是合法的，合法的定义为这一行中每一个数 (a_{i,j}cap x) 的二进制值中都有奇数个 (1)。
你需要求出对于所有的 (x)，(text{cnt}) 的取值。
(nle10^5,mle10,xle2^{20})
再次重申，(cap) 是按位与的意思。
首先我们用数学公式定义一下 (text{cnt})（因为公式复杂，所以加了 (tt large)）：
说明一下正确性。如果 (text{popcnt}(a_{i,j}cap x)) 是奇数的话，那么 ((-1)^{text{popcnt}(a_{i,j}cap x)}) 的结果就是 (-1)。最后 (1-(-1)^{text{popcnt}(a_{i,j}cap x)}) 就是 (2)，最后会被 (frac{1}{2^m}) 除去；如果 (text{popcnt}(a_{i,j}cap x)) 是偶数的话，那么 ((-1)^{text{popcnt}(a_{i,j}cap x)}) 的结果就是 (1)。最后 (1-(-1)^{text{popcnt}(a_{i,j}cap x)}) 就是 (0)，那么整行的结果都是 (0)。
然后我们发现它是可以展开的：
然后我们有一个性质：
也就是 (sum_{i=1}^ntext{popcnt}(a_icap x)) 的奇偶性与 (text{popcnt}((oplus_{i=1}^na_i)cap x)) 的相同。这点在上面的新的运算符 (circ) 的性质中有类似的体现。
容易得到：
我们发现一加一减的可以容斥，我们容斥计算 (f_i) 表示 (n) 行的所有式子中 ((-1)^i) 前面的系数和。
// num 处理到第几列 // x 当前的指数 // mu 当前的系数（+1 or -1） void dfs(ll *a, ll num, ll x, ll mu) { if(num > m) { f[x] += mu; return; } dfs(a, num + 1, x, mu); // 不加入第 num 列，系数不变 dfs(a, num + 1, x ^ a[num], -mu); } 这样我们就可以进一步化简：
我们突然发现后面这个 ((-1)^i) 取值只有两种，当 (xcap i) 是奇数时取值为 (-1)，否则为 (1)。
好了，现在我们的问题转换为了求出：
这不就是 FWT 中的异或变换吗：
综上，我们发现这题就是推式子容斥之后得到 FWT 的形式。
原题需要将输出加密：
#include <bits/stdc++.h> #define ll long long using namespace std; const ll P = 1e9 + 7; #define N 100010 #define M 20 #define K 21 ll n, m, k; ll f[1 << K]; ll a[N][M]; // num 处理到第几列 // x 当前的指数 // mu 当前的系数（+1 or -1） void dfs(ll *a, ll num, ll x, ll mu) { if(num > m) { f[x] += mu; return; } dfs(a, num + 1, x, mu); // 不加入第 num 列，系数不变 dfs(a, num + 1, x ^ a[num], -mu); } void Xor(ll *a, ll type) { for(ll x = 2; x <= (1 << k); x <<= 1) { ll z = x >> 1; for(ll i = 0; i < (1 << k); i += x) { for(ll j = 0; j < z; j ++) { (a[i + j] += a[i + j + z]) %= P; (a[i + j + z] = a[i + j] - 2 * a[i + j + z]) %= P; (a[i + j] *= type) %= P; (a[i + j + z] *= type) %= P; } } } } ll qpow(ll x, ll y) { if(y == 0) return 1; if(y % 2 == 1) return x * qpow(x, y - 1) % P; ll tmp = qpow(x, y / 2); return tmp * tmp % P; } int main() { while(scanf("%lld %lld %lld", &n, &m, &k) != EOF) { for(ll i = 0; i < (1 << k); i ++) f[i] = 0; for(ll i = 1; i <= n; i ++) { for(ll j = 1; j <= m; j ++) { scanf("%lld", &a[i][j]); } dfs(a[i], 1, 0, 1); } Xor(f, 1); ll pw = 1, inv = qpow(1 << m, P - 2), ans = 0; for(ll i = 0; i < (1 << k); i ++) { ans ^= f[i] * pw % P * inv % P; (pw *= 3) %= P; } printf("%lldn", ans); } } 「AT ABC212H」 Nim Counting
给定两个数 (N,K)，以及一个长度为 (K) 的整数数组 ((A_1,A_2,cdots, A_K))。
两个人玩 Nim 游戏。
现在通过以下方式生成一个游戏：
任意选择一个 (1le Mle N)，(M) 表示石子堆数。
对于每一堆，其石子数是 (A) 中任意一个数。
对于 (sum_{i=1}^N K^i) 种游戏，求先手获胜的游戏数，答案对 (998244353) 取模。
(nle2times10^5,Kle2^{16},a_ile2^{16})
根据玩 Nim 游戏的经验，可以发现先手获胜当且仅当 (bigoplus_{i=0}^n A_ineq 0)。
所以我们定义 dp 式子 (f_{i,j}) 表示有 (i) 个石堆，且石堆异或和为 (j) 的获胜方案数，有：
答案就是 (sum_{i=1}^nsum_{jneq0} f_{i,j})。
直接转移是朴素的，发现上面的式子刚好是 FWT 异或操作，也就是：
我们定义 (a) 是一个全是 (1) 的数组即可。
同时，我们发现其实不需要真的进行 (n) 次卷积，其实只需要将 FWT 变换过之后的结果 (A)，求出 (A+A^2+A^3+cdots+A^n) 即可。
上面的可以通过等比数列求和公式计算。
#include <bits/stdc++.h> using namespace std; #define ll long long #define P 998244353 const ll K = 1 << 20; ll n, k, ans; ll f[K]; void FWT(ll *a, ll type) { for(ll x = 2; x <= K; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < K; i += x) { for(ll j = 0; j < k; j ++) { (a[i + j] += a[i + j + k]) %= P; (a[i + j + k] = a[i + j] - 2 * a[i + j + k]) %= P; (a[i + j] *= type) %= P; (a[i + j + k] *= type) %= P; } } } } ll qpow(ll x, ll y) { if(y == 0) return 1; if(y % 2 == 1) return x * qpow(x, y - 1) % P; ll tmp = qpow(x, y / 2); return tmp * tmp % P; } int main() { scanf("%lld %lld", &n, &k); for(ll i = 1; i <= k; i ++) { ll x; scanf("%lld", &x); f[x] ++; } FWT(f, 1); for(ll i = 0; i < K; i ++) { if(f[i] == 1) f[i] = n; else { f[i] = f[i] * (qpow(f[i], n) - 1) % P * qpow(f[i] - 1, P - 2) % P; } } FWT(f, 499122177); for(ll i = 1; i < K; i ++) { (ans += f[i]) %= P; } printf("%lld", (ans % P + P) % P); } 「AT ARC100E」 Or Plus Max
给你一个长度为 (2^n) 的序列 (a)，每个(1le Kle 2^n-1)，找出最大的 (a_i+a_j)（(i cup j le K)，(0 le i < j < 2^n)）并输出。
(nle18)
就是要求 (max_{icup j=k} a_i+a_j)。
我们维护 (f_{i}) 表示 (max_{icup j=i} a_i)，(g_i=text{max2}_{icup j=i} a_i)，(text{max2}) 表示次大值。
然后就像 FWT 的或变换一样了。
#include <bits/stdc++.h> using namespace std; #define ll long long const ll N = 1 << 18; ll n; struct node { ll mx1, mx2; node(ll a = 0, ll b = 0):mx1(a), mx2(b) {} friend node operator +(const node &x, const node &y) { if(x.mx1 > y.mx1) { return node(x.mx1, max(x.mx2, y.mx1)); } return node(y.mx1, max(y.mx2, x.mx1)); } } a[N]; void FWT(node *a) { for(ll x = 2; x <= n; x <<= 1) { ll k = x >> 1; for(ll i = 0; i < n; i += x) { for(ll j = 0; j < k; j ++) { a[i + j + k] = a[i + j] + a[i + j + k]; } } } } int main() { scanf("%lld", &n); n = 1 << n; for(ll i = 0; i < n; i ++) { scanf("%lld", &a[i].mx1); } FWT(a); for(ll i = 0; i < n; i ++) { a[i].mx1 = a[i].mx1 + a[i].mx2; } ll ans = 0; for(ll i = 1; i < n; i ++) { ans = max(ans, a[i].mx1); printf("%lldn", ans); } } 「HDU 6618」 Good Numbers
定义一个正整数 (n) 是好数当且仅当 (n) 在 8 进制表示下所有的数码出现的次数为 3 的倍数（出现 0 次亦可）。
有多少个 (k) 位的 8 进制数（不含前导 0），满足这个数是好的，且是 (p) 的倍数。对 (10^9+9) 取模。
例如：当 (k=3,p=2) 时，好数有 (222,444,666) 三个。
(kle10^{18},p<8)
考虑状压 dp，设 (f_{i,s,j}) 表示第 (i) 位，(8) 种数出现次数对 (3) 取模的状压情况，以及数对 (p) 取模的结果为 (j)。
答案就是 (f_{k,0,0})。
直接暴力枚举位数转移是朴素的，瓶颈在于 (k)，考虑优化掉 (k)。
发现我们可以使用像快速幂一样的方法，也就是倍增 dp。
转移公式就是：
其中 (t) 是转移的位数，而 (oplus) 在这里是不进位三进制加法。
发现这样多了瓶颈——我们需要枚举 (s_1) 和 (s_2)。
但是我们发现这不就是 FWT 中异或的形式吗：(c_{ioplus j}gets a_ib_j)。考虑三进制 FWT 加速。下面给出 FWT 的代码，w1 是原根的一次方，w2 是原根的二次方：
void FWT(ll *a, ll type) { for (ll x = 3; x <= N; x *= 3) { ll k = x / 3; for (ll i = 0; i < N; i += x) { for (ll j = 0; j < k; j ++) { for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k]; if (type == 1) { tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; } else { tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P; } for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l]; } } } } 因为 (1e9+9) 存在原根 (2)，然后就朴素实现了，注意位矩阵：
代码：
#include <bits/stdc++.h> using namespace std; #define ll long long const ll P = 1e9 + 9; ll qpow(ll x, ll y) { if(y == 0) return 1; if(y % 2 == 1) return x * qpow(x, y - 1) % P; ll tmp = qpow(x, y / 2); return tmp * tmp % P; } const ll G = 2; const ll w1 = qpow(G, (P - 1) / 3); const ll w2 = qpow(G, (P - 1) / 3 * 2); const ll inv3 = qpow(3, P - 2); const ll N = 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3; ll n, p; ll tmp[8][N], res[8][N], one[8][N]; ll a[8][N], b[8][N]; ll pw3[8]; ll tmp1[3], tmp2[3]; void FWT(ll *a, ll type) { for (ll x = 3; x <= N; x *= 3) { ll k = x / 3; for (ll i = 0; i < N; i += x) { for (ll j = 0; j < k; j ++) { for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k]; if (type == 1) { tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; } else { tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P; } for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l]; } } } } ll base = 1; void fun(ll x) { if(x == 1) { memset(res, 0, sizeof res); memset(tmp, 0, sizeof tmp); memset(one, 0, sizeof one); for(ll i = 1; i < 8; i ++) res[i % p][pw3[i]] = 1; for(ll i = 0; i < 8; i ++) tmp[i % p][pw3[i]] = 1; for(ll i = 0; i < 8; i ++) one[i % p][pw3[i]] = 1; for (int i = 0; i < p; i ++) { FWT(tmp[i], 1); FWT(res[i], 1); FWT(one[i], 1); } base = 8 % p; return; } if(x % 2 == 1) { fun(x - 1); memset(a, 0, sizeof a); memset(b, 0, sizeof b); for (ll i = 0; i < p; i ++) { for (ll j = 0; j < p; j ++) { ll k = (i * 8 + j) % p; for (ll x = 0; x < N; x ++) (a[k][x] += tmp[i][x] * one[j][x]) %= P, (b[k][x] += res[i][x] * one[j][x]) %= P; } } memcpy(tmp, a, sizeof a); memcpy(res, b, sizeof b); (base *= 8) %= P; return; } fun(x / 2); memset(a, 0, sizeof a); memset(b, 0, sizeof b); for (ll i = 0; i < p; i ++) { for (ll j = 0; j < p; j ++) { ll k = (i * base + j) % p; for (ll x = 0; x < N; x ++) (a[k][x] += tmp[i][x] * tmp[j][x]) %= P, (b[k][x] += res[i][x] * tmp[j][x]) %= P; } } memcpy(tmp, a, sizeof a); memcpy(res, b, sizeof b); (base *= base) %= p; } int main() { pw3[0] = 1; for(ll i = 1; i <= 8; i ++) { pw3[i] = pw3[i - 1] * 3; } while(scanf("%lld %lld", &n, &p) != EOF) { fun(n); FWT(res[0], -1); printf("%lldn", res[0][0]); } } 「CF 1103E」Radix sum
给定一个长度为 (n) 的序列 (a_1,a_2,...,a_n)，对于每一个 (p in [0,n-1])，求满足下列条件的整数序列 (i_1,i_2,...,i_n) 的方案数，对 (2^{58}) 取模：
- (forall j in [1,n] , i_j in [1,n])；
- (sumlimits_{j=1}^n a_{i_j} = p)，这里的加法定义为十进制不进位加法。
(nle10^5,a_ile10^5)
我们可以想到 dp：设计状态 (f_{i,s}) 表示考虑到第 (i) 个数，当前加法状态为 (s)。因为 FWT 变换时线性的，可以先变换为 FWT 点值表示法，然后变成自己的 (n) 次幂，最后再变换回来。
上面是平凡的，但是题目给出了模数 (2^{58})。发现没有单位根，所以考虑扩域。
这里的分圆多项式 (Phi_{10}(x)=x^4-x^3+x^2-x+1)。
然而我们发现 IFWT 时，需要除去进制 (10)，然而我们发现 (10) 在 (2^{58}) 下没有逆元。实际上我们发现 (5) 在 (2^{58}) 下是有逆元的：(57646075230342349)，我们只需要再除去一个 (2) 就可以了。设已经除以了 (5) 的答案为 (x)，真正的答案为 (y)，也就是 (2^5yequiv xpmod{2^{64}})，显然，我们有 (yequiv frac{x}{2^5}pmod{2^{64-5}})，也就是 (yequiv frac{x}{2^5}pmod{2^{59}})，所以直接将最后的答案除以 (2^5) 即可。虽然出题人不知道为什么要模 (2^{58})，但再取下模即可。
然后就是平凡实现了：
#include <bits/stdc++.h> using namespace std; #define ll unsigned long long const ll P = 1ull << 58, N = 1e5 + 10; const ll m = 5, K = 10; ll inv5; ll n; ll pw[m + 1]; ll qpow(ll x, ll y) { if(y == 0) return 1; if(y % 2 == 1) return x * qpow(x, y - 1); ll tmp = qpow(x, y / 2); return tmp * tmp; } struct poly { ll a[30]; poly() {memset(a, 0, sizeof a);} ll operator [](ll x) const {return a[x];} ll& operator [](ll x) {return a[x];} friend poly operator *(const poly &x, const poly &y) { poly z; for(ll i = 0; i < K; i ++) { for(ll j = 0; j < K; j ++) { z[(i + j) % K] += x[i] * y[j]; } } return z; } friend poly operator *(const poly &x, const ll &y) { poly z; for(ll i = 0; i < K; i ++) { z[i] += x[i] * y; } return z; } friend poly operator +(const poly &x, const poly &y) { poly z; for(ll i = 0; i < K; i ++) { z[i] += x[i] + y[i]; } return z; } poly w(ll x) { poly res; for(ll i = 0; i < K; i ++) { res[(i + x) % K] += a[i]; } return res; } } T, f[N], one; poly qpow(poly x, ll y) { if(y == 0) return one; if(y % 2 == 1) return x * qpow(x, y - 1); poly tmp = qpow(x, y / 2); return tmp * tmp; } poly tmp1[30], tmp2[30]; void FWT(poly *a, ll type) { for(ll x = K; x <= pw[m]; x *= K) { ll k = x / K; for(ll i = 0; i < pw[m]; i += x) { for(ll j = 0; j < k; j ++) { for(ll l = 0; l < K; l ++) tmp1[l] = a[i + j + l * k], tmp2[l] = poly(); if(type == 1) { for(ll l = 0; l < K; l ++) { for(ll v = 0; v < K; v ++) { tmp2[l] = tmp2[l] + tmp1[v].w(l * v % K); } } for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l]; } else { for(ll l = 0; l < K; l ++) { for(ll v = 0; v < K; v ++) { tmp2[l] = tmp2[l] + tmp1[v].w((K - (l * v % K)) % K); } } for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l] * inv5; } } } } } ll mod(poly x){ ll n = 4; for(ll i = K - 1; i >= n; i --){ ll u = x[i]; for(ll j = 1; j <= n; j ++) x[i - j] -= u * T[n - j]; } ll u = x[0]; u >>= m; return u % P; } int main() { pw[0] = 1; for(ll i = 1; i <= m; i ++) pw[i] = pw[i - 1] * K; T[0] = 1, T[1] = -1, T[2] = 1, T[3] = -1, T[4] = 1; // 分圆多项式phi10 one[0] = 1; inv5 = 57646075230342349ull; scanf("%llu", &n); for(ll i = 1; i <= n; i ++) { ll x; scanf("%llu", &x); f[x][0] ++; } FWT(f, 1); for(ll i = 0; i < pw[m]; i ++) f[i] = qpow(f[i], n); FWT(f, -1); for(ll i = 0; i < n; i ++) cout<<mod(f[i])<<'n'; } 
