首页 > 其他分享 >多项式简陋入门

多项式简陋入门

时间：2022-10-13 16:55:35浏览次数：47  
标签：lim return 入门 int 多项式 pmod 简陋 inline

多项式全家桶

然而并没有多点求值，快速插值，转下降/上升幂，复合，复合逆

疯狂多项式，v我50

namespace efX_poly
{
	const int maxlen=(1<<23)+1,maxSqrt=1e5+1;
	inline int add(int x,int y,int m) { x+=y; return x>=m?x-m:x; }
	inline int sub(int x,int y,int m) { x-=y; return x<0?x+m:x; }
	inline int mul(int x,int y,int m) { ll z=1ll*x*y; return z>=m?z%m:z; }
	inline void Add(int &x,int y,int m) { x=add(x,y,m); }
	inline void Sub(int &x,int y,int m) { x=sub(x,y,m); }
	inline void Mul(int &x,int y,int m) { x=mul(x,y,m); }
	inline int qpow(int a,int x,int m)
	{
		int r=1;
		while(x)
		{
			if(x&1) Mul(r,a,m);
			Mul(a,a,m),x>>=1;
		} return r;
	} inline int _inv(int a,int m) { return qpow(a,m-2,m); }
	namespace set_mod
	{
		static vector<int> d;
		inline void init(int p)
		{
			int tmp=p-1;
			if(!d.empty()) d.clear();
			for(int i=2;1ll*i*i<=tmp;++i)
				if(tmp%i==0)
			{
				d.emplace_back(i);
				while(tmp%i==0) tmp/=i;
			}
			if(tmp>1) d.emplace_back(tmp);
		}
		inline bool check(int x,int p)
		{
			for(int v:d)
				if(qpow(x,(p-1)/v,p)==1)
					return false;
			return true;
		}
		inline int getg(int p)
		{
			init(p);
			irep(i,2,p)
			if(check(i,p))
				return i;
			return -1;
		}
	}
	int r[maxlen];
	inline void make_rev(int lim,int l)
	{
		static int his=-1;
		if(his==lim) return;
		for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
		his=lim;
	}
	inline pii getLimit(int n)
	{
		int lim=1,l=0;
		while(lim<n) lim<<=1,++l;
		return make_pair(lim,l);
	}
	
