多项式相关

我学科技?真的假的?

起因是在做斯特林数,发现没事就要佛佛塔或者呢塔塔,遂学之。

FFT

多项式乘法

现在要求两个多项式的卷积。

卷积

现在给你两个序列 ai,bi 并定义 ci

ci=p+q=iaibi

可以认为 ci 就是 ai,bi 的卷积。换句话说就是给你两个多项式 F(x)=i=0naixi,G(x)=i=0nbixi,然后 H(x)=F(x)G(x)=i=02ncixi,手玩一下二次式乘二次式,算出来的多项式就是卷积结果,多项式的 i 次式就是 ci

点值表示法

现在让你算两个多项式的乘法就可以 O(n2) 算了,但是不够快,考虑优化这个过程,但是发现难以优化。

前人的智慧指出:对任意 n 个数点(x不同)都可以唯一地插值出来一个对应的 n 次多项式,既然 H=(FG) 则对于原先的数点必定有 H(x0)=F(x0)G(x0)。确立这些 H(x) 的数点的过程应当是 O(n) 的,如果现在可以用一个比较高速的方法插值出 H(x) 就可以实现复杂度优化了。

单位根

如果你在看这段建议先学复数及其复平面上的三角表示和加减乘除于复平面上的表示。

考虑在复平面内给单位圆劈成 n 份,x 正方向向上的第一刀对应的向量 ωn 肯定满足 ωnn=1,因为转了一圈,就把这个东西叫单位根。

这个东西肯定就有一些性质比如

ω2n2k=ωnk

砍两倍次数每次拼两个上去和砍一倍每次拼一个一样...

ω2nk=ω2nk+n

因为转了一半过去。

分治法法塔(FFT)

现在有这样一个多项式

F(x)=a0+a1x+a2x2+...an1xn1

我直接给他奇偶劈两半,不妨认为 n 是偶数。

=(a0+a2x2+...an2xn2)+(a1x+a3x3+...an1xn1)

再设

F1(x)=a0+a2x+a4x2+...an2xn21

F2(x)=a1+a3x+a5x2+...+an1xn21

然后就有

F(x)=F1(x2)+xF2(x2)

那很明显这个 x 太丑了直接给他换成幂次等效于旋转的单位根。

F(ωnk)=F1(ωn2k)+ωnkF2(ωn2k)

=F1(ωn/2k)+ωnkF2(ωn/2k)

同理

F(ωnk+n/2)=F1(ωn/2k)ωnkF2(ωn/2k)

你知道这个玩意有多牛逼吗就是说 F(ωnk)F(ωnk+n/2) 是可以同时算出来的,加号改减号即可而且这个玩意跨度正好是 n/2

直接考虑一个递归分治状物然后根据6定理得到时间复杂度为 O(nlogn)

蛋是这个玩意常数大到飞起是无法通过板题的,至少复杂度对了常数优化都是后事。

先考虑这样一个问题:这样kuku算完算出来的是一堆 F(x) 啊插值的事情还没解决。

挨法法塔(IFFT)

现在要解决如何把点值表示法转换为一般的系数表示法,而且要充分利用单位根优势,要不然不如去考虑拉插求系数。

然后前人的智慧再次发力提出一个牛逼的构造,定义

ck=i=0n1yi(ωnk)i

yi 就是刚才算出来的那坨 F()

然后这个东西可以化。

ck=i=0n1yi(ωnk)i

=i=0n1(j=0n1aj(ωni)j)(ωnk)i

=j=0n1aji=0n1(ωnjk)i

然后先放在这,设 S(x)=i=0n1xi 并带入 ωnk 然后构造等比数列求和式

ωnkS(ωnk)S(ωnk)=(ωnk)n1

S(ωnk)=0

但是这是建立在 ωnkn=1 的基础上的,k=0 时原式为 n

所以这个

j=0n1aji=0n1(ωnjk)i

当且仅当 j=k 时为 n,否则为 0

进而 ak=ckn

所以流程就是:FFT,点乘计算,反着FFT(带个-1),然后除以 n,即次数。

放一下板子代码方便以后照搬

