ZigZagK的博客
[随机化+线段树维护凸包+最小圆覆盖]2022“杭电杯”中国大学生算法设计超级联赛(3)1006【Dusk Moon】题解
2022年7月27日 14:18
HDU
查看标签

题目概述

HDU7167

解题报告

被队内计算几何选手演了😡,本来写个最小圆覆盖就过了。

首先很显然 $n$ 个点的最小圆覆盖和这个 $n$ 个点求完凸包之后的最小圆覆盖是一样的。而且数据保证随机,那么凸包只有 $O(\log_2n)$ 个点,我们可以用线段树暴力重构维护区间凸包。

然后每次查询出区间凸包之后做最小圆覆盖就行了,最小圆覆盖在随机数据下复杂度是 $O(size)$ 的,$size$ 是点的个数。

如果不幸被卡常,并且你用的求凸包算法为Andrew,可以把维护整个凸包拆成维护上下凸包,这样每次重构时不需要sort,可以利用归并进行排序,实测能快一倍。

示例程序

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<cstdlib>
#include<algorithm>
#define X first
#define Y second
#define mp make_pair
using namespace std;
typedef long long LL;typedef double DB;
typedef pair<int,int> PI;typedef pair<DB,DB> PD;
typedef vector<PI> Polygon;
const int maxn=100000,maxt=maxn<<2;

int te,n,m;
PI a[maxn+5];
Polygon pg[maxt+5][2],res[2];
int top;PI p[maxn+5],stk[maxn+5];
PD P[maxn+5];

