FFT及NTT复习

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\) 二进制表示后为逆序关系


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


这个矩阵的逆矩阵就是每个元素变为倒数，再除以\(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;
}