#include<bits/stdc++.h>
#define db double
#define MAXN 10000005
using namespace std;
int n,m,pw,len=1,loc[MAXN];
const db pi=acos(-1.0);
struct comp{
	db x,y;
	comp(db a=0,db b=0){x=a,y=b;}
}a[MAXN],b[MAXN];
comp operator+(comp a,comp b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(comp a,comp b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(comp a,comp b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void FFT(comp *A,int typ){
	for(int i=0;i<len;i++)if(i<loc[i])swap(A[i],A[loc[i]]);
	for(int mid=1;mid<len;mid<<=1){
		comp wn(cos(pi/mid),typ*sin(pi/mid));
		for(int r=mid<<1,j=0;j<len;j+=r){
			comp w(1,0);
			for(int k=0;k<mid;k++,w=w*wn){
				comp x=A[j+k],y=w*A[j+k+mid];
				A[j+k]=x+y;
				A[j+k+mid]=x-y;
			}
		}
	}
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++)scanf("%lf",&a[i].x);
	for(int i=0;i<=m;i++)scanf("%lf",&b[i].x);
	while(len<=n+m)len<<=1,++pw;
	for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)<<(pw-1));
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<=len;i++)a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/len+0.5));
	return 0;
}

NTT

就是把单位根换成原根。这个东西是用来解决带取模卷积的,一般是mod998244353。

namespace POLY{
	int loc[MAXN];
	const int Gx=3,GI=332748118;
	struct Polynomial{
		int len;
		vector<int>a;
		int& operator[](int x){return a[x];}
		inline void fix(int x){
			len=x;
			a.resize(x+1);
		}
		inline void reset(){
			a.clear();len=0;
		}
		inline void NTT(int n,int typ){
			fix(n-1);
			for(int i=0;i<n;i++)if(i<loc[i])swap(a[i],a[loc[i]]);
			for(int len=1;len<n;len<<=1){
				int step=qpow(typ==1?Gx:GI,(mod-1)/(len<<1));
				for(int l=0;l<n;l+=(len<<1)){
					int w=1;
					for(int k=0;k<len;k++,w=w*step%mod){
						int x=a[l+k],y=w*a[l+len+k]%mod;
						a[l+k]=(x+y)%mod;
						a[l+len+k]=(x-y)%mod; 
					}
				}
			}
			if(typ==-1){
				int inv=qpow(n,mod-2);
				for(int i=0;i<n;i++)a[i]=a[i]*inv%mod;
			}
		}
	};
	Polynomial operator*(Polynomial F,Polynomial G){
		Polynomial res;
		int n=F.len+G.len;
		int len=1,pw=0;
		while(len<=n)len<<=1;
		res.fix(len-1);
		for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)*(len>>1));
		F.NTT(len,1);G.NTT(len,1);
		for(int i=0;i<len;i++)res[i]=F[i]*G[i]%mod;
		res.NTT(len,-1);
		res.fix(n);
		return res;
	}	
}

多项式相关的题感觉主角都不是多项式而是推式子...

法法塔二

ck=ikn1aibik

这个长得就特别卷积啊想想怎么转化成卷积形式。就是说

ck=i=kn1aibik

=i=0n1kai+kbi

这个加号就很恶心所以考虑设 ai=an1i,这个reverse一下就好了。

进而

=i=0n1kan1ikbi

然后直接开卷就行了。

乍一看和多项式没有任何关系啊,然后就要开始考虑构造了。

首先原式的上下指标都能稍微改一改,改成到 j 和从 j 开始的,不过不重要。

Ei=j=1iqj(ij)2j=inqj(ij)2

然后考虑构造 f(x)=qx,g(x)=1i2

则原式变为

=j=1if(j)g(ij)j=inf(j)g(ji)

前半部分可以直接卷,考虑后半部分的处理。

根据上一题的经验尝试替换指标和反转

j=inf(j)g(ji)

=j=0nif(j+i)g(j)

f(j)=f(nj)

=j=0nif(nij)g(j)

然后也可以卷了。

Normal

神秘综合题。

首先进行一波推式子。考虑对于一个点 x 什么时候会贡献答案。相当于对于当前的分治中心 y,当且仅当 xy 的路径上还没有点成为过分治中心时(即在子树中)会贡献答案,进而要求路径上不能选点,分治中心要选其他随意,这时有 1dis(x,y)+1 的概率对答案做出 sizy 上的 x1 贡献。

所以答案即为

i=1nj=1n1dis(i,j)+1

因为距离跑到分母上了所以答案没法简单统计...

