浅谈DFT

FFT相关

单位根

定义使得\(\omega^n=1\)为\(n\)的单位根

如果设\(\omega=\cos\theta i+\sin\theta\)

\(则\omega^n=(\cos n\theta i+\sin n\theta)\),就相当于在平面内为\(n\theta\)的角

即\(\omega=\cos\frac{2k\pi}{n}i+\sin\frac{2k\pi}{n}=e^{i\frac{2k\pi}{n}}\)

由此,\(\omega^\frac{n}{2}=e^{ik\pi}=-1\)

设\(\omega_n\)为\(n\)次单位根​,则\(\omega_n^{2k}=\omega_{n/2}^k\)

---

FFT

对于一个多项式\(f(x)\),\(FFT\)就是求出当\(x=\omega_n^0,\omega_n^1,\omega_n^2,\omega_n^3...\omega_n^{n-1}\)时\(f(x)\)的取值

考虑对\(f(x)\)进行的系数奇偶分类,令\(g(x)=a_0+a_2x+a_4x^2....\),\(h(x)=a_1+a_3x+a_5x^2....\)

则\(f(x)=g(x^2)+xh(x^2)\)

对于\(\omega^k_n,k\in[0,\dfrac{n}{2})\)

\(f(\omega^k_n)=g(\omega^{2k}_n)+\omega^k_nh(\omega^{2k}_{n})=g(\omega^{k}_{n/2})+\omega^k_nh(\omega^{k}_{n/2})\)

注意到\(k\in[0,\dfrac{n}{2})\)，则递归求解\(g,h\)时可以发现问题的范围缩小了一半

对于\(\omega^{k+\frac{2}{n}}_n,k\in[0,\dfrac{n}{2})\)

\(f(\omega^{k+\frac{2}{n}}_n)=f(-\omega^k_n)=g(\omega^{k}_{n/2})-\omega^k_nh(\omega^{k}_{n/2})\)

可以发现对于\(k\in[\frac{n}{2},n)\)，同样可以用\(h,g\)求解，同时范围缩小了一半

因此递归求解的时间复杂度为\(\Theta(n\log_2n)\)

对于\(IFFT\),根据\(IDFT\)与\(DFT\)的关系，可得就是令\(x=\frac{1}{\omega}\),同时\(f(x)=\dfrac{f(x)}{n}\)

如果将递归时的\(g,h\)描述为一棵树，其叶子节点即为原来的数二进制的反转

即位逆序置换，其实就是原来的那个序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标

具体在实现的时候可以用递推求出每个数对应的反转将递归转化为倍增

如果我们在计算时不具体表示\(g,h\)，而是用\(f\)表示，这个被称为蝶形运算优化

结合这两个举个例子

对\([0,1,2,3,4,5,6,7]\)划分,把\(h\)放左边,\(g\)放右边

\([0,1,2,3,4,5,6,7]\)

\([0,2,4,6,1,3,5,7]\)

\([0,4,2,6,1,5,3,7]\)

从第一层开始

\(f_1(0)=g(0)+\omega_n^0 h(0)=f_2(0)+\omega _n^0f_2(1),f_1(4)=g(0)-\omega_n^0 h(0)=f_2(0)-\omega _n^0f_2(1)\)

\(f_1(1)=g(1)+\omega_n^1 h(1)=f_2(2)+\omega _n^0f_2(3),f_1(5)=g(1)-\omega_n^1 h(1)=f_2(2)-\omega _n^0f_2(3)\)

注意到规律,设\(Rev_i\)为上面序列的位置\(i\)的值，\(Len\)为当前分治长度

则\(f_1(k)=f_2(k)+\omega_n^kf_2(k+Len),k\in[0,\dfrac{n}{2}),f_1(k)=f_2(k)-\omega_n^kf_2(k+Len),k\in[\frac{n}{2},n)\)

因此\(n\)为\(2\)的幂

NTT

考虑在模意义下的单位根,\(n\)依旧为\(2\)的幂,\(p=qn+1\)

