问题可以转化为:给定一个rate,z(rate)=∑xi×ci-rate*∑xi×disi,xi为一棵生成树使(∑xi×ci-rate*∑xi×disi)的值最小(下面会介绍求此生成树的方
法),则rate=(∑xi×ci-z(rate))/( ∑xi×disi),令rateNex=(∑xi×ci)/( ∑xi×disi)。
求解生成树xi使(∑xi×ci-rate*∑xi×disi)的值最小的方法:
∑xi×ci-rate*∑xi×disi=∑xi(ci+rate×disi),从而问题转化为以ci+rate×disi为边的权重,求解最小生成树,对于完全图,可使用
prim算法,其复杂度只与节点数有关。
代码:
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#include<queue>
#include<vector>
#include<cstring>
#include<algorithm>
#define q(a) ((a)*(a))
#define rep(i,a,b) for(int i=(a);i<(b);i++)
#define rev(i,a,b) for(int i=(a);i>=(b);i--)
#define clr(a,x) memset(a,x,sizeof a)
#define inf 0x3f3f3f3f
using namespace std;
const double eps=0.00001;
const int maxn=1005;
const int maxm=maxn*maxn/2;
double dis[maxn],cost[maxn][maxn],len[maxn][maxn],w[maxn][maxn];
int vis[maxn],pre[maxn];
double x[maxn],y[maxn],z[maxn];
int n,m;
bool pan(double mid)
{
for(int i=1;i<n;i++)
for(int j=i+1;j<=n;j++)
w[i][j]=w[j][i]=cost[i][j]-mid*len[i][j];
clr(vis,0);
double mst=0;
for(int i=1;i<=n;i++)
dis[i]=w[1][i],pre[i]=1;
vis[1]=1;
for(int i=1;i<n;i++)
{
int u=-1;
for(int j=1;j<=n;j++)
if(!vis[j]&&(u==-1||dis[j]<dis[u]))u=j;
mst+=w[pre[u]][u];
vis[u]=1;
for(int j=1;j<=n;j++)
if(!vis[j]&&w[u][j]<dis[j])
dis[j]=w[u][j],pre[j]=u;
}
if(mst>0)return 1;
return 0;
}
inline double getdis(int a,int b)
{
return sqrt(q(x[a]-x[b])+q(y[a]-y[b]));
}
int main()
{
while(~scanf("%d",&n)&&n)
{
for(int i=1;i<=n;i++)
scanf("%lf%lf%lf",&x[i],&y[i],&z[i]);
for(int i=1;i<=n;i++)
{
cost[i][i]=len[i][i]=0;
for(int j=i+1;j<=n;j++)
{
cost[i][j]=cost[j][i]=fabs(z[i]-z[j]);
len[i][j]=len[j][i]=getdis(i,j);
}
}
double l=0,r=100,mid;
while(fabs(l-r)>eps)
{
mid=(l+r)/2;
if(pan(mid))l=mid;
else r=mid;
}
printf("%.3f\n",mid);
}
return 0;
}
参考: