【原始对偶费用流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;
}