多项式1

FFT,NTT

多项式

多项式，\(F(x)=\sum_{i=0}^{n-1}a_ix^{i}\)，\(a_i\) 为常数。

上面的式子可以 \(O(n)\) 求出。

两个多项式相加，即把对应的系数相加即可，复杂度 \(O(n)\)。

两个多项式相乘，由于 \(x^n=x^0x^n=x^1x^{n-1}=x^2x^{n-2}…\)，所以复杂度会飙到 \(O(n^2)\)。

其实高精度本质上就是多项式的运算。

点值表示

由于 \(O(n^2)\) 实在太差劲了，考虑优化。

直接用系数运算一看就没法优化，没有前途。

发现如果有 \(n\) 个在 \(F(x)\) 定点 \(\{(x_1,y_1),(x_2,y_2),…,(x_n,y_n)\}\)，也可以确定一个 \(n-1\) 次的多项式 \(F(x)\)。如果发现不了就当结论。

可是我们要解决乘法，这玩意怎么乘？设 \(F(x),G(x)\) 的次数分别为 \(n,m\)，我们取 \(n+m+1\) 个点，暴力求值，复杂度依然 \(O(n^2)\)，甚至常数变大不少。

但是我们发现，点值在相乘这一块有着很大的优势，因为由 \((x_i,F(x_i)),(x_i,G(x_i))\) 两个点我们可以直接导出 \((x_i,F(x_i)G(x_i))\) 这个在多项式 \(F(x)\times G(x)\) 上的点。

所以我们的任务是快速求值。

单位根

默认会复数及其三角形式，不会可以自行查找，比较简单。

\(x^n=1\) 这玩意在复数域有 \(n\) 个解，且很好求。因为复数乘法是模长相乘，幅角相加，所以 \(x\) 的幅角肯定是 \(k\times\frac{2\pi}{n}\)，模长肯定是 \(1\)。

设 \(w_n\) 为幅角为正且最小的那个，那么这 \(n\) 个解为 \(w_n^0,w_n^1,w_n^2,…,w_n^{n-1}\)，由于 \(w_n^x=w_{n}^{n+x}\)，所以再往后推也没有意义。

根据三角形式的转化，我们可以求出 \(w_{n}^{k}=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}\)。

单位根有一些性质。

• \(w_{2n}^{2k}=w_{n}^{k}\)

相当于 \(\frac{2k\pi}n\) 分子分母各扩大两倍，分数值不变。

• \(w_{n}^{k+\frac n2}=-w_{n}^{k}\)

相当你在复平面上转的半圈，值变负。

快速傅里叶变换

我们认为 \(n=2^k,k\in \Bbb Z\)。相当于多项式超出的部分系数为 \(0\)。

由于点值表示法，我们也可以把 \(n\) 个单位根带进去，也能表示这个多项式。

\[F(x)=a_0x^0+a_1x^1+a_2x^2+…+a_{n-1}x^{n-1} \]

\[=(a_0x^0+a_2x^2+a_4x^4+…)+(a_1x^1+a_3x^3+a_5x^5+…) \]

\[=(a_0x^0+a_2x^2+a_4x^4+…)+x(a_1x^0+a_3x^2+a_5x^4+…) \]

令 \(F_1(x)=(a_0x^0+a_2x^1+a_4x^2+…),F_2(x)=(a_1x^0+a_3x^1+a_5x^2+…)\)，注意指数变了。

\[F(x)=F_1(x^2)+xF_2(x^2) \]

将单位根 \(w^k_n,0\le k<\frac n2\) 放进去得到。

\[F(w^k_n)=F_1(w^{2k}_n)+w^k_nF_2(w^{2k}_n) \]

\[=F_1(w^{k}_\frac n 2)+w^k_nF_2(w^{k}_\frac n 2) \]

同时我们有

\[F(w^{k+\frac n2}_n)=F_1(w^{2k+n}_n)+w^{k+\frac n2}_nF_2(w^{2k+n}_n) \]

\[=F_1(w^{k}_\frac n 2)+w^{k+\frac n2}_nF_2(w^{k}_\frac n 2) \]

发现 \(F(w^k_n)\) 与 \(F(w^{k+\frac n2}_n)\) 十分接近，求出 \(F_1(w^{k}_\frac n 2),F_2(w^{k}_\frac n 2)\) 之后能求出两个值。

