ZJT's Blog

大変に気分がいい~

[UOJ #159][清华集训 2015]多边形下海

2018-10-28 8:192018-10-28 8:24
计算几何旋转卡壳物理

题目链接

其实这题没那么毒瘤。。

首先得搞清楚题目里说的浮力势能到底是啥。首先,由于水池是无穷大的,多边形浸入水面并不会使水平面提高。所以,可以认为这个过程把排开的水都提升到了水平面的高度,因此浮力势能可以定义为将多边形排开的水提高到水平面所需的总能量,即$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();
}

查看详细内容

[BZOJ 3564] [SHOI2014] 信号增幅仪

2018-1-13 16:552018-1-13 16:55
BZOJ计算几何最小圆覆盖

传送门

也是模板题,不过要把所有点在一个方向上进行缩放,正交分解一下就行了。

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

const double eps=1e-6;

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

int n;
double rot,scale;

vec pre_gao(vec x){
    vec axis(cos(rot),sin(rot));
    vec x1=axis*dot(axis,x);
    vec x2=x-x1;
    return x1*(1./scale)+x2;
}

vec get_inter(vec s1,vec d1,vec s2,vec d2){
    double det=cross(d1,d2);
    vec t=vec(d2.y,-d2.x)*(1./det);
    double x=-dot(s1-s2,t);
    return s1+d1*x;
}

vec get_center(vec a,vec b,vec c,double &r){
    vec d=b+(c-b)*.5,e=a+(c-a)*.5;
    vec d1((c-b).y,-(c-b).x),e1((c-a).y,-(c-a).x);
    vec res=get_inter(d,d1,e,e1);
    r=max(max(dis(res,a),dis(res,b)),dis(res,c));
    return res;
}

vec gao3(int x,int y,double &r){
    vec c=a[x]+(a[y]-a[x])*.5;
    r=dis(a[x],a[y])/2;
    for(int i=1;i<x;i++){
        if(dis(a[i],c)<r+eps) continue;
        c=get_center(a[i],a[x],a[y],r);
    }
    return c;
}

vec gao2(int x,double &r){
    vec c=a[x];
    r=0;
    for(int i=1;i<x;i++){
        if(dis(a[i],c)<r+eps) continue;
        c=gao3(i,x,r);
    }
    return c;
}

vec gao1(double &r){
    vec c=a[1];
    r=0;
    for(int i=2;i<=n;i++){
        if(dis(a[i],c)<r+eps) continue;
        c=gao2(i,r);
    }
    return c;
}

int main(){
#ifdef DEBUG
    freopen("in","r",stdin);
#endif
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf%lf",&a[i].x,&a[i].y);
    scanf("%lf%lf",&rot,&scale);
    rot=rot*M_PI/180.;
    for(int i=1;i<=n;i++) a[i]=pre_gao(a[i]);
    random_shuffle(a+1,a+n+1);
    double res;
    gao1(res);
    printf("%.3lf\n",res);
}

查看详细内容

[BZOJ 1336] [Balkan2002] Alien最小圆覆盖

2018-1-13 16:422018-1-13 16:42
BZOJ计算几何最小圆覆盖

传送门

愉快的模板题

1337也是

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

const double eps=1e-6;

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

int n;

vec get_inter(vec s1,vec d1,vec s2,vec d2){
    double det=cross(d1,d2);
    vec t(d2.y/det,-d2.x/det);
    double x=dot(s2-s1,t);
    return s1+d1*x;
}

vec get_center(vec a,vec b,vec c,double &r){
    vec d=a+(b-a)*.5,d1((b-a).y,-(b-a).x);
    vec e=b+(c-b)*.5,e1((c-b).y,-(c-b).x);
    vec res=get_inter(d,d1,e,e1);
    r=dis(res,a);
    return res;
}

vec gao3(int x,int y,double &r){
    vec c=a[x]+(a[y]-a[x])*.5;
    r=dis(a[x],a[y])/2;
    for(int i=1;i<x;i++){
        if(dis(c,a[i])<r+eps) continue;
        c=get_center(a[i],a[x],a[y],r);
    }
    return c;
}

vec gao2(int x,double &r){
    vec c=a[x]; r=0;
    for(int i=1;i<x;i++){
        if(dis(c,a[i])<r+eps) continue;
        c=gao3(i,x,r);
    }
    return c;
}

vec gao1(double &r){
    vec c=a[1]; r=0;
    for(int i=2;i<=n;i++){
        if(dis(c,a[i])<r+eps) continue;
        c=gao2(i,r);
    }
    return c;
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf%lf",&a[i].x,&a[i].y);
    random_shuffle(a+1,a+n+1);
    double r;
    vec c=gao1(r);
    printf("%.10lf\n%.10lf %.10lf\n",r,c.x,c.y);
}

查看详细内容

最小圆覆盖

