多项式模板

最新推荐文章于 2022-02-09 08:25:04 发布
最新推荐文章于 2022-02-09 08:25:04 发布
阅读量128 收藏
点赞数
CC 4.0 BY-SA版权
原文链接：http://www.cnblogs.com/f321dd/p/9363493.html

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

  • der(求导)/pri(积分)：$O(n)$
  • mul/inv/sqrt/log/exp/sin/cos/tan：$O(n\log n)$
  • div/mod：$O(n\log n)$
  • 在$m$个点处求值：$O(n\log n+m\log^2\min(m,n))$
  • 以根集构造：$O(n\log^2n)$
  • 以点集构造（插值）：$O(n\log^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(n\log^2n)$求多项式gcd的算法。

2018-08-18

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

2018-08-19

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

2019-03-17

  • 增加了sin/cos/tan。

2019-06-28

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

转载于:https://www.cnblogs.com/f321dd/p/9363493.html

【多项式n连】各种多项式模板（从入门到入土）
shiveringkonnyaku的博客
03-04 1817
qwq。。。好像很久没有写博客了的样子。。。因为最近一段时间一直处在觉得自己太蒻了的焦虑当中。。。导致啥也没学进去。。。这样的。。。本来定下的肝高数和概率论的计划也因为各种各样的原因没能完成。。。而降低到了只是零散地又复习了一些数学。。。然后发现luogu上有很多多项式的板子。。。就去做了一点。。。这样子。。。然后发现好久没更博客了。。。就准备放上来一些。。。也算是水一水博客了。。。反正老年选手已...
【学习总结】生成函数题目，多项式模板
执子之手,与子偕老
04-22 822
学习资料： 王乐平、策爷冬令营讲义。 多项式技巧 牛顿迭代（泰勒展开） 这是一篇非常好的博客。特别是把所有多项式操作都用泰勒展开推导，以后就不用再死记，也不用担心推错了！from yyb 复合逆 拉格朗日反演 注意当f(x)本身不存在逆元的时候，可以求f(x) / x的逆元 否则（x / f（x））n 的第n - 1 项为0 例题 直接推生成函数 【BZOJ3625】小朋友和二叉树 from ...
多项式模板&题目整理
banchangyue3928的博客
11-27 305
注：多项式的题目，数组应开：N的最近2的整数次幂的4倍。 多项式乘法 FFT模板 时间复杂度\(O(n\log n)\)。 模板： void FFT(Z *a,int x,int K){ static int rev[N],lst; int n=(1<<x); if(n!=lst){ for(int i=0;i<n;i++)rev...
完整版kuangbin的ACM模板（新）.pdf
01-21
完整版kuangbin的ACM模板（全新）.pdf 从最简单的搜索到几何问题等模板 涵盖了所有ACM专题的内容
任意模数多项式乘法（模板）
qq_35577488的博客
03-25 633
题目地址 不想理解原理了，直接记板子。大概有三种做法。我选择记下NTT这种常数最大的。希望不会遇到卡常的吧. #include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define LL long long #define maxn 400005 using namespace std; const LL mod1=9982
大景深系统多项式相位模板的优化设计
01-11
本篇论文的作者王晋疆、韩瑞雨以及莫绪涛三人，提出了一种基于波前编码的大景深成像系统优化设计，该设计基于立方相位模板，并进一步提出了一种新的奇对称多项式相位模板，以改善成像系统在不同离焦条件下的成像性能...
多项式系列模板
Never give in.
04-23 760
FFT struct Z { double x,y; Z(double _x=0,double _y=0):x(_x),y(_y){}; }; Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);} Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);} Z operator *(Z ...
多项式模板（inv, sqrt, ln, exp）
lifehappy
02-04 1415
#include <bits/stdc++.h> using namespace std; const int mod = 998244353, inv2 = mod + 1 >> 1; namespace Quadratic_residue { struct Complex { int r, i; Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {} }; int I2; Complex op
NTT和多项式各种基础运算模板
wzf_2000的博客
06-03 725
#include #define mod 998244353 #define N 400009 #define ll long long using namespace std; ll w[2][N],inv[N]; ll A[N],B[N],C[N],D[N],E[N],F[N]; ll W,G; ll pow(ll a,ll b,ll c) { ll ans=1; for (;b;a=a
多项式总结及模板
caoyang1123的博客
03-30 613
1.求逆 设要对多项式A求逆，逆为B： 求在模ceil（n/2）的逆，将原式与之相减后平方，发现可得在模n意义下也是0了，然后两边同乘A，移项即可倍增了（注意此时A也是在模当前长度意义下的）。 最后倍增形式为B'=2*B-B*B*A （b’等于2b减b方a），边界为常数取逆元，也可发现多项式有无逆取决于常数有无逆元 2.开根 与求逆类似，先求ceil（n/2）,然后将A移到左边使等式右...
多项式细节梳理&模板（多项式）
a83229442的博客
12-19 265
基础 很久以前的多项式总结 现在的码风又变了。。。 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=...
多项式模板QAQ
cgh_Andy's Blog
03-06 493
失踪人口回归？ 似乎抄了挺多天的多项式题目代码。。留几个模板吧。。 （都不知道是不是好模板。。 反正也是自己看 最长不超过80行我去 UOJ #34. 多项式乘法 #include &lt;bits/stdc++.h&gt; #define me(a,x) memset(a,x,sizeof a) using namespace std; const int N=3e5+2,inf=1...
进军多项式（一）：多项式全家桶模板
jwg2732的专栏
02-09 4267
多项式乘法(NTT) 默认大家会了FFT…… 我们知道FFT中因为要用三角函数和复数，所以对精度要求高 并且常数比较大 那么有没有可以代替单位根的呢？ 我们思考一下在FFT中应用的单位根的性质 1:wn0=1w_n^0=1wn0​=1 保证了有初始值 2:wni+k=wniwnjw_n^{i+k}=w_n^iw_n^jwni+k​=wni​wnj​,也就是wnk=(wn1)kw_n^k=(w_n^1)^kwnk​=(wn1​)k 因此只要求出一个就可以递推 3:wnk+n/2=−wnkw_n^{k+n/2}
【模板】多项式全家桶
weixin_33912638的博客
03-21 177
FFT struct Z { double x,y; Z(double _x=0,double _y=0):x(_x),y(_y){}; }; Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);} Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);} Z operator *(Z x,Z...
[多项式ln][多项式exp][多项式求幂][生成函数][DP][FNT] BZOJ 3684: 大朋友和多叉树
Vectorxj的博客
03-06 861
SolutionSolutionSolution 把DP写成生成函数的形式。f(x)=x+∑d∈Dfd(x)f(x)=x+∑d∈Dfd(x)f(x)=x+\sum_{d\in D}f^d(x)设g(f(x))=xg(f(x))=xg\big(f(x)\big)=x，有g(f(x))g(x)==f(x)−∑d∈Dfd(x)x−∑d∈Dxdg(f(x))=f(x)−∑d∈Dfd(x)g(x)=x−∑...
C++模板单向链表实现多项式加法
在C++中，通过模板和链表结构，我们可以实现一种通用的多项式相加方法。该方法的实现不仅可以帮助理解数据结构的灵活运用，还可以深入理解模板编程的强大功能。下面，将详细探讨C++模板实现多项式相加的知识点。 ##...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值