多项式细节梳理&模板（多项式）

基础

现在的码风又变了。。。

FFT和NTT的板子

typedef complex<double> C;
const double PI=acos(-1);
void FFT(C*a,R op){
	for(R i=0;i<N;++i)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(R i=1;i<N;i<<=1){
		C wn=C(cos(PI/i),sin(PI/i)*op),w=1,t;
		for(R j=0;j<N;j+=i<<1,w=1)
			for(R k=j;k<j+i;++k,w*=wn)
				t=a[k+i]*w,a[k+i]=a[k]-t,a[k]+=t;
	}
}

const int YL=998244353;
LL Pow(LL b,R k=YL-2,LL a=1){
	for(;k;k>>=1,b=b*b%YL)
		if(k&1)a=a*b%YL;
	return a;
}
void NTT(R*a,R n,R op){
	for(R i=0;i<n;++i)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(R i=1;i<n;i<<=1){
		LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
		if(op==-1)wn=Pow(wn);
		for(R j=0;j<n;j+=i<<1,w=1)
			for(R k=j;k<j+i;++k,w=w*wn%YL)
				x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
	}
	if(op==-1){
		LL x=Pow(n);
		for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
	}
}

任意模数NTT

怎么看都觉得MTT给人感觉不好，用double巨佬说怕掉精度，用longdouble那常数想象一下～
然后就去学了三模NTT
参考Memory of Winter巨佬的题解
细节不多吧，首先要把\(998244353,1004535809,469762049\)背下来，它们都有一个原根为\(3\)。
然后手动两两合并同余方程，不要直接用CRT公式。
~~蒟蒻为了减少码量直接开二维数组了~~，实际上像Memory of Winter巨佬一样封装结构体会跑得更快。
upd:重写了一遍，二维数组又长又丑又慢qwq
洛谷P4245 【模板】任意模数NTT

#include<bits/stdc++.h>
#define LL long long
#define I inline
#define R register int
#define G if(++ip==ie)if(fread(ip=buf,1,SZ,stdin))
#define Wn(A) Pow(3,Mod((A-1)/(i<<1)*op,A-1),A)
using namespace std;
const int SZ=1<<19,N=1<<18,YL=1e9+7,A=998244353,B=1004535809,C=469762049;
char buf[SZ],*ie=buf+SZ,*ip=ie-1;
inline int in(){
	G;while(*ip<'-')G;
	R x=*ip&15;G;
	while(*ip>'-'){x*=10;x+=*ip&15;G;}
	return x;
}
int L,r[N];
I int Mod(R x,R YL){
	return x+(x>>31&YL);
}
I int Pow(LL b,R k,LL YL,LL a=1){
	for(;k;k>>=1,b=b*b%YL)
		if(k&1)a=a*b%YL;
	return a;
}
I int Inv(LL b,LL YL){
	return Pow(b%YL,YL-2,YL);
}
struct Z{
	int a,b,c;
	Z(){}
	Z(LL a):a(a%A),b(a%B),c(a%C){}
	Z(R a,R b,R c):a(a),b(b),c(c){}
	I Z operator!(){a+=a>>31&A,b+=b>>31&B,c+=c>>31&C;return*this;}
	I Z operator+(const Z&x){return!Z(a+x.a-A,b+x.b-B,c+x.c-C);}
	I Z operator-(const Z&x){return!Z(a-x.a  ,b-x.b  ,c-x.c  );}
	I Z operator*(const Z&x){return Z((LL)a*x.a%A,(LL)b*x.b%B,(LL)c*x.c%C);}
	I int CRT(LL YL){
		static LL AB=(LL)A*B%YL,I0=Inv(A,B),I1=Inv((LL)A*B,C),x;
		x=Mod(b-a,B)*I0%B*A+a;
		return(Mod(c-x%C,C)*I1%C*AB+x)%YL;
	}
}f[N],g[N];
void NTT(Z*a,R op){
	for(R i=0;i<L;++i)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(R i=1;i<L;i<<=1){
		Z wn(Wn(A),Wn(B),Wn(C)),w(1),x;
		for(R j=0;j<L;j+=i<<1,w=1)
			for(R k=j;k<j+i;++k,w=w*wn)
				x=w*a[k+i],a[k+i]=a[k]-x,a[k]=a[k]+x;
	}
	if(op==-1){
		Z w(Inv(L,A),Inv(L,B),Inv(L,C));
		for(R i=0;i<L;++i)a[i]=a[i]*w;
	}
}
int main(){
	R n=in(),m=in(),p=in();
	for(R i=0;i<=n;++i)f[i]=in();
	for(R i=0;i<=m;++i)g[i]=in();
	for(L=1;L<=n+m;L<<=1);
	for(R i=1;i<L;++i)r[i]=(r[i>>1]|L*(1&i))>>1;
	NTT(f,1);NTT(g,1);
	for(R i=0;i<L;++i)f[i]=f[i]*g[i];
	NTT(f,-1);
	for(R i=0;i<=n+m;++i)printf("%d ",f[i].CRT(p));
	return 0;
}