\(\omega_n^k(k\in[0,n)\)互不相同

\(\omega_{2n}^{2k}=\omega_n^k\)

\(\omega_{n}^{\frac{n}{2}}=-1\)

\(\omega_n^n =1\)

考虑\(\omega_n=g^q\),其中\(g\)为\(p\)的原根

则根据原根性质\(g^{0\times q},g^{1\times q},g^{2\times q},g^{3\times q},g^{4\times q}...g^{nq-q+1}\)互不相同

\(\omega_n^n=g^{qn}=g^{p-1}=1(mod\ p)\),欧拉定理

\(\omega_{2n}^{2k}(p=2n\times \frac{q}{2}+1)=g^{2\times\frac{q}{2}}=\omega_n^k\)

\(\omega_{n}^{\frac{n}{2}}=g^{q\frac{n}{2}}=\sqrt{(\omega_n^n)}=\pm 1\),而因为互不相同，则\(\omega_{n}^{\frac{n}{2}}=-1\)

由此,\(w_n=g^{\frac{p-1}{n}}\)

至于如何求原根,可以在\(p^{\frac{1}{4}}\)的范围内枚举，使得对于\(\phi(p)\)的每一个因子\(x\)都不满足\(g^{\frac{\phi(p)}{x}}=1(mod\ p)\)

而\(998244353\)的常用原根为\(3\)

#include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
const int MAXN=3e6+5;
const double Pi=acos(-1.0);
const int MOD=998244353;
const int g=3;
struct Cpx{
	double a,b;
	Cpx(){
		a=0;
		b=0;
	} 
	Cpx(double aa,double bb){
		a=aa;
		b=bb;
	}
	Cpx operator*(const Cpx x)const{
		return Cpx(a*x.a-b*x.b,b*x.a+a*x.b);
	}
	Cpx operator+(const Cpx x)const{
		return Cpx(a+x.a,b+x.b);
	}
	Cpx operator-(const Cpx x)const{
		return Cpx(a-x.a,b-x.b);
	} 
};
int Pow(int a,int b,int p)
{
	int res=1;
	int base=a;
	while(b)
	{
		if(b&1)
		{
			res=((long long)res*base)%p;
		}
		base=((long long)base*base)%p;
		b>>=1;
	}
	return res;
}
int inv(int a,int p)
{
	return Pow(a,p-2,p);
}
int Rev[MAXN];
struct Poly{
	vector<Cpx>V;
	vector<int>U;
	void NTT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(U.size()<Len)
		{
			U.push_back(0);
		}
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(U[i],U[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			int Wn=Pow(g,(MOD-1)/(l<<1),MOD);
			if(type==-1)
			{
				Wn=inv(Wn,MOD);
			}
			for(int i=0;i<Len;i+=(l<<1))
			{
				int W=1;
				for(int j=i;j<i+l;j++,W=((long long)W*Wn)%MOD)
				{
					int Xc=U[j];
					int Yc=((long long)U[j+l]*W)%MOD;
					U[j]=((long long)Xc+Yc)%MOD;
					U[j+l]=((long long)Xc-Yc+MOD)%MOD;
				}
			}
		}
		if(type==-1)
		{
			int Liv=inv(Len,MOD); 
			for(int i=0;i<Len;i++)
			{
				U[i]=((long long)U[i]*Liv)%MOD;	
			}
		}
	}
	void FFT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(V.size()<Len)
		{
			V.push_back(Cpx(0,0));
		}
		
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(V[i],V[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			Cpx Wn=Cpx(cos(Pi/l),type*sin(Pi/l));
			for(int i=0;i<Len;i+=(l<<1))
			{
				Cpx W=Cpx(1,0);
				for(int j=i;j<i+l;j++,W=W*Wn)
				{
					Cpx Xc=V[j];
					Cpx Yc=V[j+l]*W;
					V[j]=(Xc+Yc);
					V[j+l]=(Xc-Yc);
				}
			}
		}
		if(type==-1)
		{
			for(int i=0;i<Len;i++)
			{
				V[i].a/=Len;	
			}
		}
	}
};
Poly Mul_FFT(Poly A,Poly B)
{
	int N=A.V.size();
	int M=B.V.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.FFT(Lm,1);
	 B.FFT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.V[i]=A.V[i]*B.V[i];
	 }
	 A.FFT(Lm,-1);
	 while(A.V.size()>(N+M-1))
	 {
	 	A.V.pop_back();
	 }
	 return A;
}
Poly Mul_NTT(Poly A,Poly B)
{
	int N=A.U.size();
	int M=B.U.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.NTT(Lm,1);
	 B.NTT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.U[i]=((long long)A.U[i]*B.U[i])%MOD;
	 }
	 A.NTT(Lm,-1);
	 while(A.U.size()>(N+M-1))
	 {
	 	A.U.pop_back();
	 }
	 return A;
}
int n,m,x;
int main()
{
	scanf("%d %d",&n,&m);
	Poly A,B;
	A.U.clear();
	B.U.clear();
	for(int i=0;i<=n;i++)
	{
		scanf("%d",&x);
		A.U.push_back(x);	
	} 
	for(int i=0;i<=m;i++)
	{
		scanf("%d",&x);
		B.U.push_back(x);	
	} 
	A=(Mul_NTT(A,B));
	for(int i=0;i<A.U.size();i++)
	{
		printf("%d ",A.U[i]);
	}
}

---

任意模数NTT

普通\(NTT\)要求\(p=2^kq+1\)，如\(998244353\),但实际情景中模数一般无法满足

一般系数的上界为\(p^2n\),如果\(p\)很小的话实际上可以用\(FFT\)计算后取模

但\(long\ double\)的精度在\(1e11\)后一般会丢失