让我们分析一下复杂度，我们要去求 \(n\) 个点，每个规模是 \(O(n)\) 的，本来是 \(O(n^2)\)，但是我们只要求出一半的点就能推出另一半的点，不断递归，总复杂度是 \(n\log_2 n\) 的，常数挺大。

快速傅里叶逆变换

还没有结束，因为我们最终还是要系数的，这一堆点值怎么转回系数捏？由于其推导过程比较复杂，直接给出结论：

设上面我们求出的点值 \(b_k=\sum_{i=0}^{n-1}a_iw_{n}^{ki}\)，那么我们有 \(a_k=\frac1n\sum_{i=0}^{n-1}a_iw^{-ki}_n\)。

式子还是比较简单的，在逆变换时我们只要将原本的 \(w_n^i\) 变成 \(w_{n}^{-i}\) 计算，最后再乘一个 \(\frac1n\) 就行。

代码实现 1

由此我们可以初步实现 FFT。

复数的运算我们可以用 B++ 自带的 complex。

#include<bits/stdc++.h>
#define int long long
#define cpd complex<double>
using namespace std;
const int N=3e6+5;
const double pi=acos(-1.0);
void FFT(cpd *f,int n,int t){
    if(n==1)return;int m=n>>1;
    cpd a[m],b[m];
    for(int i=0;i<m;++i)
        a[i]=f[i<<1],b[i]=f[i<<1|1];
    FFT(a,m,t),FFT(b,m,t);
    cpd w(cos(2*pi/n),sin(2*pi/n)*t),z(1,0);
    for(int i=0;i<m;++i,z*=w)
        f[i]=a[i]+z*b[i],
        f[i+m]=a[i]-z*b[i];
}int n,m,k=1;
cpd f[N],g[N];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m;
    for(int i=0,x;i<=n;++i)
        cin>>x,f[i].real(x);
    for(int i=0,x;i<=m;++i)
        cin>>x,g[i].real(x);
    while(k<=n+m)k<<=1;
    FFT(f,k,1),FFT(g,k,1);
    for(int i=0;i<k;++i)
        f[i]*=g[i];
    FFT(f,k,-1);
    for(int i=0;i<=n+m;++i)
        cout<<(int)round(f[i].real()/k)<<' ';
    return 0;
}

蝴蝶优化

首先递归常数并不好，其次复数运算常数也大，并且每次我们都开了一些数组，总之跑起来很不舒服。

如果想要解决递归，我们就不能每次重开数组，我们要正确地交换一些数的位置以实现多次乘法。根据上面的做法，这显然是一个位运算的问题。

考虑我们递归的下标。

0 1 2 3 4 5 6 7
0 2 4 6 1 3 5 7
0 4 2 6 1 5 3 7
0 4 2 6 1 5 3 7

最后一层的二进制为：

000 100 010 110 001 101 011 111

看起来没有规律，但要是把每个二进制数 reverse 一下呢？

000 001 010 011 100 101 110 111

相当于我们知道了最后一层数组是按二进制下翻转从小到大排列出来，现在我们自下而上归并排序，形成我们处理的下标。现在我们只需要初始做一个排序，然后就可以自下而上非递归处理问题。可以用 sort，但是这也是常数。

我们希望知道 \(i\) 在这种特殊比较方式下的排名，考虑 \(i\) 的最后一位，如果为 \(1\) 那么就会有 \(\frac n2\) 个数在其之前，剩下的就是 i>>1 这个数的排名了。

也不知道根蝴蝶有什么关系。

代码实现 2

#include<bits/stdc++.h>
#define int long long
#define cpd complex<double>
using namespace std;
const int N=3e6+5;
const double pi=acos(-1.0);
int r[N];
void FFT(cpd *a,int n,int t){
    for(int i=0;i<n;++i)
        if(r[i]>i)swap(a[i],a[r[i]]);
    for(int d=1;d<n;d<<=1){
        cpd w(cos(pi/d),sin(pi/d)*t);
        for(int i=0;i<n;i+=d<<1){
            cpd z(1,0),x,y;
            for(int j=0;j<d;++j,z*=w)
                x=a[i+j],y=a[i+j+d]*z,
                a[i+j]=x+y,a[i+j+d]=x-y;
        }
    }
}int n,m,k=1;
cpd f[N],g[N];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m;
    for(int i=0,x;i<=n;++i)
        cin>>x,f[i].real(x);
    for(int i=0,x;i<=m;++i)
        cin>>x,g[i].real(x);
    while(k<=n+m)k<<=1;
    for(int i=0;i<=k;++i)
        r[i]=(r[i>>1]>>1)|(i&1?k>>1:0);
    FFT(f,k,1),FFT(g,k,1);
    for(int i=0;i<k;++i)
        f[i]*=g[i];
    FFT(f,k,-1);
    for(int i=0;i<=n+m;++i)
        cout<<(int)round(f[i].real()/k)<<' ';
    return 0;
}

