fft相关的复习
2024-09-06 06:21:44
任意长度卷积 CZT
就是一波推导
\[\begin{aligned}
b_i &= \sum_{j=0}^{n-1} \omega^{ij}a_j \\
&= \sum_{j=0}^{n-1} \omega^{\frac{i^2+j^2-(i-j)^2}{2}}a_j \\
&= \omega^{\frac{i^2}{2}} \sum_{j=0}^{n-1}\omega^{\frac{-(i-j)^2}{2}} a_j \omega^{j^2}
\end {aligned}
\]
b_i &= \sum_{j=0}^{n-1} \omega^{ij}a_j \\
&= \sum_{j=0}^{n-1} \omega^{\frac{i^2+j^2-(i-j)^2}{2}}a_j \\
&= \omega^{\frac{i^2}{2}} \sum_{j=0}^{n-1}\omega^{\frac{-(i-j)^2}{2}} a_j \omega^{j^2}
\end {aligned}
\]
后面是一个减法卷积,就可以扩展到2的幂次直接fft就好了。
2次dft计算卷积
考虑有两个长度为\(n = 2^k\)的序列\(a(x), b(x)\),我们要计算他们的dft。
构造序列\(p_k = a_k + ib_k, \; q_k = a_k - ib_k\),
有结论\(dft_q(k) = conj(dft_p((n - k) \mod n))\)。展开,考虑几何意义???
我们可以解出\(dft_a, dft_b\)。
再做一遍idft就可以了
拆系数fft
记\(M = \sqrt {mod}\),把\(x\)表示成\(x = a \times M + b, b < M\)。
\((a \times M + b)(c \times M + d) = ac \times M^2 + (ad + bc) \times M + bd\)
对每一项分开算,做7次dft就可以了。
套用上述介绍做法4次dft就够了。
实现上注意在idft的时候,直接把一个序列放在real,另一个放在imag,idft回来直接/N后计算贡献就好了。
以及我们可以直接在一个for里面做解出AB,reverse序列的事情。
下面是关键部分的代码。
poly realmain(poly a, poly b) {
int n = a.size(), m = b.size();
prepare(n + m - 1);
for (int i = 0; i < n; i++) A[i] = cpx(a[i] & 32767, a[i] >> 15);
for (int i = 0; i < m; i++) B[i] = cpx(b[i] & 32767, b[i] >> 15);
dft(A, fft_n); dft(B, fft_n);
for (int i = 0; i < fft_n; i++) {
int j = (fft_n - i) % fft_n;
cpx ax, ay, bx, by;
ax = (A[i] + A[j].conj()) * cpx(0.5, 0);
ay = (A[i] - A[j].conj()) * cpx(0, -0.5);
bx = (B[i] + B[j].conj()) * cpx(0.5, 0);
by = (B[i] - B[j].conj()) * cpx(0, -0.5);
C[j] = ax * bx + ay * by * cpx(0, 1.0);
D[j] = ay * bx + ax * by * cpx(0, 1.0);
}
dft(C, fft_n); dft(D, fft_n);
poly ans(n + m - 1, 0);
for (int i = 0; i < ans.size(); i++) {
lo ax = lo(C[i].x / fft_n + 0.5) % mod;
lo ay = lo(C[i].y / fft_n + 0.5) % mod;
lo bx = lo(D[i].x / fft_n + 0.5) % mod;
lo by = lo(D[i].y / fft_n + 0.5) % mod;
ans[i] = ax + ((by + bx) << 15) + (ay << 30);
ans[i] = (ans[i] % mod + mod) % mod;
}
return ans;
}
最新文章
- 关于SAP日期操作的几个函数
- Windows Server 2012 Recycle Bin corrupted
- PAT 1012. 数字分类 (20)
- Two Strings Are Anagrams
- burpsuite绕过本地javascripte上传文件
- CalParcess.php.
- OpenLayers学习记录(1)
- ruby 疑难点之—— attr_accessor attr_reader attr_writer
- ListView inside a ScrollView
- SAE 上传根目录不存在!请尝试手动创建:./Uploads/Picture/
- net core 静态文件
- apk反汇编之smali语法
- CEOI 2014 wall (最短路)
- BZOJ 1103: [POI2007]大都市meg(dfs序,树状数组)
- Linux中的cat、more、less、head、tail命令
- kubernetes job的原理
- android 与html交互java调js与js调java操作
- android active间数据传递
- 373. Find K Pairs with Smallest Sums (java,优先队列)
- HBase分布式集群部署与设计