FFT/NTT

前置芝士：复数及其加减乘运算，多项式。

FFT 是一个 将 用系数表示的多项式 转换为 其用点值表示法表示出来的形式。

注意：以下 \(n\) 皆默认是2的幂次。

补充：

复数：形如 \(a + bi\) 的数，其中 \(i^2 = -1\)， \(a,b\) 为实数。

多项式的系数表示法：简单来说，就是形如 \(\sum ^{n-1}_{i=0}a_{i}x^{i}\) 的柿子。

多项式的点值表示法：将 \(n\) 个不同的 \(x\) 带入多项式，得到 \(n\) 个不同的值，则根据这 \(n\) 个值可以确定出一个唯一的多项式，则这 \(n\) 个 \(x\) 的值及其对应的多项式的值即构成该多项式的点值表示。

复数可以在平面直角坐标系表示，其中横纵坐标表示 \(a,b\)。复数相乘法则：模长相乘，辐角相加，即 其到原点的距离相乘，其与 \(x\) 轴正半轴的夹角相加。

由于系数表示法的多项式进行乘法是 \(O(n^2)\) 的，点值表示法的多项式乘法是 \(O(n)\) 的，所以我们只要解决 系数表示转点值表示和其逆 的快速方法即可。于是我们就引出了下面的一种方法：

傅里叶说：以原点为圆心，1为半径画圆，然后将圆 \(n\) 等分，将其每个等分点对应的复数记为 \(w_n^k\)，分别带进多项式求得其点值表示。

进一步解释： \(w_n^k\) 即为 n 次等分圆从 \(x\) 轴正半轴开始，逆时针数，第 \(k\) 个等分点对应的复数。

我们记 \(w_n^1\) 为 \(n\) 次单位根，显然可知该“根”的 \(k\) 次方即 \(w_n^k\)。

这种点值表示方法即 离散傅里叶变换(DFT)。

但是显然傅里叶的方法还是 \(O(n^2)\) 的，难道没法优化了？

当然不。

我们记 \(A(b)\) 为 以 \(b\) 当 \(x\) 的一个多项式，然后目标多项式是 \(A(x)\)，则我们可以按其次数的奇偶分成两个多项式：

\(A^{[0]}(x) = a_0 + a_2 x ^2+ a_4 x^4+.....+a_{n-2}x^{n-2}\)

\(A^{[1]}(x) = a_1 + a_3 x^2+ a_5 x^4+.....+a_{n-1}x^{n-1}\)

然后可得 \(A(x) = A^{[0]}(x)+xA^{[1]}(x)\)

然后就可以得出：

\(A(w_n^k) = A^{[0]}(w_n^k) + w_n^kA^{[1]}(w_n^k)\)

\(A(w_{n}^{k+\frac{n}{2}}) = A^{[0]}(w_n^k) - w_n^kA^{[1]}(w_n^k)\)

上面的容易发现，下面的则是需要知道 \(w_{n}^{k+n}=w_n^k,w_{n}^{k+\frac{n}{2}}=-w_n^k\) 即可。

这样我们发现这两个东西只有一个符号的区别，所以求 \(A(w_n^k)\) 就转换为了求规模仅有一半的子问题 $ A^{[0]}(w_nk)$，求这个即可同时求出后面一半子问题的答案，然后合并即可。

显然需要 \(O(n)\) 进行 \(\log{n}\) 次，复杂度 \(O(n \log {n})\)

所以递归求解就不难想出啦。（代码不给了，太慢还和迭代版的码量差不多）

然后迭代版只需要求出按奇偶分开的最后顺序，并把数组按这个顺序重排一下，就可以用差不多的方法解决啦。

于是我们就解决了 FFT。

注意一下求 \(n\) 次单位根可以直接套柿子：\(w_n^1 = complex(\cos{2.0 * \pi / n} , \sin{2.0 * \pi / n});\)

接下来我们还需要 IFFT，即 FFT 的逆，用点值表示转换为系数表示，同样的道理，只是我们在求 \(n\) 次单位根的时候给后面的 \(sin\) 乘上一个-1，最后求答案的时候给每个数除以一个 \(n\)

还有一些小细节看代码：

#include<bits/stdc++.h>
const int maxn = 4e6 + 100;
using namespace std;

int n, m;
double PI = acos(-1);
int rev[maxn], li = 1, len;

struct Complex{
	double rl, i;
	
	Complex(){rl = 0, i = 0;}
	Complex(double real, double ii):rl(real), i(ii){}
}A[maxn], B[maxn];//手写函数，防卡

Complex operator + (Complex a, Complex b){return Complex(a.rl + b.rl, a.i + b.i);}
Complex operator - (Complex a, Complex b){return Complex(a.rl - b.rl, a.i - b.i);}
Complex operator * (Complex a, Complex b){return Complex(a.rl * b.rl - a.i * b.i, a.i * b.rl + a.rl * b.i);}
//这里建议根据复数定义手推一下

