【原始对偶费用流ver2.0】hdu4744

这场比赛我就不多说什么了。。。几乎整场都在写计算几何，恰好三维计算几何又是我的软肋，后面的题基本都没看。。。

建图就不多说了，挺水的，标算不是费用流，但我听说有人zkw费用流过了，于是我就写个原始对偶费用流，但是居然超了，本机测试60组极限1050+ms==。。。于是进行优化，lyp跟我说过一次反向spfa就可以dij下去了，但是自己推了一下,有一个对原图的边修改的过程.

考虑i->j,w[i][j]>0,操作应该是w[i][j]-=d[i]-d[j],因为d[i]-d[j]<=w[i][j]，所以w[i][j]-(d[i]-d[j])>=0,但是考虑反边w[j][i]-=d[j]-d[i],也就是-w[i][j]+d[i]-d[j],因为d[i]-d[j]<=w[i][j]，所以这个反边貌似还是负的，因此修改过一次后不是所想的全为正权边，从而不能dij，lyp以前看到的那篇论文找不到了，所以也没法继续探究下去，我自己也不是很理解这个算法的理论基础- -，于是另辟蹊径，用了春潇最喜欢的优先队列+spfa，瞬间600+ms过了...

原始对偶费用流对于点度比较大的图貌似确实有些慢...

#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <algorithm>#include <cmath>#include <ctime>#define sqr(x) ((x)*(x))const int oo=1073741819;using namespace std;int tail[250],d[250],v[250],flag[250],x[250],y[250],z[250],w[250],b[1050],p[250];int next[150000],sora[150000],flow[150000],cost[150000],st[150000],po[150000];int n,ss,s,t,phi,ans,sum,w_time,e,m1;void change(int x,int w){d[x]=w;for (x=((x+m1)>>1);x;x>>=1) if (d[b[x<<1]]>d[b[(x<<1)^1]]) b[x]=b[(x<<1)^1];else b[x]=b[x<<1];}bool spfa(int s,int t){int h,r,ne,na;for (int i=0;i<=t;i++) d[i]=p[i]=oo;h=r=0;change(t,0),p[t]=0,v[t]=1;//++r;for (;d[b[1]]!=oo;) {ne=b[1];//++h;change(ne,oo);//v[ne]=0;for (int i=ne;next[i];) {i=next[i],na=sora[i];if (flow[po[i]] && p[ne]+cost[po[i]]<p[na]) {p[na]=p[ne]+cost[po[i]];change(na,p[na]);/*if (!v[na]) {v[na]=1,++r;}*/}}}//for (int i=0;i<=t;i++) d[i]=p[i];if (p[s]>=oo) return 0;phi+=p[s];for (int i=1;i<=t;i++) for (int j=i;next[j];) {j=next[j],ne=sora[j];cost[j]-=p[i]-p[ne];}return 1;}int dfs(int x,int low){if (t==x) {ans+=phi*low;return low;}int sum=0,tmp,ne;flag[x]=w_time;for (int i=x;next[i];) {i=next[i],ne=sora[i];if (flow[i] && flag[ne]!=w_time && !cost[i]) {if (flow[i]<low) tmp=dfs(ne,flow[i]);else tmp=dfs(ne,low);flow[i]-=tmp,flow[po[i]]+=tmp,sum+=tmp,low-=tmp;if (!low) break;}}return sum;}void origin(){s=n+n+1,t=s+1,ss=t+1;//memset(flow,0,sizeof(flow));for (int i=1;i<=t;i++) tail[i]=i,next[i]=0;for (m1=1;m1<t;m1<<=1) ;m1--;for (int i=1;i<=t;i++) b[i+m1]=i;}void link(int x,int y,int z,int c){++ss,next[tail[x]]=ss,tail[x]=ss,sora[ss]=y,flow[ss]=z,cost[ss]=c,next[ss]=0;++ss,next[tail[y]]=ss,tail[y]=ss,sora[ss]=x,flow[ss]=0,cost[ss]=-c,next[ss]=0;po[ss]=ss-1,po[ss-1]=ss;}double dist(int i,int j){return sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])+sqr(z[i]-z[j]));}int main(){freopen("input.txt","r",stdin);freopen("output.txt","w",stdout);w_time=0;for (;scanf("%d",&n)==1;) {if (!n) break;origin();int tot=0;for (int i=1;i<=n;i++) {scanf("%d%d%d%d",&x[i],&y[i],&z[i],&w[i]);link(s,i,w[i],0);link(i+n,t,w[i],0);tot+=w[i];}for (int i=1;i<=n;i++)for (int j=i+1;j<=n;j++) {int tmp=floor(dist(i,j));link(i,j+n,oo,tmp);link(j,i+n,oo,tmp);}ans=0,sum=0;for (phi=0;spfa(s,t);)for (int tmp;w_time++,tmp=dfs(s,oo);) sum+=tmp;if (sum!=tot) printf("-1\n");else printf("%d\n",ans);}//cout<<clock()<<endl;return 0;}