再优化

我们一共要做三遍 FFT，其实只用两遍。

考虑将 \(G\) 放到 \(F\) 的虚部上，然后求 \(F^2\)，将其虚部取出来除以 \(2\) 就是答案，原因是：

\[(a+bi)^2=(a^2-b^2)+2abi \]

#include<bits/stdc++.h>
#define int long long
#define cpd complex<double>
using namespace std;
const int N=3e6+5;
const double pi=acos(-1.0);
int r[N];
void FFT(cpd *a,int n,int t){
    for(int i=0;i<n;++i)
        if(r[i]>i)swap(a[i],a[r[i]]);
    for(int d=1;d<n;d<<=1){
        cpd w(cos(pi/d),sin(pi/d)*t);
        for(int i=0;i<n;i+=d<<1){
            cpd z(1,0),x,y;
            for(int j=0;j<d;++j,z*=w)
                x=a[i+j],y=a[i+j+d]*z,
                a[i+j]=x+y,a[i+j+d]=x-y;
        }
    }
}int n,m,k=1;
cpd f[N];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n>>m;
    for(int i=0,x;i<=n;++i)
        cin>>x,f[i].real(x);
    for(int i=0,x;i<=m;++i)
        cin>>x,f[i].imag(x);
    while(k<=max(n,m)*2)k<<=1;
    for(int i=0;i<=k;++i)
        r[i]=(r[i>>1]>>1)|(i&1?k>>1:0);
    FFT(f,k,1);
    for(int i=0;i<k;++i)
        f[i]*=f[i];
    FFT(f,k,-1);
    for(int i=0;i<=n+m;++i)
        cout<<(int)round(f[i].imag()/k/2)<<' ';
    return 0;
}

NTT

在一些模意义下，单位根是可以用原根代替的，例如模数为 \(998244353\) 时有原根 \(3\)，可以理解为 \(w_{998244352}^1=3\)，并且 \(998244352=2^{23}\times119\)，这意味着对于 FFT，我们可以轻松的求出 \(w_{998244352}^{\frac{998244352}{2^k}}\)，只要 \(k\le 23\)，对于多项式乘法已经足够。

这样的好处是我们不用再使用复数计算。

类似的模数有 \(1004535809=2^{21}\times 479+1,469762049=2^{26}\times 7+1\)。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=4e6+5,mod=998244353,ivg=(mod+1)/3;
inline int qp(int x,int y){
	int ans=1;
	for(;y;y>>=1,x=x*x%mod)
		if(y&1)ans=ans*x%mod;
	return ans;
}
int r[N];
inline void FFT(int *a,int L,int type){
	for(int i=0;i<L;++i)
		if(r[i]<i)swap(a[i],a[r[i]]);
	for(int d=2;d<=L;d<<=1){
		int wn=qp(3,(mod-1)/d),k=d>>1;
		for(int i=0;i<L;i+=d){
			for(int j=0,w=1;j<k;++j,w=w*wn%mod){
				int tmp=a[i+j+k]*w%mod;
				a[i+j+k]=(a[i+j]-tmp+mod)%mod;
				a[i+j]=(a[i+j]+tmp)%mod;
			}
		}
	}
	if(type==-1){
		reverse(a+1,a+L);int inv=qp(L,mod-2);
		for(int i=0;i<L;++i)(a[i]*=inv)%=mod;
	}
}
int a[N],b[N];
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
	int n,m,L=1;
	cin>>n>>m;
	for(int i=0;i<=n;++i)cin>>a[i];
	for(int i=0;i<=m;++i)cin>>b[i];
	while(L<=n+m)L<<=1;
	for(int i=0;i<L;++i)r[i]=(i&1)*(L>>1)+(r[i>>1]>>1);
	FFT(a,L,1),FFT(b,L,1);
	for(int i=0;i<L;++i)(a[i]*=b[i])%=mod;
	FFT(a,L,-1);
	for(int i=0;i<=n+m;++i)
		cout<<a[i]<<' ';
    return 0;
}