void FFT(Complex* a, int opt){
	for(int i = 0; i < li; i++)if(i < rev[i])swap(a[i], a[rev[i]]);//按位置重排
	
 for(int dep = 1; dep <= len; dep++){
		int mn = (1 << dep);
		Complex wn = Complex(cos(2.0 * PI / mn), opt * sin(2.0 * PI / mn));//求单位根，注意逆

		for(int k = 0; k < li; k += mn){
			Complex w = Complex(1, 0);//根，后面不断累乘

			for(int j = 0; j < mn / 2; j++){
				Complex t = w * a[k + j + mn / 2];
				Complex u = a[k + j];

				a[k + j] = u + t;
				a[k + j + mn / 2] = u - t;

				w = w * wn;
           //这里就稍微根据前面的柿子感性理解一下，之所以用u t来表示是因为复数运算费时间，这样能卡出一个常数（就是蝴蝶操作）
			}
		}
	}

	if(opt == -1){
		for(int i = 0; i <= n + m; i++)a[i].rl = a[i].rl / li;
	}//特判是否是逆
}

signed main(){
	cin >> n >> m;
	for(int i = 0; i <= n; i++){cin >> A[i].rl;}
	for(int i = 0; i <= m; i++){cin >> B[i].rl;}
	
	while(li <= n + m){li = (li << 1), len++;}//求一下多项式次数li

	for(int i = 0; i < li; i++){rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));}//最后的数组位置
	FFT(A, 1);
	FFT(B, 1);//FFT搞一下
	
	for(int i = 0; i <= li; i++){A[i] = A[i] * B[i];}//乘起来
	
	FFT(A, -1);//IFFT搞回来
	
	for(int i = 0; i <= n + m; i++)cout << (int)(A[i].rl + 0.5) << " ";//输出，注意四舍五入（精度误差）
	
	exit(0);
}

在模数不恰当时，一般只能使用FFT。

NTT

我们发现，FFT的精度和double类型的运算实在过于感人，常数巨大而且可能炸精度。

FFT的单位根是在复数域内满足需要性质的一个根，那么NTT就利用在数论内满足需要性质的另一个单位根，即原根。

不妨设 \(n\) 是 2 的正整数次幂。

考虑设 \(g_n=g^{\frac{p-1}{n}}\)，则因为 \(g\) 是原根，所以有如下性质：

\[g_{n}^{\frac{n}{2}}= -1 \pmod p \]

\[g_n^{n}=1\pmod p \]

\[g_{xn}^{xk}=g_{n}^k \pmod p \]

容易看出这些都是 FFT 中单位根需要的性质。

于是 NTT 只需要修改 FFT 中单位根的部分，并且对于 INTT，我们只需把 IDFT 中共轭复数改为原根的逆元即可。

当然，对于模数 998244353，其为 \(2^{23}\times 119 + 1\)，因此对于大概 \(8e6\) 以内大小都是能做的，但是很多模数没有比较大的2的正整数次幂整除 \(p-1\)，所以一般只能取模数是 998244353。

code:

#include<bits/stdc++.h>
#define fint register int
#define ll long long
const int maxn = 4e6 + 100, mod = 998244353, g = 3;
using namespace std;
int n, m, A1[maxn], B1[maxn], ans[maxn];//A1,B1为需要进行乘法的，ans为答案 
//以下开始是完整FFT 
int rev[maxn], li = 1, len, A[maxn], B[maxn];
int qpow(int x, int k){fint res = 1; for(; k; k = (k >> 1), x = 1ll * x * x % mod)if(k & 1)res = 1ll * res * x % mod; return res;}
void NTT(int* a, int opt){
	for(fint i = 0; i < li; i++)if(i < rev[i])swap(a[i], a[rev[i]]);
    for(fint dep = 1; dep <= len; dep++){
		fint mn = (1 << dep);
		fint gn = qpow(g, (mod - 1) / mn); if(opt == -1)gn = qpow(gn, mod - 2);
		for(fint k = 0; k < li; k += mn){
			fint w = 1;
			for(int j = 0; j < mn / 2; j++){
				fint t = 1ll * w * a[k + j + mn / 2] % mod;
				fint u = a[k + j];
				a[k + j] = (u + t) % mod;
				a[k + j + mn / 2] = (u - t) % mod;
				w = 1ll * w * gn % mod;
			}
		}
	}
	if(opt == -1){
		for(fint i = 0; i <= n + m; i++)a[i] = 1ll * a[i] * qpow(li, mod - 2) % mod;;
	}
}
void work(int a[], int b[], int lena, int lenb){
	for(fint i = 0; i <= lena; i++){A[i] = a[i];}
	for(fint i = 0; i <= lenb; i++){B[i] = b[i];}
	while(li <= lena + lenb){li = (li << 1), len++;}
	for(fint i = 0; i < li; i++){rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));}
	NTT(A, 1);
	NTT(B, 1);
	for(fint i = 0; i <= li; i++){A[i] = 1ll * A[i] * B[i] % mod;}
	NTT(A, -1);
	for(fint i = 0; i <= lena + lenb; i++)ans[i] = (A[i] % mod + mod) % mod;
}
signed main(){
	cin >> n >> m;
	for(fint i = 0; i <= n; i++){cin >> A1[i];}
	for(fint i = 0; i <= m; i++){cin >> B1[i];}
	work(A1, B1, n, m);
	for(fint i = 0; i <= n + m; i++)cout << ans[i] << " ";
}