进军多项式（一）：多项式全家桶模板

多项式乘法(NTT)

默认大家会了FFT……
我们知道FFT中因为要用三角函数和复数，所以对精度要求高
并且常数比较大
那么有没有可以代替单位根的呢？
我们思考一下在FFT中应用的单位根的性质
1:$w_n^0=1$
保证了有初始值

2:$w_n^{i+k}=w_n^i{w}_n^j$,也就是$w_n^k=(w_n^1{)}^k$
因此只要求出一个就可以递推

3:$w_n^{k+n/2}=-w_n^k$
保证了计算的简洁性

4:$w_n^{0,1\dots \dots n-1}$互不相同
保证了代入了不同的数

5:$(w_{2n}^1{)}^2=w_n^1$
保证了计算式可以分治

前人们发现原根$g$也有相似的优良性质
简单介绍原根和阶
对于一个质数$p$,和一个数a，a在mod p意义下的阶就是
$$a^r\equiv 1(modp)$$
的最小整数$r$
一个重要的性质是
$a^0,a^1,a^2\dots \dots ，a^{r-1}$互不相同
证明的话假设有
$$a^x\equiv a^y(modp)$$
那么

$$a^{|x-y|}\equiv 1(modp)$$
显然$|x-y|<r$,但是根据阶的定义，r是最小的满足这个式子的数，矛盾
那么是原根呢？
一个质数的原根g是满足下列条件的数的集合

$g$在模p意义下的阶是$p-1$

当然一个质数不只有一个原根，但总有一个最小原根

根据前人的总结，发现若定义$g_n^1=g^{\frac{(p-1)}{n}}$
那么这个东西有和$w_n^1$相似的性质
我们一一来看

1:$$g_n^0=(g^{\frac{(p-1)}{n}}{)}^0=1$$
满足

2:$$g_n^{i+j}=(g^{\frac{(p-1)}{n}}{)}^{i+j}=(g^{\frac{(p-1)}{n}}{)}^i(g^{\frac{(p-1)}{n}}{)}^j=g_n^i{g}_n^j$$

3:$$g_n^{k+n/2}=(g^{\frac{(p-1)}{n}}{)}^k(g^{\frac{(p-1)}{n}}{)}^{n/2}=g^{\frac{k(p-1)}{n}}{g}^{\frac{(p-1)}{2}}$$
而我们知道$g^{p-1}\equiv 1$,因此$g^{\frac{(p-1)}{2}}\equiv +1或-1$
又因为上述阶的性质，所有g的次幂互不相同，因此$g^{\frac{(p-1)}{2}}\equiv -1$
代入上式
$$g_n^{k+n/2}=g^{\frac{k(p-1)}{n}}{g}^{\frac{(p-1)}{2}}=-g^{\frac{k(p-1)}{n}}=-g_n^k$$

4:由阶的性质可得

5：$$(g_{2n}^1{)}^2=(g^{\frac{(p-1)}{2n}}{)}^2=g^{\frac{(p-1)}{n}}=g_n^1$$
至此我们发现w需要的性质g都有
问题解决啦
……
吗？
上述式子中有一个问题就是$g^{\frac{(p-1)}{n}}$中$(p-1)$可能不是n的倍数
如果您了解FFT，就会知道这里的$n$是2的整次幂
因此，一个模数合法要求(p-1)的因子中含有大量2的幂
这样的模数不算特别多，一个最常见的就是998244353
因此，我们只需要代入模数，就可以仿照FFT的方式啦
这就是NTT

void NTT(Poly &f,int opt,int n)
{
	for(int i=0;i<n;i++)
	if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(int len=2;len<=n;len<<=1)
	{
		int l=len/2;
		int w=Pow(opt?G:IG,(mod-1)/len);
		for(int i=0;i<n;i+=len)
		{
			LL g=1;
			for(int k=i;k<i+l;k++)
			{
				LL v=g*f[k+l]%mod;
				f[k+l]=(f[k]-v+mod)%mod;
				f[k]=(f[k]+v)%mod;
				g=1ll*g*w%mod;
			}
		} 
	}
}

然后是多项式乘法，和FFT类似的

