项式

项式

关注
发信
关注(28)粉丝(399)

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

基础

很久以前的多项式总结

现在的码风又变了。。。

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

}
05-08 15:23
Powered by LMLPHP ©2024 项式 0.063801