FFT，FWT算法比对——HDU6966. I love sequences（FWT）
FFT——傅里叶变换
离散傅里叶变换
离散傅里叶变换就是将一串长度为\(n\)的序列\(X_n\)按照指定的方式进行线性变换。
约定以下的\(n\)均为\(2\)的幂次，矩阵的行号和列号均从\(0\)开始。
傅里叶变换的矩阵为\(A_n=(a_{i,j})_{n\times n}\)，其中
可以证明\(A_n\)可逆。
显然，如果朴素的进行离散傅里叶变换，复杂度将会是\(O(n^2)\)的。下面考虑分治来优化复杂度。
设\(X_n=[x_0,x_1,\dots,x_{n-1}]^T\)是进行傅里叶变换的序列，\(M_n=A_nX_n\)，\(M_n=(m_k)_{n\times1}\)，
当\(k<\frac{n}{2}\)时，经过变换后的\(M_n\)列向量的第\(k\)行和第\(k+\frac{n}{2}\)行有如下的递推式
根据上式，要对整体长度为\(n\)的序列进行傅里叶变换，就需要分别对偶数序列和奇数序列进行傅里叶变换，然后再通过上面的式子转化为整体序列的傅里叶变换。故可递归地处理，然而递归处理常数较大，有没有什么办法迭代处理呢？
为了方便表示，记
设变换\(C_n\)，其意义为：将序列的偶数项放在前半部分，将奇数项提出放在后半部分，则有
变换\(C_n\)的目的是将奇数项和偶数项分别放在两段连续的区间里
根据变换\(C_n\)的意义，矩阵\(C_n\)的第\(i\)行，第\(j\)列满足下式：
显然\(C_n\)可逆。
从而\(M_n\)向量有如下递推式：
由于\(X_n\)是任意的，所以
观察上式，\(A_n\)总是可以表示一个矩阵左乘\(C_n\)，不妨定义矩阵\(B_n\)
由于\(C_n\)可逆，所以\(B_n\)也可逆，故可设\(H_nB_n=A_n\)，则有
所以
由于\(H_1B_1=[1],B_1=[1]\)，所以\(H_1=[1]\)。
这里的\(H_n\)可以很方便地通过循环迭代实现。
设执行\(H_n\)变换所需要的工作量为\(T(n)\)，显然\(T(n)=2T(n/2)+O(n)\)，故而\(T(n)=O(n\log n)\)。
而傅里叶变换\(A\)可拆分成先进行\(B\)变换，再进行\(H\)变换的复合变换，即\(A=HB\)，所以接下来研究变换\(B\)的性质。
在二进制下我们发现\(C_n\)总是将最低位为\(0\)的放在前半部分，最低位为\(1\)的放在后半部分。而正常排列的二进制数，总是最高位为\(0\)的在前半部分，最高位为\(1\)的放在后半部分，故而我们想到了二进制逆向表示。
记：在\(\log n\)位的二进制下，\(i\)二进制表示的逆向表示为\(r(i,n)\)。
猜想：\(B_n\)变换将下标\(i\)与下标\(r(i,n)\)的数进行对换
形式化地写出猜想即为，\(B_m=(b_{i,j})\)，其中
接下来使用数学归纳法证明此猜想。
由于
\[B_2=\begin{bmatrix} B_1&0\\ 0&B_1\\ \end{bmatrix}C_2=\begin{bmatrix} 1&0\\ 0&1\\ \end{bmatrix}\begin{bmatrix} 1&0\\ 0&1\\ \end{bmatrix}=\begin{bmatrix} 1&0\\ 0&1\\ \end{bmatrix} \]所以当\(m=2\)时，猜想成立。
假设当\(m=\frac{n}{2}\)时，猜想成立，当\(m=n\)时
\[B_{n}=\begin{bmatrix} B_\frac{n}{2}&O\\ O&B_\frac{n}{2} \end{bmatrix}C_{n} \]设\(D_n=\begin{bmatrix} B_\frac{n}{2}&O\\ O&B_\frac{n}{2} \end{bmatrix}\)，由归纳假设，\(D_n\)的第\(i\)行，第\(j\)列为
\[d_{i,j}=\begin{cases} 1,&(i<\frac{n}{2}\land j=r(i,\frac{n}{2}))\lor (i\geq\frac{n}{2}\land j-\frac{n}{2}=r(i-\frac{n}{2},\frac{n}{2}))\\ 0,&\text{其他} \end{cases} \]所以
\[\begin{aligned} b_{i,j}&=\sum\limits_{k=0}^{n-1}d_{i,k}c_{k,j}\\ &=\sum\limits_{k=0}^{n-1}\begin{cases} 1,&( i<\frac{n}{2} \land j=2k \land k=r(i,\frac{n}{2})) \lor ( i\geq\frac{n}{2} \land j=2k-n+1 \land k-\frac{n}{2}=r(i-\frac{n}{2},\frac{n}{2}))\\ 0,&\text{其他} \end{cases} \end{aligned} \]由于当\(0\leq i<n\land 0\leq j <n\)时，有
\[\begin{aligned} ( &i<\frac{n}{2} \land j=2k \land k=r(i,\frac{n}{2})) \lor ( i\geq\frac{n}{2} \land j=2k-n+1 \land k-\frac{n}{2}=r(i-\frac{n}{2},\frac{n}{2}))\\ \iff& ( i<\frac{n}{2} \land j=2r(i,\frac{n}{2}) \land k=\frac{j}{2}) \lor ( i\geq\frac{n}{2} \land j=2r(i-\frac{n}{2},\frac{n}{2})+1 \land k=\frac{j+n-1}{2}))\\ \iff& ((i<\frac{n}{2} \land j=2r(i,\frac{n}{2})) \lor (i\geq\frac{n}{2} \land j=2r(i-\frac{n}{2},\frac{n}{2})+1)) \land (k=\frac{j}{2}\lor k=\frac{j+n-1}{2})\\ \iff& j=r(i,n)\land (k=\frac{j}{2}\lor k=\frac{j+n-1}{2}) \end{aligned} \]所以
\[d_{i,j}=\sum\limits_{k=0}^{n-1}\begin{cases} 1,&j=r(i,n)\land (k=\frac{j}{2}\lor k=\frac{j+n-1}{2})\\ 0,&\text{其他} \end{cases} =\begin{cases} 1,&j=r(i,n)\\ 0,&\text{其他} \end{cases}=\begin{cases} 1,&i=r(j,n)\\ 0,&\text{其他} \end{cases} \]故当\(m=n\)时，猜想也成立。
故而执行\(B_n\)变换，可以先预处理出\(r(i,n)\)，即二进制逆序，这个过程显然我们可以使用\(O(n)\)递推。
在\(\log n\)位二进制下，\(i\)的逆向可以由\(\frac{i}{2}\)的逆向得出，\(\frac{i}{2}\)给最高位空出了一位\(0\)，剩下的位数为\(i\)的低\(\log n-1\)位，所以\(\frac{i}{2}\)逆序最低位必为\(0\)，高\(\log n-1\)位的数即为\(i\)的低\(\log n-1\)位在\(\log n-1\)位二进制下的逆序，所以要让\(\frac{i}{2}\)的逆序右移一位，作为\(i\)逆序的低\(\log n-1\)位的贡献。至于\(i\)逆序最高位的贡献，由\(i\)的最低位提供，故有如下递推式。
rev[i]=(rev[i>>1]>>1)|((i&1)<<(n-1))
所以最终快速傅里叶变换的复杂度为\(O(n\log n)\)。
离散傅里叶逆变换
关于傅里叶逆变换，可以有两种方式理解，
第一种，可以验证\(A\)的逆矩阵满足
设\(E=AA^{-1}\)
\[e_{i,j}=\sum\limits_{k=0}^{n-1}a_{i,k}(a^{-1})_{k,j} =\frac{1}{n}\sum\limits_{k=0}^{n-1}\omega_n^{(i-j)k} \]若\(i=j\)，则\(e_{i,j}=\frac{1}{n}\cdot n=1\)
若\(i\neq j\)，则\(e_{i,j}=\frac{1}{n}\frac{1-\omega_n^{(i-j)n}}{1-\omega_n^{i-j}}=\frac{1}{n}\frac{1-1}{1-\omega_n^{i-j}}=0\)
所以傅里叶逆变换本质上也是一种傅里叶变换，只是将\(\omega_n\)替换成了\(\omega_n^{-1}\)，可以想见，和上述的变换过程几乎一模一样，只是需要最后将所有元素除以\(n\)。
第二种，既然\(A_n=H_nB_n\)，那么\(A_n^{-1}=B_n^{-1}H_n^{-1}\)，由于\(A_n\)是对称矩阵，所以\(A_n^{-1}\)也是对称矩阵，并且\(B_n\)的逆变换由意义可以得出就是它本身，加之\(B_n\)也是对称矩阵。
故而
由于
可以计算
故而\(H_n^{-1}\)的转置
而\(W_n^{-1}\)显然为
\(H_n\)由\(H_1\)推出一共要推\(\log n\)次，所以系数\(\frac{1}{2}\)最终乘了\(\log n\)次，故最终相当于除以\(n\)，这和第一种理解是一致的。
代码：
#include <algorithm>
#include <cmath>
using namespace std;
const double pi = acos(-1);
const int maxn = 2e5;
struct Complex {
double re, im;
Complex operator+(const Complex& rhs) const {
return {re + rhs.re, im + rhs.im};
}
Complex operator-(const Complex& rhs) const {
return {re - rhs.re, im - rhs.im};
}
Complex operator*(const Complex& rhs) const {
return {re * rhs.re - im * rhs.im, re * rhs.im + im * rhs.re};
}
Complex operator/(double rhs) const { return {re / rhs, im / rhs}; }
};
int rev[maxn];
// n需要代入2的幂次
// f = 1: fft
// f = -1： 逆fft
void FFT(Complex* a, int n, int f) {
// 执行B变换
for (int i = 0; i < n; i++)
if (i < rev[i]) swap(a[i], a[rev[i]]);
// 执行H变换
for (int i = 1; i < n; i <<= 1) {
// 执行H_{2i}变换
Complex d_omega = {cos(pi / i), f * sin(pi / i)};
for (int j = 0; j < n; j += i << 1) {
Complex omega = {1, 0};
for (int k = 0; k < i; k++) {
Complex x = a[k + j], y = omega * a[k + j + i];
a[k + j] = x + y;
a[k + j + i] = x - y;
omega = omega * d_omega;
}
}
}
if (f == -1)
for (int i = 0; i < n; i++) a[i] = a[i] / n;
}
// k: 大于等于n的最小的满足2的幂次的数
// l: log(k)
int k, l;
void init(int n) {
k = 1, l = 0;
while (k < n) k <<= 1, l++;
// 预处理二进制逆序
rev[0] = 0;
for (int i = 1; i < k; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
傅里叶变换的性质
介绍了傅里叶变换，却还没有提及它的用途，所以我们要研究它的性质。
先提出卷积的概念，设有两个序列\(X\)和\(Y\)，\(X\)和\(Y\)的卷积记作\(Z=X*Y\)，则有
记序列\(X\)和序列\(Y\)对应位置相乘为\(Z=X\cdot Y\)，则有
记对序列\(X\)的傅里叶变换为\(Z=FFT(X)\)，则有
傅里叶变换的性质为
下面证明这件事
考虑序列第\(k\)位的情况，要证明结论，就要证明
\[\left(\sum\limits_ia_{k,i}x_i\right) \left(\sum\limits_ja_{k,j}y_j\right)=\sum\limits_ta_{k,t}\sum\limits_{i+j=t}x_iy_j\\ \sum\limits_i\sum\limits_ja_{k,i}a_{k,j}x_i y_j=\sum\limits_t\sum\limits_{i+j=t}a_{k,t}x_iy_j\\ \sum\limits_i\sum\limits_ja_{k,i}a_{k,j}x_i y_j=\sum\limits_i\sum\limits_{j}a_{k,i+j}x_iy_j\\ a_{k,i}a_{k,j}=a_{k,i+j} \]显然傅里叶变换矩阵中的元素\(a_{i,j}=\omega_n^{ij}\)满足上式，故性质成立。
通过上面的讨论，我们发现这个性质最大的用处就是加速卷积计算，原本朴素计算卷积需要\(O(n^2)\)的时间复杂度，现在可以通过傅里叶变换，转化成对应乘积，将对应乘积的结果进行傅里叶逆变换，就是最终卷积的结果。傅里叶变换与傅里叶逆变换的时间复杂度都是\(O(n\log n)\)，对应乘积的复杂度是\(O(n)\)，所以最终卷积的复杂度为\(O(n\log n)\)。
FWT——沃尔什变换
通常意义上的卷积为
如果下面的式子不是\(i+j=k\)，而是\(i\oplus j=k\)，其中\(\oplus\)是按位运算的二元运算符（注意：不一定是二进制），这种卷积叫做位运算卷积。
此时，如果朴素地计算位运算卷积，其复杂度是\(O(n^2)\)的，考虑如何优化，我们可以参照FFT的思路。
傅里叶变换有良好的性质
我们能否构造一个变换使得对于位运算卷积也有这样良好的性质呢？即
仿照FFT性质的证明，要满足这样的性质，就要求变换矩阵的元素满足
由于是按位运算，我们可以按位考虑，所以可以将\(a_{k,i}\)定义为每一位\(a_{k,i_x}\)的乘积
由于每一位的取值是有限的，故而可以暴力地列出每一位的运算规则，最终通过解方程解出\(a_{k,i}\)，具体的操作，我们以常见的与、或、异或为例。我们设FWT变换的矩阵为\(A_n\)，对于二进制而言\(n\)为\(2\)的幂次。
与卷积
与运算的表为
| \(i\) | \(j\) | \(i\land j\) |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 0 |
| 1 | 0 | 0 |
| 1 | 1 | 1 |
故而
可以解得
由于需要保证变换可逆，所以我们尽量选择线性无关的解
由上面的解，可以直接得出\(n=2\)时的变换\(A_2=\begin{bmatrix}0&1\\1&1\end{bmatrix}\)，由于\(a_{k,i}\)定义为每一位\(a_{k,i_x}\)的乘积
所以我们可以得到矩阵递推式（注意：矩阵前一半列的列号二进制最高位为\(0\)，后一半列的列号二进制最高位为\(1\)）
故而可以进行递推。
其逆变换只需要求逆矩阵即可
或卷积
或运算的表为
| \(i\) | \(j\) | \(i\lor j\) |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 1 |
故而
可以解得
依样画葫芦，
得出\(n=2\)时的变换\(A_2=\begin{bmatrix}1&0\\1&1\end{bmatrix}\)，递推式为
逆递推式为
异或卷积
异或运算的表为
或运算的表为
| \(i\) | \(j\) | \(i\oplus j\) |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |
故而
可以解得
依样画葫芦，
得出\(n=2\)时的变换\(A_2=\begin{bmatrix}1&-1\\1&1\end{bmatrix}\)，递推式为
逆递推式为
// n需要代入2的幂次
// f = 1: fwt
// f = -1： 逆fwt
void FWT(int* a, int n, int f) {
for (int i = 1; i < n; i <<= 1) {
for (int j = 0; j < n; j += i << 1) {
for (int k = 0; k < i; k++) {
int x = a[j + k], y = a[j + k + i];
if (f == 1) {
// a[j + k] = ...
// a[j + k + i] = ...
} else {
// a[j + k] = ...
// a[j + k + i] = ...
}
}
}
}
}
HDU6966. I love sequences
题面：
问题描述：
给定三个序列\(a=[a_1,a_2,\dots,a_n],b=[b_1,b_2,\dots,b_n],c=[c_1,c_2,\dots,c_n]\)每一个含有\(n\)个非负整数。你需要计算
其中序列\(d=[d_{p,1},d_{p,2},\dots,d_{p,n}]\)按以下的规则计算：
对于每一个整数\(p(1\leq p\leq n)\)，每一个整数\(k(1\leq k\leq n)\)，有
其中\(\oplus\)代表一个三进制二元按位运算符，满足\(k_x=\gcd(i_x,j_x)\)。
比如\(12\oplus 7=(110)_3\oplus (021)_3=(111)_3=13\)
输入：
这个问题只有一个样例。
第一行包含一个整数\(n(n\leq2\times10^5)\)，表示序列\(a,b,c\)的长度
第二行包含\(n\)个整数，表示\(a_1,a_2,\dots,a_n\)
第三行包含\(n\)个整数，表示\(b_1,b_2,\dots,b_n\)
第四行包含\(n\)个整数，表示\(c_1,c_2,\dots,c_n\)
输出：
输出一个整数，为答案模\(10^9+7\)
样例输入：
4
2 3 4 5
1 2 3 4
9 8 7 6
样例输出：
1643486
分析：
位运算卷积考虑FWT，考虑\(f(k,i)f(k,j)=f(k,i\oplus j)\)成立的充要条件
解得
形成矩阵
其逆矩阵
代码：
#include <cstdio>
typedef long long Lint;
const int maxn = 540000;
const int mod = 1e9 + 7;
int a[maxn], b[maxn], c[maxn];
int aa[maxn], bb[maxn], cc[maxn];
void FWT(int* a, int n, int flag) {
for (int len = 1; len < n; len *= 3) {
for (int i = 0; i < n; i += len * 3) {
for (int j = 0; j < len; j++) {
int x = a[i + j], y = a[i + j + len], z = a[i + j + 2 * len];
if (flag == 1) {
a[i + j + len] = (x + z) % mod;
a[i + j + 2 * len] = ((x + y) % mod + z) % mod;
} else {
a[i + j + len] = (mod - y + z) % mod;
a[i + j + 2 * len] = (mod - x + y) % mod;
}
}
}
}
}
int cal(int n, int c) {
int tot = 1;
while (tot <= n)
tot *= 3;
for (int i = 1; i <= n; i++) {
aa[i] = a[i];
bb[i] = b[i];
}
for (int i = n + 1; i < tot; i++) {
aa[i] = bb[i] = 0;
}
FWT(aa, tot, 1), FWT(bb, tot, 1);
for (int i = 0; i < tot; i++) {
cc[i] = (Lint)aa[i] * bb[i] % mod;
}
FWT(cc, tot, -1);
int res = 0;
for (int k = 1, ck = c; k < tot; k++, ck = (Lint)ck * c % mod) {
res = (res + (Lint)cc[k] * ck % mod) % mod;
}
return res;
}
int main() {
int n;
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
scanf("%d", a + i);
}
for (int i = 1; i <= n; i++) {
scanf("%d", b + i);
}
for (int i = 1; i <= n; i++) {
scanf("%d", c + i);
}
int ans = 0;
for (int p = 1; p <= n; p++) {
ans = (ans + cal(n / p, c[p])) % mod;
}
printf("%d\n", ans);
return 0;
}
浙公网安备 33010602011771号