多项式半家桶

多项式乘法

FFT

因为有浮点数参与运算，所以可能会出现精度的问题

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const double pi=acos(-1);
const int maxn=4e6+5;
struct fs{
	double x,y;
	friend fs operator + (rg const fs& A,rg const fs& B){
		return (fs){A.x+B.x,A.y+B.y};
	}
	friend fs operator - (rg const fs& A,rg const fs& B){
		return (fs){A.x-B.x,A.y-B.y};
	}
	friend fs operator * (rg const fs& A,rg const fs& B){
		return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
	}
	friend fs operator / (rg const fs& A,rg const double& B){
		return (fs){A.x/B,A.y/B};
	}
}a[maxn],b[maxn],w[25][maxn];
int wz[maxn],n,m;
void fft(rg fs A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
				A[j+k]=x+y;
				A[j+k+len]=x-y;
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		for(rg int i=0;i<lim;i++) A[i]=A[i]/(double)lim;
	}
}
int main(){
	n=read(),m=read();
	for(rg int i=0;i<=n;i++) a[i].x=read();
	for(rg int i=0;i<=m;i++) b[i].x=read();
	rg int lim=1,bit=0;
	for(;lim<=n+m;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++){
		wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	}
	for(rg int i=0,len=1;len<lim;len<<=1,i++){
		for(rg int j=0;j<len;j++){
			w[i][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
		}
	}
	fft(a,lim,1),fft(b,lim,1);
	for(rg int i=0;i<lim;i++) a[i]=a[i]*b[i];
	fft(a,lim,-1);
	for(rg int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x+0.5));
	printf("\n");
	return 0;
}

NTT

用原根代替复数

但必须保证模数 $p = a \cdot 2^k + 1$

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=4e6+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,m,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);	
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
				A[j+k]=addmod(x,y);
				A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++){
			A[i]=mulmod(A[i],ny);
		}
	}
}
int main(){
	n=read(),m=read();
	for(rg int i=0;i<=n;i++) a[i]=read();
	for(rg int i=0;i<=m;i++) b[i]=read();
	rg int lim=1,bit=0;
	for(;lim<=n+m;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0,k=1;k<lim;k<<=1,i++){
		w[i][0]=1;
		w[i][1]=ksm(G,(mod-1)/(k<<1));
		for(rg int j=2;j<k;j++){
			w[i][j]=mulmod(w[i][j-1],w[i][1]);
		}
	}
	ntt(a,lim,1),ntt(b,lim,1);
	for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],b[i]);
	ntt(a,lim,-1);
	for(rg int i=0;i<=n+m;i++) printf("%d ",a[i]);
	printf("\n");
	return 0;
}

$FFT$ 和 $NTT$ 预处理单位根要快不少

预处理单位根有两种写法，一种是 $O(nlogn)$ 的，另一种是 $O(n)$ 的

第一种因为访问内存连续实测更快

任意模数多项式乘法（MTT）

当多项式乘法进行取模的数并不是我们喜欢的模数时，就要用到 $MTT$

实现的方法有很多，这里介绍一种只需要 $2$ 次 $DFT$ 和 $2$ 次 $IDFT$ 的优秀做法

首先，我们对原来的多项式拆系数

令$k=sqrt(mod)$，则

\[\begin{aligned} A(x)=kA_0(x)+A_1(x)\\ B(x)=kB_0(x)+B_1(x)\end{aligned} \]

即拆成了 $A(x)/k$ 和 $A(x)\%k$ 两部分

那么

\[\begin{aligned} A(x)B(x)&=(kA_0(x)+A_1(x))(kB_0(x)+B_1(x))\\ &=k^2A_0(x)B_0(x)+A_1(x)B_1(x)+kA_0(x)B_1(x)+kA_1(x)B_0(x) \end{aligned} \]

这样就有一个$4$ 次 $DFT$ 和 $4$ 次 $IDFT$ 的暴力做法