#define Poly vector<int>
#define len(x) (int)x.size()
#define turn(x,n) x.resize(n)
void Retry(int n)
{
	for(int i=0;i<n;i++)
	tr[i]=((tr[i>>1]>>1)|((i&1)?(n>>1):0));
}
Poly Mul(Poly f,Poly g)
{
	int n=max(len(f),len(g));
	turn(f,n);
	for(int i=0;i<n;i++)
	f[i]=1ll*f[i]*g[i]%mod;
	return f;
}
Poly Times(Poly f,Poly g)
{
	int n=len(f),m=len(g);
	for(m+=n-1,n=1;n<m;n<<=1);
	LL invn=Inv(n);
	turn(f,n);turn(g,n);
	Retry(n);
	NTT(f,1,n);NTT(g,1,n);
	f=Mul(f,g);
	NTT(f,0,n);f=Mul(f,invn);
	turn(f,min(m,mxlen));
	return f;
}

多项式求逆

已知多项式$F(x)$，求$G(x)$满足$G(x)F(x)\equiv 1(mod98244353)$
由于题目良心，所以模数正好是NTT模数
我们考虑用倍增的方法求
下文都用$F_{\ast }(x)$表示次数是$F(x)$一半的多项式，其余一样
那么我们假设已经倍增的求出$G_{\ast }(x)$，要求$G(x)$
那么显然
$$G(x)\equiv G_{\ast }(x)(mod{x}^{n/2})$$
移项
$$G(x)-G_{\ast }(x)\equiv 0(mod{x}^{n/2})$$
两边同时平方
$$(G(x)-G_{\ast }(x){)}^2\equiv 0(mod{x}^n)$$
$$G(x{)}^2-2G(x)G_{\ast }(x)+G_{\ast }(x{)}^2\equiv 0(mod{x}^n)$$
同时乘上$F(x)$
$$G(x)-2G_{\ast }(x)+G_{\ast }(x{)}^{2}F(x)\equiv 0(mod{x}^n)$$
再次移项
$$G(x)\equiv 2G_{\ast }(x)-G_{\ast }(x{)}^{2}F(x)(mod{x}^n)$$
$$G(x)\equiv G_{\ast }(x)(2-G_{\ast }(x)F(x))(mod{x}^n)$$
倍增计算，复杂度$O(n\log n)$

Poly Inv(Poly f)
{
	Poly g=Poly(1,Inv(f[0]));
	for(mxlen=2;mxlen<2*len(f);mxlen<<=1)
	{
		Poly h=f;
		turn(h,mxlen);
		g=Times(g,Sub(2,Times(g,h)));
	}
	turn(g,len(f));
	return g;
}

多项式求导

考虑
已知$F(x)$，求导数$F^{\prime }(x)$
$$F(x+dx)=\sum _{i=0}^n{f}_i(x+dx{)}^i$$
代入二项式定理
$$F(x+dx)=\sum _{i=0}^n{f}_i(C_i^0{x}^i+C_i^1{x}^{i-1}dx+C_i^2{x}^{i-2}d{x}^2\dots \dots )$$
$$F(x+dx)-F(x)=\sum _{i=0}^n{f}_i(C_i^1{x}^{i-1}dx+C_i^2{x}^{i-2}d{x}^2\dots \dots )$$
代入导数定义式
$$F^{\prime }(x)=\frac{F(x+dx)-F(x)}{dx}=\sum _{i=0}^n{f}_{i}i{x}^{i-1}$$
也就是
$$F^{\prime }(x)=\sum _{i=0}^{n-1}{f}_{i+1}(i+1)x^i$$

Poly Der(Poly f)
{
	for(int i=0;i<len(f)-1;i++)
	f[i]=Mul(i+1,f[i+1]);
	turn(f,len(f)-1);
	return f;
}

多项式积分

把求导逆过来就好了

Poly Int(Poly f)
{
	turn(f,len(f)+1);
	for(int i=len(f)-1;i>=1;i--)
	f[i]=Mul(f[i-1],inv[i]);
	f[0]=0;
	return f;
}

多项式牛顿迭代

