FFT
FFT和NTT是循环卷积,如果数组开小了,高位的值会加到低位上去
DFT
\[f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7\\ f(x)=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6)\\ f(x)=G(x^2)+xH(x^2)\\ \\ f(\omega_n^k)=G(\omega_{\frac{n}{2}}^k)+\omega_{n}^{k}\times H(\omega_{\frac{n}{2}}^k)\\ f(\omega_n^{k+\frac{n}{2}})=G(\omega_{\frac{n}{2}}^k)-\omega_{n}^{k}\times H(\omega_{\frac{n}{2}}^k) \]Comp tmp[MAX_N];
void DFT(Comp* f, int n, int rev) {
if (n == 1) return;//一阶单位根只有1
for (int i = 0; i < n; ++i) tmp[i] = f[i];
for (int i = 0; i < n; ++i) {
if (i & 1) f[n / 2 + i / 2] = tmp[i];
else f[i / 2] = tmp[i];
}//系数分开,偶左奇右
Comp *g = f, *h = f + n / 2;
DFT(g, n / 2, rev), DFT(h, n / 2, rev);//求g、h关于它们单位根的点值
Comp cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * rev / n));
for (int k = 0; k < n / 2;++k) {//cur=w_{n}^{k}
tmp[k] = g[k] + cur * h[k];
tmp[k + n / 2] = g[k] - cur * h[k];
cur *= step;
}
for (int i = 0; i < n; ++i) f[i] = tmp[i];
}
位逆序置换
模拟递归过程,原数组中下标为 \(i\) 的元素递归最后会到下标为 \(R(i)\) ,\(R(i)\) 与 \(i\) 二进制表示后为逆序关系
![image-20240122145951095](C:\Users\Zhouxj's Lenovo\AppData\Roaming\Typora\typora-user-images\image-20240122145951095.png)
void change(Complex y[], int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i) if (i < rev[i]) swap(y[i], y[rev[i]]);// 保证每对数只翻转一次
return;
}
IDFT
![image-20240122152513241](C:\Users\Zhouxj's Lenovo\AppData\Roaming\Typora\typora-user-images\image-20240122152513241.png)
这个矩阵的逆矩阵就是每个元素变为倒数,再除以\(2^{len}\)
\[\frac{1}{\omega_k}=e^{-i\frac{2\pi}{k}}=\cos(\frac{2\pi}{k})+\sin(-\frac{2\pi}{k}) \]/*
* len 必须是 2^k 形式
* on == 1 时是 DFT,on == -1 时是 IDFT
*/
void fft(Complex y[], int len, int on) {
change(y, len); // 位逆序置换
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h)); // wn:当前单位复根的间隔
for (int j = 0; j < len; j += h) {
// 计算当前单位复根,一开始是 1 = w^0_n,之后是以 wn 为间隔递增: w^1_n
Complex w(1, 0);
for (int k = j; k < j + h / 2; k++) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t; y[k + h / 2] = u - t;
w = w * wn;
}
}
}
// 如果是 IDFT,它的逆矩阵的每一个元素不只是原元素取倒数,还要除以长度 len。
if (on == -1) {
for (int i = 0; i < len; i++) {
y[i].x /= len;
}
}
}
NTT
\(\omega_n=g_n=G^{\frac{p-1}{n}}(G通常为3)\)
\(\frac{1}{g_n}\)就是直接求逆
const int maxn = (1 << 21) + 10,p = 998244353,_G = 3;
inline void getr(int lim){
for(ri i = 1;i < lim;++i) r[i] = (r[i >> 1] >> 1) | (i & 1 ? (lim >> 1) : 0);
}
inline int qp(int x,int k,int res = 1){
for(;k;k>>=1,x=1ll*x*x%p) if(k&1) res=1ll*res*x%p;
return res;
}
inline void ntt(int *f,int lim,int type){
getr(lim);
for(ri i = 1;i < lim;++i) if(i < r[i]) swap(f[i],f[r[i]]);
for(ri i = 2;i <= lim;i <<= 1){
ll g = gn[i];
if(type == -1) g = qp(g,p - 2);
for(ri j = 0;j < lim;j += i){
ll gi = 1;
for(ri k = j;k < j + (i >> 1);++k,gi = gi * g % p){
ll x = f[k],y = 1ll * f[k + (i >> 1)] * gi % p;
f[k] = (x + y) % p;
f[k + (i >> 1)] = (x - y + p) % p;
}
}
}
if(type == -1){
ll inv = qp(lim,p - 2);
for(ri i = 0;i < lim;++i) f[i] = f[i] * inv % p;
}
}
void mul(int *f,int *g,int n,int m){
int lim = 1;
while(lim <= n + m) lim <<= 1;
ntt(f,lim,1); ntt(g,lim,1);
for(ri i = 0;i < lim;++i) f[i] = 1ll * f[i] * g[i] % p;
ntt(f,lim,-1);
}
signed main(){
n = rd(); m = rd();
for(ri i = 2;i < maxn;i <<= 1) gn[i] = qp(_G,(p-1)/i);
for(ri i = 0;i <= n;++i) F[i] = rd();
for(ri i = 0;i <= m;++i) G[i] = rd();
mul(F,G,n,m);
for(ri i = 0;i <= n+m;++i) printf("%lld ",F[i]); puts("");//都是从低到高
return 0;
}
标签:frac,复习,int,FFT,rev,NTT,len,++,omega
From: https://www.cnblogs.com/Lumos1921/p/17980569