考虑优化，设

\[\begin{aligned} P(x)=A_0(x)+iA_1(x)\\ Q(x)=A_0(x)-iA_1(x)\end{aligned} \]

由于 $A, B$ 的虚部都为 $0$

$P,Q$ 的每一项系数都互为共轭，同样每一个点值也互为共轭

那么只需对 $P$ 做一次 $DFT$ ，就可以通过共轭 $O(n)$ 求出 $Q$ 的点值表示

具体地说，在点值表示中 $P$ 的第 $k$ 项与 $Q$ 的第 $n-k$ 项是共轭

要特判下标为 $0$ 的情况

求出了 $P$ 和 $Q$ ，我们就能用二元一次方程的知识推出 $A_0$ 和 $A_1$

对于 $B$ 也是同理

此时，我们进行了两次 $DFT$ 得到了所有多项式的点值表示

下面就是怎么把它们变成系数表示

仍然借鉴上面的思路，构造两个多项式

$P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x)$

$Q(x)=A_0(x)B_1(x)+iA_1(x)B_1(x)$

通过已知的点值求出此时 $P, Q$ 的点值

然后分别对 $P, Q$ 做 $IDFT$

由于 $A_0 B_0,A_0 B_1,A_1 B_0,A_1 B_1$这四个多项式卷起来后的系数表示中虚部一定为 $0$

那么此时 $P$ 的实部和虚部就分别为 $A_0(x) B_0(x)$和 $A_1(x) B_0(x)$

同样 $Q$ 的实部和虚部就分别为 $A_0(x) B_1(x)$ 和 $A_1(x) B_1(x)$。

这样，我们又进行了两次 $IDFT$ 得到了所有多项式的系数表示

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=4e5+5;
const double pi=acos(-1.0);
struct fs{
	double x,y;
	friend fs operator +(rg const fs& A,rg const fs& B){
		return (fs){A.x+B.x,A.y+B.y};
	}
	friend fs operator -(rg const fs& A,rg const fs& B){
		return (fs){A.x-B.x,A.y-B.y};
	}
	friend fs operator *(rg const fs& A,rg const fs& B){
		return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
	}
	friend fs operator /(rg const fs& A,rg const double& B){
		return (fs){A.x/B,A.y/B};
	}
}w[22][maxn],a0[maxn],b0[maxn],a1[maxn],b1[maxn],p[maxn],q[maxn];
const fs I=(fs){0,1};
int wz[maxn];
void fft(rg fs A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
				A[j+k]=x+y,A[j+k+len]=x-y;
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		for(rg int i=0;i<lim;i++) A[i]=A[i]/lim;
	}
}
void mtt(rg fs A0[],rg fs A1[],rg int lim,rg int len){
	for(rg int i=0;i<lim;i++) p[i].x=p[i].y=q[i].x=q[i].y=0;
	for(rg int i=0;i<=len;i++) p[i].x=A0[i].x,p[i].y=A1[i].x;
	fft(p,lim,1);
	q[0]=(fs){p[0].x,-p[0].y};
	for(rg int i=1;i<lim;i++) q[i]=(fs){p[lim-i].x,-p[lim-i].y};
	for(rg int i=0;i<lim;i++){
		A0[i]=(p[i]+q[i])/2.0;
		A1[i]=(q[i]-p[i])*I/2.0;
	}
}
int n,m,mod,blo,a[maxn],b[maxn];
long long getval(rg double val){
	if(val>=0) return (long long)(val+0.5)%mod;
	else return (long long)(val-0.5)%mod;
}
int main(){
	n=read(),m=read(),mod=read();
	blo=sqrt(mod)+1;
	for(rg int i=0;i<=n;i++) a[i]=read()%mod;
	for(rg int i=0;i<=m;i++) b[i]=read()%mod;
	rg int lim=1,bit=0;
	for(;lim<=n+m;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0;j<len;j++) w[t0][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
	}
	for(rg int i=0;i<=n;i++) a0[i].x=a[i]%blo,a1[i].x=a[i]/blo;
	for(rg int i=0;i<=m;i++) b0[i].x=b[i]%blo,b1[i].x=b[i]/blo;
	mtt(a0,a1,lim,n),mtt(b0,b1,lim,m);
	for(rg int i=0;i<lim;i++) p[i].x=p[i].y=q[i].x=q[i].y=0;
	for(rg int i=0;i<lim;i++){
		p[i]=a0[i]*b0[i]+I*a1[i]*b0[i];
		q[i]=a0[i]*b1[i]+I*a1[i]*b1[i];
	}
	fft(p,lim,-1),fft(q,lim,-1);
	for(rg int i=0;i<=n+m;i++){
		rg long long nans=getval(q[i].y)*blo*blo%mod+(getval(q[i].x)+getval(p[i].y))*blo%mod+getval(p[i].x);
		printf("%lld ",nans%mod);
	}
	printf("\n");
	return 0;
}