考虑\(p'=\lceil\sqrt p \rceil\),则对于一个数\(X=p’A+B,(A<p',B<p')\)

考虑\(XY=(p'A_x+B_x)(p'A_y+B_y)=AxAy(p')^2+p'(A_yB_x+A_xB_y)+B_xB_y\)

如果只计算\(A_xA_y\)等，则实际的上界为\(pn\)

#include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
const int MAXN=3e5+5;
const long double Pi=acos(-1.0);
const int p=32000;
struct Cpx{
	long double a,b;
	Cpx(){
		a=0;
		b=0;
	} 
	Cpx(long double aa,long double bb){
		a=aa;
		b=bb;
	}
	Cpx operator*(const Cpx x)const{
		return Cpx(a*x.a-b*x.b,b*x.a+a*x.b);
	}
	Cpx operator+(const Cpx x)const{
		return Cpx(a+x.a,b+x.b);
	}
	Cpx operator-(const Cpx x)const{
		return Cpx(a-x.a,b-x.b);
	} 
};
int Rev[MAXN*4];
struct Poly{
	vector<Cpx>V;
	void FFT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(V.size()<Len)
		{
			V.push_back(Cpx(0,0));
		}
		
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(V[i],V[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			Cpx Wn=Cpx(cos(Pi/l),type*sin(Pi/l));
			for(int i=0;i<Len;i+=(l<<1))
			{
				Cpx W=Cpx(1,0);
				for(int j=i;j<i+l;j++,W=W*Wn)
				{
					Cpx Xc=V[j];
					Cpx Yc=V[j+l]*W;
					V[j]=(Xc+Yc);
					V[j+l]=(Xc-Yc);
				}
			}
		}
		if(type==-1)
		{
			for(int i=0;i<Len;i++)
			{
				V[i].a/=Len;	
			}
		}
	}
}G,F;
Poly Mul_FFT(Poly A,Poly B)
{
	int N=A.V.size();
	int M=B.V.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.FFT(Lm,1);
	 B.FFT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.V[i]=A.V[i]*B.V[i];
	 }
	 A.FFT(Lm,-1);
	 while(A.V.size()>(N+M-1))
	 {
	 	A.V.pop_back();
	 }
	 return A;
}
int n,m,MOD;
int x;
int main()
{
	scanf("%d %d %d",&n,&m,&MOD);
	Poly A1,B1,C1,D1;
	A1.V.clear();
	B1.V.clear();
	C1.V.clear();
	D1.V.clear();
	Poly A2,B2,C2,D2;
	A2.V.clear();
	B2.V.clear();
	C2.V.clear();
	D2.V.clear();
	for(int i=0;i<=n;i++)
	{
		scanf("%d",&x);
		int ax=x/p;
		int bx=x%p;
		A1.V.push_back(Cpx(ax,0));
		B1.V.push_back(Cpx(ax,0));
		C1.V.push_back(Cpx(bx,0));
		D1.V.push_back(Cpx(bx,0));
	}
	for(int i=0;i<=m;i++)
	{
		scanf("%d",&x);
		int ax=x/p;
		int bx=x%p;
		A2.V.push_back(Cpx(ax,0));
		B2.V.push_back(Cpx(bx,0));
		C2.V.push_back(Cpx(ax,0));
		D2.V.push_back(Cpx(bx,0));
	}
	A1=Mul_FFT(A1,A2);
	B1=Mul_FFT(B1,B2);
	C1=Mul_FFT(C1,C2);
	D1=Mul_FFT(D1,D2);
	for(int i=0;i<A1.V.size();i++)
	{
		long long F1=(long long)(A1.V[i].a+0.5);
		long long F2=(long long)(B1.V[i].a+0.5);
		long long F3=(long long)(C1.V[i].a+0.5);
		long long F4=(long long)(D1.V[i].a+0.5);
		long long F=(((F1*p)%MOD*p)%MOD+(((F2+F3)%MOD)*p)%MOD+F4)%MOD;
		printf("%lld ",F);
	}				
}

任意长度FFT

\(bluestein\)算法,用于解决任意长度的\(DFT\),常用于计算循环卷积

考虑\(DFT\)的式子

\(f(\omega_n^k)=\sum\limits_{j=0}^{n-1}a_je^{^{\frac{2\pi jk}{n}i}}\)

注意到\(jk=\dfrac{-(j-k)^2}{2}+\dfrac{j^2}{2}+\dfrac{k^2}{2}\)

\(f(\omega_n^k)=e^{^{\frac{\pi k^2}{n}i}}\sum\limits_{j=0}^{n-1}a_je^{^{\frac{\pi j^2}{n}i}}e^{^{\frac{-\pi (j-k)^2}{n}i}}\)

注意到\(a_je^{^{\frac{\pi j^2}{n}i}}e^{^{\frac{-\pi (j-k)^2}{n}i}}\)是卷积的形式

\(j-k\)可能为负，要整体平移\(n\)在复原

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#define eps 1e-9
using namespace std;
const int MAXN=2e5+5;
const double Pi=acos(-1.0);
struct Cpx{
	double a,b;
	Cpx(){
		a=0;
		b=0;
	} 
	Cpx(double aa,double bb){
		a=aa;
		b=bb;
	}
	Cpx operator*(const Cpx x)const{
		return Cpx(a*x.a-b*x.b,b*x.a+a*x.b);
	}
	Cpx operator+(const Cpx x)const{
		return Cpx(a+x.a,b+x.b);
	}
	Cpx operator-(const Cpx x)const{
		return Cpx(a-x.a,b-x.b);
	} 
	Cpx operator/(const Cpx x)const{
		return Cpx((a*x.a+b*x.b)/(x.a*x.a+x.b*x.b),(b*x.a-a*x.b)/(x.a*x.a+x.b*x.b));
	}
};
int Rev[MAXN*4];
struct Poly{
	vector<Cpx>V;
	void FFT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(V.size()<Len)
		{
			V.push_back(Cpx(0,0));
		}
		
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(V[i],V[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			Cpx Wn=Cpx(cos(Pi/l),type*sin(Pi/l));
			for(int i=0;i<Len;i+=(l<<1))
			{
				Cpx W=Cpx(1,0);
				for(int j=i;j<i+l;j++,W=W*Wn)
				{
					Cpx Xc=V[j];
					Cpx Yc=V[j+l]*W;
					V[j]=(Xc+Yc);
					V[j+l]=(Xc-Yc);
				}
			}
		}
		if(type==-1)
		{
			for(int i=0;i<Len;i++)
			{
				V[i].a/=Len;
				V[i].b/=Len;	
			}
		}
	}
}G,F;
Poly Mul_FFT(Poly A,Poly B)
{
	int N=A.V.size();
	int M=B.V.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.FFT(Lm,1);
	 B.FFT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.V[i]=A.V[i]*B.V[i];
	 }
	 A.FFT(Lm,-1);
	 while(A.V.size()>(N+M-1))
	 {
	 	A.V.pop_back();
	 }
	 return A;
}
Poly FFT_Bluestein(Poly X,int R,int type)
{
	Poly A,B;
	A.V.clear();
	B.V.clear();
	int N=X.V.size();
	for(int i=0;i<=(N-1)+R;i++)
	{
		int Fx=((i-N)*((i-N)));
		A.V.push_back(Cpx(cos(Pi*Fx/N),type*-1*sin(Pi*Fx/N)));
	}
	for(int i=0;i<=N-1;i++)
	{
		B.V.push_back(Cpx(cos(Pi*i*i/N),type*sin(Pi*i*i/N))*X.V[i]);
	} 

	A=Mul_FFT(A,B);
	B.V.clear();
	for(int i=0;i<R;i++)
	{
		int Fx=((i*i));
		B.V.push_back(A.V[i+N]*Cpx(cos(Pi*Fx/N),type*sin(Pi*Fx/N)));
	}
	if(type==-1)
	{
		for(int i=0;i<R;i++)
		{
			B.V[i].a/=N;
			B.V[i].b/=N;
		}
	} 
	return B;
}
int n;
double x;
int main()
{
	scanf("%d",&n);
	Poly B,C;
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		B.V.push_back(Cpx(x,0));
	}
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&x);
		C.V.push_back(Cpx(x,0));
	}
	B=FFT_Bluestein(B,n,1);
	C=FFT_Bluestein(C,n,1);
	for(int i=0;i<n;i++)
	{
		B.V[i]=C.V[i]*B.V[i];
	}
	B=FFT_Bluestein(B,n,-1);////
	for(int i=0;i<B.V.size();i++)
	{
		printf("%.4lf\n",B.V[i].a);
	}
}

