题目链接

其实这题没那么毒瘤。。

首先得搞清楚题目里说的浮力势能到底是啥。首先,由于水池是无穷大的,多边形浸入水面并不会使水平面提高。所以,可以认为这个过程把排开的水都提升到了水平面的高度,因此浮力势能可以定义为将多边形排开的水提高到水平面所需的总能量,即$E=\iint_D -ydxdy$,其中$D$为多边形在水平面下的部分。在计算中,只需求出多边形水平面下部分的重心,即可根据重心到水平面的距离计算该积分。

为了保持平衡,重力与浮力必须相等。记$S_0$为整个多边形的面积,$S_1$为多边形在水平面下的面积,我们有$\rho S_0=S_1$,因此$S_1$是确定的,与旋转的角度无关。

类似于旋转卡壳,我们对多边形进行旋转,此时水平面最多与多边形的两条边相交。我们大胆猜测,当与水平面相交的两条边不变时,能量是单峰的。因此我们在每一段中三分角度即可,共有$O(n)$段,因此总时间复杂度为$O(n\log n)$。

实现上,我们有几个小的优化:为了保持稳定,重力与浮力须位于同一直线上,因此多边形整体的重心到水下部分的重心的连线,一定垂直于水平面。因此,我们可以把三分换为二分,每次判断这两条直线的夹角来调整即可。同时,我们可以不二分角度,而是二分其中一个端点在对应线段上的位置。这样在计算中不会涉及到三角函数,可以减小常数,提高精度。

ps:截止到2018.10.27,UOJ上这题所有AC提交的做法都是错误的,均为将角度的范围$[0,2\pi)$等分为四段,在每段中进行三分。随便构造个正五边形再加点随机偏移就卡掉了。

#include <bits/stdc++.h>
#define MAXN 200010
using namespace std;

const double eps=1e-7;
bool equal(double x,double y){ return x-y<eps && y-x<eps; }

struct vec{
    double x,y;
    vec(double _x=0,double _y=0):x(_x),y(_y){}
    vec friend operator+(vec x,vec y){ return vec(x.x+y.x,x.y+y.y); }
    vec friend operator-(vec x,vec y){ return vec(x.x-y.x,x.y-y.y); }
    vec friend operator*(vec x,double k){ return vec(x.x*k,x.y*k); }
    double friend dot(vec x,vec y){ return x.x*y.x+x.y*y.y; }
    double friend cross(vec x,vec y){ return x.x*y.y-x.y*y.x; }
}p[MAXN],cs[MAXN],c_tot;

int n;
double S0,S1,rho;
double ss[MAXN],ts[MAXN];
double ans,ans_ang;

double get_area(vec x,vec y,vec z){ return abs(cross(y-x,z-x))/2; }
vec get_cent(vec x,vec y,vec z){ return (x+y+z)*(1.0/3); }

double get_area(vec a,vec b,vec c,vec d){ return get_area(a,b,c)+get_area(a,c,d); }
vec get_cent(vec a,vec b,vec c,vec d){ return get_cent(a,b,c)*get_area(a,b,c)+get_cent(a,c,d)*get_area(a,c,d); }

double s_area(int l,int r){
    if(l<=r) return ss[r]-ss[l]-get_area(p[1],p[l],p[r]);
    else return ss[n]-s_area(r,l);
}

vec s_cent(int l,int r){
    if(l<=r) return cs[r]-cs[l]-get_cent(p[1],p[l],p[r])*get_area(p[1],p[l],p[r]);
    else return cs[n]-s_cent(r,l);
}

double trans(vec a,vec b,vec c,vec d,double x0,double y0,double x){
    vec p1=a+(b-a)*x0,p2=c+(d-c)*y0,p3=a+(b-a)*x;
    double temp=cross(p3-p1,p2-p1);
    double temp2=cross(d-c,p3-p2);
    if(equal(temp2,0)) return 1e50;
    return temp/temp2+y0;
}

