多项式的运算

多项式不难，真的不难。
代码比较丑，勿喷

前置知识讲解：
FFT/NTT:https://www.cnblogs.com/fenghaoran/p/7107608.html
二次剩余：https://zhuanlan.zhihu.com/p/166123245
微积分：https://zhuanlan.zhihu.com/p/94592123

配套题单：https://www.luogu.com.cn/training/95194

多项式乘法

FFT/NTT即可
模板：https://www.luogu.com.cn/problem/P3803
FFT：

点击查看代码

void FFT(ljq *p,int op,int n) {
	int l=log2(n);
	for(int i=0; i<n; i++)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0; i<n; i++)
		if(i<r[i])
			swap(p[i],p[r[i]]);
	for(int i=1; i<n; i*=2) {
		ljq w_n(cos(pi/i),op*sin(pi/i));
		for(int j=0; j<n; j+=i*2) {
			ljq w(1,0),t1,t2;
			for(int k=0; k<i; k++,w*=w_n)
				t1=p[j+k],t2=p[j+k+i]*w,p[k+j]=t1+t2,p[k+j+i]=t1-t2;
		}
	}
	if(op==-1)
		for(int i=0;i<=n;i++)
			p[i]/=n;
}

NTT：

点击查看代码

void NTT(int *p,int op,int n) {
	int l=log2(n);
	r[0]=0;
	for(int i=0; i<n; i++)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0; i<n; i++)
		if(i<r[i])
			swap(p[i],p[r[i]]);
	for(int i=1; i<n; i*=2) {
		int tmp=power(g,(mod-1)/(i*2));
		if(op==-1)
			tmp=power(tmp,mod-2);
		for(int j=0; j<n; j+=i*2) {
			int w=1,t1,t2;
			for(int k=0; k<i; k++,w=1ll*w*tmp%mod)
				t1=p[j+k],t2=1ll*p[j+k+i]*w%mod,p[k+j]=(0ll+t1+t2+mod)%mod,p[k+j+i]=(0ll+t1-t2+mod)%mod;
		}
	}
	if(op==-1) {
		int invn=power(n,mod-2);
		for(int i=0; i<n; i++)
			p[i]=1ll*p[i]*invn%mod;
	}
}

多项式求逆

保证(a_0 eq 0) ,有 (A(x) imes B(x) equiv 1 pmod {x^n})，求(B(x))
假设我们已经求出