另外，三模NTT不能将三个及以上的多项式一次都乘起来，因为实际值域太大，不能保证模意义下的唯一性。
正确的做法是两个两个乘。
比如，任意模数多项式求逆的核心代码
洛谷P4239 【模板】多项式求逆（加强版）

void PolyInv(Z*a,Z*b,Z*a1,R m){
    if(m==1){b[0]=Inv(a[0].b,YL);return;}
    PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,12*m);
    Init(2*m);NTT(b,1);NTT(a1,1);
    for(R i=0;i<L;++i)a1[i]=b[i]*a1[i];
    NTT(a1,-1);
    for(R i=0;i<L;++i)a1[i]=Mod(-a1[i].CRT()+2*(i==0),YL);
    NTT(a1,1);
    for(R i=0;i<L;++i)b[i]=b[i]*a1[i];
    NTT(b,-1);memset(a1,0,12*L);memset(b+m,0,12*(L-m));
    for(R i=0;i<m;++i)b[i]=b[i].CRT();
}

多项式运算

分治乘法

若干个总项数不超过\(n\)的多项式的卷积，根据合并果子的方法，可以在不超过\(O(n\log^2n)\)的时间内完成。
有时候是若干个二项式相乘，还可以写成循环的形式，好处是减少Rader排序r数组的预处理次数。
第一次尝试这种写法：洛谷CF981H K Paths

void Solve(R x){
    R n=4*e[x].size();
    for(R i=0;i<n;i+=4)
        a[i]=1,a[i+1]=s[e[x][i>>2]];
    for(R i=4;i<n;i<<=1){
        for(R j=1;j<i;++j)
            r[j]=(r[j>>1]>>1)|(1&j?i>>1:0);
        for(R j=0;j+i<n;j+=i<<1){
            R*p=a+j,*q=p+i;
            NTT(p,i,1);NTT(q,i,1);
            for(R k=0;k<i;++k)p[k]=(LL)p[k]*q[k]%YL;
            NTT(p,i,-1);
            memset(q,0,4*i);
        }
    }
}

分治FFT

某些卷积式有一些特点，直接算计算量庞大，却可以通过CDQ分治思想，考虑左边对右边的贡献。
洛谷P4721 【模板】分治 FFT

void Solve(R l,R r){
	if(l+1==r)return;
	R m=(l+r)>>1;
	Solve(l,m);R n=Init(m-l+r-l);
	memcpy(a,f+l,4*(m-l));memset(a+m-l,0,4*(n-m+l));
	memcpy(b,g  ,4*(r-l));memset(b+r-l,0,4*(n-r+l));
	NTT(a,n,1);NTT(b,n,1);
	for(R i=0;i<n;++i)a[i]=(LL)a[i]*b[i]%YL;
	NTT(a,n,-1);
	for(R i=m;i<r;++i)f[i]=Mod(f[i]+a[i-l]);
	Solve(m,r);
}

泰勒展开

将函数\(f(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^{(n)}(x_0)}{n!}(x-x_0)^n+\xi\)

牛顿迭代推导倍增求解式