参考了一下主代码手的代码，再一次进行优化，考虑到新的标号D,与原来标号的关系,D[i]=D[j]+c[i][j]-d[i]+d[j],因此其实不需要每次把整个边表的权值变掉，只要把以前的标号记下来就可以了，速度又得到了提升，自己出数据快了300ms,原题快了100ms

update 2014.4.14 还是用stl为好...当手写zkw的时候注意0只能作为废节点，保证值为无穷大，否则死循环...

#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <algorithm>#include <cmath>#include <ctime>#define sqr(x) ((x)*(x))const int oo=1073741819;using namespace std;int tail[250],d[250],v[250],flag[250],x[250],y[250],z[250],w[250],b[1050],p[250],P[250];int next[150000],sora[150000],flow[150000],cost[150000],st[150000],po[150000];int n,ss,s,t,phi,ans,sum,w_time,e,m1;void change(int x,int w){d[x]=w;for (x=((x+m1)>>1);x;x>>=1) if (d[b[x<<1]]>d[b[(x<<1)^1]]) b[x]=b[(x<<1)^1];else b[x]=b[x<<1];}bool spfa(int s,int t){int ne,na;for (int i=0;i<=t;i++) d[i]=p[i]=oo;p[t]=0;change(t,p[t]);for (;d[b[1]]!=oo;) {ne=b[1];change(ne,oo);for (int i=ne;next[i];) {i=next[i],na=sora[i];if (flow[po[i]] && p[ne]+cost[po[i]]+P[ne]<p[na]+P[na]) {p[na]=p[ne]+cost[po[i]]+(P[ne]-P[na]);change(na,p[na]);}}}if (p[s]>=oo) return 0;for (int i=1;i<=t;i++) if (p[i]<oo) P[i]+=p[i];return 1;}int dfs(int x,int low){if (t==x) {ans+=P[s]*low;return low;}int sum=0,tmp,ne;flag[x]=w_time;for (int i=x;next[i];) {i=next[i],ne=sora[i];if (flow[i] && flag[ne]!=w_time && cost[i]+P[ne]==P[x]) {if (flow[i]<low) tmp=dfs(ne,flow[i]);else tmp=dfs(ne,low);flow[i]-=tmp,flow[po[i]]+=tmp,sum+=tmp,low-=tmp;if (!low) break;}}return sum;}void origin(){s=n+n+1,t=s+1,ss=t+1;//memset(flow,0,sizeof(flow));for (int i=1;i<=t;i++) tail[i]=i,next[i]=0;for (m1=1;m1<t;m1<<=1) ;m1--;for (int i=1;i<=t;i++) b[i+m1]=i;}void link(int x,int y,int z,int c){++ss,next[tail[x]]=ss,tail[x]=ss,sora[ss]=y,flow[ss]=z,cost[ss]=c,next[ss]=0;++ss,next[tail[y]]=ss,tail[y]=ss,sora[ss]=x,flow[ss]=0,cost[ss]=-c,next[ss]=0;po[ss]=ss-1,po[ss-1]=ss;}double dist(int i,int j){return sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])+sqr(z[i]-z[j]));}int main(){freopen("input.txt","r",stdin);freopen("output.txt","w",stdout);w_time=0;for (;scanf("%d",&n)==1;) {if (!n) break;origin();int tot=0;for (int i=1;i<=n;i++) {scanf("%d%d%d%d",&x[i],&y[i],&z[i],&w[i]);link(s,i,w[i],0);link(i+n,t,w[i],0);tot+=w[i];}for (int i=1;i<=n;i++)for (int j=i+1;j<=n;j++) {int tmp=floor(dist(i,j));link(i,j+n,oo,tmp);link(j,i+n,oo,tmp);}ans=0,sum=0;for (int i=0;i<=t;i++) P[i]=p[i]=0;for (;spfa(s,t);)for (int tmp;w_time++,tmp=dfs(s,oo);) sum+=tmp;if (sum!=tot) printf("-1\n");else printf("%d\n",ans);}//printf("%d\n",clock()); return 0;}