[A(x) imes B'(x) equiv 1 pmod {x^{frac{n}{2}}} ]

[A(x) imes B'(x) equiv 1 pmod {x^{frac{n}{2}}} ]

[A(x) imes B'(x)-1 equiv 0 pmod {x^{frac{n}{2}}} ]

左右两遍同时平方(模数也要平方)

[A(x)^2B'(x)^2+1-2 imes A(x) B'(x) equiv 0 pmod {x^n} ]

[2 imes A(x) B'(x)-A(x)^2B'(x)^2 equiv 1 pmod {x^n} ]

[A(x){2B'(x)-A(x)B'(x)^2} equiv 1 pmod {x^n} ]

那么此时的逆元就变成了(2B'(x)-A(x)B'(x)^2)
倍增上去做即可。
时间复杂度：(O(nlogn))
模板：https://www.luogu.com.cn/problem/P4238

点击查看代码

void inv(int *f,int *p,int m) {
	if(m==1)
		return p[0]=power(f[0],mod-2),void();
	inv(f,p,(m+1)>>1);
	int n;
	for(n=1; n<=m; n*=2);
	n*=2;
	for(int i=0; i<m; i++)
		c[i]=f[i];
	for(int i=m; i<n; i++)
		c[i]=0;
	NTT(c,1,n);
	NTT(p,1,n);
	for(int i=0; i<n; i++)
		p[i]=1ll*(2ll-1ll*c[i]*p[i]%mod+mod)%mod*p[i]%mod;
	NTT(p,-1,n);
	for(int i=m; i<n; i++)
		p[i]=0;
}

多项式开根

有(B(x)^2=A(x) pmod {x^n})，求(B(x))
同样，我们使用倍增法。

[B'(x)^2 equiv A(x) pmod{x^{frac{n}{2}}} ]

[B'(x)^2-A(x) equiv 0 pmod {x^{frac{n}{2}}} ]

[[B'(x)^2-A(x)]^2 equiv 0 pmod {x^n} ]

[[B'(x)^2+A(x)]^2 equiv 4 imes B'(x)^2 imes A(x) pmod {x^n} ]

[[frac{B'(x)^2+A(x)}{2 B'(x)}]^2 equiv A(x) pmod {x^n} ]

[B(x)=frac{B'(x)^2+A(x)}{2 B'(x)} ]

对于初值，求个二次剩余就好了。
二次剩余:https://zhuanlan.zhihu.com/p/166123245
时间复杂度：(O(nlogn))
模板：https://www.luogu.com.cn/problem/P5205
二次剩余版：https://www.luogu.com.cn/problem/P5277

点击查看代码

void Sqrt(int *f,int *p,int m){
	if(m==1)
		return p[0]=getx2(f[0]),void();
	int n;
	for(n=1;n<=m;n*=2);n*=2;
	for(int i=0;i<n;i++)
		p[i]=0;
	Sqrt(f,p,(m+1)>>1);
	for(int i=0;i<=m;i++)
		d[i]=f[i];
	for(int i=m;i<n;i++)
		d[i]=0;
	for(int i=0;i<n;i++)
		iv[i]=0;
	inv(p,iv,m);
	NTT(d,1,n);NTT(p,1,n);NTT(iv,1,n);
	for(int i=0;i<n;i++)
		p[i]=1ll*(0ll+p[i]+1ll*d[i]*iv[i]%mod)%mod*power(2,mod-2)%mod;
	NTT(p,-1,n);
	for(int i=0;i<n;i++)
		p[i]=1ll*p[i];
	for(int i=m;i<n;i++)
		p[i]=0;
}

多项式除法|取模

求(f(x))除以(g(x))的商(Q(x))与余数(R(x))，其中(deg_{R(x)}<deg_{f(x)})。

[f(x)=g(x)Q(x) +R(x) ]

设(deg_{f(x)}=n,deg_{g(x)}=m,deg_{Q(x)}=n-m,deg_{R(x)}) 小于 (m),这里假设(deg_{R(x)}=m-1)
我们将(x^{-1})带入(f(x))得(f(x^{-1}))
设(f^T(x)=x^nf(x^{-1}))，实质上就是反转系数。
我们把(x^{-1})带入原式。

[f(x^{-1})=g(x^{-1})Q(x^{-1})+R(x^{-1}) ]

左右两边同时乘上(x^n)

[x^nf(x^{-1})=x^ng(x^{-1})Q(x^{-1})+x^nR(x^{-1}) ]

[f^T(x)=g^T(x)Q^T(x)+x^{n-m+1}R^T(x) ]

将这个式子放到(pmod x^{n-m+1})意义下，我们发现(x^{n-m+1}R^T(x))被消去了，而(Q_T(x))项的系数为(n-m)小于(n-m+1)，不会受到影响，于是就有
(f^T(x)=g^T(x)Q^T(x) pmod {x^{n-m+1}})
根据多项式求逆可以求出(Q_T(x))，然后求出(Q(x))，如果要求(R(x))回代一下就好了。
时间复杂度：(O(nlogn))
模板：https://www.luogu.com.cn/problem/P4512

点击查看代码

void div(int *F,int *G,int *Q,int *R,int n,int m) {
	int sum=n+m,len=1;
	for(len=1; len<=sum; len*=2);len*=2;
	turn(F,ft,n);turn(G,gt,m);
	for(int i=n-m+1; i<=len; i++)
		ft[i]=0,gt[i]=0;
	for(int i=0; i<=n; i++)
		d[i]=0,iv[i]=0;
	inv(gt,iv,n-m+1);
	NTT(ft,1,len);NTT(iv,1,len);
	for(int i=0; i<len; i++)
		d[i]=1ll*ft[i]*iv[i]%mod;
	NTT(d,-1,len);
	for(int i=n-m+1; i<=len; i++)
		d[i]=0;
	turn(d,Q,n-m);
	for(int i=0; i<=len; i++)
		ft[i]=F[i];
	for(int i=0; i<=len; i++)
		gt[i]=G[i];
	NTT(ft,1,len);NTT(gt,1,len);NTT(Q,1,len);
	for(int i=0; i<=len; i++)
		R[i]=(0ll+ft[i]-1ll*gt[i]*Q[i]%mod+mod)%mod;
	NTT(Q,-1,len);NTT(R,-1,len);
	for(int i=m; i<=len; i++)
		R[i]=0;
	for(int i=n-m+1; i<=len; i++)
		Q[i]=0;
}

前置知识:求导 微积分

https://zhuanlan.zhihu.com/p/94592123

多项式求导

我们可以对多项式的每一项进行求导，然后根据导数的加法运算求出多项式的导数。
考虑单项：((a_ix^i)'=a_i(x^i)'=a_iix^{i-1})
因此多项式(A(x)=a_0+a_1x^1+a_2x^2+...+a_nx^n)的导数(A'(x)=a_1 imes 1+(a_2 imes 2)x+(a_3 imes 3)x^2+...+(a_n imes n)x^{n-1})

点击查看代码

void Der(int *f,int *p,int n) {
	for(int i=0; i<n; i++)
		p[i]=1ll*f[i+1]*(i+1)%mod;
	p[n]=0;
}

多项式积分

多项式积分及求(f(x))的原函数(F(x))

[F(x)=int f(x)dx=sum_{i=1}^{n+1}frac{a_{i-1}}{i}x^i ]

点击查看代码

void Int(int *f,int *p,int n) {
	p[0]=0;
	for(int i=1; i<=n+1; i++)
		p[i]=1ll*f[i-1]*power(i,mod-2)%mod;
}

多项式ln

保证(a_0=1)求(B(x)=ln A(x) pmod {x^n})
首先两边同时求导(B'(x) equiv A'(x)frac{1}{A(x)} pmod {x^n})
求出(B'(x))，然后积分求出(B(x))
模板:https://www.luogu.com.cn/problem/P4725

点击查看代码

void Ln(int *f,int *p,int n) {
	int len=1; for(len=1; len<=n; len*=2); len*=2;
	for(int i=0; i<=len; i++)
		ft[i]=gt[i]=0;
	Der(f,ft,n); inv(f,gt,n+1);
	NTT(ft,1,len); NTT(gt,1,len);
	for(int i=0; i<=len; i++)
		ft[i]=1ll*ft[i]*gt[i]%mod;
	NTT(ft,-1,len);
	Int(ft,p,n-1);
}

多项式的泰勒展开式

设一个多项式(f(x))为(n)次多项式

[f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_nx^n ]

对(f(x))进行(n)次微分。

[f'(x)=1 imes a_1+ 2 imes a_2 x+3 imes a_3 x^2+...+n imes a_n x^{n-1} ]

[f''(x)= 1 imes 2 imes a_2 + 2 imes 3 imes a_3 x+3 imes 4 imes a_4 x^2+...+(n-1) imes n imes a_n x^{n-2} ]

[f'''(x)=1 imes 2 imes 3 imes a_3 +2 imes 3 imes 4 imes a_4 x+3 imes 4 imes 5 imes a_5 x^2+...+(n-2) imes (n-1) imes n imes a_n imes x^{n-3}]

[... ]

[f^{(n)}(x)= 1 imes 2 imes 3 imes ... imes n imes a_n ]

将(x=0)带入上述各式可得

[a_0=f(0),a_1=frac{f'(0)}{1!},a_2=frac{f''(0)}{2!},a_3=frac{f'''(0)}{3!},...,a_n=frac{f^{(n)}(0)}{n!} ]

带回原式可得

[f(x)=f(0)+frac{f'(0)}{1!}x+frac{f''(0)}{2!}x^2+frac{f'''(0)}{3!}x^3+...+frac{f^{(n)}(0)}{n!}x^n ]

将((x-x_0))代入可得

[f(x)=f(x_0)+frac{f'(x_0)}{1!}(x-x_0)+frac{f''(x_0)}{2!}(x-x_0)^2+frac{f'''(x_0)}{3!}(x-x_0)^3+...+frac{f^{(n)}(x_0)}{n!}(x-x_0)^n ]

多项式牛顿迭代

已知(g(x)),求(f(x) pmod {x^n})满足(g(f(x)) equiv 0 pmod {x^n})
(g(f(x)) equiv 0 pmod {x})的解需要提前求出
假设已经知道了(g(f(x)) equiv 0 pmod {x^{frac{n}{2}}})的解(f_0(x))
将(g(f(x)))在(g(f(x)))处泰勒展开得

[sum_{i=0}^{+ infty}frac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i equiv 0 pmod {x^n} ]

由于(f(x)-f_0(x))的最低非零此项的系数最低为(lceil{frac{n}{2}} ceil)

[forall i in [2,+infty],(f(x)-f_0(x))^i equiv 0 pmod {x^n} ]

故有

[g(f_0(x))+g'(f_0(x)) imes (f(x)-f_0(x)) equiv 0 pmod {x^n} ]

[f(x) equiv f_0(x) - frac{g(f_0(x))}{g'(f_0(x))} pmod {x^n} ]

多项式exp

保证(a_0=0),求(f(x) equiv e^{A(x)} pmod {x^n})((a_0=0))
两边同时取自然对数得

[ln f(x) equiv A(x) pmod {x^n} ]

[ln f(x) - A(x) equiv 0 pmod {x^n} ]

构造(g(f(x))=ln f(x) -A(x))，(g'(f(x))=frac{1}{f(x)})牛顿迭代得

[f(x) equiv f_0(x)-frac{g(f_0(x)}{g'(f_0(x))} equiv f_0(x)(1 - ln f_0(x) + A(x)) pmod{x^n} ]

当(n=1)时，(f(x)=e^0=1)
时间复杂度为(O(nlogn))
模板:https://www.luogu.com.cn/problem/P4726

点击查看代码

void Exp(int *f,int *p,int n) {
	if(n==1)
		return p[0]=1,void();
	int len=1;for(len=1; len<=n; len*=2);len*=2;
	for(int i=0; i<=len; i++)
		d[i]=0,p[i]=0;
	Exp(f,p,(n+1)>>1);
	for(int i=0; i<len; i++)
		gt[i]=0;
	Ln(p,gt,n);
	for(int i=n; i<len; i++)
		d[i]=0;
	for(int i=0; i<n; i++)
		d[i]=f[i];
	NTT(p,1,len);NTT(gt,1,len);NTT(d,1,len);
	for(int i=0; i<len; i++)
		p[i]=1ll*(1ll-gt[i]+d[i]+mod)%mod*1ll*p[i]%mod;
	NTT(p,-1,len);
	for(int i=n; i<len; i++)
		p[i]=0;
}

多项式k次幂

保证(a_0=1) , (计算)A^k(x)( 如果直接用快速幂来做的话时间复杂度为)O(nlognlogk)( 我们可以将幂的形式转换为ln和exp )A^k(x)=e{ln A^k(x)}=e{k imes ln A(x)}( 时间复杂度为)O(nlogn)$

加强

如果(a_0 eq 1)呢？
分类讨论一下。
如果(a_0=0)我们可以将多项式平移，将其化成(A(x)=x^tB(x))
如果(a_0 eq 1)且(a_0 eq 0)，我们可以将常数项除掉,(A(x)=a_0 C(x))
然后就有(A(x)=x^t[a_0 C(x)])

[A^k(x)={x^t[a_0 C(x)]}^k=x^{tk} a_0^k C^k(x) ]

再加强

如果(k leq 10^{10^5})呢？
已知(A^k(x)={x^t[a_0 C(x)]}^k=x^{tk} a_0^k C^k(x))
(x^{tk})好解决，如果(tk)很大，答案直接(A^k(x)=0)就好了。
(a_0^k)的解决上需要用到费马小定理(a^{p-1} equiv 1pmod p)
所以(a_0^k=a_0^{k mod (p-1)})
那么(C^k(x))呢
对于多项式(C(x))我们有(C^p(x) equiv 1 pmod p)
这个应该可以用多项式定理进行证明。
然后(C^k(x) equiv C^{k mod p}(x) pmod p)
模板:https://www.luogu.com.cn/problem/P5245

点击查看代码

void Power(int *f,int *p,int n,int k1,int k2,int ov)
{
	int num0=0;
	for(int i=0;i<n;i++)
	{
		if(f[i]==0)
			num0++;
		else
			break;
	}
	if(f[0]==0&&ov)
		return ;
	if(1ll*num0*k1>=1ll*n)
		return ;
	for(int i=0;i<n;i++)
		e[i]=f[i+num0];
	int nop=e[0],ivnop=power(nop,mod-2),pwnop=power(nop,k2);
	for(int i=0;i<n;i++)
		e[i]=1ll*e[i]*ivnop%mod;
	Ln(e,p,n);
	for(int i=0;i<n;i++)
		e[i]=1ll*k1*p[i]%mod;
	for(int i=0;i<n;i++)
		tt[i]=0;
	Exp(e,tt,n);
	for(int i=k1*num0;i<n;i++)
		p[i]=1ll*tt[i-num0*k1]*pwnop%mod;
	for(int i=0;i<num0*k1;i++)
		p[i]=0;
}

多项式三角函数

求出(sin f(x),cos f(x), an f(x) pmod {x^n})
根据欧拉公式(e^{ix}=cos x+isin x)以及(e^{-ix}=cos x -i sin x)
我们可以得到三角函数的另一种表达式

[sin x =frac{e^{ix}-e^{-ix}}{2i} ]

[cos x =frac{e^{ix}+e^{-ix}}{2} ]

带入得

[sin f(x)=frac{e^{i f(x)}-e^{-i f(x)}}{2i} ]

[cos f(x) =frac{e^{i f(x)}+e^{-i f(x)}}{2} ]

求出(sin f(x)) 和 (cos f(x)) 后根据 ( an f(x)=frac{sin f(x)}{cos f(x)}) 求出 ( an f(x))
时间复杂度:O(nlogn)
模板:https://www.luogu.com.cn/problem/P5264

点击查看代码

void Sin(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=1ll*f[i]*I%mod,c2[i]=-c1[i];
	Exp(c1,d1,n);Exp(c2,d2,n);
	for(int i=0;i<n;i++)
		p[i]=(d1[i]-d2[i]+mod)%mod*1ll*power(2*I,mod-2)%mod;
}
void Cos(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=1ll*f[i]*I%mod,c2[i]=-c1[i];
	Exp(c1,d1,n);Exp(c2,d2,n);
	for(int i=0;i<n;i++)
		p[i]=(d1[i]+d2[i])%mod*1ll*power(2,mod-2)%mod;
}
void Tan(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d2[i]=0;
	Sin(f,c1,n);Cos(f,c2,n);
	inv(c2,d2,n);
	int len=1;for(len=1;len<=n;len*=2);len*=2;
	NTT(c1,1,len);NTT(d2,1,len);
	for(int i=0;i<len;i++)
		p[i]=1ll*c1[i]*d2[i]%mod;
	NTT(p,-1,len);
	for(int i=n;i<len;i++)
		p[i]=0;
}

多项式反三角函数

求出(arcsin f(x),arccos f(x),arctan f(x) pmod {x^n})
我们对它们求导再积分

[(arcsin x)'=frac{1}{sqrt{1-x^2}} ]

[arcsin x= int frac{1}{sqrt{1-x^2}}dx ]

[(arccos x)'=-frac{1}{sqrt{1-x^2}} ]

[arccos x=-int frac{1}{sqrt{1-x^2}}dx ]

[(arctan x)'=frac{1}{x^2+1} ]

[arctan x =int frac{1}{x^2+1}dx ]

代入(f(x))得

[arcsin f(x)= int frac{f'(x)}{sqrt{1-f^2(x)}}dx ]

[arccos f(x)=-int frac{f'(x)}{sqrt{1-f^2(x)}}dx ]

[arctan f(x) =int frac{f'(x)}{f^2(x)+1}dx ]

模板:https://www.luogu.com.cn/problem/P5265

点击查看代码

void Arcsin(int *f,int *p,int n) {
	int len=1;for(len=1; len<=n; len*=2);len*=2;
	for(int i=0; i<len; i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0; i<n; i++)
		c1[i]=f[i];
	Der(c1,c2,n);NTT(c1,1,len);
	for(int i=0; i<len; i++)
		c1[i]=1ll*c1[i]*c1[i]%mod;
	NTT(c1,-1,len);
	for(int i=n; i<len; i++)
		c1[i]=0;
	c1[0]=(1-c1[0]+mod)%mod;
	for(int i=1; i<n; i++)
		c1[i]=(mod-c1[i])%mod;
	Sqrt(c1,d1,n);
	for(int i=n; i<len; i++)
		d1[i]=0;
	inv(d1,d2,n);
	for(int i=n; i<len; i++)
		d2[i]=0;
	NTT(c2,1,len);NTT(d2,1,len);
	for(int i=0; i<len; i++)
		d1[i]=1ll*c2[i]*d2[i]%mod;
	NTT(d1,-1,len);
	for(int i=n; i<len; i++)
		d1[i]=0;
	Int(d1,p,n);
}
void Arccos(int *f,int *p,int n) {
	Arcsin(f,p,n);
	for(int i=0; i<n; i++)
		p[i]=mod-p[i];
}
void Arctan(int *f,int *p,int n) {
	int len=1;
	for(len=1; len<=n; len*=2);
	len*=2;
	for(int i=0; i<len; i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0; i<n; i++)
		c1[i]=f[i];
	Der(c1,c2,n);NTT(c1,1,len);
	for(int i=0; i<len; i++)
		c1[i]=1ll*c1[i]*c1[i]%mod;
	NTT(c1,-1,len);
	c1[0]=(c1[0]+1)%mod;
	for(int i=n; i<len; i++)
		c1[i]=0;
	inv(c1,d1,n);NTT(d1,1,len);NTT(c2,1,len);
	for(int i=0; i<len; i++)
		d2[i]=1ll*c2[i]*d1[i]%mod;
	NTT(d2,-1,len);
	for(int i=n; i<len; i++)
		d2[i]=0;
	Int(d2,p,n);
}

提醒：多项式代码非常玄学，较为难调，调代码的过程中请注意多项式长度问题和数组清零问题。

多项式大集合:

点击查看代码

#include<bits/stdc++.h>
#define mod 1004535809
#define g 3
#define I 86583718
#define maxn 1000000
using namespace std;
int r[maxn],c[maxn],d[maxn],e[maxn],ft[maxn],gt[maxn],iv[maxn],tt[maxn];
int c1[maxn],c2[maxn],d1[maxn],d2[maxn];
inline int read() {
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9') {
		if(ch=='-')
			f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9') {
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*f;
}
int power(int x,long long y) {
	int ans=1,lazy=x;
	while(y) {
		if(y&1)
			ans=1ll*ans*lazy%mod;
		lazy=1ll*lazy*lazy%mod;
		y>>=1;
	}
	return ans;
}
int ff(int x)
{
	int ans=1;
	for(int i=1;i<=x;i++)
		ans=1ll*ans*i%mod;
}
int calc(int r1,int r2,int y,int q)
{
	int ans1=1,ans2=0,last1=r1,last2=r2;
	while(y)
	{
		int t1=last1,t2=last2,p1=ans1,p2=ans2;
		if(y&1)
			ans1=(1ll*p1*t1%mod+1ll*p2*t2%mod*q%mod)%mod,ans2=(1ll*p1*t2%mod+1ll*p2*t1%mod)%mod;
		last1=(1ll*t1*t1%mod+1ll*q*t2%mod*t2%mod)%mod,last2=(1ll*t1*t2%mod+1ll*t2*t1%mod)%mod;
		y>>=1;
	}
	return ans1;
}
int finda(int n)
{
	for(int i=0;;i++)
	{
		int t=(1ll*i*i%mod-n+mod)%mod;
		if(power(t,(mod-1)/2)==mod-1)
			return i;
	}
}	
int getx2(int n)
{
	if(n==0)
		return 0;
	int a=finda(n);
	int ans=calc(a,1,(mod+1)/2,(1ll*a*a%mod-n+mod)%mod);
	return min(ans,mod-ans);
}
void NTT(int *p,int op,int n) {
	int l=log2(n);
	r[0]=0;
	for(int i=0; i<n; i++)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0; i<n; i++)
		if(i<r[i])
			swap(p[i],p[r[i]]);
	for(int i=1; i<n; i*=2) {
		int tmp=power(g,(mod-1)/(i*2));
		if(op==-1)
			tmp=power(tmp,mod-2);
		for(int j=0; j<n; j+=i*2) {
			int w=1,t1,t2;
			for(int k=0; k<i; k++,w=1ll*w*tmp%mod)
				t1=p[j+k],t2=1ll*p[j+k+i]*w%mod,p[k+j]=(0ll+t1+t2+mod)%mod,p[k+j+i]=(0ll+t1-t2+mod)%mod;
		}
	}
	if(op==-1) {
		int invn=power(n,mod-2);
		for(int i=0; i<n; i++)
			p[i]=1ll*p[i]*invn%mod;
	}
}
void inv(int *f,int *p,int m) {
	if(m==1)
		return p[0]=power(f[0],mod-2),void();
	inv(f,p,(m+1)>>1);
	int n;
	for(n=1; n<=m; n*=2);
	n*=2;
	for(int i=0; i<m; i++)
		c[i]=f[i];
	for(int i=m; i<n; i++)
		c[i]=0;
	NTT(c,1,n);
	NTT(p,1,n);
	for(int i=0; i<n; i++)
		p[i]=1ll*(2ll-1ll*c[i]*p[i]%mod+mod)%mod*p[i]%mod;
	NTT(p,-1,n);
	for(int i=m; i<n; i++)
		p[i]=0;
}
void Sqrt(int *f,int *p,int m) {
	if(m==1)
		return p[0]=getx2(f[0]),void();
	int n;
	for(n=1; n<=m; n*=2);
	n*=2;
	for(int i=0; i<n; i++)
		p[i]=0;
	Sqrt(f,p,(m+1)>>1);
	for(int i=0; i<=m; i++)
		d[i]=f[i];
	for(int i=m; i<n; i++)
		d[i]=0;
	for(int i=0; i<n; i++)
		iv[i]=0;
	inv(p,iv,m);
	NTT(d,1,n);
	NTT(p,1,n);
	NTT(iv,1,n);
	for(int i=0; i<n; i++)
		p[i]=1ll*(0ll+p[i]+1ll*d[i]*iv[i]%mod)%mod*power(2,mod-2)%mod;
	NTT(p,-1,n);
	for(int i=0; i<n; i++)
		p[i]=1ll*p[i];
	for(int i=m; i<n; i++)
		p[i]=0;
}
void turn(int *f,int *p,int n) {
	for(int i=0; i<=n; i++)
		p[i]=f[n-i];
}
void div(int *F,int *G,int *Q,int *R,int n,int m) {
	int sum=n+m,len=1;
	for(len=1; len<=sum; len*=2);
	len*=2;
	turn(F,ft,n);
	turn(G,gt,m);
	for(int i=n-m+1; i<=len; i++)
		ft[i]=0,gt[i]=0;
	for(int i=0; i<=n; i++)
		d[i]=0,iv[i]=0;
	inv(gt,iv,n-m+1);
	NTT(ft,1,len);
	NTT(iv,1,len);
	for(int i=0; i<len; i++)
		d[i]=1ll*ft[i]*iv[i]%mod;
	NTT(d,-1,len);
	for(int i=n-m+1; i<=len; i++)
		d[i]=0;
	turn(d,Q,n-m);
	for(int i=0; i<=len; i++)
		ft[i]=F[i];
	for(int i=0; i<=len; i++)
		gt[i]=G[i];
	NTT(ft,1,len);
	NTT(gt,1,len);
	NTT(Q,1,len);
	for(int i=0; i<=len; i++)
		R[i]=(0ll+ft[i]-1ll*gt[i]*Q[i]%mod+mod)%mod;
	NTT(Q,-1,len);
	NTT(R,-1,len);
	for(int i=m; i<=len; i++)
		R[i]=0;
	for(int i=n-m+1; i<=len; i++)
		Q[i]=0;
}
void Der(int *f,int *p,int n) {
	for(int i=0; i<n; i++)
		p[i]=1ll*f[i+1]*(i+1)%mod;
	p[n]=0;
}
void Int(int *f,int *p,int n) {
	p[0]=0;
	for(int i=1; i<=n+1; i++)
		p[i]=1ll*f[i-1]*power(i,mod-2)%mod;
}
void Ln(int *f,int *p,int n) {
	int len=1;for(len=1; len<=n; len*=2);len*=2;
	for(int i=0; i<len; i++)
		ft[i]=gt[i]=0;
	Der(f,ft,n);
	inv(f,gt,n+1);
	NTT(ft,1,len);
	NTT(gt,1,len);
	for(int i=0; i<len; i++)
		ft[i]=1ll*ft[i]*gt[i]%mod;
	NTT(ft,-1,len);
	for(int i=n; i<len; i++)
		ft[i]=0;
	Int(ft,p,n-1);
	for(int i=n; i<len; i++)
		p[i]=0;
}
void Exp(int *f,int *p,int n) {
	if(n==1)
		return p[0]=1,void();
	int len=1;for(len=1; len<=n; len*=2);len*=2;
	for(int i=0; i<=len; i++)
		d[i]=0,p[i]=0,gt[i]=0;
	Exp(f,p,(n+1)>>1);
	for(int i=0; i<len; i++)
		gt[i]=0;
	Ln(p,gt,n);
	for(int i=n; i<len; i++)
		d[i]=0;
	for(int i=0; i<n; i++)
		d[i]=f[i];
	NTT(p,1,len);NTT(gt,1,len);NTT(d,1,len);
	for(int i=0; i<len; i++)
		p[i]=1ll*(1ll-gt[i]+d[i]+mod)%mod*1ll*p[i]%mod;
	NTT(p,-1,len);
	for(int i=n; i<len; i++)
		p[i]=0;
}
void Power(int *f,int *p,int n,int k1,int k2,int ov)
{
	int num0=0;
	for(int i=0;i<n;i++)
	{
		if(f[i]==0)
			num0++;
		else
			break;
	}
	if(f[0]==0&&ov)
		return ;
	if(1ll*num0*k1>=1ll*n)
		return ;
	for(int i=0;i<n;i++)
		e[i]=f[i+num0];
	int nop=e[0],ivnop=power(nop,mod-2),pwnop=power(nop,k2);
	for(int i=0;i<n;i++)
		e[i]=1ll*e[i]*ivnop%mod;
	Ln(e,p,n);
	for(int i=0;i<n;i++)
		e[i]=1ll*k1*p[i]%mod;
	for(int i=0;i<n;i++)
		tt[i]=0;
	Exp(e,tt,n);
	for(int i=k1*num0;i<n;i++)
		p[i]=1ll*tt[i-num0*k1]*pwnop%mod;
	for(int i=0;i<num0*k1;i++)
		p[i]=0;
}
void Sin(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=1ll*f[i]*I%mod,c2[i]=-c1[i];
	Exp(c1,d1,n);Exp(c2,d2,n);
	for(int i=0;i<n;i++)
		p[i]=(d1[i]-d2[i]+mod)%mod*1ll*power(2*I,mod-2)%mod;
}
void Cos(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=1ll*f[i]*I%mod,c2[i]=-c1[i];
	Exp(c1,d1,n);Exp(c2,d2,n);
	for(int i=0;i<n;i++)
		p[i]=(d1[i]+d2[i])%mod*1ll*power(2,mod-2)%mod;
}
void Tan(int *f,int *p,int n)
{
	for(int i=0;i<n*4;i++)
		c1[i]=c2[i]=d2[i]=0;
	Sin(f,c1,n);Cos(f,c2,n);
	inv(c2,d2,n);
	int len=1;for(len=1;len<=n;len*=2);len*=2;
	NTT(c1,1,len);NTT(d2,1,len);
	for(int i=0;i<len;i++)
		p[i]=1ll*c1[i]*d2[i]%mod;
	NTT(p,-1,len);
	for(int i=n;i<len;i++)
		p[i]=0;
}
void Arcsin(int *f,int *p,int n)
{
	int len=1;for(len=1;len<=n;len*=2);len*=2;
	for(int i=0;i<len;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=f[i];
	Der(c1,c2,n);
	NTT(c1,1,len);
	for(int i=0;i<len;i++)
		c1[i]=1ll*c1[i]*c1[i]%mod;
	NTT(c1,-1,len);
	for(int i=n;i<len;i++)
		c1[i]=0;
	c1[0]=(1-c1[0]+mod)%mod;
	for(int i=1;i<n;i++)
		c1[i]=(mod-c1[i])%mod;
	Sqrt(c1,d1,n);
	for(int i=n;i<len;i++)
		d1[i]=0;
	inv(d1,d2,n);
	for(int i=n;i<len;i++)
		d2[i]=0;
	NTT(c2,1,len);NTT(d2,1,len);
	for(int i=0;i<len;i++)
		d1[i]=1ll*c2[i]*d2[i]%mod;
	NTT(d1,-1,len);
	for(int i=n;i<len;i++)
		d1[i]=0;
	Int(d1,p,n);
}
void Arccos(int *f,int *p,int n)
{
	Arcsin(f,p,n);
	for(int i=0;i<n;i++)
		p[i]=mod-p[i];
}
void Arctan(int *f,int *p,int n)
{
	int len=1;for(len=1;len<=n;len*=2);len*=2;
	for(int i=0;i<len;i++)
		c1[i]=c2[i]=d1[i]=d2[i]=0;
	for(int i=0;i<n;i++)
		c1[i]=f[i];
	Der(c1,c2,n);
	NTT(c1,1,len);
	for(int i=0;i<len;i++)
		c1[i]=1ll*c1[i]*c1[i]%mod;
	NTT(c1,-1,len);
	c1[0]=(c1[0]+1)%mod;
	for(int i=n;i<len;i++)
		c1[i]=0;
	inv(c1,d1,n);
	NTT(d1,1,len);NTT(c2,1,len);
	for(int i=0;i<len;i++)
		d2[i]=1ll*c2[i]*d1[i]%mod;
	NTT(d2,-1,len);
	for(int i=n;i<len;i++)
		d2[i]=0;
	Int(d2,p,n);
}
int main() {
	return 0;
}

例题与练习

例题：(ZJOI2014)力

我们来推式子。

[E_j=sum_{i=1}^{j-1} frac{q_i}{(i-j)^2}-sum_{i=j+1}^{n} frac{q_i}{(i-j)^2} ]

前后求和加上第j项不影响答案

[E_j=sum_{i=1}^{j} frac{q_i}{(i-j)^2}-sum_{i=j}^{n} frac{q_i}{(i-j)^2} ]

我们设(T_i=frac{1}{i^2})

[E_j=sum_{i=1}^{j} q_i imes T_{i-j}-sum_{i=j}^{n} q_i imes T_{i-j} ]

可以看到，前半部分已经化成一个卷积的形式，后半部分还要进一步转化
设(P_j=sum_{i=1}^{j} q_i imes T_{i-j},Q_j=sum_{i=j}^{n} q_i imes T_{i-j})

[E_j=P_j-Q_j ]

接下来重点转化(Q_j)

[Q_j=sum_{i=0}^{n-j} T_i imes q_{i+j} ]

我们设(x=n-j),(q'_j=q_{n-j})

[Q_j=sum_{i=0}^x T_i q'_{x-i} ]

于是这边也化成一个卷积的形式了

点击查看代码

int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		cin>>q[i];
		a1[i]=q[i],a2[i]=1.0/i/i,a3[n-i]=q[i];
	}
	for(len=1;len<=n;len*=2);len*=2;
	FFT(a1,1,len);FFT(a2,1,len);FFT(a3,1,len);
	for(int i=0;i<len;i++)
		P[i]=a1[i]*a2[i],Q[i]=a2[i]*a3[i];
	FFT(P,-1,len);FFT(Q,-1,len);
	for(int i=1;i<=n;i++)
		printf("%.3lf
",P[i].real()-Q[n-i].real());
	
	return 0;	
}

例题：城市规划

简要题意：求(n)个点的有标号无向连通图的数量。
这是一个计数Dp题，它的弱化版我们以前是做过的。
我们设(f_i)表示(i)个点的无向连通图的数量，(g_i)表示(i)个点的无向图的数量。

[g_i=sum_{i=1}^nf_i imes C_{n-1}^{i-1} imes g_{n-i} ]

将(g_i)和组合数带入得

[2^{frac{n(n-1)}{2}}=sum_{i=1}^nf_i imes frac{(n-1)!}{(n-i)!(i-1)!} imes 2^{frac{(n-i)(n-i-1)}{2}} ]

[frac{2^{frac{n(n-1)}{2}}}{(n-1)!}= sum_{i-1}^n frac{f_i}{(i-1)!} imes frac{2^{frac{(n-i)(n-i-1)}{2}}}{(n-i)!} ]

观察到右边是一个卷积的形式。
我们设

[F= sum_{i=1}^{+ infty} frac{f_i}{(i-1)!} x^i ]

[G= sum_{i=1}^{+ infty} frac{2^{frac{i(i-1)}{2}}}{(i-1)!} x^i ]

[H= sum_{i=1}^{+ infty} frac{2^frac{i(i-1)}{2}}{i!} x^i ]

于是就有

[G=F imes H ]

[F=G imes H^{-1} ]

多项式求逆即可。
时间复杂度:(O(nlogn))

点击查看代码

int main() {
	scanf("%d",&n);
	fact[0]=1;
	for(int i=1;i<=n;i++)
		fact[i]=1ll*fact[i-1]*i%mod;
	t3[0]=1;
	for(int i=1;i<=n;i++)
		t2[i]=1ll*power(2,1ll*i*(i-1)/2)*power(fact[i-1],mod-2)%mod;
	for(int i=1;i<=n;i++)
		t3[i]=1ll*power(2,1ll*i*(i-1)/2)*power(fact[i],mod-2)%mod;
	inv(t3,t4,n);
	int len=1;for(len=1;len<=n;len*=2);len*=2;
	NTT(t2,1,len);NTT(t4,1,len);
	for(int i=0;i<len;i++)
		t1[i]=1ll*t2[i]*t4[i]%mod;
	NTT(t1,-1,len);
	cout<<1ll*t1[n]*fact[n-1]%mod<<endl;

	return 0;
}

这题也可以利用多项式ln的方法求解，留给学有余力的同学思考。

例题：The Child and Binary Tree

设(g_i)表示(i)是否在权值集合中(在为1，不在为0)，(f_i)表示权值为(i)的神犇二叉树的个数。

[f_n=sum_{i=1}^ng_isum_{j=1}^{n-i}f_j imes f_{n-i-j} ]

[f_n=sum_{i=1}^nsum_{j=1}^{n-i}g_if_j imes f_{n-i-j} ]

不难发现这是三个多项式卷积的形式
我们设(G)为(g)的生成函数(即(G=sum_{i=1}^n g_i x^i))

[F=F^2G+1 ]

加一是因为有空树。
我们来解一下这个方程。

[F^2G-F+1=0 ]

[F^2G^2-FG+G=0 ]

[(GF-frac{1}{2})^2=frac{1}{4}-G ]

[GF=frac{1}{2} pm sqrt{frac{1}{4}-G} ]

此时我们很懵B，不知道选正号还是负号。
我们知道，左边的式子，由于(G)的常数项为(0)，所以(GF)的常数项一定为(0)。
所以右边一坨东西的常数一定为(0)，根号里面的东西开根之后为(frac{1}{2})，所以应该取负号。

[GF=frac{1-sqrt{1-4G}}{2} ]

[GF=frac{2G}{sqrt{1-4G}+1} ]

[G(F-frac{2}{sqrt{1-4G}+1})=0 ]

由于(G eq 0)

[F=frac{2}{sqrt{1-4G}+1} ]

多项式求逆+多项式开根即可。

点击查看代码

int main() {
	scanf("%d%d",&m,&n);
	for(int i=1;i<=m;i++)
	{
		scanf("%d",&x);
		a1[x]=mod-4; 
	}
	a1[0]++;
	Sqrt(a1,a2,n+1);
	a2[0]=(a2[0]+1)%mod; 
	inv(a2,a3,n+1);
	for(int i=1;i<=n;i++)
		printf("%d
",2ll*a3[i]%mod);
	
	return 0;
}

例题：[差分与前缀和](https://www.luogu.com.cn/problem/P5488

我们设(A)为(a)序列的生成函数。
求前缀和也就是给多项式(A)乘上(1+x+x^2+x^3+...=frac{1}{1-x})
求(k)维前缀和也就是乘上(frac{1}{(1-x)^k})
求差分也就是给多项式(A)乘上((1-x))，求(k)维差分也就是乘上((1-x)^k)
多项式快速幂+多项式求逆即可。

点击查看代码

int main() {
	scanf("%d",&n);
	scanf("%s",k+1);
	scanf("%d",&op);
	int len=strlen(k+1);
	for(int i=1;i<=len;i++)
	{
		if(1ll*k1*10>=mod)
			ok=1;
		k1=(10ll*k1+k[i]-'0')%mod,k2=(10ll*k2+k[i]-'0')%(mod-1);
	}
	b[0]=1,b[1]=mod-1;
	for(int i=0;i<n;i++)
		scanf("%d",&a[i]);
	Power(b,A,n,k1,k2,ok);
	inv(A,B,n);
	for(len=1;len<=n;len*=2);len*=2;
	NTT(a,1,len);NTT(A,1,len);NTT(B,1,len);
	for(int i=0;i<len;i++)
		A1[i]=1ll*a[i]*A[i]%mod,A2[i]=1ll*a[i]*B[i]%mod;
	NTT(A1,-1,len);NTT(A2,-1,len);
	if(op==0)
	{
		for(int i=0;i<n;i++)
			printf("%d ",A2[i]);
		puts("");
	}
	else
	{
		for(int i=0;i<n;i++)
			printf("%d ",A1[i]);
		puts("");
	}	
	
	return 0;
}

练习题：难度不大，可以一做
https://www.luogu.com.cn/problem/P5641
提示：卷积+NTT
[AH2017/HNOI2017]礼物
提示：推式子+卷积
有标号二分图计数
提示：EGF+多项式开根

思考题：题目较毒瘤，留给学(xian)有(zhe)余(mei)力(shi)的同学解决(不建议做)
[HNOI2019]白兔之舞
解法：单位根反演+多项式+MTT
遗忘的集合
解法：生成函数+多项式ln+莫比乌斯反演
[WC2019]数树
解法：prufer+树形dp+多项式