zoukankan      html  css  js  c++  java
  • 洛谷 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;
    }
    
  • 相关阅读:
    51nod 1113 矩阵快速幂 如题目
    poj Raising Modulo Numbers 快速幂模板(取膜)
    bzoj 1503: [NOI2004]郁闷的出纳员 平衡树
    codevs 1063 合并果子 优先队列相关
    bzoj 3224: Tyvj 1728 普通平衡树 Treap模版
    快排模板
    hdu 4353 统计点在三角形内的个数
    hdu 3264 圆的交+二分
    hdu 3685 多边形重心+凸包
    hdu 3992 AC自动机上的高斯消元求期望
  • 原文地址:https://www.cnblogs.com/ET2006/p/luogu-P3781.html
Copyright © 2011-2022 走看看