#define EOLN(x) ((x)==10 || (x)==13 || (x)==EOF)
inline char readc(){
    static char buf[100000],*l=buf,*r=buf;
    return l==r && (r=(l=buf)+fread(buf,1,100000,stdin),l==r)?EOF:*l++;
}
template<typename T> int readi(T &x){
    T tot=0;char ch=readc(),lst='+';
    while (!isdigit(ch)) {if (ch==EOF) return EOF;lst=ch;ch=readc();}
    while (isdigit(ch)) tot=(tot<<3)+(tot<<1)+(ch^48),ch=readc();
    lst=='-'?x=-tot:x=tot;return EOLN(ch);
}
template<typename T> inline T operator - (const T &a,const T &b) {return mp(a.X-b.X,a.Y-b.Y);}
LL Cross(const PI &a,const PI &b) {return (LL)a.X*b.Y-(LL)a.Y*b.X;}
void MergeUp(Polygon &a,const Polygon &b){
    if (a.empty()) {a=b;return;}
    int n=0,A=a.size(),B=b.size();PI now;
    for (int i=0,j=0;i<A || j<B;){
        if (j==B || i<A && a[i]<b[j]) now=a[i++]; else now=b[j++];
        if (!n || now!=p[n]) p[++n]=now;
    }
    top=0;
    for (int i=1;i<=n;i++){
        while (top>1 && Cross(p[i]-stk[top-1],stk[top]-stk[top-1])<=0) top--;
        stk[++top]=p[i];
    }
    a.clear();
    for (int i=1;i<=top;i++) a.push_back(stk[i]);
}
void MergeDown(Polygon &a,const Polygon &b){
    if (a.empty()) {a=b;return;}
    int n=0,A=a.size(),B=b.size();PI now;
    for (int i=0,j=0;i<A || j<B;){
        if (j==B || i<A && a[i]>b[j]) now=a[i++]; else now=b[j++];
        if (!n || now!=p[n]) p[++n]=now;
    }
    top=0;
    for (int i=1;i<=n;i++){
        while (top>1 && Cross(p[i]-stk[top-1],stk[top]-stk[top-1])<=0) top--;
        stk[++top]=p[i];
    }
    a.clear();
    for (int i=1;i<=top;i++) a.push_back(stk[i]);
}
void Build(int L,int R,int p=1){
    if (L==R) {pg[p][0].resize(1);pg[p][1].resize(1);pg[p][0][0]=pg[p][1][0]=a[L];return;}
    int mid=L+(R-L>>1);
    Build(L,mid,p<<1);Build(mid+1,R,p<<1|1);
    pg[p][0]=pg[p<<1][0];MergeUp(pg[p][0],pg[p<<1|1][0]);
    pg[p][1]=pg[p<<1][1];MergeDown(pg[p][1],pg[p<<1|1][1]);
//    printf("[%d,%d]\n",L,R);
//    for (auto x:pg[p][0]) printf("%d %d\n",x.X,x.Y);puts("");
//    for (auto x:pg[p][1]) printf("%d %d\n",x.X,x.Y);puts("");
//    puts("");
}
void Update(int pos,int l=1,int r=n,int p=1){
    if (l==r) {pg[p][0][0]=pg[p][1][0]=a[l];return;}
    int mid=l+(r-l>>1);
    pos<=mid?Update(pos,l,mid,p<<1):Update(pos,mid+1,r,p<<1|1);
    pg[p][0]=pg[p<<1][0];MergeUp(pg[p][0],pg[p<<1|1][0]);
    pg[p][1]=pg[p<<1][1];MergeDown(pg[p][1],pg[p<<1|1][1]);
}
void Ask(int L,int R,int l=1,int r=n,int p=1){
    if (L==l && r==R) return MergeUp(res[0],pg[p][0]),MergeDown(res[1],pg[p][1]);
    int mid=l+(r-l>>1);
    if (R<=mid) Ask(L,R,l,mid,p<<1); else if (L>mid) Ask(L,R,mid+1,r,p<<1|1);
    else Ask(L,mid,l,mid,p<<1),Ask(mid+1,R,mid+1,r,p<<1|1);
}
inline DB sqr(DB x) {return x*x;}
DB Norm(PD a) {return sqr(a.X)+sqr(a.Y);}
DB Dis(PD a,PD b) {return sqrt(Norm(a-b));}
LL R(PI a,PI b) {return ceil(Dis(a,b)/2)+0.5;}
PD Center(PD x,PD y,PD z){
    DB a1=y.X-x.X,b1=y.Y-x.Y,a2=z.X-x.X,b2=z.Y-x.Y;
    DB c1=(sqr(a1)+sqr(b1))/2,c2=(sqr(a2)+sqr(b2))/2,d=a1*b2-a2*b1;
    return mp(x.X+(c1*b2-c2*b1)/d,x.Y+(a1*c2-a2*c1)/d);
}
DB Solve(PD *a,int n){
    random_shuffle(a,a+n);
    PD O=a[0];DB R=0;
    for (int i=1;i<n;i++)
        if (Dis(a[i],O)>R){
            O=a[i];R=0;
            for (int j=0;j<i;j++)
                if (Dis(a[j],O)>R){
                    O=mp((a[i].X+a[j].X)/2,(a[i].Y+a[j].Y)/2);R=Dis(a[i],O);
                    for (int k=0;k<j;k++) if (Dis(a[k],O)>R) O=Center(a[i],a[j],a[k]),R=Dis(a[i],O);
                }
        }
    return R;
}
int main(){
    for (readi(te);te;te--){
        readi(n);readi(m);
        for (int i=1;i<=n;i++) readi(a[i].X),readi(a[i].Y);
        Build(1,n);
        for (int t=1,tp,x,y,z;t<=m;t++){
            readi(tp);readi(x);readi(y);
            if (tp==1) readi(z),a[x]=mp(y,z),Update(x); else {
                res[0].clear();res[1].clear();Ask(x,y);
                int si=0;PI fst,lst;
                for (auto x:res[0]) P[si++]=x;
                fst=P[0];lst=P[si-1];
                for (auto x:res[1]) if (x!=fst && x!=lst) P[si++]=x;
                if (si==1) {puts("0");continue;}
                if (si==2) {printf("%lld\n",R(P[0],P[1]));continue;}
                printf("%lld\n",LL(ceil(Solve(P,si))+0.5));
            }
        }
    }
    return 0;
}
版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!