牛顿迭代法

泰勒展开

\(f(x)=\sum\limits_{i=0}^{+\infin}\dfrac{f^{(i)}(a)}{i!}(x-a)^i\),称为\(f(x)\)在\(a\)处的展开,\(f^{(i)}(a)\)是对\(f\)求\(i\)阶导

我们要解决的问题是\(g(f(x))\equiv0(mod\ x^n)\),相当于抹去\(n\)次以后的时满足\(g\)的\(f(x)\),\(f(x)\)可以看作是自变量

对于\(g(f(x))\equiv0(mod\ x^1)\)，一般是可以直接算的

考虑倍增，现已求得\(g(f(x))\equiv0(mod\ x^{\frac{n}{2}})\)时的\(f_0(x)\)

目标是\(g(f(x))\equiv0(mod\ x^{{n}})\)时的\(f(x)\)

实际上\(f_0(x)\equiv f(x)(mod\ x^{\frac{n}{2}})\)

令\(f(x)=f_0(x)+p(x)x^{\frac{n}{2}}\)

对\(g(f(x))\)在\(f_0(x)\)处泰勒展开

\(g(f(x))=\sum\limits_{i=0}^{+\infin}\dfrac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i=\sum\limits_{i=0}^{+\infin}\dfrac{g^{(i)}(f_0(x))}{i!}(p(x)x^{\frac{n}{2}})^i=0(mod\ x^n)\)

注意到\(i\ge2\)时式子\(mod\ x^n=0\)，因此