	int inv[maxlen];
	int tmp[maxlen],tmpb[maxlen];
	inline void view_arr(int *a,int lim) { cerr<<"## "; for(int i=0;i<lim;++i) {cerr<<a[i]<<' ';} cerr_out(); }
	class poly
	{
	public:
		const int p,g,gi;
	private:
		int G[2][maxSqrt+1],GI[2][maxSqrt+1];
		inline void init()
		{
			memset(inv,0,sizeof inv);//can't use it any more
			G[0][0]=GI[0][0]=1;
			for(int i=1;i<=maxSqrt;++i)
			{
				G[0][i]=mul(G[0][i-1],g,p);
				GI[0][i]=mul(GI[0][i-1],gi,p);
			}
			G[1][0]=GI[1][0]=1;
			for(int i=1;i<=maxSqrt;++i)
			{
				G[1][i]=mul(G[1][i-1],G[0][maxSqrt],p);
				GI[1][i]=mul(GI[1][i-1],GI[0][maxSqrt],p);
			}
		}
		inline int very_quick_pow(bool o,int x)
		{
			if(o) return mul(G[1][x/maxSqrt],G[0][x%maxSqrt],p);
			return mul(GI[1][x/maxSqrt],GI[0][x%maxSqrt],p);
		}
	public:
		inline void NTT(int *a,int lim,int op)
		{
			for(int i=0;i<lim;++i) if(i<r[i]) swap(a[i],a[r[i]]);
			for(int mid=1,wn;mid<lim;mid<<=1)
			{
				wn=very_quick_pow(op==1,(p-1)/(mid<<1));
				for(int j=0,w;j<lim;j+=(mid<<1))
				{
					w=1;
					for(int k=0,x,y;k<mid;++k,w=mul(w,wn,p))
						x=a[j+k],
						y=mul(w,a[j+k+mid],p),
						a[j+k]=add(x,y,p),
						a[j+k+mid]=sub(x,y,p);
				}
			}
			if(op==-1)
			{
				int iv=_inv(lim,p);
				for(int i=0;i<lim;++i) Mul(a[i],iv,p);
			}
		}
		inline void clear(int *a,int st,int ed) { for(int i=st;i<ed;++i) a[i]=0; }
		inline void arrSub(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Sub(a[i],b[i],p); }
		inline void arrSub(int *a,int c,int lim) { for(int i=0;i<lim;++i) Sub(a[i],c,p); }
		inline void arrAdd(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Add(a[i],b[i],p); }
		inline void arrAdd(int *a,int c,int lim) { for(int i=0;i<lim;++i) Add(a[i],c,p); }
		inline void Reverse(int *a,int *b,int n) { for(int i=0;i<n;++i) b[i]=a[n-i-1]; }
		inline void arrMul(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Mul(a[i],b[i],p); }
		inline void arrMul(int *a,int c,int lim) { for(int i=0;i<lim;++i) Mul(a[i],c,p); }
		inline void Inv(int *a,int *b,int n)
		{
			static int tmp[maxlen];//开成全局就gg了
			if(n==1) { b[0]=_inv(a[0],p); return; }
			Inv(a,b,(n+1)>>1);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memcpy(tmp,a,sizeof(int)*n);
			for(int i=n;i<lim;++i) tmp[i]=0;
			NTT(tmp,lim,1),NTT(b,lim,1);
			for(int i=0;i<lim;++i) b[i]=mul(sub(2,mul(b[i],tmp[i],p),p),b[i],p);
			NTT(b,lim,-1);//tmp就不用管了，以后用的时候会刷新的
			for(int i=n;i<lim;++i) b[i]=0;
		}
		inline void Sqrt(int *a,int *b,int n)
		{
			if(n==1) { b[0]=1; return; }
			if(!inv[2]) inv[2]=_inv(2,p);
			Sqrt(a,b,(n+1)>>1);
			memset(tmp,0,sizeof(int)*(n<<1));
			Inv(b,tmp,n);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memcpy(tmpb,a,sizeof(int)*n);
			clear(tmpb,n,lim);
			NTT(tmpb,lim,1);
			NTT(tmp,lim,1);
			NTT(b,lim,1);
			for(int i=0;i<lim;++i)
				Add(b[i],mul(tmpb[i],tmp[i],p),p),
				Mul(b[i],inv[2],p);
			NTT(b,lim,-1);
			clear(b,n,lim);
		}//be very very careful!
		inline void Int(int *a,int *b,int n)
		{
			for(int i=n-1;i;--i)
			{
				if(!inv[i]) inv[i]=_inv(i,p);
				b[i]=mul(a[i-1],inv[i],p);
			} b[0]=0;
		}
		inline void Der(int *a,int *b,int n)
		{
			for(int i=1;i<n;++i) b[i-1]=mul(a[i],i,p); 
			b[n-1]=0;//attention to the bound
		}
		inline void Ln(int *a,int *b,int n)
		{
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memset(tmp,0,sizeof(int)*lim);
			memset(tmpb,0,sizeof(int)*lim);
			memcpy(tmp,a,sizeof(int)*n);
			Der(a,tmp,n);
			Inv(a,tmpb,n);
			NTT(tmp,lim,1),NTT(tmpb,lim,1);
			for(int i=0;i<lim;++i) tmp[i]=mul(tmp[i],tmpb[i],p);
			NTT(tmp,lim,-1);
			Int(tmp,tmp,n);
			memcpy(b,tmp,sizeof(int)*n);
		}
		inline void Exp(int *a,int *b,int n)
		{
			if(n==1) { b[0]=1; return; }
			Exp(a,b,(n+1)>>1);
			Ln(b,tmp,n);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			for(int i=0;i<n;++i) tmp[i]=sub(a[i],tmp[i],p);
			clear(tmp,n,lim);
			clear(b,n,lim);
			tmp[0]++;
			NTT(tmp,lim,1);
			NTT(b,lim,1);
			arrMul(b,tmp,lim);
			NTT(b,lim,-1);
			clear(b,n,lim);
		}
		inline void Pow(int *a,int k,int *b,int n)
		{
			static int tmp[maxlen];
			memset(tmp,0,sizeof(int)*(n<<1));
			Ln(a,tmp,n);
			arrMul(tmp,k,n);
			Exp(tmp,b,n);
		}//it's not totally safe
		inline void Mod(int *a,int *b,int *P,int *Q,int n,int m)
		{
			--n,--m;//强行钦定使用左闭右开区间
			static int ar[maxlen],br[maxlen];
			memset(ar,0,sizeof(int)*(n));
			memset(br,0,sizeof(int)*(n));
			Reverse(a,ar,n+1);
			Reverse(b,br,m+1);
			irep(i,n-m+1,m) br[i]=0;
			Inv(br,P,n-m+1);
			pii lm=getLimit(2*n-m+1);//ar*p=br^-1
			make_rev(lm.first,lm.second);
			NTT(ar,lm.first,1),
			NTT(P,lm.first,1);
			arrMul(ar,P,lm.first);
			NTT(ar,lm.first,-1);
			Reverse(ar,P,n-m+1);//所以p不需要ntt回来
			lm=getLimit(n+1);
			make_rev(lm.first,lm.second);
			clear(P,n-m+1,lm.first);
			memcpy(ar,b,sizeof(int)*lm.first);
			NTT(P,lm.first,1),
			NTT(ar,lm.first,1);
			arrMul(ar,P,lm.first);
			NTT(P,lm.first,-1);//p is needed
			NTT(ar,lm.first,-1);
			for(int i=0;i<lm.first;++i) Q[i]=sub(a[i],ar[i],p);
			clear(Q,m,lm.first);
		}//[a]=n,[b]=m,a%b=p...q
		poly():p(998244353),g(3),gi(332748118) { init(); }
		poly(int mod):p(mod),g(set_mod::getg(mod)),gi(_inv(g,mod)) { init(); }
	};
}

