国家集训队2018论文集《最小方差生成树》命题报告 学习笔记/生成树
国家集训队$2018$论文集《最小方差生成树》命题报告 学习笔记/生成树
这几天出的第$4$道论文题,总的来说,代码量适中$(500+)$,远远小于前天代码$(1300+),$论文可读性强,知识点基本掌握,是一道不可多得的不可做题,痛斥出题人$!$
题意大概是这样,给定一张图,求出一个价值最小的生成树
生成树的代价计算方式为生成树边权的方差
拿到这道题就去推式子,想搞出一个能够直接用式子表示权值直接生成树,结果失败,然后考虑计算方差最小,肯定是选的越接近越好,继续思考无果
这道题在考场上可以拿到$12+8$的子任务分数
暴力解法
这个$trick$就很妙.就是说每次枚举平均数,然后做一遍最小生成树,是能保证最优的还是最优
因为方差在平均值等于$u$时候取到最小,每个都是这样,感性理解的话,无论怎么偏,都不会更优
毕竟是论文题,我们需要挖掘算法的本质问题,我们现在的问题是,不一定是选最大的边好,也不一定是选最小的边好,那选相邻的一部分是不会更差的,显然离散程度会比选旁边的小
显然这里不能随便排序,一个边的边权因为和平均值有关,所以这个排序是不确定的
z这里观察一下,什么时候排名会变,平均值经过两数平均值的时候会改变,那么一共$m^2$个端点,那么任意两个端点内部是不会发生变化的,直接$m^2$个最小生成树就好了,复杂度$m^3log(m)$,这个过程是可以优化的,可以考虑用$LCT$调整最小生成树,那就用新的小边去替换大边,具体实现不是很会,大概也会了,无非是看边的排名有没有变化,然后加进去看一下是否更优,只不过使用的$LCT$维护形态,$O(m^2logm)$常数巨大...
引出另一道题$Topcoder \ SRM\ 611\ Ega\ itarianism2$
$N$个点的完全图求最小生成树,边多的时候,用$prim$是更优的,$O(n^6)[n<=20]$
那么对于原题,其实还能继续优化,我们在这已经有了一个$O(m^3logm)$的优秀算法(在这里不考虑$LCT$解法,毕竟常数好大,$splay$常数都$20+$了,还是过不了),我们考虑一下,是不是还有些东西没有充分利用,的确有一部分,区间之间转移的是不是不用每次重新跑最小生成树
这里给出结论$:$
这个结论挺强的,就是说一条边如果存在于生成树中,一定会出现在一段连续的区间的生成树里面,感性想一想,一个边什么时候会消失.他的排名发生变化了,或者其余的边发生变化了,导致被替代,放到二次函数就很显然了
引理
在确定$u$的情况下,如果边$i$在所有权值比他小的边都放进去之后两端点联通,那么就不在生成树里面,否则在生成树里面
这样的话,我们可以把函数图像全都画出来,然后找出$l_i,r_i$,怎么找
我们只需要找到平均值区间让他在里面就好了,枚举每条边,然后把边权按平均值排个序,表示变化点,然后左端点就是找一个最小的平均值使得这条边恰好不用放进去,就是把所有小的边都放进去,并查集维护连通性就好了,小的边怎么找,就是把左函数单独提取出来,找到所有交点,从大到小排,看加到哪个就不用这条边了,那么就找到了$l$,至于$r$是以一样的,为啥是从大到小排序,显然的,越近的话排名越小
我们现在得到了每个边左右区间,这个时候区间就少了,只有$m$个,然后计算每个区间的生成树就好$tql$,复杂度$O(m^2logm)$
复杂度瓶颈在于求$l_i,r_i,$可以考虑从小到大加边,这样是没问题的,以为这样保证了上面要求的按平均值排序,那么这个每个边该不该加进去就好判断了,加边然后$LCT$求最大生成树就好了,还是说这里是大常数
这个可以用一个结论优化一部分,边$i,j$,如果最大生成树$u_i,v_i$的最短边$j$,那么$l_i=r_j=(w_i+w_j)/2$,感性想一想函数图像就很显然了
中间有一个关于乘积最小生成树优化,然后没用就寄了
上面已经解决了$T=1$,那么$T=2$其实是对第一问扩展罢了
**写到这我已经自闭了**
首先我们求出了最小方差生成树,暴力的按照$T=1$的话,可以考虑在$[l_i,r_i]$外找个最小值,当然会$TLE$,当然正解只是对这个思路优化了一下
我们只需要求$n$条边的答案就好了
删掉了生成树的一条边,$l_i$变化的边一定在右边,$r_i$变化的边一定在左边
那么要求出变化的$l,r$,还是从小到大插入,替换掉的最小边就是$l_i=r_j=(w_i+w_j)/2$,复杂度$O(n^3/n^2logn+nm)$,高精度常数炸掉
继续考虑优化,我们中间有一步直接得到了$m$个区间
每次修改我们改变的边只有$n$个,那么$l_i,r_i$改变的也是$n$个,不知道为什么...
$Ans=\sigma^2\times N^2=s_2\times N-s_1^2$
由于我们范围已经确定了,那么这个东西就搞一个差分就好了,对于更改的位置呢
把原序列变化成了$n$个区间,我们这个东西的变化值有了,就是这一段$s$会发生的变化就知道了,我们就是求着每个区间的最小值就好了
$s_1$是和,$s_2$是平方和
$Ans=(s_{i2}+\det s_2)\times N-(s_{i1}^2+2\times s_{i1}\times\det s_1+\det^2s_1)$
这是个斜率优化式子,线段树维护凸壳然后指针查询就好了
粘一下$60pts$代码,最后那个维护很麻烦,大坑未补.
现在$450$行,写完貌似得破$600$了...
#include<bits/stdc++.h> #define base 1000000000 #define int long long #define MAXM 100005 #define MAXN 305 using namespace std; struct node { int x,y,w; int id; }rd[MAXM]; bool cmp(node a,node b) { return a.w<b.w; } int n,m,T,vis0[MAXM],L0[MAXM],R0[MAXM]; struct Num{ int p,len; int a[5]; Num(int k=0){ p=len=0; memset(a,0,sizeof(a)); if (k<0)len=-1; else{ while (k){ a[len++]=k%base; k/=base; } } } }Check,S1,S2,ans,z1[MAXM],z2[MAXM],dS1[MAXM<<1],dS2[MAXM<<1],Ans[MAXM]; namespace IO{ int num[100]; int x; char c; int read(){ x=0,c=getchar(); while ((c<'0')||(c>'9'))c=getchar(); while ((c>='0')&&(c<='9')){ x=x*10+c-'0'; c=getchar(); } return x; } void write(Num x,char c='\0'){ if (x.len<0){ putchar('-'); putchar('1'); putchar(c); return; } for(int i=0;i<x.len;i++) for(int j=0;j<9;j++){ num[++num[0]]=x.a[i]%10; x.a[i]/=10; } while ((num[0])&&(!num[num[0]]))num[0]--; if (!num[0])putchar('0'); while (num[0])putchar(num[num[0]--]+'0'); putchar(c); } }; vector<pair<long long,int> >v0; namespace LCT { int fa[MAXM<<1],val[MAXM<<1],vis[MAXM],Max[MAXM<<1],Min[MAXM<<1],st[MAXM<<1],ch[MAXM<<1][2]; int tag[MAXM<<1]; bool Nrt(int x) { int f=fa[x]; if(ch[f][0]==x||ch[f][1]==x) return true; return false; } //一个点的所有儿子里面只有一个实边 //splay维护的是一棵深度递增的节点 //显然是左边的时候上边的,右边的是下面的 //那么显然的维护出来的是这个点到最深点的一条链 void push_up(int x) { Min[x]=min(Min[ch[x][0]],Min[ch[x][1]]); Max[x]=max(Max[ch[x][0]],Max[ch[x][1]]); if(x>n) { Min[x]=min(Min[x],x-n); Max[x]=max(Max[x],x-n); } //其实这里只有边才有权值 //那么边权就是大小了 } void push_poz(int x) { swap(ch[x][0],ch[x][1]); tag[x]^=1; } void push_down(int x) { if(tag[x]) { if(ch[x][0]) push_poz(ch[x][0]); if(ch[x][1]) push_poz(ch[x][1]); tag[x]^=1; } } bool get_which(int x) { return ch[fa[x]][1]==x; } void rotate(int x) { int f=fa[x],g_fa=fa[f],k=(ch[f][1]==x),w=ch[x][!k]; if(Nrt(f)) ch[g_fa][ch[g_fa][1]==f]=x; ch[x][!k]=f; ch[f][k]=w; if(w) fa[w]=f; f[fa]=x; fa[x]=g_fa; push_up(f); push_up(x); } void splay(int x) { int y=x,z=0; st[++z]=y; while(Nrt(y)) st[++z]=y=fa[y]; while(z) push_down(st[z--]); while(Nrt(x)) { y=fa[x];z=fa[y]; if(Nrt(y)) { rotate(get_which(x)==get_which(y)?y:x); } rotate(x); } push_up(x); } void access(int x) { for(int y=0;x;x=fa[y=x]) { splay(x); ch[x][1]=y; push_up(x); } } void mkrt(int x) { access(x); splay(x); push_poz(x); } int Fdrt(int x) { access(x); splay(x); while(ch[x][0]) { push_down(x); x=ch[x][0]; } splay(x); return x; } void spilt(int x,int y) { mkrt(x); access(y); splay(y); } void link(int x,int y) { mkrt(x); if(Fdrt(y)!=x) { fa[x]=y; } } void cut(int x,int y) { mkrt(x); if(Fdrt(y)==x&&fa[y]==x&&!ch[y][0]) { fa[y]=ch[x][1]=0; push_up(x); } } int query_min(int x,int y) { mkrt(x); if(Fdrt(y)!=x) return 0; return Min[x]; } int query_max(int x,int y) { mkrt(x); if(Fdrt(y)!=x) return 0; return Max[x]; } int add_min(int id) { int pos=query_min(rd[id].y,id+n); //反正是加边,一开始先加一半的生成树不会变化 //显然不会加半个边,那么现在加另一半就好了 //查另一半,pos就是更新完最小生成树的位置 //如果有的话就更新就好了 //由于维护的是最大生成树 //感觉是最大生成树其实是最小的位置 //反正都需要删掉hhh // cout<<"pos: "<<pos<<endl; if(pos) { vis[0]--,vis[pos]=0; cut(rd[pos].y,pos+n); } vis[0]++,vis[id]=1; link(rd[id].y,id+n); return pos; } int add_max(int id) { int pos=query_max(rd[id].y,id+n); if(pos) { vis[0]--,vis[pos]=0; cut(rd[pos].y,pos+n); } vis[0]++,vis[id]=1; link(rd[id].y,id+n); return pos; } void Init() { Min[0]=0x3f3f3f3f,Max[0]=0; //这就是LCT最小生成树的过程 for(int i=1;i<=n+m;i++) { fa[i]=tag[i]=ch[i][0]=ch[i][1]=0; push_up(i); } for(int i=1;i<=m;i++) { link(rd[i].x,i+n); } } void clear() { for(int i=1;i<=m;i++) { if(vis[i]) { cut(rd[i].y,i+n); vis[0]--,vis[i]=0; } } //至今尚未知晓vis何用 //vis应该是这条边在不在最小生成树里面 } } namespace Clac { int cmp(Num x,Num y){ if (x.len!=y.len){ if (x.len<y.len)return -1; return 1; } for(int i=x.len-1;i>=0;i--) if (x.a[i]!=y.a[i]){ if (x.a[i]<y.a[i])return -1; return 1; } return 0; } Num min(Num x,Num y){ if (cmp(x,y)<0)return x; return y; } Num add(Num x,Num y){ Num ans; if (x.p==y.p){ ans.p=x.p,ans.len=max(x.len,y.len); for(int i=0;i<ans.len;i++){ ans.a[i]+=x.a[i]+y.a[i]; if (ans.a[i]>=base){ ans.a[i]-=base; ans.a[i+1]++; } } if ((ans.len<5)&&(ans.a[ans.len]))ans.len++; return ans; } if (cmp(x,y)<0)swap(x,y); ans.p=x.p,ans.len=x.len; for(int i=0;i<ans.len;i++){ ans.a[i]+=x.a[i]-y.a[i]; if (ans.a[i]<0){ ans.a[i]+=base; ans.a[i+1]--; } } while ((ans.len)&&(!ans.a[ans.len-1]))ans.len--; return ans; } Num dec(Num x,Num y){ y.p^=1; return add(x,y); } Num mul(Num x,Num y){ Num ans; ans.p=(x.p^y.p),ans.len=x.len+y.len-1; for(int i=0;i<x.len;i++) for(int j=0;j<y.len;j++){ ans.a[i+j]+=x.a[i]*y.a[j]; ans.a[i+j+1]+=ans.a[i+j]/base; ans.a[i+j]%=base; } if ((ans.len<5)&&(ans.a[ans.len]))ans.len++; return ans; } };//-------------------封装高精度 void get_LR() { memset(L0,0,sizeof(L0)); memset(R0,0,sizeof(R0)); LCT::clear(); for(int i=1;i<=m;i++) { L0[i]=LCT::add_min(i); // cout<<"L0: "<<L0[i]<<endl; if(L0[i]) R0[L0[i]]=i; } } void get_v() { v0.clear(); for(int i=1;i<=m;i++) { if(!L0[i]) { v0.push_back(make_pair(0,i)); } else { v0.push_back(make_pair((rd[L0[i]].w+rd[i].w<<1)+1,i)); } if(R0[i]) { v0.push_back(make_pair((rd[i].w+rd[R0[i]].w<<1)+1,-i)); } } } void clac() { //果不其然,这里要排序 //这是就是统计贡献了 sort(v0.begin(),v0.end()); int tot=0; S1=S2=0,ans.len=10; for(int i=0;i<v0.size();i++) { int pos=abs(v0[i].second); if(v0[i].second>0) { tot++; S1=Clac::add(S1,z1[pos]),S2=Clac::add(S2,z2[pos]); } else { tot--; S1=Clac::dec(S1,z1[pos]),S2=Clac::dec(S2,z2[pos]); } if(v0[i].first!=v0[i+1].first) { if(!v0[i].first) Check=S1; ans=Clac::min(ans,Clac::dec(S2,Clac::mul(S1,S1))); } } // cout<<"tot: "<<tot<<endl; if(tot!=n-1) { ans=-1; return ; } S1=S2=0; memset(vis0,0,sizeof(vis0)); for(int i=0;i<v0.size();i++) { int pos=abs(v0[i].second); vis0[pos]^=1; if(v0[i].second>0) { S1=Clac::add(S1,z1[pos]),S2=Clac::add(S2,z2[pos]); } else { S1=Clac::dec(S1,z1[pos]),S2=Clac::dec(S2,z2[pos]); } if(v0[i].first!=v0[i+1].first&&(!Clac::cmp(ans,Clac::dec(S2,Clac::mul(S1,S1)))))break; } } signed main() { cin>>n>>m>>T; for(int i=1;i<=m;i++) { cin>>rd[i].x>>rd[i].y>>rd[i].w; rd[i].id=i; } sort(rd+1,rd+1+m,cmp); LCT::Init(); for(int i=1;i<=m;i++) { z1[i]=rd[i].w; z2[i]=Clac::mul(n-1,Clac::mul(z1[i],z1[i])); // IO::write(z1[i],'\n'); // IO::write(z2[i],'\n'); } get_LR(); // cout<<"VIS: "<<LCT::vis[0]<<endl; if(LCT::vis[0]!=n-1) { if(T==1) { puts("-1"); } else { for(int i=1;i<=m;i++) { puts("-1"); } } return 0; } get_v(),clac(); if(T==1) { IO::write(ans,'\n'); return 0; } for(int i=1;i<=m;i++) { if(!vis0[i]) Ans[rd[i].id]=ans; } }
机房dalao教我压行...
#include<bits/stdc++.h> #define base 1000000000 #define int long long #define MAXM 100005 #define MAXN 305 using namespace std; struct node{int x,y,w,id;}rd[MAXM]; bool cmp(node a,node b){return a.w<b.w;} int n,m,T,vis0[MAXM],L0[MAXM],R0[MAXM]; struct Num{int p,len;int a[5];Num(int k=0){p=len=0;memset(a,0,sizeof(a));if (k<0)len=-1;else{while (k){ a[len++]=k%base;k/=base;}}}}Check,S1,S2,ans,z1[MAXM],z2[MAXM],dS1[MAXM<<1],dS2[MAXM<<1],Ans[MAXM]; namespace IO{ int num[100],x; char c; int read(){x=0,c=getchar();while ((c<'0')||(c>'9'))c=getchar();while ((c>='0')&&(c<='9')){x=x*10+c-'0';c=getchar();}return x;} void write(Num x,char c='\0'){if (x.len<0){putchar('-');putchar('1');putchar(c);return;} for(int i=0;i<x.len;i++)for(int j=0;j<9;j++){num[++num[0]]=x.a[i]%10;x.a[i]/=10;} while ((num[0])&&(!num[num[0]]))num[0]--;if (!num[0])putchar('0');while (num[0])putchar(num[num[0]--]+'0');putchar(c);} }; vector<pair<long long,int> >v0; namespace LCT { int fa[MAXM<<1],val[MAXM<<1],vis[MAXM],Max[MAXM<<1],Min[MAXM<<1],st[MAXM<<1],ch[MAXM<<1][2],tag[MAXM<<1]; bool Nrt(int x){int f=fa[x];if(ch[f][0]==x||ch[f][1]==x) return true; return false;} void push_up(int x){Min[x]=min(Min[ch[x][0]],Min[ch[x][1]]);Max[x]=max(Max[ch[x][0]],Max[ch[x][1]]);if(x>n){Min[x]=min(Min[x],x-n);Max[x]=max(Max[x],x-n);}} void push_poz(int x){swap(ch[x][0],ch[x][1]);tag[x]^=1;} void push_down(int x){if(tag[x]){if(ch[x][0]) push_poz(ch[x][0]);if(ch[x][1]) push_poz(ch[x][1]);tag[x]^=1;}} bool get_which(int x){return ch[fa[x]][1]==x;} void rotate(int x){int f=fa[x],g_fa=fa[f],k=(ch[f][1]==x),w=ch[x][!k];if(Nrt(f)) ch[g_fa][ch[g_fa][1]==f]=x;ch[x][!k]=f;ch[f][k]=w;if(w) fa[w]=f;f[fa]=x;fa[x]=g_fa;push_up(f);push_up(x);} void splay(int x){int y=x,z=0;st[++z]=y;while(Nrt(y)) st[++z]=y=fa[y];while(z) push_down(st[z--]);while(Nrt(x)){y=fa[x];z=fa[y];if(Nrt(y)){rotate(get_which(x)==get_which(y)?y:x);}rotate(x);}push_up(x);} void access(int x){for(int y=0;x;x=fa[y=x]){splay(x);ch[x][1]=y;push_up(x);}} void mkrt(int x){access(x);splay(x);push_poz(x);} int Fdrt(int x){access(x);splay(x);while(ch[x][0]){push_down(x);x=ch[x][0];}splay(x);return x;} void spilt(int x,int y){mkrt(x);access(y);splay(y);} void link(int x,int y){mkrt(x);if(Fdrt(y)!=x){fa[x]=y;}} void cut(int x,int y){mkrt(x);if(Fdrt(y)==x&&fa[y]==x&&!ch[y][0]){fa[y]=ch[x][1]=0;push_up(x);}} int query_min(int x,int y){mkrt(x);if(Fdrt(y)!=x) return 0;return Min[x];} int query_max(int x,int y){mkrt(x);if(Fdrt(y)!=x) return 0;return Max[x];} int add_min(int id){int pos=query_min(rd[id].y,id+n);if(pos){vis[0]--,vis[pos]=0;cut(rd[pos].y,pos+n);}vis[0]++,vis[id]=1;link(rd[id].y,id+n);return pos;} int add_max(int id){int pos=query_max(rd[id].y,id+n);if(pos){vis[0]--,vis[pos]=0;cut(rd[pos].y,pos+n);}vis[0]++,vis[id]=1;link(rd[id].y,id+n);return pos;} void Init(){Min[0]=0x3f3f3f3f,Max[0]=0;for(int i=1;i<=n+m;i++){fa[i]=tag[i]=ch[i][0]=ch[i][1]=0;push_up(i);}for(int i=1;i<=m;i++){link(rd[i].x,i+n);}} void clear(){for(int i=1;i<=m;i++){if(vis[i]){cut(rd[i].y,i+n);vis[0]--,vis[i]=0;}}} } namespace Clac{ int cmp(Num x,Num y){if (x.len!=y.len){if (x.len<y.len)return -1;return 1;}for(int i=x.len-1;i>=0;i--)if (x.a[i]!=y.a[i]){if (x.a[i]<y.a[i])return -1;return 1;}return 0;} Num min(Num x,Num y){if (cmp(x,y)<0)return x;return y;} Num add(Num x,Num y){Num ans;if (x.p==y.p){ans.p=x.p,ans.len=max(x.len,y.len);for(int i=0;i<ans.len;i++){ans.a[i]+=x.a[i]+y.a[i];if (ans.a[i]>=base){ans.a[i]-=base;ans.a[i+1]++;}}if ((ans.len<5)&&(ans.a[ans.len]))ans.len++;return ans;}if (cmp(x,y)<0)swap(x,y);ans.p=x.p,ans.len=x.len;for(int i=0;i<ans.len;i++){ans.a[i]+=x.a[i]-y.a[i];if (ans.a[i]<0){ans.a[i]+=base;ans.a[i+1]--;}}while ((ans.len)&&(!ans.a[ans.len-1]))ans.len--;return ans;} Num dec(Num x,Num y){y.p^=1;return add(x,y);} Num mul(Num x,Num y){Num ans;ans.p=(x.p^y.p),ans.len=x.len+y.len-1;for(int i=0;i<x.len;i++)for(int j=0;j<y.len;j++){ans.a[i+j]+=x.a[i]*y.a[j];ans.a[i+j+1]+=ans.a[i+j]/base;ans.a[i+j]%=base;}if ((ans.len<5)&&(ans.a[ans.len]))ans.len++;return ans;} };//-------------------封装高精度 void get_LR(){memset(L0,0,sizeof(L0));memset(R0,0,sizeof(R0));LCT::clear();for(int i=1;i<=m;i++){L0[i]=LCT::add_min(i);if(L0[i]) R0[L0[i]]=i;}} void get_v(){v0.clear();for(int i=1;i<=m;i++){if(!L0[i]){v0.push_back(make_pair(0,i));}else{v0.push_back(make_pair((rd[L0[i]].w+rd[i].w<<1)+1,i));}if(R0[i]){v0.push_back(make_pair((rd[i].w+rd[R0[i]].w<<1)+1,-i));}}} void clac() { sort(v0.begin(),v0.end());int tot=0;S1=S2=0,ans.len=10; for(int i=0;i<v0.size();i++){int pos=abs(v0[i].second);if(v0[i].second>0){tot++;S1=Clac::add(S1,z1[pos]),S2=Clac::add(S2,z2[pos]);}else{tot--;S1=Clac::dec(S1,z1[pos]),S2=Clac::dec(S2,z2[pos]);}if(v0[i].first!=v0[i+1].first){if(!v0[i].first) Check=S1;ans=Clac::min(ans,Clac::dec(S2,Clac::mul(S1,S1)));}} if(tot!=n-1){ans=-1;return ;} S1=S2=0;memset(vis0,0,sizeof(vis0)); for(int i=0;i<v0.size();i++){int pos=abs(v0[i].second);vis0[pos]^=1;if(v0[i].second>0){S1=Clac::add(S1,z1[pos]),S2=Clac::add(S2,z2[pos]);}else{S1=Clac::dec(S1,z1[pos]),S2=Clac::dec(S2,z2[pos]);}if(v0[i].first!=v0[i+1].first&&(!Clac::cmp(ans,Clac::dec(S2,Clac::mul(S1,S1)))))break;} } signed main(){ cin>>n>>m>>T; for(int i=1;i<=m;i++){cin>>rd[i].x>>rd[i].y>>rd[i].w;rd[i].id=i;}sort(rd+1,rd+1+m,cmp);LCT::Init(); for(int i=1;i<=m;i++){z1[i]=rd[i].w;z2[i]=Clac::mul(n-1,Clac::mul(z1[i],z1[i]));}get_LR(); if(LCT::vis[0]!=n-1){if(T==1){puts("-1");}else{for(int i=1;i<=m;i++){puts("-1");}}return 0;} get_v(),clac();if(T==1){IO::write(ans,'\n');return 0;} }