\(g(f(x))=g(f_0(x))+g'(f_0(x))p(x)x^{\frac{n}{2}}=0(mod\ x^n)\)

\(p(x)x^{\frac{n}{2}}=-\dfrac{g(f_0(x))}{g'(f_0(x))}(mod\ x^n)\)

即\(f(x)=f_0(x)-\dfrac{g(f_0(x))}{g'(f_0(x))}(mod\ x^n)\)

发现可以用\(f_0(x)\)推出\(f(x)\)

\(eg1.\)多项式求逆

设原函数为\(h(x)\)

构造\(g(f(x))=h(x)-f^{-1}(x)\)

则\(f(x)=f_0(x)-\dfrac{h(x)-f_0^{-1}(x)}{f_0^{-2}(x)}=2f_0(x)-h(x)f_0^2(x)(mod\ x^n)\)

上式只有乘法因此可以直接多次计算，但要注意每次选取的长度,注意要拓展到\(2*n\)同时\(h\)不能太长

Poly Inverse(Poly A,int N)
{
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(inv(A.U[0],MOD));
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=((long long)Fn.U[j]*(2-((long long)Fn.U[j]*H.U[j])%MOD+MOD)%MOD)%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}

\(eg2.\)多项式取对数

即求\(\ln(h(x))\equiv f(x)(mod\ x^n)\)

同时取导数\(\dfrac{h'(x)}{h(x)}=f'(x)(mod\ x^n)\)

再积分\(\int\dfrac{h'(x)}{h(x)}dx=f(x)(mod\ x^n)\)

对于多项式求导与积分均可在\(\Theta(n)\)解决

注意\([x_0]h(x)=1\)时才能在模域内取对数，而先导再积分会丢失常数项

Poly Der(Poly x)
{
	Poly Nex;
	Nex.U.clear();
	for(int i=1;i<x.U.size();i++){
		Nex.U.push_back(((long long)i*x.U[i])%MOD);
	}
	return Nex;
}
Poly Ing(Poly x)
{
	Poly Nex;
	Nex.U.clear();
	Nex.U.push_back(0);
	for(int i=0;i<x.U.size();i++)
	{
		Nex.U.push_back(((long long)x.U[i]*inv(i+1,MOD))%MOD);
	}
	return Nex;
}
Poly Ln(Poly x,int N)
{
	Poly ex=Der(x);
	Poly ey=Inverse(x,N);
	ex=Mul_NTT(ex,ey);
	ex=Ing(ex);
	while(ex.U.size()>N)
	{
		ex.U.pop_back();
	}	
	return ex;
}

\(eg3.\)多项式\(exp\)

构造\(g(f(x))=h(x)-\ln(f(x))\)

则\(f(x)=f_0(x)+(h(x)-ln(f_0(x)))f_0(x)=f_0(x)(1+h(x)-ln(f_0(x)))(mod\ x^n)\)

关于长度

实际上在计算的时候每个多项式都应该取到\(n\),只不过\(f_0(x)\)长度只有\(\frac{n}{2}\),但要注意如果长度取得长,乘法的结果的长度也要对应,否则就算成循环卷积了

Poly Exp(Poly A,int N)
{
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(1);
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		Poly Fln=Ln(Fn,l);
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		Fln.NTT(Lm+1,1);
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=((long long)Fn.U[j]*(((long long)H.U[j]+1-Fln.U[j]+MOD)%MOD))%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}

\(eg4.\)多项式带余除法

定义\(f_R(x)\)为\(f(x)\)反转系数后的结果

发现\(f_R(x)=x^nf(\frac{1}{x})\)

对于\(f(x)=q(x)g(x)+r(x)\)

设\(f(x)\)的次数\(n\),\(g(x)\)的次数\(m\),则\(q(x)\)的次数为\(n-m\),\(r(x)\)的次数小于\(m\)

\(f(\frac{1}{x})=q(\frac{1}{x})g(\frac{1}{x})+r(\frac{1}{x})\)

\(x^nf(\frac{1}{x})=x^{n-m}q(\frac{1}{x})x^ng(\frac{1}{x})+x^{n-m+1}x^{m-1}r(\frac{1}{x})=f_R(x)=q_R(x)g_R(x)+x^{n-m+1}x^{m-1}r(\frac{1}{x})\)

\(f_R(x)=q_R(x)g_R(x)(mod\ x^{n-m+1})\)

\(q_R(x)=\dfrac{f_R(x)}{g_R(x)}(mod\ x^{n-m+1})\)

注意到\(q\)的最高次数为\(n-m\),因此这就是\(q_R(x)\)的真实值,且次数一定为\(n-m\)

而\(r(x)=f(x)-q(x)g(x)\)

Poly Div(Poly F,Poly G)
{
	int N=F.U.size()+1;
	int M=G.U.size()+1;
	reverse(F.U.begin(),F.U.end());
	reverse(G.U.begin(),G.U.end());
	Poly Mx=Inverse(G,N-M+1);
	F=Mul_NTT(F,Mx);
	Poly Nx;
	Nx.U.clear();
	for(int i=0;i<N-M+1;i++)
	{
		Nx.U.push_back(F.U[i]);
	}
	reverse(Nx.U.begin(),Nx.U.end());
	return Nx;
}
Poly Mod(Poly F,Poly G)
{
	Poly Gx=Div(F,G);
	Gx=Mul_NTT(Gx,G);
	Poly Nex;
	Nex.U.clear();
	int M=G.U.size();
	for(int i=0;i<M;i++)
	{
		Nex.U.push_back(((long long)F.U[i]-Gx.U[i]+MOD)%MOD);
	}
	return Nex;
};

\(eg5.\)多项式开根

构造函数\(g(f(x))=h(x)-f^2(x)\)

\(f(x)=f_0(x)+\frac{h(x)-f_0^2(x)}{2f_0(x)}=\frac{f_0(x)}{2}+\frac{h(x)}{2f_0(x)}\)

Poly Sqrt(Poly A,int N)
{
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(1);
	int inv2=(MOD-MOD/2);
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		Poly Fln=Inverse(Fn,l);
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		Fln.NTT(Lm+1,1);
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=(((long long)Fn.U[j]*inv2)%MOD+(((long long)inv2*Fln.U[j])%MOD*H.U[j])%MOD)%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}

eg6.多项式快速幂

\(G(x)=F^k(x)\)

同时取对数\(\ln G(x)=k\ln F(x)\)

再取\(exp\)

\(G(x)=exp(k\ln F(x))\)

注意\([x^0]F(x)\not=1\)时是无法取\(\ln\)的

找到第一个\(t\),\([x^t]F(x)\not=0\)

\(G(x)=exp(k\ln F(x))=exp(k\ln \dfrac{F(x)}{a_tx^t}a_tx^t)=exp(k\ln \dfrac{F(x)}{a_tx^t})a_t^kx^{kt}\)

#include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
const int MAXN=1e5+5;
const int MOD=998244353;
const int g=3;
const long double Pi=acos(-1.0);
const int p=32000;
int Rev[MAXN*4];
struct Cpx{
	long double a,b;
	Cpx(){
		a=0;
		b=0;
	} 
	Cpx(long double aa,long double bb){
		a=aa;
		b=bb;
	}
	Cpx operator*(const Cpx x)const{
		return Cpx(a*x.a-b*x.b,b*x.a+a*x.b);
	}
	Cpx operator+(const Cpx x)const{
		return Cpx(a+x.a,b+x.b);
	}
	Cpx operator-(const Cpx x)const{
		return Cpx(a-x.a,b-x.b);
	} 
};
int Pow(int a,int b,int pff){
	int res=1;
	int base=a;
	while(b)
	{
		if(b&1)
		{
			res=((long long)res*base)%pff;
		}
		base=((long long)base*base)%pff;
		b>>=1;
	}
	return res;
}
int inv(int a,int pff){
	return Pow(a,pff-2,pff);
}
struct Poly{
	vector<int>U;
	vector<Cpx>V;
	int size()
	{
		return U.size();
	}
	void push_back(int x)
	{
		U.push_back(x);
		return;
	}
	void clear()
	{
		U.clear();
		return;
	}
	void FFT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(V.size()<Len)
		{
			V.push_back(Cpx(0,0));
		}
		
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(V[i],V[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			Cpx Wn=Cpx(cos(Pi/l),type*sin(Pi/l));
			for(int i=0;i<Len;i+=(l<<1))
			{
				Cpx W=Cpx(1,0);
				for(int j=i;j<i+l;j++,W=W*Wn)
				{
					Cpx Xc=V[j];
					Cpx Yc=V[j+l]*W;
					V[j]=(Xc+Yc);
					V[j+l]=(Xc-Yc);
				}
			}
		}
		if(type==-1)
		{
			for(int i=0;i<Len;i++)
			{
				V[i].a/=Len;	
			}
		}
	}
	void NTT(int Limit,int type)
	{
		int Len=(1<<Limit);
		for(int i=0;i<Len;i++)
		{
			Rev[i]=((Rev[i>>1]>>1)|((i&1)<<(Limit-1)));
		}
		
		while(U.size()<Len)
		{
			U.push_back(0);
		}
		for(int i=0;i<Len;i++)
		{
			if(i<Rev[i])
			{
				swap(U[i],U[Rev[i]]);
			}
		}
		for(int l=1;l<Len;l<<=1)
		{
			int Wn=Pow(g,(MOD-1)/(l<<1),MOD);
			if(type==-1)
			{
				Wn=inv(Wn,MOD);
			}
			for(int i=0;i<Len;i+=(l<<1))
			{
				int W=1;
				for(int j=i;j<i+l;j++,W=((long long)W*Wn)%MOD)
				{
					int Xc=U[j];
					int Yc=((long long)U[j+l]*W)%MOD;
					U[j]=((long long)Xc+Yc)%MOD;
					U[j+l]=((long long)Xc-Yc+MOD)%MOD;
				}
			}
		}
		if(type==-1)
		{
			int Liv=inv(Len,MOD); 
			for(int i=0;i<Len;i++)
			{
				U[i]=((long long)U[i]*Liv)%MOD;	
			}
		}
	}
};
Poly Mul_FFT(Poly A,Poly B){
	int N=A.V.size();
	int M=B.V.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.FFT(Lm,1);
	 B.FFT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.V[i]=A.V[i]*B.V[i];
	 }
	 A.FFT(Lm,-1);
	 while(A.V.size()>(N+M-1))
	 {
	 	A.V.pop_back();
	 }
	 return A;
}
Poly Mul_NTT(Poly A,Poly B){
	int N=A.U.size();
	int M=B.U.size();
	int nox=1;
	int Lm=0;
	while(nox<=(N+M-2))
	{
		nox<<=1;
		Lm++;
	 } 
	 A.NTT(Lm,1);
	 B.NTT(Lm,1);
	 for(int i=0;i<nox;i++)
	 {
	 	A.U[i]=((long long)A.U[i]*B.U[i])%MOD;
	 }
	 A.NTT(Lm,-1);
	 while(A.U.size()>(N+M-1))
	 {
	 	A.U.pop_back();
	 }
	 return A;
}
Poly Mul_MTT(Poly A,Poly B,int P){
	int n=A.U.size()-1;
	int m=B.U.size()-1;
	Poly A1,B1,C1,D1;
	A1.V.clear();
	B1.V.clear();
	C1.V.clear();
	D1.V.clear();
	Poly A2,B2,C2,D2;
	A2.V.clear();
	B2.V.clear();
	C2.V.clear();
	D2.V.clear();
	for(int i=0;i<=n;i++)
	{
		int x=A.U[i];
		int ax=x/p;
		int bx=x%p;
		A1.V.push_back(Cpx(ax,0));
		B1.V.push_back(Cpx(ax,0));
		C1.V.push_back(Cpx(bx,0));
		D1.V.push_back(Cpx(bx,0));
	}
	for(int i=0;i<=m;i++)
	{
		int x=B.U[i];
		int ax=x/p;
		int bx=x%p;
		A2.V.push_back(Cpx(ax,0));
		B2.V.push_back(Cpx(bx,0));
		C2.V.push_back(Cpx(ax,0));
		D2.V.push_back(Cpx(bx,0));
	}
	A1=Mul_FFT(A1,A2);
	B1=Mul_FFT(B1,B2);
	C1=Mul_FFT(C1,C2);
	D1=Mul_FFT(D1,D2);
	Poly Polat;
	Polat.U.clear();
	for(int i=0;i<A1.V.size();i++)
	{
		long long F1=(long long)(A1.V[i].a+0.5);
		long long F2=(long long)(B1.V[i].a+0.5);
		long long F3=(long long)(C1.V[i].a+0.5);
		long long F4=(long long)(D1.V[i].a+0.5);
		long long F=(((F1*p)%P*p)%P+(((F2+F3)%P)*p)%P+F4)%P;
		Polat.U.push_back(F);
	}
	return Polat;
}
Poly Inverse(Poly A,int N){
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(inv(A.U[0],MOD));
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=((long long)Fn.U[j]*(2-((long long)Fn.U[j]*H.U[j])%MOD+MOD)%MOD)%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}
Poly Der(Poly x){
	Poly Nex;
	Nex.U.clear();
	for(int i=1;i<x.U.size();i++){
		Nex.U.push_back(((long long)i*x.U[i])%MOD);
	}
	return Nex;
}
Poly Ing(Poly x){
	Poly Nex;
	Nex.U.clear();
	Nex.U.push_back(0);
	for(int i=0;i<x.U.size();i++)
	{
		Nex.U.push_back(((long long)x.U[i]*inv(i+1,MOD))%MOD);
	}
	return Nex;
}
Poly Ln(Poly x,int N){
	Poly ex=Der(x);
	Poly ey=Inverse(x,N);
	ex=Mul_NTT(ex,ey);
	ex=Ing(ex);
	while(ex.U.size()>N)
	{
		ex.U.pop_back();
	}	
	return ex;
}
Poly Exp(Poly A,int N){
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(1);
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		Poly Fln=Ln(Fn,l);
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		Fln.NTT(Lm+1,1);
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=((long long)Fn.U[j]*(((long long)H.U[j]+1-Fln.U[j]+MOD)%MOD))%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}
Poly Div(Poly F,Poly G){
	int N=F.U.size()+1;
	int M=G.U.size()+1;
	reverse(F.U.begin(),F.U.end());
	reverse(G.U.begin(),G.U.end());
	Poly Mx=Inverse(G,N-M+1);
	F=Mul_NTT(F,Mx);
	Poly Nx;
	Nx.U.clear();
	for(int i=0;i<N-M+1;i++)
	{
		Nx.U.push_back(F.U[i]);
	}
	reverse(Nx.U.begin(),Nx.U.end());
	return Nx;
}
pair<Poly,Poly> Mod(Poly F,Poly G){
	Poly Gx=Div(F,G);
	Gx=Mul_NTT(Gx,G);
	Poly Nex;
	Nex.U.clear();
	int M=G.U.size();
	for(int i=0;i<M;i++)
	{
		Nex.U.push_back(((long long)F.U[i]-Gx.U[i]+MOD)%MOD);
	}
	return make_pair(Gx,Nex);
};
Poly Sqrt(Poly A,int N){
	Poly Fn;
	Fn.U.clear();
	Fn.U.push_back(1);
	int inv2=(MOD-MOD/2);
	if(N==1)
	{
		return Fn;
	}
	for(int l=2,Lm=1;(l>>1)<N;l<<=1,Lm++)
	{
		Poly H;
		H.U.clear();
		for(int j=0;j<l;j++)
		{
			if(j<A.U.size())
			{
				H.U.push_back(A.U[j]);
			}
			else
			{
				H.U.push_back(0);
			}
		}	
		Poly Fln=Inverse(Fn,l);
		H.NTT(Lm+1,1);
		Fn.NTT(Lm+1,1);
		Fln.NTT(Lm+1,1);
		for(int j=0;j<l*2;j++)
		{
			Fn.U[j]=(((long long)Fn.U[j]*inv2)%MOD+(((long long)inv2*Fln.U[j])%MOD*H.U[j])%MOD)%MOD;
		}
		Fn.NTT(Lm+1,-1);
		while(Fn.U.size()>l)
		{
			Fn.U.pop_back();
		}
	}
	while(Fn.U.size()>N)
	{
		Fn.U.pop_back();
	}
	return Fn;
}
int Bostan_Mori(Poly F,Poly G,int N){
//	printf("fic\n");
	if(F.size()>=G.size())
	{
		pair<Poly,Poly>DJNB=Mod(F,G);
		int Fsd=0;
		if(N<DJNB.first.size())
		{
			Fsd=DJNB.first.U[N];
		}
		return ((long long)Bostan_Mori(DJNB.second,G,N)+(Fsd))%MOD;
	}
	else
	{
		while(N)
		{
	//		printf("fic\n");
			Poly GT=G;
			for(int i=1;i<GT.size();i+=2)
			{
				GT.U[i]=(MOD-GT.U[i]);
			}
			Poly V=Mul_NTT(G,GT);
			Poly U=Mul_NTT(F,GT);
			F.U.clear();
			G.U.clear();
			for(int i=N&1;i<U.size();i+=2)
			{
				F.U.push_back(U.U[i]);
			}
			for(int i=0;i<V.size();i+=2)
			{
				G.U.push_back(V.U[i]);
			}
			N>>=1;
		}
		if(F.size()==0)
		{
			return 0;
		}
		return ((long long)F.U[0]*inv(G.U[0],MOD))%MOD;
	}
}
int Homo_recursion(Poly A,Poly R,int n){
	int k=R.size();
	Poly S=Mul_NTT(A,R);
	S.U.resize(k-1);
	return Bostan_Mori(S,R,n);
}
Poly BM(Poly A){
	int N=A.size();
	Poly R,Rp;
	int Pa,Pap;
	R.clear();
	R.push_back(1);
	Poly Nex;
	Nex.clear();
	int P=-1;
	for(int i=0;i<N;i++)
	{
		Nex.push_back(A.U[i]);
		Pa=0;
		for(int j=0;j<R.U.size();j++)
		{
			Pa=((long long)Pa+((long long)Nex.U[(i)-j]*R.U[j])%MOD)%MOD;
		}
		if(Pa==0)
		{
			continue;
		}
		else
		{
			if(P==-1)
			{
				R.U.resize(i+2,0);
				R.U[0]=1;
				Rp.U.clear();
				Rp.U.push_back(1);
				Pap=Pa;
				P=i;
			}
			else
			{
				int IDSY=((long long)Pa*inv(Pap,MOD))%MOD;
				if(Rp.size()+i-P>R.size())
				{
					Poly WORNMDESBTIR;
					WORNMDESBTIR.clear();
					WORNMDESBTIR.U.resize(Rp.size()+i-P,0);
					for(int j=0;j<Rp.size();j++)
					{
						WORNMDESBTIR.U[j+i-P]=(MOD-((long long)Rp.U[j]*IDSY)%MOD);
					}
					for(int j=0;j<R.size();j++)
					{
						WORNMDESBTIR.U[j]=((long long)WORNMDESBTIR.U[j]+R.U[j])%MOD;
					}
					
					Rp=R;	
					Pap=Pa;
					P=i;
					R=WORNMDESBTIR;
				}
				else
				{
					for(int j=0;j<Rp.size();j++)
					{
						R.U[j+i-P]=((long long)R.U[j+i-P]-((long long)Rp.U[j]*IDSY)%MOD+MOD)%MOD;
					}
				}		
			}
			
		}
	}
	
	return R; 
}
Poly Pow_Poly(Poly A,string k,int n)
{
	Poly Nex;
	Nex.clear();
	int t=n+1;
	for(int i=0;i<A.size();i++)
	{
		if(A.U[i]!=0)
		{
			t=i;
			break;
		}
	}
	int Rv=inv(A.U[t],MOD);
	for(int i=t;i<A.size();i++)
	{
		Nex.U.push_back(((long long)A.U[i]*Rv)%MOD);
	}
	Nex=Ln(Nex,n);
	int K1=0,K2=0;
	int K3=0;
	for(int i=0;i<k.size();i++)
	{
		if(K3<=n)
		{
			K3=(K3*10+(k[i]-'0'));
		}
		K1=((long long)K1*10+(k[i]-'0'))%MOD;
		K2=((long long)K2*10+(k[i]-'0'))%(MOD-1);
	}
//	printf("%d %d %d?\n",K1,K2,K3);
	if((long long)t*K3>n||(t==(n+1)))
	{
		Poly E;
		E.U.resize(n,0);
	//	printf("guck");
		return E;
	}
	for(int i=0;i<Nex.size();i++)
	{
		Nex.U[i]=((long long)Nex.U[i]*K1)%MOD;
	} 
	Nex=Exp(Nex,n);
	int Exc=Pow(A.U[t],K2,MOD);
	for(int i=0;i<Nex.size();i++)
	{
		Nex.U[i]=((long long)Nex.U[i]*Exc)%MOD;	
	} 
	Poly Rxf;
	Rxf.U.resize(n,0);
	int Rf=t*K3;
	//printf("%d??\n",Rf);
	for(int i=0;i<Nex.size();i++)
	{
		if(i+Rf>=n)
		{
			break;
		}
		Rxf.U[i+Rf]=Nex.U[i];
		
	}
	return Rxf;
}
signed main(){
}