(bzoj1337 || 洛谷P1742 最小圆覆盖 )|| (bzoj2823 || 洛谷P2533 [AHOI2012]信号塔)
用随机增量法。讲解:https://blog.csdn.net/jokerwyt/article/details/79221345
设点集A的最小覆盖圆为g(A)
可以发现:
1.g(A)是唯一的
2.g(A)可以由<=3个点唯一确定
①由一个点确定(所有点重合时)
②由两个点确定(圆直径必定为两点连线)
③由三个点确定(圆必定是过三点的唯一圆)
把A中在g(A)上的点叫做关键点
可以发现:若点b不在g(A)内(b不属于A),则b是g(A+b)的关键点
证明:如果b不是g(A+b)的关键点,即b在g(A+b)内,那么显然g(A+b)=g(A),b在g(A)内,与题意矛盾
记点l,点l+1,点l+2,..,点r构成的集合为l..r
现在题目就是要求g(1..i)
设g(A,B)表示点集A的最小覆盖圆,并且点集B中所有点都在这个圆上
设Cir(x,y)表示一个圆心在点x,半径为y的圆
设Cir(a,b,c)表示由a,b,c三点确定的圆
g(1..i)
=g(1..i-1)(如果i属于g(1..i-1))
=g(1..i-1,i)(如果i不属于g(1..i-1))
g(空集,i)=Cir(i,0)
g(1..p,i)
=g(1..p-1,i)(如果p属于g(1..p-1,i))
=g(1..p-1,{i,p})(如果p不属于g(1..p-1,i))
g(空集,{i,j})=Cir((i,j中点),dis(i,j)/2)
g(1..q,{i,p})
=g(1..q-1,{i,p})(如果q属于g(1..q,{i,p}))
=g(1..q-1,{i,p,q})(如果q不属于g(1..q,{i,p}))
g(任意集合,{i,j,k})=Cir(i,j,k)
可以把以上过程改写成循环(以上相当于给出了从g(1..i-1)转移到g(1..i)的方法)
(以下代码交洛谷)
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<cmath> 6 #include<ctime> 7 using namespace std; 8 #define fi first 9 #define se second 10 #define mp make_pair 11 #define pb push_back 12 typedef long long ll; 13 typedef unsigned long long ull; 14 typedef pair<int,int> pii; 15 namespace X 16 { 17 const double eps=1e-10; 18 struct Point 19 { 20 double x,y; 21 Point(double x=0,double y=0):x(x),y(y){} 22 }; 23 typedef Point Vec; 24 Vec operator+(const Vec& a,const Vec& b) 25 { 26 return Vec(a.x+b.x,a.y+b.y); 27 } 28 Vec operator-(const Vec& a,const Vec& b) 29 { 30 return Vec(a.x-b.x,a.y-b.y); 31 } 32 Vec operator*(double a,const Vec& b) 33 { 34 return Vec(a*b.x,a*b.y); 35 } 36 Vec operator*(const Vec& a,double b) 37 { 38 return Vec(b*a.x,b*a.y); 39 } 40 Vec operator/(const Vec& a,double b) 41 { 42 return Vec(a.x/b,a.y/b); 43 } 44 int dcmp(double x) 45 //正为1,负为-1,0为0 46 { 47 if(fabs(x)<eps) return 0; 48 return x<0?-1:1; 49 } 50 bool operator==(const Vec& a,const Vec& b) 51 //判向量相等 52 { 53 return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; 54 } 55 double dot(const Vec& a,const Vec& b) 56 //点积 57 { 58 return a.x*b.x+a.y*b.y; 59 } 60 double cross(const Vec& a,const Vec& b) 61 //叉积 62 { 63 return a.x*b.y-a.y*b.x; 64 } 65 double len(const Vec& x) 66 //向量长度 67 { 68 return sqrt(dot(x,x)); 69 } 70 double angle(const Vec& a,const Vec& b) 71 //夹角,0~180° 72 { 73 return acos(dot(a,b)/len(a)/len(b)); 74 } 75 double angle1(const Vec& a,const Vec& b) 76 //夹角,带方向,a到b逆时针为正,顺时针为负,共线为0 77 { 78 return acos(dot(a,b)/len(a)/len(b))*(dcmp(cross(a,b))); 79 } 80 Vec rotate(const Vec& a,double rad) 81 //旋转,正为逆时针,负为顺时针 82 { 83 return Vec(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); 84 } 85 Vec normal(const Vec& x) 86 //单位法线,左转90°后除以自身长度 87 { 88 double l=len(x); 89 return Vec(-x.y/l,x.x/l); 90 } 91 Vec normal1(const Vec& x) 92 //法线,不归一化 93 { 94 return Vec(-x.y,x.x); 95 } 96 Point lineline(const Point& p,const Vec& v,const Point& q,const Vec& w) 97 //直线交点,GetLineIntersection 98 { 99 return p+v*cross(w,p-q)/cross(v,w); 100 } 101 double ptline(const Point& p,const Point& a,const Point& b) 102 //point_to_line,点到直线距离,DistanceToLine 103 { 104 Vec v1=b-a,v2=p-a; 105 return fabs(cross(v1,v2)/len(v1)); 106 } 107 double ptseg(const Point& p,const Point& a,const Point& b) 108 //point_to_segment,点到线段距离,DistanceToSegment 109 { 110 if(a==b) return len(p-a); 111 //Vec v1=b-a,v2=a-p,v3=p-b; 112 Vec v1=b-a,v2=p-a,v3=p-b; 113 if(dcmp(dot(v1,v2))<0) return len(v2); 114 else if(dcmp(dot(v1,v3))>0) return len(v3); 115 else return fabs(cross(v1,v2)/len(v1)); 116 } 117 double area2(const Point& a,const Point& b,const Point& c) 118 //叉积对应平行四边形的面积 119 { 120 return cross(b-a,c-a); 121 } 122 double thrarea(const Point& a,const Point& b,const Point& c) 123 //三角形面积,绝对值 124 { 125 return fabs(cross(b-a,c-a)/2); 126 } 127 bool ifpar(const Vec& a,const Vec& b) 128 //ifParallel 129 //是否共线/平行 130 { 131 return dcmp(cross(a,b))==0; 132 } 133 Point pointline(const Point& p,const Point& a,const Vec& v) 134 //点在直线上投影,GetLineProjection 135 { 136 return a+v*(dot(p-a,v)/dot(v,v)); 137 } 138 double sqr(double x){return x*x;} 139 double dis(const Point &x,const Point &y) 140 { 141 return sqrt(sqr(x.x-y.x)+sqr(x.y-y.y)); 142 } 143 }; 144 using namespace X; 145 Point p[100010]; 146 int n; 147 Point a; 148 double r; 149 void make(const Point &x,const Point &y) 150 { 151 Vec t=(y-x)/2; 152 a=x+t; 153 r=len(t); 154 } 155 void make(const Point &x,const Point &y,const Point &z) 156 { 157 Vec t1=(y-x)/2,t2=(z-y)/2; 158 Point x1=x+t1,y1=y+t2; 159 t1=normal1(t1);t2=normal1(t2); 160 a=lineline(x1,t1,y1,t2); 161 r=dis(a,x); 162 } 163 bool ok(const Point &x) 164 { 165 return dis(x,a)<=r; 166 } 167 int main() 168 { 169 int i,j,k; 170 srand(time(0)); 171 scanf("%d",&n); 172 for(i=1;i<=n;i++) 173 scanf("%lf%lf",&p[i].x,&p[i].y); 174 random_shuffle(p+1,p+n+1); 175 a=p[1];r=0; 176 for(i=2;i<=n;i++) 177 if(!ok(p[i])) 178 { 179 a=p[i];r=0; 180 for(j=1;j<i;j++) 181 if(!ok(p[j])) 182 { 183 make(p[i],p[j]); 184 for(k=1;k<j;k++) 185 if(!ok(p[k])) 186 { 187 make(p[i],p[j],p[k]); 188 } 189 } 190 } 191 printf("%.10f\n",r); 192 printf("%.10f %.10f\n",a.x,a.y); 193 return 0; 194 }
此题同上面题,改一下输出格式就一样
bzoj数据稍微有点卡精度,把存储半径改成存储半径的平方就可以A了
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<cmath> 6 #include<ctime> 7 using namespace std; 8 #define fi first 9 #define se second 10 #define mp make_pair 11 #define pb push_back 12 typedef long long ll; 13 typedef unsigned long long ull; 14 typedef pair<int,int> pii; 15 namespace X 16 { 17 const double eps=1e-6; 18 struct Point 19 { 20 double x,y; 21 Point(double x=0,double y=0):x(x),y(y){} 22 }; 23 typedef Point Vec; 24 Vec operator+(const Vec& a,const Vec& b) 25 { 26 return Vec(a.x+b.x,a.y+b.y); 27 } 28 Vec operator-(const Vec& a,const Vec& b) 29 { 30 return Vec(a.x-b.x,a.y-b.y); 31 } 32 Vec operator*(double a,const Vec& b) 33 { 34 return Vec(a*b.x,a*b.y); 35 } 36 Vec operator*(const Vec& a,double b) 37 { 38 return Vec(b*a.x,b*a.y); 39 } 40 Vec operator/(const Vec& a,double b) 41 { 42 return Vec(a.x/b,a.y/b); 43 } 44 int dcmp(double x) 45 //正为1,负为-1,0为0 46 { 47 if(fabs(x)<eps) return 0; 48 return x<0?-1:1; 49 } 50 bool operator==(const Vec& a,const Vec& b) 51 //判向量相等 52 { 53 return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; 54 } 55 double dot(const Vec& a,const Vec& b) 56 //点积 57 { 58 return a.x*b.x+a.y*b.y; 59 } 60 double cross(const Vec& a,const Vec& b) 61 //叉积 62 { 63 return a.x*b.y-a.y*b.x; 64 } 65 double len(const Vec& x) 66 //向量长度 67 { 68 return sqrt(dot(x,x)); 69 } 70 double angle(const Vec& a,const Vec& b) 71 //夹角,0~180° 72 { 73 return acos(dot(a,b)/len(a)/len(b)); 74 } 75 double angle1(const Vec& a,const Vec& b) 76 //夹角,带方向,a到b逆时针为正,顺时针为负,共线为0 77 { 78 return acos(dot(a,b)/len(a)/len(b))*(dcmp(cross(a,b))); 79 } 80 Vec rotate(const Vec& a,double rad) 81 //旋转,正为逆时针,负为顺时针 82 { 83 return Vec(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); 84 } 85 Vec normal(const Vec& x) 86 //单位法线,左转90°后除以自身长度 87 { 88 double l=len(x); 89 return Vec(-x.y/l,x.x/l); 90 } 91 Vec normal1(const Vec& x) 92 //法线,不归一化 93 { 94 return Vec(-x.y,x.x); 95 } 96 Point lineline(const Point& p,const Vec& v,const Point& q,const Vec& w) 97 //直线交点,GetLineIntersection 98 { 99 return p+v*cross(w,p-q)/cross(v,w); 100 } 101 double ptline(const Point& p,const Point& a,const Point& b) 102 //point_to_line,点到直线距离,DistanceToLine 103 { 104 Vec v1=b-a,v2=p-a; 105 return fabs(cross(v1,v2)/len(v1)); 106 } 107 double ptseg(const Point& p,const Point& a,const Point& b) 108 //point_to_segment,点到线段距离,DistanceToSegment 109 { 110 if(a==b) return len(p-a); 111 //Vec v1=b-a,v2=a-p,v3=p-b; 112 Vec v1=b-a,v2=p-a,v3=p-b; 113 if(dcmp(dot(v1,v2))<0) return len(v2); 114 else if(dcmp(dot(v1,v3))>0) return len(v3); 115 else return fabs(cross(v1,v2)/len(v1)); 116 } 117 double area2(const Point& a,const Point& b,const Point& c) 118 //叉积对应平行四边形的面积 119 { 120 return cross(b-a,c-a); 121 } 122 double thrarea(const Point& a,const Point& b,const Point& c) 123 //三角形面积,绝对值 124 { 125 return fabs(cross(b-a,c-a)/2); 126 } 127 bool ifpar(const Vec& a,const Vec& b) 128 //ifParallel 129 //是否共线/平行 130 { 131 return dcmp(cross(a,b))==0; 132 } 133 Point pointline(const Point& p,const Point& a,const Vec& v) 134 //点在直线上投影,GetLineProjection 135 { 136 return a+v*(dot(p-a,v)/dot(v,v)); 137 } 138 double sqr(double x){return x*x;} 139 double dis2(const Point &x,const Point &y) 140 { 141 return sqr(x.x-y.x)+sqr(x.y-y.y); 142 } 143 double len2(const Vec &x){return dot(x,x);} 144 }; 145 using namespace X; 146 Point p[1000100]; 147 int n; 148 Point a; 149 double r2; 150 void make(const Point &x,const Point &y) 151 { 152 Vec t=(y-x)/2; 153 a=x+t; 154 r2=len2(t); 155 } 156 void make(const Point &x,const Point &y,const Point &z) 157 { 158 Vec t1=(y-x)/2,t2=(z-y)/2; 159 Point x1=x+t1,y1=y+t2; 160 t1=normal1(t1);t2=normal1(t2); 161 a=lineline(x1,t1,y1,t2); 162 r2=dis2(a,x); 163 } 164 bool ok(const Point &x) 165 { 166 return dcmp(dis2(x,a)-r2)<=0; 167 } 168 int main() 169 { 170 int i,j,k; 171 srand(time(0)); 172 scanf("%d",&n); 173 for(i=1;i<=n;i++) 174 scanf("%lf%lf",&p[i].x,&p[i].y); 175 random_shuffle(p+1,p+n+1); 176 a=p[1];r2=0; 177 for(i=2;i<=n;i++) 178 if(!ok(p[i])) 179 { 180 a=p[i];r2=0; 181 for(j=1;j<i;j++) 182 if(!ok(p[j])) 183 { 184 make(p[i],p[j]); 185 for(k=1;k<j;k++) 186 if(!ok(p[k])) 187 { 188 make(p[i],p[j],p[k]); 189 } 190 } 191 } 192 printf("%.2f %.2f ",a.x,a.y); 193 printf("%.2f\n",sqrt(r2)); 194 return 0; 195 }