多项式泰勒展开

$F(x)=\sum_{n=0}^{\infty}\frac{F^{(n)}(x_0)(x-x_0)^n}{n!}(F^{(n)}(x)为F(x)的n阶导)$

怎样更好地理解并记忆泰勒展开式

多项式牛顿迭代

设 $F(G(x)) \equiv 0 \pmod{x^n}$

已知 $F(G_0(x)) \equiv 0 \pmod{x^\frac{n}{2}}$

对于 $F(G(x))$ 进行泰勒展开

$F(G(x))=\sum_{n=0}^{\infty}\frac{F^{(n)}(G_0(x))(G(x)-G_0(x))^n}{n!} \pmod{x^n}$

因为 $G(x)-G_0(x)$ 在 $\frac{n}{2}$ 次项之前的次数都是相同的，所以作差就变成了 $0$

平方之后 $n$ 次项之前都变为了 $0$，在模 $x^n$ 次方意义下是 $0$

之后再做乘法还是 $0$

所以只需要展开前两项就可以了

即 $F(G_0(x))+F'(G_0(x))(G(x)-G_0(x)) \equiv 0 \pmod{x^n}$

所以 $G(x)=G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))}$

多项式乘法逆

两种方法

方法一

设 $F(x) G(x) \equiv 1 \pmod {x^n}$

已知 $F(x) G_0(x) \equiv 1 \pmod {x^{\frac{n}{2}}}$

两式相减，可以得到 $F(x)G(x)-F(x)G_0(x) \equiv 0 \pmod {x^{\frac{n}{2}}}$

$F(x)(G(x)-G_0(x)) \equiv 0 \pmod {x^{\frac{n}{2}}}$

左边的部分肯定不为 $0$

所以一定有 $G(x)-G_0(x) \equiv 0 \pmod {x^{\frac{n}{2}}}$

两边同时平方

$G(x)^2+G_0(x)^2-2G(x)G_0(x) \equiv 0 \pmod {x^n}$

两边同时乘 $F(x)$

得到 $G(x)+F(x)G_0(x)^2-2G_0(x) \equiv 0 \pmod {x^n}$

即 $G(x)=2G_0(x)-F(x)G_0(x)^2 \equiv 0 \pmod {x^n}$

方法二

使用牛顿迭代

令 $H(x)=\frac{1}{x}-F(x)$

已知 $H(G_0(x)) \equiv 0 \pmod{x^\frac{n}{2}}$

求 $H(G(x)) \equiv 0 \pmod{x^n}$

套用牛顿迭代公式

