G56 最小圆覆盖 随机增量法

视频链接：https://www.bilibili.com/video/BV1NL411X7xe/

1. Luogu P1742 最小圆覆盖

题意：给出 n 个点，让你画一个最小的包含所有点的圆

增量法：

1. 假设圆 O 是前 𝑖−1 个点的最小覆盖圆，加入第 𝑖 个点，如果在圆内或边上则什么也不做。否则，新得到的最小覆盖圆肯定经过第 𝑖 个点
2. 然后以第 𝑖 个点为基础（半径为0），加入第 𝑗 个点（𝑗<𝑖），若第 𝑗 个点在圆外，则最小覆盖圆必经过 𝑖,𝑗 两个点
3. 然后以 𝑖,𝑗 两点为基础，加入第 𝑘 个点（𝑘<𝑗<𝑖），若第 𝑘 个点在圆外，则最小覆盖圆必经过 𝑖,𝑗,𝑘 三个点
4. 遍历完所有点之后，所得到的圆就是覆盖所有点得最小圆

随机化：可以避免最坏情况的出现，时间为 O(n)

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

const int N=100010;
const double PI=acos(-1);
int n;
struct Point{double x,y;} p[N];
struct Circle{Point p; double r;} C;

Point operator+(Point a, Point b){
  return {a.x+b.x, a.y+b.y};
}
Point operator-(Point a, Point b){
  return {a.x-b.x, a.y-b.y};
}
Point operator*(Point a, double t){
  return {a.x*t, a.y*t};
}
Point operator/(Point a, double t){
  return {a.x/t, a.y/t};
}
double operator*(Point a, Point b){
  return a.x*b.y-a.y*b.x;
}
Point rotate(Point a, double b){
  return {a.x*cos(b)-a.y*sin(b),a.x*sin(b)+a.y*cos(b)};
}
double dis(Point a, Point b){
  return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point cross(Point a,Point u,Point b,Point v){ 
  double t=(a-b)*v/(v*u);
  return a+u*t; //直线交点
}
pair<Point, Point> midperp(Point a, Point b){
  return {(a+b)/2, rotate(b-a,PI/2)}; //中垂线
}
Circle cover(Point a, Point b){
  return {(a+b)/2, dis(a,b)/2}; //覆盖两点的圆
}
Circle cover(Point a, Point b, Point c){
  auto u=midperp(a, b), v=midperp(a, c);
  auto p=cross(u.x, u.y, v.x, v.y);
  return {p, dis(p, a)}; //覆盖三点的圆
}
void increment(){ //增量法
  C={p[1], 0};
  for(int i=2; i<=n; i++){
    if(C.r<dis(C.p, p[i])){
      C={p[i], 0}; //一点圆
      for(int j=1; j<i; j++){
        if(C.r<dis(C.p, p[j])){
          C=cover(p[i], p[j]); //两点圆
          for(int k=1; k<j; k++)
            if(C.r<dis(C.p, p[k]))
              C=cover(p[i],p[j],p[k]); //三点圆
        }
      }
    }
  }  
}
int main(){
  scanf("%d", &n);
  for(int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
  random_shuffle(p+1, p+n+1); //随机化
  increment(); //增量法
  printf("%.10lf\n%.10lf %.10lf\n", C.r, C.p.x, C.p.y);
  return 0;
}

2. Luogu P4288 [SHOI2014]信号增幅仪

题意：给出 n 个点，长轴方向和扩大倍数，求一个覆盖所有点的最小椭圆

思路：

1. 把坐标系逆时针旋转 a 度，让 x 轴沿长轴方向，沿 x 轴把椭圆压缩成圆。 这等价于所有点都顺时针旋转 a 度，再把所有点的 x 坐标变为原来的 1/b

2. 最后跑一个最小圆覆盖即可

时间：O(n)

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

const int N=50010;
const double PI=acos(-1);
int n;
struct Point{double x,y;} p[N];
struct Circle{Point p; double r;} C;

Point operator+(Point a, Point b){
  return {a.x+b.x, a.y+b.y};
}
Point operator-(Point a, Point b){
  return {a.x-b.x, a.y-b.y};
}
Point operator*(Point a, double t){
  return {a.x*t, a.y*t};
}
Point operator/(Point a, double t){
  return {a.x/t, a.y/t};
}
double operator*(Point a, Point b){
  return a.x*b.y-a.y*b.x;
}
Point rotate(Point a, double b){
  return {a.x*cos(b)-a.y*sin(b),a.x*sin(b)+a.y*cos(b)};
}
double dis(Point a, Point b){
  return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point cross(Point a,Point u,Point b,Point v){ 
  double t=(a-b)*v/(v*u);
  return a+u*t; //直线交点
}
pair<Point, Point> midperp(Point a, Point b){
  return {(a+b)/2, rotate(b-a,PI/2)}; //中垂线
}
Circle cover(Point a, Point b){
  return {(a+b)/2, dis(a,b)/2}; //覆盖两点的圆
}
Circle cover(Point a, Point b, Point c){
  auto u=midperp(a, b), v=midperp(a, c);
  auto p=cross(u.x, u.y, v.x, v.y);
  return {p, dis(p, a)}; //覆盖三点的圆
}
void increment(){ //增量法
  C={p[1], 0};
  for(int i=2; i<=n; i++){
    if(C.r<dis(C.p, p[i])){
      C={p[i], 0}; //一点圆
      for(int j=1; j<i; j++){
        if(C.r<dis(C.p, p[j])){
          C=cover(p[i], p[j]); //两点圆
          for(int k=1; k<j; k++)
            if(C.r<dis(C.p, p[k]))
              C=cover(p[i],p[j],p[k]); //三点圆
        }
      }
    }
  }  
}
int main(){
  scanf("%d", &n);
  for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
  double a, b;
  scanf("%lf%lf", &a, &b);
  for(int i=1; i<=n; i++){
    p[i]=rotate(p[i], -a/180*PI);
    p[i].x/=b; //椭圆压缩成圆
  }
  random_shuffle(p+1, p+n+1); //随机化
  increment(); //增量法
  printf("%.3lf\n", C.r);
}

练习

 HDU3007 Buried memory