我的多项式全家桶

第一个自己实现的多项式全家桶。打比赛时终于有板子了!

代码是从之前学转置原理的那篇博客里升级来的。

但是功能很残?效率很逊?接口很怪?评价是能用就行。

目前封装了:乘、逆、除(取模)、开方(无二次剩余,不会 qwq)、对数、指数、多点求值、快速插值、常系数齐次线性递推、Berlekamp-Massey 算法。

为什么突然开始多项式?原因是在尝试改鸡的时候发现了 sjy 的究极神秘题解,感觉 BM 算法竟然有意想不到的用处,于是学了 BM 算法,然后递归到了常系数其次线性递推,然后递归到了除法,然后升级了多项式板子。

不过话说 BM 算法其实与多项式没有很大关系不是吗?

#include <cstdio>
#include <vector>
#include <cassert>
#include <algorithm>
#define lc (p<<1)
#define rc (p<<1|1)
#define ALL(p) p.begin(),p.end()
#define IL inline
// #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<22,stdin)),p1==p2?EOF:*p1++)
using namespace std;
char buf[1<<22],*p1=buf,*p2=buf;
int read(){
	char c=getchar();int x=0;bool f=0;
	while(c<48||c>57) f|=(c=='-'),c=getchar();
	do x=(x<<1)+(x<<3)+(c^48),c=getchar();
	while(c>=48&&c<=57);
	if(f) return -x;
	return x;
}
typedef vector<int> vi;
const int P=998244353;
typedef long long ll;
IL int qp(int a,int b=P-2){
	int res=1;
	while(b){
		if(b&1) res=(ll)res*a%P;
		a=(ll)a*a%P;b>>=1;
	}
	return res;
}
int len,ilen,bt;
int rev[1<<20],cw[1<<20|1];
IL void init(int _len){ // mod x^len
	len=1;bt=-1;
	while(len<_len) len<<=1,++bt;
	int w=qp(3,(P-1)>>(bt+1));
	cw[0]=cw[len]=1;
	for(int i=1;i<len;++i){
		cw[i]=(ll)cw[i-1]*w%P;
		rev[i]=(rev[i>>1]>>1)|((i&1)<<bt);
	}
	ilen=qp(len);
}
struct poly{
	vi f;
	poly():f(){}
	poly(int Len):f(Len){}
	poly(initializer_list<int> Init):f(Init){}
	IL int& operator[](const int &x){return f[x];}
	IL const int& operator[](const int &x)const{return f[x];}
	IL void NTT(){
		f.resize(len,0);
		for(int i=1;i<len;++i) if(rev[i]<i) swap(f[rev[i]],f[i]);
		for(int i=1,tt=len>>1;i<len;i<<=1,tt>>=1)
			for(int j=0;j<len;j+=(i<<1))
				for(int k=j,t=0;k<(j|i);++k,t+=tt){
					int x=f[k],y=(ll)f[k|i]*cw[t]%P;
					if((f[k]=x+y)>=P) f[k]-=P;
					if((f[k|i]=x-y)<0) f[k|i]+=P;
				}
	}
	IL void INTT(){
		for(int i=1;i<len;++i) if(rev[i]<i) swap(f[rev[i]],f[i]);
		for(int i=1,tt=len>>1;i<len;i<<=1,tt>>=1)
			for(int j=0;j<len;j+=(i<<1))
				for(int k=j,t=len;k<(j|i);++k,t-=tt){
					int x=f[k],y=(ll)f[k|i]*cw[t]%P;
					if((f[k]=x+y)>=P) f[k]-=P;
					if((f[k|i]=x-y)<0) f[k|i]+=P;
				}
		for(int i=0;i<len;++i) f[i]=(ll)f[i]*ilen%P;
	}
	IL int size(){return f.size();}
	IL void reduce(){while(!f.empty()&&!f.back()) f.pop_back();}
	IL void rever(){reverse(ALL(f));}
	IL void show(){
		for(int x:f) printf("%d ",x);
		putchar('\n');
	}
	IL void trunc(int lim){
		if(lim<int(f.size())) f.erase(f.begin()+lim,f.end());
	}
	IL poly inv(int lim){ // mod x^lim
		assert(f[0]);
		poly cur({qp(f[0])});
		for(int t=1;t<lim;t<<=1){
			poly ff(t<<2);
			copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
			init(t<<2);ff.NTT();cur.NTT();
			poly tmp(len);
			for(int i=0;i<len;++i){
				tmp[i]=(2ll-(ll)cur[i]*ff[i])%P*cur[i]%P;
				if(tmp[i]<0) tmp[i]+=P;
			}
			tmp.INTT();
			cur.f.swap(tmp.f);
			cur.trunc(t<<1);
		}
		cur.trunc(lim);
		return cur;
	}
	friend poly operator+(poly A,poly B){
		int mx=max(A.size(),B.size());
		A.f.resize(mx,0);B.f.resize(mx,0);
		poly C(mx);
		for(int i=0;i<mx;++i){
			C[i]=A[i]+B[i];
			if(C[i]>=P) C[i]-=P;
		}
		return C;
	}
	friend poly operator-(poly A,poly B){
		int mx=max(A.size(),B.size());
		A.f.resize(mx,0);B.f.resize(mx,0);
		poly C(mx);
		for(int i=0;i<mx;++i){
			C[i]=A[i]-B[i];
			if(C[i]<0) C[i]+=P;
		}
		C.reduce();
		return C;
	}
	friend poly operator*(poly A,poly B){
		init(A.size()+B.size()-1);A.NTT();B.NTT();
		poly C(len);
		for(int i=0;i<len;++i) C[i]=(ll)A[i]*B[i]%P;
		C.INTT();C.reduce();
		return C;
	}
	poly operator*(int X){
		int n=size();
		poly F(n);
		for(int i=0;i<n;++i) F[i]=(ll)f[i]*X%P;
		return F;
	}
	friend poly operator%(poly F,poly G){
		int n=F.size()-1,m=G.size()-1;
		if(n<m) return F;
		F.rever();G.rever();
		poly Q=G.inv(n-m+1);
		Q.prod(Q,F);Q.f.resize(n-m+1,0);
		F.rever();G.rever();Q.rever();
		poly R=F-Q*G;
		R.reduce();
		return R;
	}
	IL void prod(poly A,poly B){
		init(A.size()+B.size()-1);A.NTT();B.NTT();
		f.resize(len);
		for(int i=0;i<len;++i) f[i]=(ll)A[i]*B[i]%P;
		INTT();reduce();
	}
	IL void prodT(poly A,poly B){
		int an=A.size()-1,bn=B.size()-1;
		reverse(ALL(B.f));prod(A,B);
		for(int i=0;i<=an;++i) f[i]=f[i+bn];
		trunc(an+1);
	}
	IL int calc(int t){
		int pw=1,res=0;
		for(int x:f){
			res=(res+(ll)pw*x)%P;
			pw=(ll)pw*t%P;
		}
		return res;
	}
	IL poly deriv(){
		int n=f.size();
		poly D(n-1);
		for(int i=1;i<n;++i) D[i-1]=(ll)f[i]*i%P;
		return D;
	}
	IL poly integ(){
		int n=f.size();
		poly D(n+1),inv(n+1);
		for(int i=1;i<=n;++i){
			if(i>1) inv[i]=(ll)inv[P%i]*(P-P/i)%P;
			else inv[1]=1;
			D[i]=(ll)f[i-1]*inv[i]%P;
		}
		return D;
	}
	IL poly getln(int lim){ // mod x^lim
		assert(f[0]==1);
		poly F=deriv()*inv(lim);
		F.f.resize(lim-1,0);
		return F.integ();
	}
	IL poly getexp(int lim){ // mod x^lim
		assert(f[0]==0);
		poly F({1});
		for(int t=1;t<lim;t<<=1){
			poly ff(t<<1);
			copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
			F=F*(poly({1})-F.getln(t<<1)+ff);
			F.f.resize(t<<1,0);
		}
		F.trunc(lim);
		return F;
	}
	IL poly sqrt(int lim){ // without quadratic residue , mod x^lim
		assert(f[0]==1);
		poly F({1});
		for(int t=1;t<lim;t<<=1){
			poly ff(t<<1);
			copy(f.begin(),f.begin()+min(t<<1,size()),ff.f.begin());
			F=(F+F.inv(t<<1)*ff)*((P+1)>>1);
			F.f.resize(t<<1,0);
		}
		F.trunc(lim);
		return F;
	}
}F;
namespace extend{
	const int N=200003;
	int n;
	int px[N],py[N],res[N];
	namespace multieval{
		poly G[N<<2],H[N<<2];
		void calc(int p,int l,int r){
			if(l==r){G[p]=poly({1,(P-px[r])%P});return;}
			int mid=(l+r)>>1;
			calc(lc,l,mid);
			calc(rc,mid+1,r);
			G[p].prod(G[lc],G[rc]);
		}
		void solve(int p,int l,int r){
			H[p].trunc(r-l+1);
			if(l==r){res[r]=H[p][0];return;}
			int mid=(l+r)>>1;
			G[p].f.clear();G[p].f.shrink_to_fit();
			H[lc].prodT(H[p],G[rc]);
			H[rc].prodT(H[p],G[lc]);
			H[p].f.clear();H[p].f.shrink_to_fit();
			solve(lc,l,mid);
			solve(rc,mid+1,r);
		}
		void sol(){
			// init n, px, py first
			calc(1,1,n);
			H[1].prodT(F,G[1].inv(n+1));
			solve(1,1,n);
		}
	}
	namespace interplot{
		poly T[N<<2];
		void eval(int p,int l,int r){
			if(l==r){T[p]=poly({(P-px[r])%P,1});return;}
			int mid=(l+r)>>1;
			eval(lc,l,mid);
			eval(rc,mid+1,r);
			T[p].prod(T[lc],T[rc]);
		}
		poly proc(int p,int l,int r){
			if(l==r) return poly({py[r]});
			int mid=(l+r)>>1;
			poly L=proc(lc,l,mid);
			poly R=proc(rc,mid+1,r);
			L.prod(L,T[rc]);
			R.prod(R,T[lc]);
			return L+R;
		}
		poly solve(){
			// init n, px, py first
			eval(1,1,n);
			F=T[1].deriv();
			multieval::sol();
			for(int i=1;i<=n;++i) py[i]=(ll)py[i]*qp(res[i])%P;
			return proc(1,1,n);
		}
	}
}
namespace berlekamp_massey{
	const int N=200003;
	int n,k,kk,las,pos,ex;
	int s[N],ss[N],a[N],tmp[N];
	int solve(){
		n=read();ex=read();
		for(int i=0;i<n;++i) a[i]=read();
		bool fl=1;
		for(int i=0;i<n;++i){
			int cur=a[i];
			for(int j=1;j<=k;++j) cur=(cur-(ll)a[i-j]*s[j])%P;
			if(cur<0) cur+=P;
			if(!cur) continue;
			if(fl){
				k=i+1;kk=0;pos=i;las=cur;fl=0;
				continue;
			}
			int nk=max(k,i-pos+kk);
			int val=(ll)cur*qp(las)%P;
			for(int j=1;j<i-pos;++j) tmp[j]=0;
			tmp[i-pos]=val;
			for(int j=1;j<=kk;++j) tmp[i-pos+j]=(ll)(P-val)*ss[j]%P;
			if(k-i<kk-pos){
				kk=k;pos=i;las=cur;
				for(int j=1;j<=k;++j) ss[j]=s[j];
			}
			k=nk;
			for(int j=1;j<=k;++j){
				s[j]+=tmp[j];
				if(s[j]>=P) s[j]-=P;
			}
		}
		poly M(k+1),F(2),G({1});
		M[k]=1;
		for(int i=1;i<=k;++i)
			if(s[i]) M[k-i]=P-s[i];
		F[1]=1;
		while(ex){
			if(ex&1) G=G*F%M;
			F=F*F%M;ex>>=1;
		}
		G.f.resize(k,-1);
		int ans=0;
		for(int i=0;i<k;++i) ans=(ans+(ll)G[i]*a[i])%P;
		if(ans<0) ans+=P;
		return ans;
	}
}
int main(){
	int n=read();
	poly F(n);
	for(int i=0;i<n;++i) F[i]=read();
	F=F.sqrt(n);
	F.show();
	return 0;
}
posted @ 2024-01-22 08:42  yyyyxh  阅读(65)  评论(0编辑  收藏  举报