多项式模板

这里给出了一个多项式/形式幂级数的类的实现。由于仅考虑形式幂级数的前n项，即$!mod{x^n}$下的等价类，故其形式与多项式相同，因而在类的实现上没有对两者作区分。没有返回值的方法会将结果用于更新自身。对于所有用于形式幂级数的方法，参数n代表$!mod{x^n}$，且保证调用后前n项均可访问。三角函数使用了要求$sqrt{-1}$存在的实现。下表说明了各方法的复杂度。

• der(求导)/pri(积分)：$O(n)$
• mul/inv/sqrt/log/exp/sin/cos/tan：$O(nlog n)$
• div/mod：$O(nlog n)$
• 在$m$个点处求值：$O(nlog n+mlog^2min(m,n))$
• 以根集构造：$O(nlog^2n)$
• 以点集构造（插值）：$O(nlog^2n)$

注意，对于某些算式，相比直接调用以上方法进行计算，使用针对性的实现能使效率大幅提高。

---

#include<algorithm>
#include<vector>
#define RAN(a)a.begin(),a.end()
using namespace std;
typedef unsigned long long u64;
typedef unsigned u32;
namespace num{
	const u32 p=998244353;
	const u32 g=3;
	inline u32 imod(u32 a){
		return a<p?a:a-p;
	}
	inline u32 ipow(u32 a,u32 n){
		u32 s=1;
		for(;n;n>>=1){
			if(n&1)
				s=(u64)s*a%p;
			a=(u64)a*a%p;
		}
		return s;
	}
	class inv_t{
	public:
		inv_t():f(1,1){}
		u32 operator[](int n){
			int m=f.size();
			if(m<n){
				f.resize(n);
				for(int i=m+1;i<=n;++i)
					f[i-1]=(u64)(p-p/i)*f[p%i-1]%p;
			}
			return f[n-1];
		}
		u32 operator()(u32 a)const{
			return ipow(a,p-2);
		}
	private:
		vector<u32>f;
	}inv;
}
using namespace num;
class poly{
public:
	vector<u32>a;
	poly(){}
	explicit poly(int n):a(n){}
	u32&operator[](int i){
		return a[i];
	}
	const u32&operator[](int i)const{
		return a[i];
	}
	int size()const{
		return a.size();
	}
	void swap(poly&b){
		a.swap(b.a);
	}
	void der(){
		fix();
		if(size()){
			for(int i=1;i<size();++i)
				a[i-1]=(u64)i*a[i]%p;
			a.pop_back();
		}
	}
	void pri(){
		fix();
		shl();
		for(int i=size()-1;i>0;--i)
			a[i]=(u64)a[i]*num::inv[i]%p;
	}
	static int len(int n){
		while(n^n&-n)
			n+=n&-n;
		return n;
	}
	void fft(int n,bool f){
		a.resize(n);
		if(n<=1)
			return;
		for(int i=0,j=0;i<n;++i){
			if(i<j)
				std::swap(a[i],a[j]);
			int k=n>>1;
			while((j^=k)<k)
				k>>=1;
		}
		vector<u32>w(n/2);
		w[0]=1;
		for(int i=1;i<n;i<<=1){
			for(int j=i/2-1;~j;--j)
				w[j<<1]=w[j];
			int m=(p-1)/i/2;
			u64 s=ipow(g,f?p-1-m:m);
			for(int j=1;j<i;j+=2)
				w[j]=s*w[j-1]%p;
			for(int j=0;j<n;j+=i<<1){
				u32*b=&a[0]+j,*c=b+i;
				for(int k=0;k<i;++k){
					u32 v=(u64)w[k]*c[k]%p;
					c[k]=imod(b[k]+p-v);
					b[k]=imod(b[k]+v);
				}
			}
		}
	}
	void dft(int n){
		fft(n,0);
	}
	void idft(){
		int n=size();
		fft(n,1);
		u64 f=num::inv(n);
		for(int i=0;i<n;++i)
			a[i]=f*a[i]%p;
	}
	void fix(){
		while(size()&&!a.back())
			a.pop_back();
	}
	void mul(poly b){
		fix();
		b.fix();
		int n=len(size()+b.size()-1);
		dft(n);
		b.dft(n);
		for(int i=0;i<n;++i)
			a[i]=(u64)a[i]*b[i]%p;
		idft();
		fix();
	}
	void mod(int n){
		a.resize(n);
	}
	void inv(int n){
		int m=len(n);
		mod(m);
		vector<u32>b(1,num::inv(a[0]));
		a.swap(b);
		for(int i=2;i<=m;i<<=1){
			int l=i<<1;
			poly c(l);
			for(int j=0;j<i;++j)
				c[j]=b[j];
			c.dft(l);
			dft(l);
			for(int j=0;j<l;++j)
				a[j]=a[j]*(2+p-(u64)a[j]*c[j]%p)%p;
			idft();
			mod(i);
		}
		mod(n);
	}
	void sqrt(int n){
		int m=len(n);
		mod(m);
		vector<u32>b(1,1);
		a.swap(b);
		u64 w=(p+1)/2;
		for(int i=2;i<=m;i<<=1){
			poly c(i);
			for(int j=0;j<i;++j)
				c[j]=b[j];
			vector<u32>t=a;
			inv(i);
			mul(c);
			mod(i);
			for(int j=0;j<i>>1;++j)
				a[j]=w*(a[j]+t[j])%p;
			for(int j=i>>1;j<i;++j)
				a[j]=w*a[j]%p;
		}
		mod(n);
	}
	void log(int n){
		mod(n);
		poly b=*this;
		der();
		b.inv(n-1);
		mul(b);
		mod(n-1);
		pri();
		mod(n);
	}
	void exp(int n){
		int m=len(n);
		mod(m);
		vector<u32>b(1,1);
		a.swap(b);
		for(int i=2;i<=m;i<<=1){
			poly c=*this;
			log(i);
			for(int j=0;j<i;++j)
				a[j]=imod(b[j]+p-a[j]);
			++a[0]%=p;
			mul(c);
			mod(i);
		}
		mod(n);
	}
	void sin(int n){
		mod(n);
		u64 w=ipow(g,(p-1)/4);
		for(int i=0;i<n;++i)
			a[i]=a[i]*w%p;
		exp(n);
		poly b=*this;
		b.inv(n);
		w=(p+1)/2*(p-w)%p;
		for(int i=0;i<n;++i)
			a[i]=(a[i]+p-b[i])*w%p;
	}
	void cos(int n){
		mod(n);
		u64 w=ipow(g,(p-1)/4);
		for(int i=0;i<n;++i)
			a[i]=a[i]*w%p;
		exp(n);
		poly b=*this;
		b.inv(n);
		w=(p+1)/2;
		for(int i=0;i<n;++i)
			a[i]=(a[i]+b[i])*w%p;
	}
	void tan(int n){
		mod(n);
		u64 w=ipow(g,(p-1)/4)*2;
		for(int i=0;i<n;++i)
			a[i]=a[i]*w%p;
		exp(n);
		++a[0]%=p;
		inv(n);
		for(int i=0;i<n;++i)
			a[i]=a[i]*w%p;
		(a[0]+=p-w/2)%=p;
	}
	poly mod_bf(poly b){
		fix();
		b.fix();
		int n=size();
		int m=b.size();
		if(n<m)
			return poly();
		poly q(n-m+1);
		u64 w=num::inv(b[m-1]);
		for(;n>=m;--n)
			if(u64 s=a[n-1]*w%p){
				q[n-m]=s;
				for(int i=n-1;i>=n-m;--i)
					a[i]=(a[i]+(p-s)*b[i-n+m])%p;
			}
		fix();
		return q;
	}
	void gcd(poly b){
		fix();
		b.fix();
		while(b.size()){
			mod_bf(b);
			swap(b);
		}
	}
	void div(poly b){
		fix();
		b.fix();
		int n=size()-b.size()+1;
		if(n<=0){
			a.clear();
			return;
		}
		reverse(RAN(a));
		mod(n);
		reverse(RAN(b.a));
		b.inv(n);
		mul(b);
		mod(n);
		reverse(RAN(a));
		fix();
	}
	void mod(poly b){
		fix();
		b.fix();
		int m=b.size();
		if(size()>=m){
			poly c=*this;
			div(b);
			mul(b);
			mod(m-1);
			for(int i=0;i<m-1;++i)
				a[i]=imod(c[i]+p-a[i]);
			fix();
		}
	}
	u32 val(u32 x)const{
		u32 s=0;
		for(int i=size()-1;~i;--i)
			s=((u64)s*x+a[i])%p;
		return s;
	}
	u32 operator()(u32 x)const{
		return val(x);
	}
	template<class ite1,class ite2>
	void val(int n,ite1 x,ite2 y)const;
	template<class ite>
	poly(int n,ite x);
	template<class ite1,class ite2>
	poly(int n,ite1 x,ite2 y);
	class seg;
private:
	void shl(){
		a.insert(a.begin(),0);
	}
	void mod(poly b,const poly&f);
	template<class ite>
	static poly gen(int n,ite a);
	template<class ite>
	static poly gen(int x,int y,ite&a);
	template<class ite>
	static poly gen_bf(int n,ite&a);
};
void poly::mod(poly b,const poly&f){
	fix();
	b.fix();
	int m=b.size();
	if(size()>=m){
		poly c=*this;
		div(b);
		dft(f.size());
		for(int i=0;i<f.size();++i)
			a[i]=(u64)a[i]*f[i]%p;
		idft();
		for(int i=0;i<m-1;++i)
			a[i]=imod(c[i]+p-a[i]);
		mod(m-1);
		fix();
	}
}
template<class ite>
poly::poly(int n,ite x){
	*this=gen(n,x);
}
template<class ite>
poly poly::gen(int n,ite a){
	return gen(0,n,a);
}
template<class ite>
poly poly::gen(int x,int y,ite&a){
	if(y-x<=200){
		return gen_bf(y-x,a);
	}else{
		int m=x+y>>1;
		poly b=gen(x,m,a);
		b.mul(gen(m,y,a));
		return b;
	}
}
template<class ite>
poly poly::gen_bf(int n,ite&a){
	poly f(1);
	f[0]=1;
	for(int i=0;i<n;++i){
		f.shl();
		u64 s=p-*a++;
		for(int j=0;j<=i;++j)
			f[j]=(f[j]+s*f[j+1])%p;
	}
	return f;
}
class poly::seg{
public:
	template<class ite>
	seg(int n,const ite&x):a(n){
		ite t=x;
		for(int i=0;i<n;++i)
			a[i]=*t++;
		dfs(0,n,1);
	}
	poly gen()const{
		return t[1].a;
	}
	template<class ite>
	poly gen(const ite&b)const{
		ite t=b;
		return gen(t);
	}
	template<class ite>
	void val(const ite&b,const poly&f)const{
		ite t=b;
		val(t,f);
	}
	template<class ite>
	seg(int n,ite&x):a(n){
		for(int i=0;i<n;++i)
			a[i]=*x++;
		dfs(0,n,1);
	}
	template<class ite>
	poly gen(ite&b)const{
		poly f=t[1].a;
		f.der();
		return gen(0,a.size(),1,b,f);
	}
	template<class ite>
	void val(ite&b,const poly&f)const{
		val(0,a.size(),1,b,f);
	}
private:
	struct pair{
		poly a,b;
	};
	vector<pair>t;
	vector<u32>a;
	void dfs(int x,int y,int k){
		if(t.size()<=k)
			t.resize(k+1);
		if(y-x<=200){
			t[k].a=poly::gen(y-x,a.begin()+x);
		}else{
			int m=x+y>>1,i=k<<1,j=i^1;
			dfs(x,m,i);
			dfs(m,y,j);
			int n=len(y-x+1);
			t[i].b=t[i].a;
			t[i].b.dft(n);
			t[j].b=t[j].a;
			t[j].b.dft(n);
			t[k].a=t[i].b;
			for(int l=0;l<n;++l)
				t[k].a[l]=(u64)t[k].a[l]*t[j].b[l]%p;
			t[k].a.idft();
		}
	}
	template<class ite>
	void val(int x,int y,int k,ite&b,poly f)const{
		if(k!=1)
			f.mod(t[k].a,t[k].b);
		else
			f.mod(t[k].a);
		if(y-x<=200){
			for(int i=0;i<y-x;++i)
				*b++=f(a[x+i]);
		}else{
			int m=x+y>>1,i=k<<1,j=i^1;
			val(x,m,i,b,f);
			val(m,y,j,b,f);
		}
	}
	template<class ite>
	poly gen(int x,int y,int k,ite&b,poly f)const{
		if(k!=1)
			f.mod(t[k].a,t[k].b);
		else
			f.mod(t[k].a);
		if(y-x<=200){
			vector<u64>c(y-x);
			for(int i=0;i<y-x;++i){
				u64 n=a[x+i];
				u64 m=*b++;
				m=m*num::inv(f(n))%p;
				u64 s=0;
				for(int j=y-x;j>0;--j){
					s=(t[k].a[j]+s*n)%p;
					c[j-1]+=s*m%p;
				}
			}
			poly d(y-x);
			for(int i=0;i<y-x;++i)
				d[i]=c[i]%p;
			return d;
		}else{
			int m=x+y>>1;
			int i=k<<1;
			int j=i^1;
			poly c=gen(x,m,i,b,f);
			poly d=gen(m,y,j,b,f);
			int n=len(y-x+1);
			c.dft(n);
			d.dft(n);
			for(int l=0;l<n;++l)
				c[l]=((u64)c[l]*t[j].b[l]+(u64)d[l]*t[i].b[l])%p;
			c.idft();
			return c;
		}
	}
};
template<class ite1,class ite2>
poly::poly(int n,ite1 x,ite2 y){
	*this=seg(n,x).gen(y);
}
template<class ite1,class ite2>
void poly::val(int n,ite1 x,ite2 y)const{
	int m=size();
	if(min(m,n)<=100)
		for(int i=0;i<n;++i)
			*y++=val(*x++);
	else
		for(int l=n;;l/=2)
			if(l*2<=m){
				for(int i=0;i<n;i+=l)
					seg(min(l,n-i),x).val(y,*this);
				break;
			}
}

---

2018-08-14

• 计划今后增加$O(nlog^2n)$求多项式gcd的算法。

2018-08-18

• 增加了$O(n^2)$的gcd，未保证结果为首一多项式。
• 分离了基础部分和点值有关部分。

2018-08-19

• 为避免影响编译器优化，删去了在运行时确定模数的功能。

2019-03-17

• 增加了sin/cos/tan。

2019-06-28

• 修复了除法中非预期地降低效率的问题。