已知多项式函数\(F(x)\)，用倍增法求解\(B\)使得\(F(B)\equiv0(\mod x^{2^k})\)。
假设已经求出\(B,F(B)\equiv0(\mod x^{2^k})\)，现在求\(B_1,F(B_1)\equiv0(\mod x^{2^{k+1}})\)。
\(F(B_1)\)在\(x=B\)处只展开一项得到
\(F(B_1)=F(B)+F'(B)(B_1-B)\equiv0(\mod x^{2^{k+1}})\)
\(B_1=B-\frac{F(B)}{F'(B)}\)，带入即可。
接下来的\(A\)为已知多项式，B为待求多项式。

求逆

\(F(B)=AB-1=0\)
\(B_1=B-\frac{AB-1}{A}=B-B(AB-1)=2B-AB^2\)

开方

\(F(B)=B^2-A=0\)
\(B_1=B-\frac{B^2-A}{2B}=\frac{B^2+A}{2B}\)

对数

没必要迭代。
\(B=\ln A\)
\(B=\int\frac{A'}{A}\)

指数

\(B=e^A\)
\(\ln B-A=0\)
\(B_1=B-\frac{\ln B-A}{\frac1B}=B(1+A-\ln B)\)

快速幂

\(B=A^k=e^{k\ln A}\)
这常数看着就吓人。。。

除法&取余

https://www.luogu.org/problemnew/solution/P4512

代码模板

写的时候保证了一定的稳定性（数组清空问题），因此牺牲了一丁点效率。
参数中，a为已知，b为待求，a1、b1为辅助数组。
传入前需自己确保b、a1、b1为空。
返回后保证a维持原状，b存好答案，a1、b1为空。
还有，~~线性递推，多点求值，插值~~待填坑

const int YL=998244353;
int Inv[N],r[N];//在外面初始化逆元
inline int Mod(const int x){
	return x>=YL?x-YL:x;
}
LL Pow(LL b,R k=YL-2,LL a=1){
	for(;k;k>>=1,b=b*b%YL)
		if(k&1)a=a*b%YL;
	return a;
}
int Init(R m){
    R n=1;while(n<m<<1)n<<=1;
    for(R i=0;i<n;++i)r[i]=(r[i>>1]|(1&i)*n)>>1;
    return n;
}
void NTT(R*a,R n,R op){
	for(R i=0;i<n;++i)
		if(i<r[i])swap(a[i],a[r[i]]);
	for(R i=1;i<n;i<<=1){
		LL wn=Pow(3,(YL-1)/(i<<1)),w=1,x;
		if(op==-1)wn=Pow(wn);
		for(R j=0;j<n;j+=i<<1,w=1)
			for(R k=j;k<j+i;++k,w=w*wn%YL)
				x=w*a[k+i]%YL,a[k+i]=Mod(a[k]-x+YL),a[k]=Mod(a[k]+x);
	}
	if(op==-1){
		LL x=Pow(n);
		for(R i=0;i<n;++i)a[i]=a[i]*x%YL;
	}
}
void PolyRev(R*a,R*b,R n){
	for(R i=0;i<n;++i)b[i]=a[n-i-1];
}
void PolyDer(R*a,R*b,R n){
	for(R i=1;i<n;++i)b[i-1]=(LL)a[i]*i%YL;
	if(a==b)a[n-1]=0;
}
void PolyInt(R*a,R*b,R n){
	for(R i=n;i;--i)b[i]=(LL)a[i-1]*Inv[i]%YL;
	if(a==b)a[0]=0;
}
void PolyInv(R*a,R*b,R*a1,R m){
	if(m==1){b[0]=Pow(a[0]);return;}
	PolyInv(a,b,a1,(m+1)>>1);memcpy(a1,a,4*m);
	R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
	for(R i=0;i<n;++i)b[i]=(YL+2-(LL)b[i]*a1[i]%YL)*b[i]%YL;
	NTT(b,n,-1);memset(a1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolySqrt(R*a,R*b,R*a1,R*b1,R m){
	if(m==1){b[0]=sqrt(a[0]);return;}
	PolySqrt(a,b,a1,b1,(m+1)>>1);PolyInv(b,b1,a1,m);memcpy(a1,a,4*m);
	R n=Init(m);NTT(a1,n,1);NTT(b1,n,1);
	for(R i=0;i<n;++i)a1[i]=(LL)a1[i]*b1[i]%YL;
	NTT(a1,n,-1);
	for(R i=0;i<m;++i)b[i]=(LL)(b[i]+a1[i])*((YL+1)>>1)%YL;
	memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyLn(R*a,R*b,R*a1,R m){
	PolyInv(a,b,a1,m);PolyDer(a,a1,m);
	R n=Init(m);NTT(b,n,1);NTT(a1,n,1);
	for(R i=0;i<n;++i)b[i]=(LL)b[i]*a1[i]%YL;
	NTT(b,n,-1);memset(a1,0,4*n);PolyInt(b,b,m);
}
void PolyExp(R*a,R*b,R*a1,R*b1,R m){
	if(m==1){b[0]=1;return;}
	PolyExp(a,b,a1,b1,(m+1)>>1);PolyLn(b,b1,a1,m);memcpy(a1,a,4*m);
	R n=Init(m);NTT(b,n,1);NTT(a1,n,1);NTT(b1,n,1);
	for(R i=0;i<n;++i)b[i]=(LL)b[i]*(YL+1+a1[i]-b1[i])%YL;
	NTT(b,n,-1);memset(a1,0,4*n);memset(b1,0,4*n);memset(b+m,0,4*(n-m));
}
void PolyDiv(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
	PolyRev(a,a1,m);PolyInv(a1,b,A1,M-m+1);PolyRev(A,A1,M);
	R n=Init(M);NTT(b,n,1);NTT(A1,n,1);
	for(R i=0;i<n;++i)b[i]=(LL)b[i]*A1[i]%YL;
	NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
	reverse(b,b+M-m+1);memset(b+M-m+1,0,4*(n-M+m-1));
}
void PolyMod(R*A,R*a,R*b,R*A1,R*a1,R M,R m){
	PolyDiv(A,a,b,A1,a1,M,m);memcpy(A1,A,4*M);memcpy(a1,a,4*m);
	R n=Init(M);NTT(b,n,1);NTT(A1,n,1);NTT(a1,n,1);
	for(R i=0;i<n;++i)b[i]=Mod(A1[i]-(LL)b[i]*a1[i]%YL+YL);
	NTT(b,n,-1);memset(A1,0,4*n);memset(a1,0,4*n);
}

posted @ 2018-12-19 21:03 Flash_Hu 阅读(784) 评论(1) 编辑收藏举报

刷新页面返回顶部

FlashHu

——蝌