\(\text{poly}\) 提供了一个变幻模数的接口，不过使得常数变大了一倍，如果没有必要写任意模数的话应该拆开

函数的规范是 传入参数，传出参数，数组长度

比如求 \(g(x)=f^-1(x)\pmod {x^n}\) 就是 Inv(f,g,n)，不难注意到这样写我们的数组是左闭右开的，事实上，我们要求所有函数这样（所以取模里面硬减了一下来）

更详细的约定是

  • NTT(a,b,c)，a为需要变幻的多项式，b表示次数，c为 \(-1\) 时进行逆变换，其他时候进行正变幻

  • Inv(a,b,c)，a为需要求逆的多项式，b表示逆元（在 \(\mathbb{Z}_p\) 的意义下），c为次数，之后默认最后一个参数表示数列长度（也意味着是在 \(\bmod {x^n}\) 的意义下）

  • clear(a,b,c)，清空a在 \([l,r)\) 的值

  • arrAdd(a,b,c)，如果b为数列，则将对应位相加，如果b为常数，则a的每一位加上b

  • arrSub(a,b,c)，同上，加法变为减法

  • arrMul(a,b,c)，同上，加法变为乘法

  • Sqrt(a,b,c)，求解 \(b\equiv\sqrt a\pmod{x^n}\)

  • Der(a,b,c)，求解 \(b=a^\prime\)

  • Int(a,b,c)，求解 \(b=\int a\)，常数项置 \(0\)

  • Ln(a,b,c)，求解 \(b=\ln a\)，应当有 \(a_0=1\)

  • Exp(a,b,c)，求解 \(b=e^x\)，应当有 \(a_0=0\)

  • Pow(a,k,b,c)，求解 \(b=a^k\)，应当有 \(a_0=1\)

  • Mod(a,b,p,q,n,m)，a的长度为 \(n\)，b的长度为 \(m\)，求解 \(c=\lfloor\dfrac{a}{b}\rfloor,d=a\bmod b\)，就是小学学的大除法

原理都不难，主要就是多项式牛顿迭代

多项式牛迭是为了解决下面这个问题

\[G(F(x))\equiv 0\pmod{x^n} \]

我们期望用 \(G(x)\) 表示或者求出 \(F(x)\)，当然，它们都是多项式

考虑一个倍增的做法

我们假设已经求出了 \(F_0(x)\) 满足