2018-1-13 16:372018-1-15 22:58
计算几何最小圆覆盖

突然发现自己连这种经典模型都没做过,果然我还是太弱了qwq

最小圆覆盖就是用一个圆去覆盖平面上的一个点集,要求这个圆尽可能小。这个问题可以用随机增量法解决。

这个算法的流程是这样的:先把所有点random_shuffle一下,然后从前往后一个个加点。假设我们已经求出了$1\sim i-1$号点的最小圆覆盖(记为圆$C_{i-1}$),现在要求$C_i$。如果$i$号点在这个圆内(或圆上),则$C_i=C_{i-1}$。否则,我们可以确定$i$号点一定在$C_i$的边界上,且$C_i$的边界上最多只有$3$个点。如果出现了这种情况,我们就求一个子问题:覆盖了前$i-1$个点,且点$i$在边界上的最小的圆。

这个子问题的求法也类似,还是一个个加点。设当前加的点为$j$号点,我们已经得到了包含$1\sim j-1$且边界经过$i$号点的最小圆$C'_{j-1}$,现在要求$C'_j$。如果$j$号点在$C'_{j-1}$内,就可以跟之前一样令$C'_j=C'_{j-1}$。否则,我们可以确定$j$号点和$i$号点都在$C'_j$的边界上,然后又转化成了另一个子问题:覆盖了前$j-1$个点,且点$i$和$j$都在边界上的最小的圆。

这个子问题的求法也类似,还是一个个加点(怎么感觉这句话见过一次)。设当前加的点是$k$号点,我们已经得到了包含$1\sim k-1$且边界经过$j,i$号点的最小圆$C''_{k-1}$,现在要求$C''_k$。如果$k$号点在$C''_{k-1}$内,就可以跟之前一样令$C''_k=C''_{k-1}$。否则,我们可以确定$i,j,k$三个点都在$C''_k$的边界上。但是这时候就不用转化子问题了!因为三个点已经能确定一个圆了。。。

我们来分析一下这个算法的复杂度。记$T_1(n),T_2(n),T_3(n)$分别为三个子问题的复杂度。首先显然$T_3(n)=O(n)$,因为最多会求$n$次三点确定一个圆。由于点的顺序是随机的,且$C'_n$的边界上只会有$O(1)$个点,所以$n$号点在$C'_n$边界上的概率是$O(\frac 1n)$的。由此可得:

$$ \begin{aligned} T_2(n)&=T_2(n-1)+O(\frac 1n)T_3(n)\\&=T_2(n-1)+O(\frac 1n)O(n)\\&=T_2(n-1)+O(1)\\&=O(n) \end{aligned} $$

同理$T_1(n)=T_1(n-1)+O(\frac 1n)T_2(n)=O(n)$。所以整个算法的期望时间复杂度是线性的。

关于怎么三点确定一个圆,可以随便找两条边,然后算中垂线交点就行了(其实就是三角形的外心)。

具体实现的话可以看这道例题

查看详细内容

集训队作业之ARC082-E

2017-11-29 14:22017-11-29 14:2
集训队作业AtCoder计算几何

传送门

显然题目问的那坨东西能转化成有多少点集的子集能圈出一个非退化的凸包(即面积大于$0$),所以把多点共线的情况减掉就行了。枚举直线上的第一个点和第二个点,统计出有多少点在这条线上,然后随便乘一乘加一加减一减就行了。注意不要把同一条直线算重复了。

#include <cstdio>
#include <cstring>
#define MAXN 210
#define LL long long

const LL P=998244353;

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

struct vec2{
    int x,y;
    vec2(int _x=0,int _y=0):x(_x),y(_y){}
    friend vec2 operator-(vec2 x,vec2 y){ return vec2(x.x-y.x,x.y-y.y); }
    friend int dot(vec2 x,vec2 y){ return x.x*y.x+x.y*y.y; }
    friend int cross(vec2 x,vec2 y){ return x.x*y.y-x.y*y.x; }
}p[MAXN];

int n;

int main(){
#ifdef DEBUG
    freopen("E.in","r",stdin);
#endif
    LL inv2=getPow(2,P-2);
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%d%d",&p[i].x,&p[i].y);
    LL ans=getPow(2,n);
    LL res=0;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if(i^j){
                vec2 d=p[j]-p[i];
                int cnt=2;
                for(int k=1;k<=n;k++)
                    if(cross(p[k]-p[j],d)==0 && k!=i && k!=j){
                        if(dot(p[k]-p[j],d)<0){
                            cnt=0;
                            break;
                        }
                        cnt++;
                    }
                LL temp=getPow(2,cnt)-cnt-1;
                res=(res+temp)%P;
            }
    res=res*inv2%P;
    res=(res+n+1)%P;
    ans-=res;
    ans=(ans%P+P)%P;
    printf("%lld\n",ans);
}

查看详细内容