ZJOI 2014 星系调查（推导）

题意

https://loj.ac/problem/2201

思路

说白了就是一条路径上有 (n) 个二维坐标，求一条直线使得所有点到此直线的距离和最小。

设这条直线为 (y=kx+b) ，距离和为 (delta)。

[delta=sum{(kx_i-y_i+b)^2over k^2+1} ]

[delta={k^2sum x_i^2+sum y_i^2+nb^2-2ksum x_iy_i+2kbsum x_i-2bsum y_iover k^2+1} ]

由于 (k,b) 的取值互不影响，我们先假设 (k) 是一个常量，变形如下

[delta={nb^2+(2ksum x_i-2sum y_i)b+k^2sum x_i^2-2ksum x_iy_i+sum y_i^2over k^2+1} ]

发现就是一个关于 (b) 的二次方程，开口朝上，顶点处取最小值，即

[b=-{2ksum x_i-2sum y_iover 2n}={sum y_iover n}-k{sum x_iover n} ]

设 (displaystylear x={sum x_iover n},displaystylear y={sum y_iover n})

得出 (b=ar y-kar x)

代入并化简成关于 (k) 的式子得到

[delta={{(sum x_i^2-2ar xsum x_i+nar x^2)k^2+(-2sum x_iy_i+2ar y sum x_i+2ar xsum y_i-2nar xar y)k+(sum y_i^2-2ar ysum y_i+nar y^2)}over k^2+1} ]

令

[egin{array}{} A&=sum x_i^2-2ar xsum x_i+nar x^2\ &=sum x_i^2-{(sum x_i)^2over n}\ B&=-2sum x_iy_i+2ar y sum x_i+2ar xsum y_i-2nar xar y\ &=-2sum x_iy_i+2{sum x_isum y_iover n}\ C&=sum y_i^2-2ar ysum y_i+nar y^2\ &=sum y_i^2-{(sum y_i)^2over n} end{array} ]

那么

[delta={Ak^2+Bk+Cover k^2+1} ]

把 (k) 当作主元化简得

[(A-delta)k^2+Bk+C-delta=0 ]

这个二次方程的 (Delta) 为 (B^2-4(A-delta)(C-delta))

[B^2-4(A-delta)(C-delta)ge 0\ -4delta ^2+4(A+C)delta-4AC+B^2 ge 0 ]

解得 (displaystyledelta={A+Cpm sqrt{A^2-2AC+B^2+C^2}over 2})

根号前取负号即可达到最小值.

所以维护六元组 ((sum x_i,sum y_i,sum x_i^2,sum y_i^2,sum x_iy_i,sum1)) 即可求出 (A,B,C) ，求出 (delta) 的最小值。

树的情况只用维护到根路径上的信息即可。

基环树的情况类似，断开一条环边，再枚举可行路径（至多两条）即可。

代码

