[Luogu3803] 多项式乘法 [FFT/NTT]

题意
给定 n 次多项式 F(x) 和 m 次多项式 G(x)
求：F(x) ⊗ G(x) 。 n, m ≤ 10⁶

---

链接
https://www.luogu.org/problemnew/show/P3803

---

FFT
长度为 $2^n$ 的情况下
按照下标奇偶分类
显然有 $A(x)=A_1(x^2)+x{A}_2(x^2)$
代入 ${\omega }_n^k=e^{\frac{2i\pi k}{n}}$
$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})$
$A({\omega }_n^{k+n/2})=A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k})=A(-{\omega }_n^k)$
（实际上我们有 ${\omega }_n^{k+n/2}=-{\omega }_n^k$
逆变换只需要把单位根用其共轭复数替换，并将求出来的结果除以 n 。
（ DFT 是多点求值，可以表示成系数矩阵乘上范德蒙德矩阵得到点值）
（ IDFT 是反过来，那就是点值乘上范德蒙德矩阵的逆矩阵得到系数矩阵）
（然后求一求求一求搞一搞搞一搞就能得到关系了 -> Click Here
逆变换的具体推导过程也可以 -> Click Here
实现的时候我们把原数组下标二进制翻转就能直接得到递归到最底层的下标
利用它从最底层一层层向上合并即可

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<cctype>
#include<ctime>
using namespace std;
#define R register
const int MAXN = 2097152+16;
int N, M, Limit = 1, Rev[MAXN], Log;
struct Complex
{
	double r, i;
	inline Complex(const double& a = 0, const double& b = 0)
	{
		r = a;
		i = b;
	}
	inline Complex operator + (const Complex& o) const
	{
		return Complex(r + o.r, i + o.i);
	}
	inline Complex operator - (const Complex& o) const
	{
		return Complex(r - o.r, i - o.i);
	}
	inline void operator += (const Complex& o)
	{
		r += o.r, i += o.i;
	}
	inline void operator -= (const Complex& o)
	{
		r -= o.r, i -= o.i;
	}
	inline Complex operator * (const Complex& o) const
	{
		return Complex(r * o.r - i * o.i, r * o.i + i * o.r);
	}
	inline void operator *= (const Complex& o)
	{
		*this = *this * o;
	}
}A[MAXN], B[MAXN], Wn[MAXN][2];
void FFT(Complex* A, const bool& Type = 0)
{
	static Complex t;
	for (R int i = 0; i < Limit; ++i) if(Rev[i] > i) swap(A[i], A[Rev[i]]);
	for (R int Mid = 1; Mid < Limit; Mid <<= 1)
	{
		for (R int Len = Mid << 1, Pos = 0; Pos < Limit; Pos += Len)
		{
			for (R int Delta = 0; Delta < Mid; ++Delta)
			{
				t = Wn[Limit/Len*Delta][Type] * A[Pos + Mid + Delta];
				A[Pos + Mid + Delta] = A[Pos + Delta] - t;
				A[Pos + Delta] += t;
			}
		}
	}
}
int main()
{
	scanf("%d%d", &N, &M);
	for (R int i = 0; i <= N; ++i) scanf("%lf", &A[i].r);
	for (R int i = 0; i <= M; ++i) scanf("%lf", &B[i].r);
	while (Limit <= N + M) Limit <<= 1, ++Log;
	for (R int i = 0; i < Limit; ++i) Rev[i] = (Rev[i>>1]>>1) | ((i&1)<<(Log - 1));
	const double tpi = 2.0 * acos(-1.0) / Limit;
	register double tmp;
	for (R int i = 0; i < Limit; ++i) tmp = tpi * i, Wn[i][0] = Complex(cos(tmp), sin(tmp)), Wn[i][1] = Complex(cos(tmp), -sin(tmp));
	FFT(A); FFT(B);
	for (R int i = 0; i < Limit; ++i) A[i] *= B[i];
	FFT(A, 1);
	for (R int i = 0; i <= N + M; ++i) printf("%d ", (int)(A[i].r/Limit+0.5));
	return 0;
}

---

NTT
单位复根实现 FFT 的重要性质：
1. ${\omega }_n^k(0\le k<n)$ 互不相同 → 点值合法
敏感的话容易想到 $g^r(\mathrm{m}\mathrm{o}\mathrm{d}n)(0\le r<n)$ 也是如此