已知函数$F(x)$,求函数$G(x)$使得$F(G(x))=0$
直接给出公式吧
$$G(x)=G_{\ast }(x)-\frac{F(G_{\ast }(x))}{F^{\prime }(G_{\ast }(x))}$$

多项式ln

已知函数$F(x)$,求$G(x)=ln(F(x))$
$$G(x)=ln(F(x))$$
两遍同时求导
$$G^{\prime }(x)=(ln(F(x)){)}^{\prime }$$
根据复合函数求导的链式法则
$$(F(G(x)){)}^{\prime }=F^{\prime }(G(x))G^{\prime }(x)$$
而我们知道
$$l{n}^{\prime }(x)=\frac{1}{x}$$
右边的式子也就是
$$G^{\prime }(x)=\frac{F^{\prime }(x)}{F(x)}$$
这个时候就可以算了
然后再积分回去就可以得到原函数了

Poly Ln(Poly f)
{
	return Int(Times(Der(f),Inv(f)));
}

多项式exp

已知函数$F(x)$,求$G(x)=e^{F(x)}$
两边同时取$ln$
$$lnG(x)=F(x)$$
移项得
$$lnG(x)-F(x)=0$$
设$H(G(x))=lnG(x)-F(x)$
也就是
$$H(G(x))=0$$
很符合牛顿迭代的式子
代入公式
$$G(x)=G_{\ast }(x)-\frac{H(G_{\ast }(x))}{H^{\prime }(G_{\ast }(x))}$$
因为$F(x)$是常数，所以不会对导数有贡献
$H^{\prime }(G(x))=\frac{1}{G(x)}$
回带
$$G(x)=G_{\ast }(x)-G_{\ast }(x)H(G_{\ast }(x))$$
$$G(x)=G_{\ast }(x)-G_{\ast }(x)(ln(G_{\ast }(x))-F(x))$$
$$G(x)=G_{\ast }(x)(1-ln(G_{\ast }(x))+F(x))$$
倍增求解就行啦

Poly Exp(Poly f)
{
	Poly g=Poly(1,1);
	for(mxlen=2;mxlen<4*len(f);mxlen<<=1)
	{
		Poly A=f;
		turn(A,mxlen);
		g=Times(g,Add(Sub(1,Ln(g)),A));
	}
	turn(g,len(f));
	return g;
}

多项式快速幂

已知多项式$F(x)$,求$G(x)=F(x{)}^k$
我们知道
$$e^{ln(x)}=x$$
所以
$$x^k=(e^{ln(x)}{)}^k$$
$$x^k=e^{kln(x)})$$
我们先对$F(x)$取$ln$,然后乘以k，再exp回去就行啦
这里我们也可以知道这个k完全可以是mod p同余的，因此即使指数很大也没事啦

Poly Pow(Poly f,int k)
{
	f=Ln(f);
	f=Mul(f,k);
	f=Exp(f);
	return f;
}

多项式开根

已知函数$F(x)$,求$G(x{)}^2=F(x)$
方法一：调用多项式快速幂，不过常数巨大
方法二：使用牛顿迭代
因为
$$G(x{)}^2=F(x)$$
$$G(x{)}^2-F(x)=0$$
设$H(G(x))=G(x{)}^2-F(x)$
求导得$H^{\prime }(G(x))=2G(x)$
代入牛顿迭代
$$G(x)=G_{\ast }(x)-\frac{H(G_{\ast }(x))}{H^{\prime }(G_{\ast }(x))}$$
$$G(x)=G_{\ast }(x)-\frac{G_{\ast }(x{)}^2-F(x)}{2G_{\ast }(x)}$$
$$G(x)=\frac{G_{\ast }(x{)}^2+F(x)}{2G_{\ast }(x)}$$
倍增求解即可

Poly Sqrt(Poly f)
{
	Poly g=Poly(1,1);
	for(mxlen=2;mxlen<4*len(f);mxlen<<=1)
	{
		Poly A=f;turn(A,mxlen);
		g=Times(Add(Times(g,g),A),Mul(Inv(g),inv[2]));
	} 
	turn(g,len(f));
	return g;
}

多项式带余除法