\[G(F_0(x))\equiv 0\pmod{x^{\frac{n}{2}}} \]

显然有

\[F_0(x)\equiv F(x)\pmod{x^{\frac n2}} \]

把 \(G(F(x))\) 看做关于 \(F(x)\) 的函数，在 \(F_0(x)\) 处展开

可得

\[\begin{align*} G(F(x))&=G(F_0(x))+\frac{G^\prime(F_0(x))}{1!}(F(x)-F_0(x))+\cdots\\ &=\sum_{i=0}^\infty\frac{G^{(i)}(F(x))}{i!}(F(x)-F_0(x))^i \end{align*} \]

我们只关心 \(F(x)\) 在 \(\bmod x^n\) 的意义下的值，而显然有 \(F(x)-F_0(x)\equiv0\pmod{x^\frac{n}{2}}\)

也就是说，\((F(x)-F_0(x))^2\equiv0\pmod {x^n}\)

于是只需要保留泰勒展开的前两项，即有

\[F(x)=F_0(x)-\frac{G(F_0(x)}{G^\prime(F_0(x))} \]

最后 \(n=1\) 的边界一般都可以转化为数论问题，然后就可以倍增了

复杂度一般都是 \(O(n\log n)\)，用主定理看一眼，都是 \(T(n)=O(n\log n)+T(\dfrac{n}{2})\)，当然，还有一些这里没有的操作，有些就不是这个复杂度，如果你写半在线卷积或者分治 \(\text{NTT}\) 的话也不是

符号化方法

咕咕咕咕

标签：lim,return,入门,int,多项式,pmod,简陋,inline
From： https://www.cnblogs.com/efX-bk/p/poly.html

相关文章

  • DOS出初入门学习
    打开CMD方式开始+系统+命令提示符Windows+R打开运行，输入cmd打开控制台（推荐使用）在任意位置按住Shift+鼠标右键点击，在此处打开Powershall窗口资源管理器的地址栏......
  • Salesforce入门课程，清理非活跃用户这8点必须要注意！
     非活跃用户确实给Salesforce管理员带来了诸多问题。创建用户后，可以将其设为非活跃状态，但永远不能将其删除。随着时间的推移，大多数组织的非活跃用户会远多于活跃用户。......
  • C++入门（上）
    ......
  • 7天入门小程序开发 | 07-小程序用户登录及多账号管理
            回顾下之前的内容，我们能够设计小程序页面，能够实现用户与数据库之间的交互，也已经使用了云开发中的云存储、数据库、云函数，能够从前到后实现简单小程序的开发......
  • 3-Air724UG模块(4G全网通GPRS开发)-下载DTU固件和入门使用
    <p><iframename="ifd"src="https://mnifdv.cn/resource/cnblogs/LearnAir724UG"frameborder="0"scrolling="auto"width="100%"height="1500"></iframe></p> 说明......
  • WEB简介与HTTP入门
    一、Web简介1、什么是Web学习Web安全当然要简单的了解什么是Web，Web与生活息息相关，上个网站浏览新闻，看个视频等其中涉及到几个基本的点。从通信，会接触到URL，到协议，......
  • 软件测试入门学习
    caohongxing的博客​软件测试2小时入门​​​https://study.163.com/course/courseMain.htm?courseId=1004794006&trace_c_p_k2=debbdb37dde34011af67c8e4f996a17a​​......
  • Markdown入门
    Markdown轻松入门刚注册的博客园来记录小白成长笔记，第一天先来打个卡吧！Markdown是一种语法，这种语法需要Typora编辑器结合，有好多编辑器都支持这种语法，例如印象笔记、博客......
  • 7.函数入门
    函数函数的作⽤函数的使⽤步骤函数的参数作⽤函数的返回值作⽤函数的说明⽂档函数嵌套函数是组织好的，可重复使用的，用来实现单一，或相关联功能的代码段。函数能......
  • 【每周CV论文推荐】 深度学习人脸检测入门必读文章
    欢迎来到《每周CV论文推荐》。在这个专栏里，还是本着有三AI一贯的原则，专注于让大家能够系统性完成学习，所以我们推荐的文章也必定是同一主题的。人脸图像是整个图像处理领域里......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99