P3803 【模板】多项式乘法（FFT）

(color{#0066ff}{题目描述})

给定一个n次多项式F(x)，和一个m次多项式G(x)。

请求出F(x)和G(x)的卷积。

(color{#0066ff}{输入格式})

第一行2个正整数n,m。

接下来一行n+1个数字，从低到高表示F(x)的系数。

接下来一行m+1个数字，从低到高表示G(x)的系数

(color{#0066ff}{输出格式})

一行n+m+1个数字，从低到高表示F(x)∗G(x)的系数。

(color{#0066ff}{输入样例})

1 2
1 2
1 2 1

(color{#0066ff}{输出样例})

1 4 5 2

(color{#0066ff}{数据范围与提示})

保证输入中的系数大于等于 0 且小于等于9。

对于100%的数据：(n, m leq {10}^6) , 共计20个数据点，2s。

数据有一定梯度。

空间限制：256MB

(color{#0066ff}{题解})

对于两个多项式相乘，显然暴力是(O(n^2))的

如何优化呢？

我们知道，n+1个点可以唯一确定一个n次多项式

那么我们对于(A(x)*B(x)=C(x))

可以拆成这样

(left{egin{matrix}(a_1,b_1) \ (a_2,b_2) \(a_3,b_3) \ . \ . \ . \(a_{n+1},b_{n+1)}end{matrix} ight} * left{egin{matrix}(c_1,d_1) \ (c_2,d_2) \(c_3,d_3) \ . \ . \ . \(c_{n+1},d_{n+1)}end{matrix} ight} = left{egin{matrix}(a_1*c_1,b_1*d_1） \ (a_2*c_2,b_2*d_2) \(a_3*c_3,b_3*d_3) \ . \ . \ . \(a_{n+1}*c_{n+1},b_{n+1}*d_{n+1})end{matrix} ight})

上面这个点值表达式的运算显然是(O(n))的

我们现在要把A和B转成点值表达式，然后乘过去，最后再转换回来

现在考虑转成点值表达式

比如(A(x)=a_0+a_1x+a_2x^2)

那么点值表达式就是

((x_1,a_0+a_1x_1+a_2x_1^2)(x_2,a_0+a_1x_2+a_2x_2^2)(x_3,x_0+a_1x_3+a_2x_3^2))

但是带入求y用秦九韶是(O(n))的，而且要带入n个点，所以就(O(n^2))了

现在考虑优化这个过程

引入一个东西，单位复数根

定义(omega_n^n=1，omega_n)有n个

对于复数坐标系的两个点（可以理解为两个向量）

两个单位向量A，B，模长为(a=b=1)，角度为(alpha,eta)

欧拉公式：(e^{ix}=cos x+i*sin x)

所以二者相乘即为(ae^{ialpha}*be^{ieta}=abe^{i(alpha+eta)}=Ae^{iTheta})

也就是说，复数相乘集合意义，模长相乘，极角相加

因此，(w_n^n=1)，则(nTheta\%2pi=0)

(Theta=frac{2pi}{n})

定义(omega_0=(1,0),omega_1=omega_n^1)

则(omega_n^k=omega_0*omega_1^k)

而且有(omega_n^{k+frac{n}{2}}=-omega_n^k)，因为关于原点对称

还有(omega_{2n}^{2k}=omega_n^k),n份取k份等同于2n份取2k份

回到多项式(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})

我们硬凑一下，使得(A(x))共有(2^k)次方项（不足后面补系数0）

然后开始分治,把(A(x))拆成(A_0(x))和(A_1(x))

(A_0(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{frac{n-2}{2}})

(A_1(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{frac{n-2}{2}})

不难发现(A(x)=A_0(x^2)+x*A_1(x^2))

我们把(omega_n^k)带入

(A(omega_n^k)=A_0(omega_n^{2k})+omega_n^k*A_1(omega_n^{2k})=A_0(omega_{frac{n}{2}}^k)+omega_n^k*A_1(omega_{frac{n}{2}}^k))

(A(omega_n^{k+frac{n}{2}})=A_0(-omega_n^{k})-omega_n^k*A_1(-omega_n^{k})=A_0(omega_{frac{n}{2}}^k)-omega_n^k*A_1(omega_{frac{n}{2}}^k))

两式子差距只在正负号！

所以求出了(A_0,A_1)，相加即为第一个，相减即为第二个

也就是说，对于当前序列，我们求出了一半，另一半也出来了

那么我们去找一找对应关系

({a_0} ||||| {a_1} ||||| {a_2} ||||| {a_3} ||||| {a_4} ||||| {a_5} ||||| {a_6} ||||| {a_7})

({000} ||| {001} ||| {010} ||| {011} ||| {100} ||| {101} ||| {110} ||| {111})

(|||||||||||{a_0,a_2,a_4,a_6} ||||||||||||||||||||||||| {a_1,a_3,a_5,a_7}||||||||||||)

(||||||{a_0,a_4} |||||||||| {a_2,a_6} ||||||||||||| {a_1,a_5} |||||||||| {a_3,a_7} |||||)

({a_0} ||||| {a_4} ||||| {a_2} ||||| {a_6} ||||| {a_1} ||||| {a_5} ||||| {a_3} ||||| {a_7})

({000} ||| {100} ||| {010} ||| {110} ||| {001} ||| {101} ||| {011} ||| {111})

卧槽，这是二进制翻转啊

我们用(r_i)代表数i翻转后是几

则(r[i]=r[i>>1]>>1|(i&1)*(len<<1))

要求当前位置的翻转后的数，把当前最后一位删去

那么当前的数的r已经求出来了

将这个数翻转（注意，二进制位数固定，比如长度为8，那么(00000001 o 10000000)）,再<<1，

这样现在的数除了最高位其他都是当前的翻转数了

这时只要考虑原来最后一位是不是1即可

通过FFT，我们求出了A(x)的点值表达式

B(x)同理，然(A(x)*B(x)=C(x))，(O(n))求

现在的问题是怎么转回系数表达式

定理

(left{egin{matrix} y_0\y_1\y_2\.\.\.\y_{n-1}end{matrix} ight}=left{egin{matrix} 1 & 1 & 1 & ...& 1 \1 & omega_n & omega_n^2 & ... & omega_n^{n-1}\1 & omega_n^2 & omega_n^4 & ... & omega_n^{2(n-1)}\. & . & . & . & .\. & . & . & . & .\. & . & . & . & .\1 & omega_n^{n-1} & omega_n^{2(n-1)} & ... & omega_n^{(n-1)(n-1)}end{matrix} ight}*left{egin{matrix} a_0\a_1\a_2\.\.\.\a_{n-1}end{matrix} ight})

左面是点值表达式，右面是系数表达式

现在已知(Y=W*A)

则(A=Y*W^{-1})

于是要矩阵求逆

不过这个矩阵又有个定理

它的逆矩阵为：指数全取负，再所有数/n之后的矩阵

那么就简单了

#include<bits/stdc++.h>
#define LL long long
LL in() {
	char ch; LL x = 0, f = 1;
	while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
	for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
	return x * f;
}
const double pi = acos(-1);
const int maxn = 3e6 + 41;
struct node {
	double x, y;
	node(double x = 0, double y = 0): x(x), y(y) {}
	friend node operator + (const node &a, const node &b) { return node(a.x + b.x, a.y + b.y); }
	friend node operator - (const node &a, const node &b) { return node(a.x - b.x, a.y - b.y); }
	friend node operator * (const node &a, const node &b) { return node(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
	friend node operator / (const node &a, const double &b) { return node(a.x / b, a.y / b); }
}A[maxn], B[maxn], C[maxn];
int len, n, m, r[maxn];
void FFT(node *D, int flag) {
	for(int i = 0; i < len; i++) if(i < r[i]) std::swap(D[i], D[r[i]]);
	for(int l = 1; l < len; l <<= 1) {
		node w0(cos(pi / l), flag * sin(pi / l));
		for(int i = 0; i < len; i += (l << 1)) {
			node w(1, 0), *a0 = D + i, *a1 = D + i + l;
			for(int k = 0; k < l; k++, a0++, a1++, w = w * w0) {
				node tmp = *a1 * w;
				*a1 = *a0 - tmp;
				*a0 = *a0 + tmp;
			}
		}
	}
	if(!(~flag)) for(int i = 0; i < len; i++) D[i] = D[i] / len;
}
int main() {
	n = in(), m = in();
	for(len = 1; len <= n + m; len <<= 1);
	for(int i = 0; i <= n; i++) A[i] = in();
	for(int i = 0; i <= m; i++) B[i] = in();
	for(int i = 1; i < len; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) * (len >> 1));
	FFT(A, 1), FFT(B, 1);
	for(int i = 0; i < len; i++) C[i] = A[i] * B[i];
	FFT(C, -1);
	for(int i = 0; i <= n + m; i++) printf("%d%c", (int)round(C[i].x), i == n + m? '
' : ' ');
	return 0;
}

NTT模数写法

#include<bits/stdc++.h>
#define LL long long
LL in() {
    char ch; LL x = 0, f = 1;
    while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
    for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
    return x * f;
}
using std::vector;
const int mod = 998244353;
const int maxn = 3e6 + 10;
LL ksm(LL x, LL y) {
    LL re = 1LL;
    while(y) {
        if(y & 1) re = re * x % mod;
        x = x * x % mod;
        y >>= 1;
    }	
    return re;
}
        
void FNTT(vector<int> &A, int len, int flag) { 
	A.resize(len);
	int *r = new int[maxn];
	r[0] = 0;
	for(int i = 0; i < len; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) * (len >> 1));
	for(int i = 0; i < len; i++) if(i < r[i]) std::swap(A[i], A[r[i]]);
    for(int l = 1; l < len; l <<= 1) {
		int w0 = ksm(3, (mod - 1) / (l << 1));
		for(int i = 0; i < len; i += (l << 1)) {
			int w = 1, a0 = i, a1 = l + i;
			for(int k = 0; k < l; k++, a0++, a1++, w = 1LL * w * w0 % mod) {
                int tmp = 1LL * A[a1] * w % mod;
                A[a1] = ((A[a0] - tmp) % mod + mod) % mod;
                A[a0] = (A[a0] + tmp) % mod;
			}
		}
	}
	if(!(~flag)) {
		std::reverse(A.begin() + 1, A.end());
		int inv = ksm(len, mod - 2);
		for(int i = 0; i < len; i++) A[i] = 1LL * inv * A[i] % mod;
	}
	delete []r;
}

vector<int> operator * (vector<int> A, vector<int> B) {
    int tot = A.size() + B.size() - 1;
    int len = 1;
    while(len <= tot) len <<= 1;
    FNTT(A, len, 1);
    FNTT(B, len, 1);
    vector<int> ans;
    ans.resize(len);
    for(int i = 0; i < len; i++) ans[i] = 1LL * A[i] * B[i] % mod;
    FNTT(ans, len, -1);
    ans.resize(tot);
    return ans;
}
signed main() {
    int n = in(), m = in();
    vector<int> A, B, C;
    for(int i = 0; i <= n; i++) A.push_back(in());
    for(int i = 0; i <= m; i++) B.push_back(in());
    C = A * B;
    for(int i = 0; i <= n + m; i++) printf("%d%c", C[i], i == n + m? '
' : ' ');
    return 0;
}