$\begin {aligned} G(x) &= G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))} \\ &= G_0(x)-\frac{\frac{1}{G_0(x)}-F(x)}{-\frac{1}{G_0(x)^2}} \\ &= 2G_0(x)-F(x)G_0(x)^2 \end {aligned}$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=262145,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(B,lim,1),ntt(finv,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	for(rg int i=0;i<n;i++) a[i]=read();
	rg int lim=1;
	for(;lim<n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	getinv(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

多项式开根

也可以用两种方法推导

方法一

设 $G(x)^2 \equiv F(x) \pmod {x^n}$

已知 $G_0(x)^2 \equiv F(x)\pmod {x^{\frac{n}{2}}}$

两式相减，可以得到 $G(x)^2-G_0(x)^2 \equiv 0 \pmod {x^{\frac{n}{2}}}$

$(G(x)+G_0(x))(G(x)-G_0(x)) \equiv 0 \pmod {x^{\frac{n}{2}}}$

左边的部分肯定不为 $0$

所以一定有 $G(x)-G_0(x) \equiv 0 \pmod {x^{\frac{n}{2}}}$

两边同时平方

$G(x)^2+G_0(x)^2-2G(x)G_0(x) \equiv 0 \pmod {x^n}$

得到 $F(x)+G_0(x)^2-2G(x)G_0(x) \equiv 0 \pmod {x^n}$

即 $G(x) \equiv \frac{F(x)+G_0(x)^2}{2G_0(x)} \pmod {x^n}$

方法二

使用牛顿迭代

令 $H(x)=x^2-F(x)$

已知 $H(G_0(x)) \equiv 0 \pmod{x^\frac{n}{2}}$

求 $H(G(x)) \equiv 0 \pmod{x^n}$

套用牛顿迭代公式

$\begin {aligned} G(x) &= G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))} \\ &= G_0(x)-\frac{G_0(x)^2-F(x)}{2G_0(x)} \\ &= \frac{F(x)+G_0(x)^2}{2G_0(x)} \end {aligned}$

代码

普通版代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=262145,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		memset(B,0,maxn<<2);
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(B,lim,1),ntt(finv,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqr1[maxn],fsqr2[maxn];
void getsqr(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=1;
		return;
	}
	getsqr(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	getinv(B,fsqr1,deg);
	for(rg int i=0;i<deg;i++) fsqr2[i]=A[i];
	for(rg int i=deg;i<lim;i++) fsqr2[i]=0;
	ntt(B,lim,1),ntt(fsqr1,lim,1),ntt(fsqr2,lim,1);
	rg int ny=ksm(2,mod-2);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqr2[i],mulmod(B[i],B[i])),mulmod(ny,fsqr1[i]));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	for(rg int i=0;i<n;i++) a[i]=read();
	rg int lim=1;
	for(;lim<n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	getsqr(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

加强版代码(二次剩余)

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=262145,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		memset(B,0,maxn<<2);
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(B,lim,1),ntt(finv,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqr1[maxn],fsqr2[maxn],O;
struct fs{
	int x,y;
	fs(){}
	fs(rg int aa,rg int bb){
		x=aa,y=bb;
	}
	friend fs operator *(const fs& A,const fs& B){
		return fs(addmod(mulmod(A.x,B.x),mulmod(O,mulmod(A.y,B.y))),addmod(mulmod(A.x,B.y),mulmod(A.y,B.x)));
	}
};
fs fpow(rg fs ds,rg int zs){
	rg fs nans=fs(1,0);
	while(zs){
		if(zs&1) nans=nans*ds;
		ds=ds*ds;
		zs>>=1;
	}
	return nans;
}
int getans(rg int val){
	rg int now=0;
	while(1){
		now=rand()%mod+1;
		O=delmod(mulmod(now,now),val);
		if(ksm(O,(mod-1)>>1)==mod-1) break;
	}
	fs nans=fpow(fs(now,1),(mod+1)>>1);
	return std::min(nans.x,mod-nans.x);
}
void getsqr(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=getans(A[0]);
		return;
	}
	getsqr(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	getinv(B,fsqr1,deg);
	for(rg int i=0;i<deg;i++) fsqr2[i]=A[i];
	for(rg int i=deg;i<lim;i++) fsqr2[i]=0;
	ntt(B,lim,1),ntt(fsqr1,lim,1),ntt(fsqr2,lim,1);
	rg int ny=ksm(2,mod-2);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqr2[i],mulmod(B[i],B[i])),mulmod(ny,fsqr1[i]));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	for(rg int i=0;i<n;i++) a[i]=read();
	rg int lim=1;
	for(;lim<n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	getsqr(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

多项式求导

\[F(x)=\sum^{n}_{i=0}a_ix^i \]

\[F(x+dx)=\sum^{n}_{i=0}a_i(x+dx)^i \]

\[F(x+dx)=\sum^{n}_{i=0}a_i(\binom{i}{0}x^i+\binom{i}{1}x^{i-1}dx+\binom{i}{2}x^{i-2}dx^2+...+\binom{i}{i-1}xdx^{i-1}+\binom{i}{i}dx^i) \]

\[F(x+dx)-F(x)=\sum^{n}_{i=1}a_i(\binom{i}{1}x^{i-1}dx+\binom{i}{2}x^{i-2}dx^2+...+\binom{i}{i-1}xdx^{i-1}+\binom{i}{i}dx^i) \]

$F'(x)=\frac{F(x+dx)-F(x)}{dx}=\sum^{n}_{i=1}a_i\binom{i}{1}x^{i-1}$

因为 $dx$ 趋进于无穷小，所以后面的几项都变成了 $0$

$F'(x)=\sum^{n}_{i=0}a_{i + 1}(i+1)x^i$

很容易发现，其实就是求导之前的多项式系数整体左移一位再乘上一个 $i+1$，得到的就是它的导数。

代码

void getdao(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i-1]=mulmod(A[i],i);
	}
	B[deg-1]=0;
}

多项式积分

求导的逆运算

代码

void getfen(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i]=mulmod(A[i-1],ksm(i,mod-2));
	}
	B[0]=0;
}