已知$F(x)$,$G(x)$,求$H(x),R(x)$满足，$F(x)=H(x)G(x)+R(x)$
设
$F_r(x)=x^{n}F(\frac{1}{x})$
我们看这个$F_r(x)$究竟是什么
$$F_r(x)=x^{n}F(x^{-1})=\sum _{i=0}^n{f}_i{x}^n{x}^{-i}=\sum _{i=0}^n{f}_i{x}^{n-i}=\sum _{i=0}^n{f}_{n-i}{x}^i$$
也就是将$F(x)$系数翻转
我们观察式子
$$F(x)=H(x)G(x)+R(x)$$
带入$\frac{1}{x}$
$$F(\frac{1}{x})=H(\frac{1}{x})G(\frac{1}{x})+R(\frac{1}{x})$$
同乘$x^n$
$$x^{n}F(\frac{1}{x})=x^{n-m}H(\frac{1}{x})x^{m}G(\frac{1}{x})+x^{n-m+1}{x}^{m-1}R(\frac{1}{x})$$
$$F_r(x)=H_r(x)G_r(x)+x^{n-m+1}{R}_r(x)$$
同模$x^{n+m-1}$
$$F_r(x)\equiv H_r(x)G_r(x)$$
$$H_r(x)\equiv \frac{G_r(x)}{F_r(x)}$$
多项式求逆算一下，再翻转系数

Poly Divide(Poly f,Poly g,int n,int m)
{
	Poly Tf,Tg;
	turn(Tf,n+1);
	turn(Tg,n+1);
	for(int i=0;i<=n;i++)
	Tf[n-i]=f[i];
	for(int i=0;i<=m;i++)
	Tg[m-i]=g[i];
	for(int i=n-m+2;i<=m;i++) Tg[i]=0;
	Tg=Inv(Tg);
	Tf=Times(Tf,Tg);
	Reverse(Tf,n-m+1);
	for(int i=n-m+1;i<=n;i++)
	Tf[i]=0;
	return Tf;	
}

$R(x)=F(x)-G(x)H(x)$

Poly Mod(Poly f,Poly g,int n,int m)
{
	Poly h=Divide(f,g,n,m);
	Poly res=Sub(f,Times(g,h));
	return res; 
}