假设 $p=qn+1,{\omega }_n=g^q$ 且 $n$ 是 $2$ 的幂。 ←所以 ntt 也跟 ftt 一样，需要长度为 $2^n$
NTT 在 $(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 意义下进行。

2. ${\omega }_{2n}^{2k}={\omega }_n^k$ → 分治①
${\omega }_n=g^q$ 那么 ${\omega }_{2n}=g^{q/2}$ 因为 $p=qn+1=\frac{q}{2}\cdot 2n+1$
那么 ${\omega }_{2n}^{2k}=g^{kq}={\omega }_n^k$

3. ${\omega }_n^{k+n/2}=-{\omega }_n^k$
我们有 ${\omega }_n^{k+n/2}=g^{(k+n/2)q}=g^{kq}\cdot g^{nq/2}$
$g^{nq/2}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 是否成立呢？
显然 $g^{nq}=g^{p-1}=g^{\varphi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
$g^{nq/2}=\pm 1$ 且 $g^0=1̸=g^{nq/2}$ 所以 $g^{nq/2}=-1$
那么 ${\omega }_n^{k+n/2\equiv -{\omega }_n^k}(\mathrm{m}\mathrm{o}\mathrm{d}p)$

4. $k̸=0,\sum _{r=0}^{n-1}{\omega }_n^{kr}=0$
卜难证啊 一模一样

具体的 -> Click Here

然后换一换单位根就可以快乐 ntt
只要数据不会爆模数就可以啦？

$998244353=119\ast 2^{23}+1,g=3$
$1004535809,g=3$
忘记了或者遇到没见过的 ntt 模数可以手打一个找原根？

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<cctype>
#include<ctime>
using namespace std;
const int MAXN = 2929229;
const long long p = 998244353;
const long long g = 3;
int n, m, Lim = 1, Log;
long long A[MAXN], B[MAXN], Wn[MAXN][2];
int Rev[MAXN];
long long qpowm(long long a, long long b)
{
	long long ret = 1;
	a %= p;
	while (b)
	{
		if (b & 1)
		{
			ret *= a;
			ret %= p;
		}
		a *= a;
		a %= p;
		b >>= 1;
	}
	return ret;
}
inline long long Adjust(long long x)
{
	x = (x % p) + p;
	return (x>=p)?(x-p):x;
}
void NTT(long long *a, const bool& Type = 0)
{
	for (register int i = 0; i < Lim; ++i) if (Rev[i] > i) swap(a[i], a[Rev[i]]);
	register long long t;
	for (register int Mid = 1; Mid < Lim; Mid <<= 1)
	{
		for (register int Len = Mid << 1, Pos = 0, qwq = Lim / Len; Pos < Lim; Pos += Len)
		{
			for (register int Sub = 0; Sub < Mid; ++Sub)
			{
				t = Wn[qwq*Sub][Type] * a[Pos + Mid + Sub];
				a[Pos + Mid + Sub] = Adjust(a[Pos + Sub] - t);
				a[Pos + Sub] = (a[Pos + Sub] + t) % p;
//				cout << a[Pos+Mid+Sub] << " " << a[Pos+Mid] << endl;
			}
		}
	}
}
int main()
{
	scanf("%d%d", &n, &m);
	for (register int i = 0; i <= n; ++i) scanf("%lld", &A[i]);
	for (register int i = 0; i <= m; ++i) scanf("%lld", &B[i]);
	n += m;
	while (Lim <= n) Lim <<= 1, ++Log;
	for (register int i = 0; i < Lim; ++i) Rev[i] = (Rev[i>>1]>>1)|((i&1)<<(Log-1));
	const long long q = (p - 1) / Lim;
	register long long t = qpowm(g, q);
	const long long G = qpowm(qpowm(g, p-2), q);
	for (register int i = 0; i < Lim; ++i)
	{
		Wn[i][0] = (!i) ? 1 : (Wn[i-1][0] * t % p);
		Wn[i][1] = (!i) ? 1 : (Wn[i-1][1] * G % p);
//		cout << Wn[i][0] << " " << Wn[i][1] << endl;
	}
	NTT(A); NTT(B);
	for (register int i = 0; i < Lim; ++i) A[i] = A[i] * B[i] % p;
	NTT(A, 1);
	t = qpowm(Lim, p - 2);
	for (register int i = 0; i <= n; ++i) printf("%lld ", A[i] * t % p);
	return 0;
}

---

（无关）
才知道f社出新作了
火星了 两 个 月（震声）
我自尽罢

我去撒库拉抹油了.jpg