多项式对数函数（多项式 ln）

对两边同时求导

$G(x)=F(A(x)),F(x)=ln(x)$

$G'(x)=F'(A(x))A'(x)=\frac{A'(x)}{A(x)}$

所以对 $A$ 求导求个逆再积一下分就算出了 $G$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=262145,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(finv,lim,1),ntt(B,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(finv[i],B[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i-1]=mulmod(i,A[i]);
	}
	B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i]=mulmod(A[i-1],ksm(i,mod-2));
	}
	B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<lim;i++) fln1[i]=fln2[i]=0;
	getinv(A,fln1,deg);
	getdao(A,fln2,deg);
	ntt(fln1,lim,1),ntt(fln2,lim,1);
	for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
	ntt(fln1,lim,-1);
	getfen(fln1,B,deg);
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	for(rg int i=0;i<n;i++) a[i]=read();
	rg int lim=1;
	for(;lim<n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][1],w[t0][i-1]);
	}
	getln(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

多项式指数函数（多项式 exp）

求 $G(x) \equiv e^{F(x)} \pmod{x^n}$

两边同时取对数

有 $ln(G(x)) \equiv F(x) \pmod{x^n}$

令 $H(x)=ln(x)-F(x)$

已知 $H(G_0(x)) \equiv 0 \pmod{x^\frac{n}{2}}$

求 $H(G(x)) \equiv 0 \pmod{x^n}$

套用牛顿迭代公式

$\begin {aligned} G(x) &= G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))} \\ &= G_0(x)-\frac{ln(G_0(x))-F(x)}{\frac{1}{G_0(x)}} \\ &= G_0(x)(1-\ln G_0(x)+F(x)) \end {aligned}$

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=262145,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(finv,lim,1),ntt(B,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(finv[i],B[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i-1]=mulmod(i,A[i]);
	}
	B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i]=mulmod(A[i-1],ksm(i,mod-2));
	}
	B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<lim;i++) fln1[i]=fln2[i]=0;
	getinv(A,fln1,deg);
	getdao(A,fln2,deg);
	ntt(fln1,lim,1),ntt(fln2,lim,1);
	for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
	ntt(fln1,lim,-1);
	getfen(fln1,B,deg);
}
int fexp1[maxn],fexp2[maxn];
void getexp(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=1;
		return;
	}
	getexp(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	getln(B,fexp1,deg);
	for(rg int i=0;i<deg;i++) fexp2[i]=A[i];
	for(rg int i=deg;i<lim;i++) fexp2[i]=0;
	ntt(B,lim,1),ntt(fexp1,lim,1),ntt(fexp2,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],addmod(delmod(1,fexp1[i]),fexp2[i]));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	for(rg int i=0;i<n;i++) a[i]=read();
	rg int lim=1;
	for(;lim<n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][1],w[t0][i-1]);
	}
	getexp(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

多项式快速幂

$G(x)=(F(x))^k\pmod{x^n}$

取下对数

$\ln G(x)=k\ln F(x)\pmod{x^n} $

再取一下指数

$G(x)=e^{k\ln F(x)}\pmod{x^n} $

代码

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);	
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
				A[j+k]=addmod(x,y);
				A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++){
			A[i]=mulmod(A[i],ny);
		}
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(B,lim,1),ntt(finv,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i-1]=mulmod(A[i],i);
	}
	B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
	for(rg int i=1;i<deg;i++){
		B[i]=mulmod(A[i-1],ksm(i,mod-2));
	}
	B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
	memset(fln1,0,maxn<<2);
	memset(fln2,0,maxn<<2);
	getdao(A,fln1,deg);
	getinv(A,fln2,deg);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	ntt(fln1,lim,1),ntt(fln2,lim,1);
	for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
	ntt(fln1,lim,-1);
	getfen(fln1,B,deg);
}
int fexp1[maxn],fexp2[maxn];
void getexp(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=1;
		return;
	}
	getexp(A,B,(deg+1)>>1);
	getln(B,fexp1,deg);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) fexp2[i]=A[i];
	for(rg int i=deg;i<lim;i++) fexp2[i]=0;
	ntt(fexp1,lim,1),ntt(fexp2,lim,1),ntt(B,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],addmod(delmod(1,fexp1[i]),fexp2[i]));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fmi[maxn],nowk;