考虑点分治,对于当前分治中心 u 会得到若干距离不同的路径,如果知道不同距离路径的个数就可以简单计算了。

但是点分治不同子树枝杈间可以靠 u 两两合并路径,发现这一过程恰好是多项式卷积的形式(adisxmergebdisyabdisx+y)。但是还没完,如果有多个子树直接就假掉了...前两个子树卷积完再来一个就要卷积一个跨三个子树的路径了,显然是错的。

转而考虑容斥,用在 u 中任意选两个的方案减去所有单个子树中任选两个的不合法方案就是当前分治中心的贡献,发现正好可以变为 u 子树自己卷积减去所有它子树内自己卷积的结果,用fft处理。

然后分析复杂度,因为每次卷积的多项式大小和点分树子树大小一致所以复杂度是 T(n)=O(nlogn)+2(T(n2)+O(n2logn))=O(nlog2n) 的。

#include<bits/stdc++.h>
#define db double
#define MAXN 65005
using namespace std;
int n,pw,len=1,loc[MAXN];
db ans;
const db pi=acos(-1.0);
struct comp{
	db x,y;
	comp(db a=0,db b=0){x=a,y=b;}
}a[MAXN],f[MAXN],b[MAXN];
db P1[MAXN],P2[MAXN];
comp operator+(comp a,comp b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(comp a,comp b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(comp a,comp b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
struct Polynomial{
	int len;
	vector<comp>a;
	comp& operator[](int x){return a[x];}
	inline void fix(int x){
		len=x;
		a.resize(x+1);
	}
	inline void reset(){
		a.clear();len=0;
	}
	inline void FFT(int n,db typ){
		fix(n-1);
		for(int i=0;i<n;i++)if(i<loc[i])swap(a[i],a[loc[i]]);
		for(int mid=1;mid<n;mid<<=1){
			comp wn(cos(pi/mid),typ*sin(pi/mid));
			for(int r=mid<<1,j=0;j<n;j+=r){
				comp w(1,0);
				for(int k=0;k<mid;k++,w=w*wn){
					comp x=a[j+k],y=w*a[j+k+mid];
					a[j+k]=x+y;
					a[j+k+mid]=x-y;
				}
			}
		}
		if(typ==-1)for(int i=0;i<n;i++)a[i].x/=(db)n;
	}
};
Polynomial operator*(Polynomial F,Polynomial G){
	Polynomial res;
	int n=F.len+G.len;
	int len=1,pw=0;
	while(len<=n)len<<=1;
	res.fix(len-1);
	for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)*(len>>1));
	F.FFT(len,1);G.FFT(len,1);
	for(int i=0;i<len;i++)res[i]=F[i]*G[i];
	res.FFT(len,-1);
	res.fix(n);
	return res;
}
struct node{
	int v,nxt;
}edge[MAXN<<1];
int h[MAXN],tmp;
inline void add(int u,int v){
	edge[++tmp]=(node){v,h[u]};
	h[u]=tmp;
}
int siz[MAXN],dp[MAXN],sum,root;
bool vis[MAXN];
inline void getrt(int u,int fa){
	siz[u]=1;
	dp[u]=0;
	for(int i=h[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(v==fa||vis[v])continue;
		getrt(v,u);
		siz[u]+=siz[v];
		dp[u]=max(dp[u],siz[v]);
	}
	dp[u]=max(dp[u],sum-dp[u]);
	if(dp[u]<dp[root])root=u;
}
inline void getdis(int u,int fa,int d,Polynomial &F){
	++F[d].x;
	for(int i=h[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(v==fa||vis[v])continue;
		getdis(v,u,d+1,F);
	}
}
inline void work(int u){
	for(int i=h[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(vis[v])continue;
		Polynomial son;
		son.fix(siz[v]);
		getdis(v,u,1,son);
		son=son*son;
		for(int i=0;i<=son.len;i++)ans-=son[i].x/(i+1);
	}
	Polynomial son;
	son.fix(siz[u]);
	getdis(u,0,0,son);
	son=son*son;
	for(int i=0;i<=son.len;i++)ans+=son[i].x/(i+1);
}
inline void solve(int u){
	vis[u]=1;
	//printf("solving nowu=%d\n",u);
	work(u);
	for(int i=h[u];i;i=edge[i].nxt){
		int v=edge[i].v;
		if(vis[v])continue;
		sum=siz[v],root=0;
		getrt(v,u);
		solve(root);
	}
}
signed main(){
	scanf("%d",&n);
	for(int i=1,u,v;i<n;i++){
		scanf("%d%d",&u,&v);
		++u,++v;
		add(u,v);
		add(v,u);
	}	
	root=0,sum=dp[root]=n;
	getrt(1,0);
	solve(root);
	printf("%.4f",ans);
	return 0;
}

BZOJ3771

发现丢失斧头的价值求和可以用多项式卷积的指数表示(恰好斧头的价值都是唯一的),系数则为方案。

考虑去重,分讨丢了一把,两把和三把斧头的情况。

对于丢了一把斧头的情况则为给定的 ai

对于两把斧头的情况:因为不能重复地丢一把斧头,设输入给定的数据表示的多项式为 F,令 G 为仅重复选两次表示出的多项式(即将F的所有次数乘二)则贡献为 FFG

对于丢三把斧头的情况,这时候不合法的情况可以为重复选三把斧头或者重复选两把斧头,保留上文 F,G 的基础上设 H 是仅重复选三次得到的多项式。总方案为 FFF,减去选两个一样且第三个不一样的方案 FGH,但是这样的情况在三个f中会重复出现三次所以乘三,最后减去选三个一样的,贡献为 FFF3(FGH)H

这个题自排列是不合法的所以丢多于一个的给方案除个自排列就行。

有点卡精,建议统计答案全用double最后输出的时候再转int。

求和

众所周知第二类斯特林数有通项公式

S(n,m)=1m!i=0m(1)i(mi)(mi)n=i=0m(1)i(mi)ni!(mi)!

考虑转化原式为

f(n)=j=0n2jj!i=jnS(i,j)

=j=0n2jj!i=jnk=0j(1)k(jk)ik!(jk)!

=j=02jj!k=0j(1)ki=0n(jk)ik!(jk)!

考虑拆分后面的分式为卷积形式,设

F(x)=(1)xx!,G(x)=i=0nxix!=xn+11(x1)x!

Ans=j=0n2jj!k=0jF(k)G(jk)

预处理一下开卷就完事了。

然后可以顺手捏掉第二类斯特林数行。

万径人踪灭

挺有意思的。

答案相当于是回文子序列个数-回文串个数,后者好说主要看回文子序列怎么求。

两个字符回文的前提是相等且关于轴线对称,所有满足对称的位置对 (x1,x2) 都有 x1+x2=20,恰好满足卷积的下标构成。

考虑分别处理出所有 ab 位置的多项式并自己卷积,则卷积后的第 i 项系数 ai 就应该是关于 i/2 位置对称的字符对个数,对答案的贡献就应是 2ai2

礼物

设现在往上又给了 x

Ans=i=1n(aibi+x)2

=i=1n(ai2+bi2)+nx2+2(ab)xi=1naibi

发现除了 i=1naibi 都可以 O(m) 求出。

根据以往的经验把式子变成 ai=ani

i=1naibi=i=0nanibi

这样做的原因是对手环的旋转等效于将任一元素倍长然后旋转 k 就变为

ck=i=0n+kan+kibi

所以翻转并倍长一个序列然后卷积求 [n,2n) 的系数极大值即可。

BZOJ5093

对单点的贡献:相当于枚举度数计算,剩下 n1 个点随便连边即可,n 个点就是式子乘 n

Ans=ni=0n1(n1i)ik2(n12)

=n2(n12)j=0kS(k,j)i=0n1(n1i)ij

=n2(n12)j=0kS(k,j)i=0n1(n1jij)nj

=n2(n12)j=0knj2nj1S(k,j)

然后就是处理一行斯特林数

S(i,j)=k=0j(1)k(jk)ik!(jk)!

=k=0jF(k)G(jk),F(x)=(1)xx!,G(x)=xix!

ntt卷一卷就行。

posted @   Cl41Mi5deeD  阅读(11)  评论(1编辑  收藏  举报
相关博文:
阅读排行:
· 全程不用写代码,我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· 记一次.NET内存居高不下排查解决与启示
· MongoDB 8.0这个新功能碉堡了,比商业数据库还牛
· 白话解读 Dapr 1.15:你的「微服务管家」又秀新绝活了
点击右上角即可分享
微信分享提示