#include<bits/stdc++.h>
#define FOR(i,x,y) for(int i=(x),i##END=(y);i<=i##END;++i)
#define DOR(i,x,y) for(int i=(x),i##END=(y);i>=i##END;--i)
template<typename T,typename _T>inline bool chk_min(T &x,const _T y){return y<x?x=y,1:0;}
template<typename T,typename _T>inline bool chk_max(T &x,const _T y){return x<y?x=y,1:0;}
typedef long long ll;
const int N=(int)1e5+5;
template<const int maxn,const int maxm,typename T>struct Linked_list
{
	int head[maxn],nxt[maxm],tot;T to[maxm];
	Linked_list(){clear();}
	T &operator [](const int x){return to[x];}
	void clear(){memset(head,-1,sizeof(head)),tot=0;}
	void add(int u,T v){to[++tot]=v,nxt[tot]=head[u],head[u]=tot;}
	#define EOR(i,G,u) for(int i=G.head[u];~i;i=G.nxt[i])
};
struct DisjointSet
{
	int fa[N];
	void init(int n){FOR(i,1,n)fa[i]=i;}
	int get_fa(int x){return x==fa[x]?x:fa[x]=get_fa(fa[x]);}
	void merge(int x,int y)
	{
		x=get_fa(x),y=get_fa(y);
		if(x==y)return;
		fa[x]=y;
	}
};
struct hexa
{
	int a,b,c,d,e,f;
	hexa(int _a=0,int _b=0,int _c=0,int _d=0,int _e=0,int _f=0)
	{
		a=_a,b=_b,c=_c,d=_d,e=_e,f=_f;
	}
	hexa operator +(const hexa &_)const
	{
		return hexa(a+_.a,b+_.b,c+_.c,d+_.d,e+_.e,f+_.f);
	}
	hexa operator -(const hexa &_)const
	{
		return hexa(a-_.a,b-_.b,c-_.c,d-_.d,e-_.e,f-_.f);
	}
	hexa add(int x,int y)
	{
		return hexa(a+x,b+y,c+x*x,d+y*y,e+x*y,f+1);
	}
	double solve()
	{
		double A=c-1.0*a*a/f,B=-2*e+2.0*a*b/f,C=d-1.0*b*b/f;
		return (A+C-sqrt(A*A-2*A*C+B*B+C*C))/2;
	}
};
hexa sum[N];
DisjointSet D;
Linked_list<N,N<<1,int>G;
int dep[N],fa[N],sz[N],son[N],top[N];
int X[N],Y[N];
int n,m,q;
int U,V;

void dfs(int u,int f,int d)
{
	dep[u]=d,fa[u]=f,sz[u]=1,son[u]=0;
	sum[u]=sum[f].add(X[u],Y[u]);
	EOR(i,G,u)
	{
		int v=G[i];
		if(v==f)continue;
		dfs(v,u,d+1);
		sz[u]+=sz[v];
		if(sz[v]>sz[son[u]])son[u]=v;
	}
}
void hld(int u,int f,int tp)
{
	top[u]=tp;
	if(son[u])hld(son[u],u,tp);
	EOR(i,G,u)
	{
		int v=G[i];
		if(v==f||v==son[u])continue;
		hld(v,u,v);
	}
}
int get_lca(int x,int y)
{
	while(top[x]!=top[y])
	{
		if(dep[top[x]]<dep[top[y]])std::swap(x,y);
		x=fa[top[x]];
	}
	return dep[x]<dep[y]?x:y;
}
int get_dis(int x,int y)
{
	return dep[x]+dep[y]-2*dep[get_lca(x,y)];
}
hexa get_val(int x,int y)
{
	int lca=get_lca(x,y);
	return (sum[x]-sum[lca]+sum[y]-sum[lca]).add(X[lca],Y[lca]);
}
bool intersect(int a,int b,int c,int d)
{
	int lca1=get_lca(a,b),lca2=get_lca(c,d);
	return 	get_dis(a,lca2)+get_dis(lca2,b)==get_dis(a,b)||
			get_dis(c,lca1)+get_dis(lca1,d)==get_dis(c,d);
}

int main()
{
	scanf("%d%d",&n,&m);
	FOR(i,1,n)scanf("%d%d",&X[i],&Y[i]);
	D.init(n);
	FOR(i,1,m)
	{
		int u,v;
		scanf("%d%d",&u,&v);
		if(D.get_fa(u)==D.get_fa(v))U=u,V=v;
		else
		{
			G.add(u,v),G.add(v,u);
			D.merge(u,v);
		}
	}
	dfs(1,0,1);
	hld(1,0,1);
	scanf("%d",&q);
	while(q--)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		if(n==m-1)printf("%.5lf
",get_val(x,y).solve());
		else
		{
			double res=get_val(x,y).solve();
			if(!intersect(x,U,y,V))chk_min(res,(get_val(x,U)+get_val(y,V)).solve());
			if(!intersect(x,V,y,U))chk_min(res,(get_val(x,V)+get_val(y,U)).solve());
			printf("%.5lf
",res);
		}
	}
	return 0;
}