void getmi(rg int A[],rg int B[],rg int deg){
	getln(A,fmi,n);
	for(rg int i=0;i<n;i++) fmi[i]=mulmod(fmi[i],nowk);
	getexp(fmi,B,n);
}
char s[maxn];
int main(){
	for(rg int i=0,k=1;k<maxn;k<<=1,i++){
		w[i][0]=1;
		w[i][1]=ksm(G,(mod-1)/(k<<1));
		for(rg int j=2;j<k;j++){
			w[i][j]=mulmod(w[i][j-1],w[i][1]);
		}
	}
	n=read();
	scanf("%s",s+1);
	rg int len=strlen(s+1);
	for(rg int i=1;i<=len;i++){
		nowk=addmod(s[i]-'0',mulmod(nowk,10));
	}
	for(rg int i=0;i<n;i++) a[i]=read();
	getmi(a,b,n);
	for(rg int i=0;i<n;i++) printf("%d ",b[i]);
	printf("\n");
	return 0;
}

多项式除法

$F_R(x)$ 等于 $F(x)$ 系数翻转得到的多项式

$F(x)=Q(x)G(x)+R(x)$

$F(\frac{1}{x})=Q(\frac{1}{x})G(\frac{1}{x})+R(\frac{1}{x})$

$x^n F(\frac{1}{x})=x^{n-m} Q(\frac{1}{x}) x^m G(\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}R_R(x)$

$F_R(x) \equiv Q_R(x)G_R(x)+x^{n-m+1}R_R(x)\pmod {x^{n-m+1}}$

$F_R(x) \equiv Q_R(x)G_R(x)\pmod {x^{n-m+1}}$

