洛谷 P3781

洛谷题面传送门

SDOI 2017 R2 D1 T3,nb tea %%%

讲个笑话,最近我在学动态 dp,wjz 在学 FWT,而我们刚好在同一天做到了这道题,而这道题刚好又是 FWT+动态 dp

首先考虑怎样暴力计算答案,我们记 (dp_{u,j}) 表示以 (u) 为根的子树中有多少个连通块包含 (u) 且权值的异或和为 (j),初始 (dp_{u,val_u}=1),每次遍历 (u) 的一个子树 (v) 就对这个子树就对这两个子树的 (dp) 做一个合并,即 (dp_{u,x}leftarrow dp_{u,x}+sumlimits_{y=0}^{m-1}dp_{u,y} imes dp_{v,xoplus y}),最终答案即为 (sumlimits_{u}dp_{u,k})。正确性显然,时间复杂度 (mathcal O(qnm^2)),可以通过前四个测试点。

考虑优化,首先一个非常明显的优化是,DP 转移方程式长得一脸 xor 卷积的样子,如果我们记 (f*g) 表示 (f,g) 两个集合幂级数的 FWTxor,那么上述式子可以改写为 (dp_u=dp_u+dp_u*dp_v=dp_u*(dp_v+1)),因此考虑将所有 (dp_u) 都变为 ( ext{FWT}(dp_u)),那么 (dp_u) 的初始值就变为 (dp_{u,i}= ext{FWT}(E_{val_u})_i),其中 (E_i) 为满足 (f_i=1,f_j=0(j e i)) 的集合幂级数 (f),这个可以通过预处理所有 ( ext{FWT}(E_i)) 实现 (mathcal O(m)) 初始化。转移操作可以根据 FWT 那一套理论变成 (dp_{u,i}leftarrow dp_{u,i} imes(dp_{v,i}+1)),这样即可实现 (mathcal O(m)) 转移。然后每次操作完了之后再 IFWT 回来即可,时间复杂度降到了 (mathcal O(qnm+qmlog m)),还是只能通过前四个测试点(

注意到上述 (dp) 对于不带修改的情况是 efficient enough 的,但是带上修改就直接萎掉了,因此考虑擅长处理修改操作的动态 (dp) 来解决这个问题,按照动态 (dp) 的套路我们将树剖成一条条重链,(dp) 分为轻儿子和重儿子处理,我们记 (dpl_{u,i}=sumlimits_{vin ext{lightson}(u)}(dp_{v,i}+1)),那么记 (w=wson_u),则有 (dp_{u,i}=dpl_{u,i} imes ext{FWT}(E_{val_u})_i imes (dp_{w,i}+1))

但是光记录一个 (dp) 值是远远不够的,因为最终我们要求的是整棵子树中 (dp_{u,k}) 的值之和,所有我们不得不再额外记录 (sum_{u,i}) 表示子树中所有点的 (dp_{u,i}) 之和,那么有 (sum_{u,i}=sumlimits_{vin ext{son}(u)}sum_{v,i}+dp_{u,i}),按照套路我们还是记 (suml_{u,i}=sumlimits_{vin ext{lightson}(u)}sum_{v,i}),那么有 (sum_{u,i}=sum_{w,i}+dp_{u,i}+suml_{u,i}=sum_{w,i}+dpl_{u,i} imes ext{FWT}(E_{val_u})_i imes(dp_{w,i}+1)+suml_{u,i})

考虑将这东西写成矩阵的形式,那么有:

[egin{bmatrix}dp_{u}&sum_{u}&1end{bmatrix}=egin{bmatrix}dp_{w}&sum_{w}&1end{bmatrix} imes egin{bmatrix} dpl_u imes ext{FWT}(E_{val_u})&dpl_u imes ext{FWT}(E_{val_u})&0\ 0&1&0\ dpl_u imes ext{FWT}(E_{val_u})&dpl_u imes ext{FWT}(E_{val_u})+suml_u&1 end{bmatrix} ]

其中 (f imes g) 就对应项相乘好了,(f+g) 也同理。

(A_u=egin{bmatrix} dpl_u imes ext{FWT}(E_{val_u})&dpl_u imes ext{FWT}(E_{val_u})&0\ 0&1&0\ dpl_u imes ext{FWT}(E_{val_u})&dpl_u imes ext{FWT}(E_{val_u})+suml_u&1 end{bmatrix}),那么对于一个点 (u) 而言,记它到重链底经过的节点依次是 (u=v_1,v_2,cdots,v_k),那么有

[egin{bmatrix}dp_{u}&sum_{u}&1end{bmatrix}=egin{bmatrix}0&0&1end{bmatrix} imesprodlimits_{i=k}^1A_{v_i} ]

这个可以树链剖分+线段树维护。

修改操作就按照动态 (dp) 的套路不断跳重链并撤销原来的 (dp_{top_u})(dpl_{fa[top_u]})(suml_{fa[top_u]}) 的影响并加入新的贡献即可,时间复杂度 (qlog^2nm+qmlog m),LOJ 上可以通过,而洛谷上由于某两位毒瘤提供的毒瘤卡树剖的数据,只能获得 (80) 分的好成绩。

说起来轻巧,实现起来一堆细节需要注意:

  1. 直接矩阵乘法会多 (27) 的常数,导致 TLE,因此需要按照套路进行优化,注意到这个 (3 imes 3) 的矩阵中只有四个位置是有用的,因此可以只维护这四个位置的值,即 (egin{bmatrix}a&b&0\0&1&0\c&d&1end{bmatrix}),那么有 (egin{bmatrix}a_1&b_1&0\0&1&0\c_1&d_1&1end{bmatrix} imes egin{bmatrix}a_2&b_2&0\0&1&0\c_2&d_2&1end{bmatrix}=egin{bmatrix}a_1a_2&a_1b_2+b_1&0\0&1&0\a_2c_1+c_2&b_2c_1+d_1+d_2&1end{bmatrix}),这样常数可以降到 (4)
  2. 注意线段树 pushup 的顺序。
  3. 在撤销原来的贡献时会出现除以 (0) 的情况,因此可以将 (dpl_{u,i}) 存成一个个结构体,每个结构体用 (x imes 0^y) 表示一个数,每次乘以 (0) 时令 (y) 加一,除以 (0) 则令 (y) 减一,这样可以避免这个问题(u1s1 蒟蒻是第一次遇到这个套路呢,大佬不喜勿喷)
  4. 注意计算新加入的贡献时是计算线段树上 ([dfn[top[x]]],dfn[bot[top[x]]]) 内矩阵的乘积,而不是 ([dfn[x]],dfn[bot[top[x]]]),蒟蒻因为这个错误调了 1h,心态炸裂。

码了 212 行……

const int MAXN=3e4;
const int MAXV=1<<7;
const int MOD=1e4+7;
const int INV2=5004;
int n,m,val[MAXN+5],inv[MOD+4];
void getinv(){
	for(int i=(inv[0]=inv[1]=1)+1;i<MOD;i++) inv[i]=inv[MOD%i]*(MOD-MOD/i)%MOD;
}
void FWTxor(int *a,int len,int type){
	for(int i=2;i<=len;i<<=1)
		for(int j=0;j<len;j+=i)
			for(int k=0;k<(i>>1);k++){
				int X=a[j+k],Y=a[(i>>1)+j+k];
				if(~type) a[j+k]=(X+Y)%MOD,a[(i>>1)+j+k]=(X-Y+MOD)%MOD;
				else a[j+k]=(X+Y)*INV2%MOD,a[(i>>1)+j+k]=(X-Y+MOD)*INV2%MOD;
			}
}
struct num0{//number expressed as x*0^y
	int x,y;
	num0(int v=1){(!v)?(y=x=1):(x=v,y=0);}
	num0 operator *(const int &rhs){
		(!rhs)?(++y):(x=x*rhs%MOD);
		return *this;
	}
	num0 operator /(const int &rhs){
		(!rhs)?(--y):(x=x*inv[rhs]%MOD);
		return *this;
	}
	int num(){return y?0:x;}
};
struct poly{
	int a[MAXV+5];
	poly(){memset(a,0,sizeof(a));}
	poly(int x){for(int i=0;i<m;i++) a[i]=x;}
	poly operator +(poly rhs) const{
		poly res;
		for(int i=0;i<m;i++) res.a[i]=(a[i]+rhs.a[i])%MOD;
		return res;
	}
	poly operator *(poly rhs) const{
		poly res(1);
		for(int i=0;i<m;i++) res.a[i]=a[i]*rhs.a[i]%MOD;
		return res;
	}
	void FWT(){FWTxor(a,m,1);}
	void IFWT(){FWTxor(a,m,-1);}
} e[MAXV+5];
int hd[MAXN+5],to[MAXN*2+5],nxt[MAXN*2+5],ec=0;
void adde(int u,int v){to[++ec]=v;nxt[ec]=hd[u];hd[u]=ec;}
int siz[MAXN+5],fa[MAXN+5],dep[MAXN+5],wson[MAXN+5];
int top[MAXN+5],dfn[MAXN+5],tim=0,rid[MAXN+5];
int bot[MAXN+5];
void dfs1(int x=1,int f=0){
	siz[x]=1;fa[x]=f;
	for(int e=hd[x];e;e=nxt[e]){
		int y=to[e];if(y==f) continue;
		dep[y]=dep[x]+1;dfs1(y,x);siz[x]+=siz[y];
		if(siz[y]>siz[wson[x]]) wson[x]=y;
	}
}
void dfs2(int x=1,int tp=1){
	top[x]=tp;rid[dfn[x]=++tim]=x;if(wson[x]) dfs2(wson[x],tp);
	for(int e=hd[x];e;e=nxt[e]){
		int y=to[e];if(y==wson[x]||y==fa[x]) continue;
		dfs2(y,y);
	}
}
int f[MAXN+5][MAXV+5],sum[MAXN+5][MAXV+5],suml[MAXN+5][MAXV+5];
num0 fl[MAXN+5][MAXV+5];
void dfs3(int x=1){
	for(int i=0;i<m;i++) f[x][i]=e[val[x]].a[i],fl[x][i]=1;
	for(int e=hd[x];e;e=nxt[e]){
		int y=to[e];if(y==fa[x]) continue;dfs3(y);
		for(int i=0;i<m;i++) f[x][i]=f[x][i]*(f[y][i]+1)%MOD;
		for(int i=0;i<m;i++) sum[x][i]=(sum[x][i]+sum[y][i])%MOD;
	} for(int i=0;i<m;i++) sum[x][i]=(sum[x][i]+f[x][i])%MOD;
}
void dfs4(int x=1){
	if(wson[x]) dfs4(wson[x]);
	for(int e=hd[x];e;e=nxt[e]){
		int y=to[e];if(y==fa[x]||y==wson[x]) continue;dfs4(y);
		for(int i=0;i<m;i++) fl[x][i]=fl[x][i]*((f[y][i]+1)%MOD);
		for(int i=0;i<m;i++) suml[x][i]=(suml[x][i]+sum[y][i])%MOD;
	}
}
struct mat{
	poly a,b,c,d;
	mat(){}
	mat operator *(const mat &rhs){
		mat res;res.a=a*rhs.a;res.b=b+a*rhs.b;
		res.c=rhs.a*c+rhs.c;res.d=rhs.b*c+d+rhs.d;
		return res;
	}
};
void print(mat x){
	for(int i=0;i<m;i++) printf("%d%c",x.a.a[i]," 
"[i==m-1]);
	for(int i=0;i<m;i++) printf("%d%c",x.b.a[i]," 
"[i==m-1]);
	for(int i=0;i<m;i++) printf("%d%c",x.c.a[i]," 
"[i==m-1]);
	for(int i=0;i<m;i++) printf("%d%c",x.d.a[i]," 
"[i==m-1]);
}
mat get(int x){
	mat res;
	for(int i=0;i<m;i++) res.a.a[i]=res.b.a[i]=res.c.a[i]=fl[x][i].num()*e[val[x]].a[i]%MOD;
	for(int i=0;i<m;i++) res.d.a[i]=(res.a.a[i]+suml[x][i])%MOD;
	return res;
}
struct node{int l,r;mat v;} s[MAXN*4+5];
void pushup(int k){s[k].v=s[k<<1|1].v*s[k<<1].v;}
void build(int k,int l,int r){
	s[k].l=l;s[k].r=r;if(l==r) return s[k].v=get(rid[l]),void();
	int mid=l+r>>1;build(k<<1,l,mid);build(k<<1|1,mid+1,r);pushup(k);
}
mat query(int k,int l,int r){
	if(l<=s[k].l&&s[k].r<=r) return s[k].v;
	int mid=s[k].l+s[k].r>>1;
	if(r<=mid) return query(k<<1,l,r);
	else if(l>mid) return query(k<<1|1,l,r);
	else return query(k<<1|1,mid+1,r)*query(k<<1,l,mid);
}
void modify(int k,int p){
	if(s[k].l==s[k].r) return s[k].v=get(rid[p]),void();
	int mid=s[k].l+s[k].r>>1;
	if(p<=mid) modify(k<<1,p);else modify(k<<1|1,p);
	pushup(k);
}
void change(int x){
	while(x){
		if(fa[top[x]]){
			mat res=query(1,dfn[top[x]],dfn[bot[top[x]]]);
//			printf("%d
",fa[top[x]]);
//			for(int i=0;i<m;i++) printf("{%d,%d}%c",fl[fa[top[x]]][i].x,fl[fa[top[x]]][i].y," 
"[i==m-1]);
//			for(int i=0;i<m;i++) printf("%d%c",(res.c.a[i]+1)%MOD," 
"[i==m-1]);
//			print(get(fa[top[x]]));
			for(int i=0;i<m;i++) fl[fa[top[x]]][i]=fl[fa[top[x]]][i]/((res.c.a[i]+1)%MOD);
			for(int i=0;i<m;i++) suml[fa[top[x]]][i]=(suml[fa[top[x]]][i]-res.d.a[i]+MOD)%MOD;
		} modify(1,dfn[x]);
		if(fa[top[x]]){
			mat res=query(1,dfn[top[x]],dfn[bot[top[x]]]);
//			print(res);
			for(int i=0;i<m;i++) fl[fa[top[x]]][i]=fl[fa[top[x]]][i]*((res.c.a[i]+1)%MOD);
			for(int i=0;i<m;i++) suml[fa[top[x]]][i]=(suml[fa[top[x]]][i]+res.d.a[i])%MOD;
//			for(int i=0;i<m;i++) printf("{%d,%d}%c",fl[fa[top[x]]][i].x,fl[fa[top[x]]][i].y," 
"[i==m-1]);
//			for(int i=0;i<m;i++) printf("%d%c",(res.c.a[i]+1)%MOD," 
"[i==m-1]);
//			print(get(fa[top[x]]));
		} x=fa[top[x]];
	}
}
int main(){
	scanf("%d%d",&n,&m);getinv();
	for(int i=0;i<m;i++) e[i].a[i]=1,e[i].FWT();
//	for(int i=0;i<m;i++) for(int j=0;j<m;j++) printf("%d%c",e[i].a[j]," 
"[j==m-1]);
	for(int i=1;i<=n;i++) scanf("%d",&val[i]);
	for(int i=1,u,v;i<n;i++) scanf("%d%d",&u,&v),adde(u,v),adde(v,u);
	dfs1();dfs2();dfs3();dfs4();build(1,1,n);
//	for(int i=1;i<=n;i++) for(int j=0;j<m;j++) printf("%d%c",f[i][j]," 
"[j==m-1]);
//	for(int i=1;i<=n;i++) for(int j=0;j<m;j++) printf("%d%c",sum[i][j]," 
"[j==m-1]);
//	for(int i=1;i<=n;i++) for(int j=0;j<m;j++) printf("%d%c",fl[i][j].num()," 
"[j==m-1]);
//	for(int i=1;i<=n;i++) for(int j=0;j<m;j++) printf("%d%c",suml[i][j]," 
"[j==m-1]);
//	for(int i=0;i<m;i++) printf("%d%c",t.d.a[i]," 
"[i==m-1]);
	for(int i=1;i<=n;i++) if(top[i]==i){
		int cur=i;while(wson[cur]) cur=wson[cur];
		bot[i]=cur;
	} int qu;scanf("%d",&qu);
	while(qu--){
		char opt[9];scanf("%s",opt+1);
		if(opt[1]=='C'){
			int x,v;scanf("%d%d",&x,&v);
			val[x]=v;change(x);
		} else {
			int k;scanf("%d",&k);
			mat res=query(1,dfn[1],dfn[bot[1]]);
//			for(int i=0;i<m;i++) printf("%d%c",res.d.a[i]," 
"[i==m-1]);
			res.d.IFWT();
			printf("%d
",res.d.a[k]);
		}
	}
	return 0;
}
原文地址:https://www.cnblogs.com/ET2006/p/luogu-P3781.html