void init(){
    for(int i=3;i<=n;i++){
        ts[i]=get_area(p[i],p[i-1],p[1]);
        ss[i]=ss[i-1]+ts[i];
        cs[i]=cs[i-1]+get_cent(p[i],p[i-1],p[1])*ts[i];
    }
    c_tot=cs[n]*(1.0/ss[n]);
}

double sgn(int l,int r,double x,double y){
    vec p1=p[l]+(p[l+1]-p[l])*x,p2=p[r]+(p[r+1]-p[r])*y;
    double a1=s_area(l%n+1,r),a2=get_area(p1,p[l+1],p[r],p2);
    vec c1=s_cent(l%n+1,r),c2=get_cent(p1,p[l+1],p[r],p2);
    double at=a1+a2;
    //assert(equal(at,S1));
    vec ct=(c1+c2)*(1.0/at);
    return dot(ct-c_tot,p2-p1);
}

double get_ang(vec p1,vec p2){
    double res=M_PI*2-atan2(p2.y-p1.y,p2.x-p1.x);
    if(res<-eps) res+=M_PI*2;
    if(res>M_PI*2-eps) res-=M_PI*2;
    return res;
}

void update_ans(int l,int r,double x,double y){
    vec p1=p[l]+(p[l+1]-p[l])*x,p2=p[r]+(p[r+1]-p[r])*y;
    double a1=s_area(l%n+1,r),a2=get_area(p1,p[l+1],p[r],p2);
    vec c1=s_cent(l%n+1,r),c2=get_cent(p1,p[l+1],p[r],p2);
    double at=a1+a2;
    vec ct=(c1+c2)*(1.0/at);
    double len=sqrt(dot(p1-p2,p1-p2));
    double h1=cross(ct-p1,p2-p1)*(1.0/len),h2=cross(p2-p1,c_tot-p1)*(1.0/len);
    double s=h1*S1+h2*S0*rho;
    if(s<ans){
        ans=s;
        ans_ang=get_ang(p1,p2);
    }
}

void calc(int l,int r,double x0,double x1,double y0){
    double y1=trans(p[l],p[l+1],p[r],p[r+1],x0,y0,x1);
    double sgn0=sgn(l,r,x0,y0),sgn1=sgn(l,r,x1,y1);
    if(equal(sgn0,0)) update_ans(l,r,x0,y0);
    else if(equal(sgn1,0)) update_ans(l,r,x1,y1);
    else if(sgn0*sgn1<0){
        while(x1-x0>1e-8){
            double xmid=(x0+x1)/2,ymid=trans(p[l],p[l+1],p[r],p[r+1],x0,y0,xmid);
            double ts=sgn(l,r,xmid,ymid);
            if(sgn0*ts<0) x1=xmid,y1=ymid,sgn1=ts;
            else x0=xmid,y0=ymid,sgn0=ts;
        }
        update_ans(l,r,x0,y0);
    }
}

void solve(){
    scanf("%d%lf",&n,&rho);
    for(int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
    p[0]=p[n]; p[n+1]=p[1];
    init();
    S0=ss[n]; S1=S0*rho;
    int l=1,r;
    for(r=1;r<=n && ss[r+1]<S1+eps;r++);
    double kl=0,kr=(S1-ss[r])/(ss[r+1]-ss[r]);
    ans=1e50;
    while(1){
        double tl=trans(p[r],p[r+1],p[l],p[l+1],kr,kl,1),tr=trans(p[l],p[l+1],p[r],p[r+1],kl,kr,1);
        if(tl<1+eps){
            calc(l,r,kl,tl,kr);
            r=r%n+1;
            kr=0;
            kl=tl;
        }else{
            //assert(tr<1+eps);
            calc(l,r,kl,1,kr);
            l=l%n+1;
            kl=0;
            kr=tr;
            if(l==1) break;
        }
    }
    printf("%.10f\n",ans_ang);
}

int main(){
#ifdef DEBUG
    freopen("159.in","r",stdin);
#endif
    int T;
    scanf("%d",&T);
    while(T--) solve();
}