$Q_R(x) \equiv F_R(x)G_R^{-1}(x)\pmod {x^{n-m+1}}$

求一遍 $G_R$ 的逆，然后就可以利用多项式乘法求出 $Q$

根据

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

直接计算即可

代码

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);	
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
				A[j+k]=addmod(x,y);
				A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++){
			A[i]=mulmod(A[i],ny);
		}
	}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
	if(deg==1){
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	ntt(B,lim,1),ntt(finv,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int m,fdiv1[maxn],fdiv2[maxn],ansdiv1[maxn],ansdiv2[maxn];
void getdiv(rg int A[],rg int B[],rg int deg1,rg int deg2){
	std::reverse(A,A+deg1);
	std::reverse(B,B+deg2);
	for(rg int i=0;i<deg1-deg2+1;i++) fdiv2[i]=B[i];
	getinv(fdiv2,fdiv1,deg1-deg2+1);
	rg int lim=1,bit=0;
	for(;lim<(deg1-deg2+1)<<1;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0;i<deg1-deg2+1;i++) fdiv2[i]=A[i];
	ntt(fdiv1,lim,1),ntt(fdiv2,lim,1);
	for(rg int i=0;i<lim;i++) fdiv1[i]=mulmod(fdiv1[i],fdiv2[i]);
	ntt(fdiv1,lim,-1);
	for(rg int i=deg1-deg2+1;i<lim;i++) fdiv1[i]=0;
	memcpy(ansdiv1,fdiv1,sizeof(fdiv1));
	std::reverse(ansdiv1,ansdiv1+deg1-deg2+1);
	for(;lim<(deg2-1)<<1;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	std::reverse(A,A+deg1);
	std::reverse(B,B+deg2);
	for(rg int i=0;i<deg2-1;i++) fdiv1[i]=ansdiv1[i],fdiv2[i]=B[i];
	for(rg int i=deg2-1;i<lim;i++) fdiv1[i]=fdiv2[i]=0;
	ntt(fdiv1,lim,1),ntt(fdiv2,lim,1);
	for(rg int i=0;i<lim;i++) fdiv1[i]=mulmod(fdiv1[i],fdiv2[i]);
	ntt(fdiv1,lim,-1);
	for(rg int i=0;i<lim;i++) fdiv1[i]=delmod(A[i],fdiv1[i]);
	for(rg int i=deg2-1;i<lim;i++) fdiv1[i]=0;
	memcpy(ansdiv2,fdiv1,sizeof(fdiv1));
}
int main(){
	for(rg int i=0,k=1;k<maxn;k<<=1,i++){
		w[i][0]=1;
		w[i][1]=ksm(G,(mod-1)/(k<<1));
		for(rg int j=2;j<k;j++){
			w[i][j]=mulmod(w[i][j-1],w[i][1]);
		}
	}
	n=read(),m=read();
	for(rg int i=0;i<=n;i++) a[i]=read();
	for(rg int i=0;i<=m;i++) b[i]=read();
	getdiv(a,b,n+1,m+1);
	for(rg int i=0;i<n-m+1;i++) printf("%d ",ansdiv1[i]);
	printf("\n");
	for(rg int i=0;i<m;i++) printf("%d ",ansdiv2[i]);
	printf("\n");
	return 0;
}

posted @ 2021-04-07 21:35 liuchanglc 阅读(207) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

liuchanglc

谁终将声震人间,必长久深自缄默;谁终将点燃闪电,必长久如云漂泊

多项式半家桶

多项式乘法

FFT

代码

NTT

代码

任意模数多项式乘法（MTT）

代码

多项式泰勒展开

多项式牛顿迭代

多项式乘法逆

方法一

方法二

代码

多项式开根

方法一

方法二

代码

多项式求导

代码

多项式积分

代码

多项式对数函数（多项式 ln）

代码

多项式指数函数（多项式 exp）

代码

多项式快速幂

代码

多项式除法

代码

公告