ZJT's Blog

Keep your determination.

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

2018-1-13 11:552018-1-13 11: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 11:422018-1-13 11: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 11:372018-1-15 17: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 9:22017-11-29 9: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);
}

查看详细内容