完整代码

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
#define int long long 
const int mod = 998244353;
const int N = 1e6+7;
int Add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int Sub(int x,int y){return x-y<0?x-y+mod:x-y;}
int Mul(int x,int y){return 1ll*x*y%mod;}
int Pow(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return res;
}
int Inv(int x){return Pow(x,mod-2);}
int inv[N];
void Init()
{
	inv[1]=1;
	for(int i=2;i<N;i++)
	inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
#define Poly vector<int>
#define len(x) (int)x.size()
#define turn(x,n) x.resize(n)
Poly Sub(int x,Poly f)
{
	for(int i=0;i<len(f);i++)
	f[i]=mod-f[i];
	f[0]=Add(f[0],x);
	return f;
}
Poly Add(Poly f,int x)
{
	f[0]=Add(f[0],x);
	return f;
}
Poly Sub(Poly f,int x)
{
	f[0]=Sub(f[0],x);
	return f;
}
Poly Mul(Poly f,int x)
{
	for(int i=0;i<len(f);i++)
	f[i]=Mul(f[i],x);
	return f;
}
Poly Add(Poly f,Poly g)
{
	int n=max(len(f),len(g));
	turn(f,n);
	for(int i=0;i<n;i++)
	f[i]=Add(f[i],g[i]);
	return f; 
}
Poly Sub(Poly f,Poly g)
{
	int n=max(len(f),len(g));
	turn(f,n);
	for(int i=0;i<n;i++)
	f[i]=Sub(f[i],g[i]);
	return f;
}
const int G=3,IG=Pow(G,mod-2);
int tr[N];
void Retry(int n)
{
	for(int i=0;i<n;i++)
	tr[i]=((tr[i>>1]>>1)|((i&1)?(n>>1):0));
}
void NTT(Poly &f,int opt,int n)
{
	for(int i=0;i<n;i++)
	if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(int len=2;len<=n;len<<=1)
	{
		int l=len/2;
		int w=Pow(opt?G:IG,(mod-1)/len);
		for(int i=0;i<n;i+=len)
		{
			LL g=1;
			for(int k=i;k<i+l;k++)
			{
				LL v=g*f[k+l]%mod;
				f[k+l]=(f[k]-v+mod)%mod;
				f[k]=(f[k]+v)%mod;
				g=1ll*g*w%mod;
			}
		} 
	}
}
Poly Mul(Poly f,Poly g)
{
	int n=max(len(f),len(g));
	turn(f,n);
	for(int i=0;i<n;i++)
	f[i]=1ll*f[i]*g[i]%mod;
	return f;
}
int mxlen=100000000;
Poly Times(Poly f,Poly g)
{
	int n=len(f),m=len(g);
	for(m+=n-1,n=1;n<m;n<<=1);
	LL invn=Inv(n);
	turn(f,n);turn(g,n);
	Retry(n);
	NTT(f,1,n);NTT(g,1,n);
	f=Mul(f,g);
	NTT(f,0,n);f=Mul(f,invn);
	turn(f,min(m,mxlen));
	return f;
}
Poly Inv(Poly f)
{
	Poly g=Poly(1,Inv(f[0]));
	for(mxlen=2;mxlen<2*len(f);mxlen<<=1)
	{
		Poly h=f;
		turn(h,mxlen);
		g=Times(g,Sub(2,Times(g,h)));
	}
	turn(g,len(f));
	return g;
}
Poly Der(Poly f)
{
	for(int i=0;i<len(f)-1;i++)
	f[i]=Mul(i+1,f[i+1]);
	turn(f,len(f)-1);
	return f;
}
Poly Int(Poly f)
{
	turn(f,len(f)+1);
	for(int i=len(f)-1;i>=1;i--)
	f[i]=Mul(f[i-1],inv[i]);
	f[0]=0;
	return f;
}
Poly Ln(Poly f)
{
	return Int(Times(Der(f),Inv(f)));
}
Poly Exp(Poly f)
{
	Poly g=Poly(1,1);
	for(mxlen=2;mxlen<4*len(f);mxlen<<=1)
	{
		Poly A=f;
		turn(A,mxlen);
		g=Times(g,Add(Sub(1,Ln(g)),A));
	}
	turn(g,len(f));
	return g;
}
Poly Pow(Poly f,int k)
{
	f=Ln(f);
	f=Mul(f,k);
	f=Exp(f);
	return f;
}
Poly Sqrt(Poly f)
{
	Poly g=Poly(1,1);
	for(mxlen=2;mxlen<4*len(f);mxlen<<=1)
	{
		Poly A=f;turn(A,mxlen);
		g=Times(Add(Times(g,g),A),Mul(Inv(g),2));
	} 
	turn(g,len(f));
	return g;
}
void Put(Poly f,int n)
{
	for(int i=0;i<n;i++)
	printf("%lld ",f[i]);
	printf("\n");
}
char s[N];
LL Get()
{
	LL num=0;
	scanf("%s",s+1);
	LL m=strlen(s+1);
	for(int i=1;i<=m;i++)
	{
		num=num*10ll%mod;
		num=(num+s[i]-'0')%mod;
	}
	return num;
}
void Reverse(Poly &f,int len)
{
	Poly g=f;
	for(int i=0;i<len;i++)
	f[i]=g[len-1-i];
}
Poly Divide(Poly f,Poly g,int n,int m)
{
	Poly Tf,Tg;
	turn(Tf,n+1);
	turn(Tg,n+1);
	for(int i=0;i<=n;i++)
	Tf[n-i]=f[i];
	for(int i=0;i<=m;i++)
	Tg[m-i]=g[i];
	for(int i=n-m+2;i<=m;i++) Tg[i]=0;
	Tg=Inv(Tg);
	Tf=Times(Tf,Tg);
	Reverse(Tf,n-m+1);
	for(int i=n-m+1;i<=n;i++)
	Tf[i]=0;
	return Tf;	
}
Poly Mod(Poly f,Poly g,int n,int m)
{
	Poly h=Divide(f,g,n,m);
	Poly res=Sub(f,Times(g,h));
	return res; 
}
signed main()
{
	return 0;
}