「ZJOI2014」星系调查
「ZJOI2014」星系调查
本题核心在于快速求XPs 的线性假设相斥度。
点\((x1,y1)\)到直线\(y=kx+b\)的距离的平方为\(\displaystyle {(kx1+b-y1)^2}\over {k^2+1}\)。
那么 XPs 的相斥度为\(\displaystyle \sum_{i \in 路径上的点} { {(kx_i+b-y_i)^2}\over {k^2+1}}\)。
将式子拆开：\(\displaystyle \sum_{i \in 路径上的点} {{{x_i}^2 k^2+2x_i kb-2x_iy_i k-2y_ik+b^2+{y_i}^2}\over {k^2+1}}\)
可以发现各项的系数很容易得出（路径求和）。
那么问题即求\(\displaystyle {a1 k^2 + b1 kb +c1 b^2 + d1 b + e1 k+ f1}\over {k^2+1}\)的最小值。
把和\(b\)有关的项提出来。
\[{c1 *({b+ {(b1*k+d1) \over 2*c1 }})^2 - {(b1*k+d1)^2 \over 4*c1} + a1 k^2 + e1 k+ f1}\over {k^2+1}
\]
显然\(\displaystyle ({b+ {(b1*k+d1) \over 2*c1 }})^2 =0\)最优。
然后拆开\(\displaystyle {(b1*k+d1)^2 \over 4*c1}\)，得到新的\(a1,e1,f1\)。
在化简一下为\(\displaystyle a1+ {e1 k+f1-a1 \over {k^2+1}}\)，令\(f1=f1-a1\)。
若\(e1=0\)，答案就是\(a1+min(f1,0)\)。
若\(e1!=0\)，\(\displaystyle {e1 k+f1-a1 \over {k^2+1}} = {\displaystyle e1 \over \displaystyle {k^2+1 \over (k + {f1 \over e1} )}}\)。
令\(f1={f1 \over e1}\)。
然后下面那个直接拆掉，\(\displaystyle k - f1 + {(f1*f1+1) \over (k+f1)}\)用均值不等式求一下即可。
当然还要分正负来讨论：\(\displaystyle min({e1 \over (2*\sqrt {1+f1*f1}-2*f1)},{e1 \over (-2*\sqrt {1+f1*f1}-2*f1)})\)。
然后就好了。
如果是树，那么直接路径求和，如果是基环树，那么走环的时候分两类，不走环一类，共三类。
接下来就是基本操作了。
#include<bits/stdc++.h>
#define rep(q,a,b) for(int q=a,q##_end_=b;q<=q##_end_;++q)
#define dep(q,a,b) for(int q=a,q##_end_=b;q>=q##_end_;--q)
#define mem(a,b) memset(a,b,sizeof a )
#define debug(a) cerr<<#a<<' '<<a<<"___"<<endl
using namespace std;
void in(int &r) {
static char c;
r=0;
while(c=getchar(),c<48);
do r=(r<<1)+(r<<3)+(c^48);
while(c=getchar(),c>47);
}
#define double long double
bool cur1;
const int mn=100005;
int head[mn],ne[mn<<1],to[mn<<1],cnt1;
#define link(a,b) link_edge(a,b),link_edge(b,a)
#define link_edge(a,b) to[++cnt1]=b,ne[cnt1]=head[a],head[a]=cnt1
#define travel(x) for(int q(head[x]);q;q=ne[q])
int _x[mn],_y[mn],n,m;
double get_the_val(double a1,double b1,double c1,double d1,double e1,double f1){
//a1 k^2 + b1 kb +c1 b^2 + d1 b + e1 k+ f1
//c1 (b+(b1*k+d1)/(2*c1))^2 - (b1*k+d1)^2 / (4*c1)
a1+=-b1*b1/(4*c1);
f1+=-d1*d1/(4*c1);
e1+=-b1*d1/(2*c1);
b1=0;
c1=0;
d1=0;
f1-=a1;
if(!e1)return a1+min(f1,(double)0);
//(a1 k^2 + e1 k+ f1) /over (k^2+1)
//a1 + {(e1 k + f1)\over (k^2+1)}
//e1 * 1/ ((k^2+1) /(k + f1/e1 ))
f1/=e1;
// k - f1 + (f1*f1+1)/(k+f1)
// 2* sqrt(f1*f1+1) -2*f
return a1+min(e1/(2*sqrt(1+f1*f1)-2*f1),e1/(-2*sqrt(1+f1*f1)-2*f1));
}
int si[mn],fa[mn],H[mn],son[mn],top[mn],high;
int LCA(int a,int b){
while(top[a]!=top[b])H[top[a]]>H[top[b]]?a=fa[top[a]]:b=fa[top[b]];
return H[a]<H[b]?a:b;
}
int sum_mul_x[mn],sum_x_y[mn],sum_y[mn],sum_x[mn],sum_mul_y[mn];
bool mark_in_loop[mn];
int in_which_node[mn],rt;
void dfs(int f,int x){
sum_mul_x[x]=_x[x]*_x[x]+sum_mul_x[f];
sum_mul_y[x]=_y[x]*_y[x]+sum_mul_y[f];
sum_x_y[x]=_x[x]*_y[x]+sum_x_y[f];
sum_y[x]=_y[x]+sum_y[f];
sum_x[x]=_x[x]+sum_x[f];
in_which_node[x]=rt;
H[x]=++high;
si[x]=1,fa[x]=f;
travel(x)if(to[q]!=f&&!mark_in_loop[to[q]]){
dfs(x,to[q]);
si[x]+=si[to[q]];
if(si[to[q]]>si[son[x]])son[x]=to[q];
}
--high;
}
void redfs(int f,int x,int tp){
top[x]=tp;
if(son[x])redfs(x,son[x],tp);
travel(x)if(!mark_in_loop[to[q]]&&to[q]!=f&&to[q]!=son[x])redfs(x,to[q],to[q]);
}
namespace part_1{
void solve(){
dfs(0,1);
redfs(0,1,1);
int Q,a,b;
in(Q);
while(Q--){
in(a),in(b);
int lca=LCA(a,b);
int f_lca=fa[lca];
//a1 k^2 + b1 kb +c1 b^2 + d1 b + e1 k+ f1
//y^2 + b^2 -2 by - 2 yx k +2 x kb +x^2 k^2
double a1=sum_mul_x[a]+sum_mul_x[b]-sum_mul_x[lca]-sum_mul_x[f_lca];
double b1=(sum_x[a]+sum_x[b]-sum_x[lca]-sum_x[f_lca])<<1;
double c1=H[a]+H[b]-H[lca]-H[f_lca];
double d1=-(sum_y[a]+sum_y[b]-sum_y[lca]-sum_y[f_lca])*2;
double e1=-(sum_x_y[a]+sum_x_y[b]-sum_x_y[lca]-sum_x_y[f_lca])*2;
double f1=sum_mul_y[a]+sum_mul_y[b]-sum_mul_y[lca]-sum_mul_y[f_lca];
printf("%.5Lf\n",get_the_val(a1,b1,c1,d1,e1,f1));
}
}
}
namespace part_2{
int loop[mn],loop_len;
int loop_sum_mul_x[mn],loop_sum_x_y[mn],loop_sum_y[mn],loop_sum_x[mn],loop_sum_mul_y[mn];
int last,mark[mn];
void find_loop(int f,int x){
if(mark[x]){
last=x;
return;
}
mark[x]=1;
travel(x)if(to[q]!=f){
find_loop(x,to[q]);
if(last!=-1){
loop[++loop_len]=x;
if(x==last)last=-1;
return;
}
}
}
int LCA(int a,int b){
while(top[a]!=top[b])H[top[a]]>H[top[b]]?a=fa[top[a]]:b=fa[top[b]];
return H[a]<H[b]?a:b;
}
int loop_mp_id[mn];
void solve(){
last=-1,find_loop(0,1);
rep(w,1,loop_len){
mark_in_loop[loop[w]]=1;
loop_mp_id[loop[w]]=w;
}
rep(q,1,n)if(mark_in_loop[q])rt=q,dfs(0,q),redfs(0,q,q);
int v1=0,v2=0,v3=0,v4=0,v5=0;
rep(w,1,loop_len){
int x=loop[w];
v1+=_x[x]*_x[x];
v2+=_x[x]*_y[x];
v3+=_y[x];
v4+=_x[x];
v5+=_y[x]*_y[x];
loop_sum_mul_x[w]=v1;
loop_sum_x_y[w]=v2;
loop_sum_y[w]=v3;
loop_sum_x[w]=v4;
loop_sum_mul_y[w]=v5;
}
int Q,a,b;
int a2,b2,c2,d2,e2,f2,a1,b1,c1,d1,e1,f1;
int l,r;
in(Q);
while(Q--){
in(a),in(b);
if(in_which_node[a]==in_which_node[b]){
int lca=LCA(a,b);
int f_lca=fa[lca];
a1=sum_mul_x[a]+sum_mul_x[b]-sum_mul_x[lca]-sum_mul_x[f_lca];
b1=(sum_x[a]+sum_x[b]-sum_x[lca]-sum_x[f_lca])<<1;
c1=H[a]+H[b]-H[lca]-H[f_lca];
d1=-(sum_y[a]+sum_y[b]-sum_y[lca]-sum_y[f_lca])*2;
e1=-(sum_x_y[a]+sum_x_y[b]-sum_x_y[lca]-sum_x_y[f_lca])*2;
f1=sum_mul_y[a]+sum_mul_y[b]-sum_mul_y[lca]-sum_mul_y[f_lca];
printf("%.5Lf\n",get_the_val(a1,b1,c1,d1,e1,f1));
}else{
double ans=1e18;
l=loop_mp_id[in_which_node[a]],r=loop_mp_id[in_which_node[b]];
if(l>r)swap(l,r),swap(a,b);
++l,--r;
a2=sum_mul_x[a]+sum_mul_x[b];
b2=(sum_x[a]+sum_x[b])<<1;
c2=H[a]+H[b];
d2=-(sum_y[a]+sum_y[b])<<1;
e2=-(sum_x_y[a]+sum_x_y[b])<<1;
f2=sum_mul_y[a]+sum_mul_y[b];
a1=loop_sum_mul_x[r]-loop_sum_mul_x[l-1];
b1=(loop_sum_x[r]-loop_sum_x[l-1])<<1;
c1=r-l+1;
d1=-(loop_sum_y[r]-loop_sum_y[l-1])<<1;
e1=-(loop_sum_x_y[r]-loop_sum_x_y[l-1])<<1;
f1=loop_sum_mul_y[r]-loop_sum_mul_y[l-1];
ans=min(ans,get_the_val(a1+a2,b1+b2,c1+c2,d1+d2,e1+e2,f1+f2));
a1=loop_sum_mul_x[loop_len]-loop_sum_mul_x[r+1]+loop_sum_mul_x[l-2];
b1=(loop_sum_x[loop_len]-loop_sum_x[r+1]+loop_sum_x[l-2])<<1;
c1=loop_len-(r+1)+l-2;
d1=-(loop_sum_y[loop_len]-loop_sum_y[r+1]+loop_sum_y[l-2])<<1;
e1=-(loop_sum_x_y[loop_len]-loop_sum_x_y[r+1]+loop_sum_x_y[l-2])<<1;
f1=loop_sum_mul_y[loop_len]-loop_sum_mul_y[r+1]+loop_sum_mul_y[l-2];
ans=min(ans,get_the_val(a1+a2,b1+b2,c1+c2,d1+d2,e1+e2,f1+f2));
printf("%.5Lf\n",ans);
}
}
}
}
bool cur2;
int main(){
// cerr<<(&cur2-&cur1)/1024.0/1024.0<<endl;
freopen("inv.in","r",stdin);
freopen("inv.out","w",stdout);
int a,b;
in(n),in(m);
rep(q,1,n)in(_x[q]),in(_y[q]);
rep(q,1,m)in(a),in(b),link(a,b);
if(m==n-1)part_1::solve();
else part_2::solve();
return 0;
}
