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

多项式乘法(NTT)

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

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

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

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

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

前人们发现原根$g$也有相似的优良性质
简单介绍原根和阶
对于一个质数$p$,和一个数a，a在mod p意义下的阶就是

\[a^r\equiv1(mod\;p) \]

的最小整数$r$
一个重要的性质是
$a^0,a^1,a^2……，a^{r-1}$互不相同
证明的话假设有

\[a^x\equiv a^y(mod\;p) \]

那么

\[a^{|x-y|} \equiv 1(mod\;p) \]

显然$|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{n}})^0=1$$
满足

2:$$g_n^{i+j}=(g{n}})^{i+j}=(g{n}})^i(g{n}})^j=g_{n}ig_n^j $$

3:$$g_n^{k+n/2}=(g{n}})^k(g{n}})^{n/2}=g{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}¹⁾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(mod\; 98244353)$
由于题目良心，所以模数正好是NTT模数
我们考虑用倍增的方法求
下文都用$F_*(x)$表示次数是$F(x)$一半的多项式，其余一样
那么我们假设已经倍增的求出$G_*(x)$，要求$G(x)$
那么显然

\[G(x)\equiv G_*(x)(mod\;x^{n/2}) \]

移项

\[G(x)-G_*(x)\equiv 0(mod\;x^{n/2}) \]

两边同时平方

\[(G(x)-G_*(x))^2\equiv 0(mod\;x^n) \]

\[G(x)^2-2G(x)G_*(x)+G_*(x)^2\equiv 0(mod\;x^n) \]

同时乘上$F(x)$

\[G(x)-2G_*(x)+G_*(x)^2F(x)\equiv 0(mod\;x^n) \]

再次移项

\[G(x)\equiv 2G_*(x)-G_*(x)^2F(x)(mod\;x^n) \]

\[G(x)\equiv G_*(x)(2-G_*(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'(x)$

\[F(x+dx)=\sum_{i=0}^nf_i(x+dx)^i \]

代入二项式定理

\[F(x+dx)=\sum_{i=0}^nf_i(C_i^0x^i+C_i^1x^{i-1}dx+C_i^2x^{i-2}dx^2……) \]

\[F(x+dx)-F(x)=\sum_{i=0}^nf_i(C_i^1x^{i-1}dx+C_i^2x^{i-2}dx^2……) \]

代入导数定义式

\[F'(x)=\frac{F(x+dx)-F(x)}{dx}=\sum_{i=0}^nf_iix^{i-1} \]

也就是

\[F'(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_*(x)-\frac{F(G_*(x))}{F'(G_*(x))} \]

多项式ln

已知函数$F(x)$,求$G(x)=ln(F(x))$

\[G(x)=ln(F(x)) \]

两遍同时求导

\[G'(x)=(ln(F(x)))' \]

根据复合函数求导的链式法则

\[(F(G(x)))'=F'(G(x))G'(x) \]

而我们知道

\[ln'(x)=\frac{1}{x} \]

右边的式子也就是

\[G'(x)=\frac{F'(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_*(x)-\frac{H(G_*(x))}{H'(G_*(x))} \]

因为$F(x)$是常数，所以不会对导数有贡献
$H'(G(x))=\frac{1}{G(x)}$
回带

\[G(x)=G_*(x)-G_*(x)H(G_*(x)) \]

\[G(x)=G_*(x)-G_*(x)(ln(G_*(x)-F(x))) \]

\[G(x)=G_*(x)(1-ln(G_*(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'(G(x))=2G(x)$
代入牛顿迭代

\[G(x)=G_*(x)-\frac{H(G_*(x))}{H'(G_*(x))} \]

\[G(x)=G_*(x)-\frac{G_*(x)^2-F(x)}{2G_*(x)} \]

\[G(x)=\frac{G_*(x)^2+F(x)}{2G_*(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),2));
	} 
	turn(g,len(f));
	return g;
}

多项式带余除法

已知$F(x)$,$G(x)$,求$H(x),R(x)$满足，$F(x)=H(x)G(x)+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;	
}

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;
}

posted @ 2022-02-09 08:26 Larunatrecy 阅读(58) 评论(0) 收藏举报

刷新页面返回顶部

Loading

月华如弦~天明将近~~

我总是在寻找着月光，却总是无法如愿以偿。