很好一題目，使我的最小圓覆蓋旋轉。
先假設 \(p = 1\)。這是最簡單的情況。這個時候我們就得到了一個裸的最小圓覆蓋。
當 \(p \not= 1\)，但是 \(a = 0\) 的時候。圓就變成了對稱軸與座標軸平行的橢圓，運用高中知識仿射一下，又回到了最小圓覆蓋。
在一般的情況下，我們先透過座標的旋轉回到第二種情況，再進行仿射，又是最小圓覆蓋。
目前為止實際上這個題目已經做完了。如果你還不會做，那麼你一定遇到了下面這個問題：
什麼是最小圓覆蓋
模板題傳送門
最小圓覆蓋問題是指：在二維平面上有一堆點，求一個最小半徑的圓，能夠將所有點全部都包住。
為解決這個問題，我們首先引入三個性質：
性質 1：最小覆蓋圓是唯一的
證明：假設存在兩個圓，它們的都是符合要求的最小圓。設它們的半徑為 \(r\)。那麼所有的點就必定落在下圖中的藍色區域：
而我們以兩圓的兩個交點為直徑再做一個圓（如下圖中紫色圓），顯然它包含了該區域中的任意點，並且其直徑為原本的兩圓的一條弦，必然小於原本的兩圓的直徑：
那麼，藍色圓不為最小圓，與假設衝突。
性質 2：對於一個點集，其最小覆蓋圓只為兩種情況之一：1. 圓由該點集內三個或者三個以上的點共同確定；2. 該點集內的兩個點確定圓的直徑
證明：記滿足該兩種情況為命題 \(p\)，是某個點集的最小覆蓋圓為命題 \(q\)。
假如點集內沒有點落在圓上，那麼顯然圓保持圓心不變，一直收縮。必然能夠到某一狀態，使得至少一個點落在圓上。
假如點集內只有一個點落在圓上，保持該點始終在圓上收縮半徑，也必定在某一狀態，使得另一點落在圓上。
假如點集內有兩點落在圓上，但兩點不構成直徑。讓圓心在兩點確定的線段的中垂線上移動，必然存在某一狀態，使得除了這兩點以外的某一點也落在圓上
綜上，\(\neg p\Rightarrow \neg q\)，即 \(p \Rightarrow q\)。
對於一個只有兩個點的點集，顯然條件二為最小覆蓋圓。
對於一個等邊三角形的三個頂點構成的點集，顯然條件一為最小覆蓋圓。
綜上，\(\exists\) 圓 \(O\) 滿足 \(p\)，使得 \(q\) 成立。
性質 3：對於一個已知點集 \(S\)，已知其最小覆蓋圓為圓 \(O_1\)。向點集內新增一點 \(P\)，若 \(P\) 在圓外，則新點集的最小覆蓋圓過點 \(P\)。
證明：對於一個新圓 \(O_2\)，其包含了點集 \(S\) 中的所有點和 \(P\)，且 \(P\) 在圓 \(O_2\) 內。
若 \(O_2\) 上只有兩個點。 \(S\) 中任意兩點間的距離必然小於圓 \(O_1\) 的直徑，那麼這兩點不可能構成 \(O_2\) 的直徑，與性質 2 衝突。
若 \(O_2\) 上有三個或三個以上的點，如圖（深藍色為 \(O_1\)，天青色為 \(O_2\)），\(S\) 中的點必然被包含在淺藍色區域內：
類似於性質 1 的證明，圓 \(O_1\) 不為點集 \(S\) 的最小覆蓋圓，與題設衝突。
於是得證。
引入這三個性質之後，我們就可以得到一個玄而又玄的演算法：隨機增量演算法。
首先將所有的點隨機化（隨機化的作用在之後），記該點集為 \(S\)，\(S_i\) 代表點集中的第 \(i\) 個點。
接著初始化一個半徑為 \(0\)，圓心為 \(S_1\) 的圓。
接著遍歷所有點，將當前列舉到的點丟入已知點集中，並維護一個新的最小覆蓋圓，以此類推。
問題在於，當我們已知前 \(i- 1\) 個點的最小覆蓋圓時，如何確定加入第 \(i\) 個點後的最小覆蓋圓。
若 \(S_i\) 在圓內，那麼直接繼續迴圈。
若 \(S_i\) 在圓外，由性質 3，新的最小覆蓋圓過點 \(S_i\)。此時我們已知圓上一點，不足以確定一個圓。
於是我們又初始化一個半徑為 \(0\)，圓心為 \(S_i\) 的圓。從 \(i - 1\) 開始向前列舉 \(j\)，此時性質 3 依然適用。
若 \(S_j\) 在圓內，那麼直接繼續迴圈。
若 \(S_j\) 在圓外，由性質 3，新的最小覆蓋圓過點 \(S_j\)。此時我們已知圓上兩點，依舊不夠。
於是我們又初始化一個以 \(S_i\) 和 \(S_j\) 為直徑的圓，從 \(1 \sim j - 1\) 列舉 \(k\)，依舊利用性質 \(3\)。
若 \(S_k\) 在圓內，那麼直接繼續迴圈。
若 \(S_k\) 在圓外，由性質 3，新的最小覆蓋圓過點 \(S_k\)。此時三點確定一個圓。
當 \(k\) 迴圈結束時，我們就得到了加入第 \(i\) 個點後的最小覆蓋圓。
核心程式碼如下：
random_shuffle(p + 1, p + n + 1);// 隨機化陣列
// 結構體 (coor){x 座標，y 座標}，陣列 p 是該型別
// 結構體 (circle){圓心，半徑}
c = (circle){p[1], 0};
for(int i = 2; i <= n; i++){
// double_cmp(double a, double b)：比較兩個浮點數 a 和 b 的大小
if(double_cmp(c.r, get_dis(c.o, p[i])) == -1){
c = (circle){p[i], 0};
for(int j = 1; j < i; j++){
if(double_cmp(c.r, get_dis(c.o, p[j])) == -1){
// get_dis(coor a, coor b)：兩點之間距離
c = (circle){(p[i] + p[j]) / 2, get_dis(p[i], p[j]) / 2};
for(int k = 1; k < j; k++){
// get_circle(coor a, coor b, coor c)：三點確定圓
if(double_cmp(c.r, get_dis(c.o, p[k])) == -1) c = get_circle(p[i], p[j], p[k]);
}
}
}
}
}
好現在有人問，這玩意兒不是 \(O(n^3)\) 的嗎？
這只是理論最壞時間複雜度。然而實際上，每次在進入下一層迴圈時，都會判定一次當前點是否在當前圓外。可以證明，當資料隨機時，我們只有 \(\frac{3}{n}\) 的機率遇到一個在圓外的點。這也是隨機化處理的原因。在隨機化處理後，時間複雜度就變為了 \(O(n\times\frac{3}{n}\times n\times\frac{3}{n}\times n)\)，即 \(O(n)\)。
完整程式碼如下
#include<bits/stdc++.h>
#define MAXN 100010
#define eps 1e-9
#define pi acos(-1)
using namespace std;
struct coor{ double x, y; };
struct circle{ coor o; double r; };
coor operator + (coor a, coor b){ return (coor){a.x + b.x, a.y + b.y}; }
coor operator - (coor a, coor b){ return (coor){a.x - b.x, a.y - b.y}; }
coor operator * (coor a, double b){ return (coor){a.x * b, a.y * b}; }
coor operator / (coor a, double b){ return (coor){a.x / b, a.y / b}; }
double operator * (coor a, coor b){ return a.x * b.y - a.y * b.x; }
int double_cmp(double a, double b){
if(fabs(a - b) < eps) return 0;
if(a < b) return -1;
return 1;
}
// 向量旋轉 b，單位為弧度
coor rotate(coor a, double b){ return (coor){a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)}; }
// 兩點間距離
double get_dis(coor a, coor b){
double dx = a.x - b.x;
double dy = a.y - b.y;
return sqrt(dx * dx + dy * dy);
}
// 兩直線交點，p 和 v，q 和 w 分別確定兩條直線
coor inter_section(coor p, coor v, coor q, coor w){
coor u = p - q;
double t = (w * u) / (v * w);
return p + v * t;
}
// 求中垂線，用兩點表示一條直線
pair<coor, coor> bisector(coor a, coor b){
coor p = (a + b) / 2;
coor q = rotate(b - a, pi/2);
return make_pair(p, q);
}
// 三點定圓
circle get_circle(coor a, coor b, coor c){
pair<coor, coor> x = bisector(a, b), y = bisector(a, c);
coor o = inter_section(x.first, x.second, y.first, y.second);
double r = get_dis(o, a);
return (circle){o, r};
}
int n;
circle c;
coor p[MAXN];
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);
c = (circle){p[1], 0};
for(int i = 2; i <= n; i++){
if(double_cmp(c.r, get_dis(c.o, p[i])) == -1){
c = (circle){p[i], 0};
for(int j = 1; j < i; j++){
if(double_cmp(c.r, get_dis(c.o, p[j])) == -1){
c = (circle){(p[i] + p[j]) / 2, get_dis(p[i], p[j]) / 2};
for(int k = 1; k < j; k++){
if(double_cmp(c.r, get_dis(c.o, p[k])) == -1) c = get_circle(p[i], p[j], p[k]);
}
}
}
}
}
printf("%.10lf\n",c.r);
printf("%.10f %.10f\n",c.o.x, c.o.y);
}
回到本題
好的以上是最小圓覆蓋，根據一開始的思路，我們很容易就能夠將上面那份板子改為本題的答案。
#include<bits/stdc++.h>
#define MAXN 100010
#define eps 1e-9
#define pi acos(-1)
using namespace std;
struct coor{ double x, y; };
struct circle{ coor o; double r; };
coor operator + (coor a, coor b){ return (coor){a.x + b.x, a.y + b.y}; }
coor operator - (coor a, coor b){ return (coor){a.x - b.x, a.y - b.y}; }
coor operator * (coor a, double b){ return (coor){a.x * b, a.y * b}; }
coor operator / (coor a, double b){ return (coor){a.x / b, a.y / b}; }
double operator * (coor a, coor b){ return a.x * b.y - a.y * b.x; }
int double_cmp(double a, double b){
if(fabs(a - b) < eps) return 0;
if(a < b) return -1;
return 1;
}
coor rotate(coor a, double b){ return (coor){a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)}; }
double get_dis(coor a, coor b){
double dx = a.x - b.x;
double dy = a.y - b.y;
return sqrt(dx * dx + dy * dy);
}
coor inter_section(coor p, coor v, coor q, coor w){
coor u = p - q;
double t = (w * u) / (v * w);
return p + v * t;
}
pair<coor, coor> bisector(coor a, coor b){
coor p = (a + b) / 2;
coor q = rotate(b - a, pi/2);
return make_pair(p, q);
}
circle get_circle(coor a, coor b, coor c){
pair<coor, coor> x = bisector(a, b), y = bisector(a, c);
coor o = inter_section(x.first, x.second, y.first, y.second);
double r = get_dis(o, a);
return (circle){o, r};
}
int n, a, t;
double b;
circle c;
coor p[MAXN];
int main(){
scanf("%d", &n);
for(int i = 1; i <= n; i++) scanf("%lf%lf",&p[i].x,&p[i].y);
scanf("%d%d",&a,&t);
b = a / 180.0 * pi;
// 旋轉
for(int i = 1; i <= n; i++) p[i] = rotate(p[i], b);
// 仿射
for(int i = 1; i <= n; i++) p[i].x /= (double)t;
random_shuffle(p + 1, p + n + 1);
c = (circle){p[1], 0};
for(int i = 2; i <= n; i++){
if(double_cmp(c.r, get_dis(c.o, p[i])) == -1){
c = (circle){p[i], 0};
for(int j = 1; j < i; j++){
if(double_cmp(c.r, get_dis(c.o, p[j])) == -1){
c = (circle){(p[i] + p[j]) / 2, get_dis(p[i], p[j]) / 2};
for(int k = 1; k < j; k++){
if(double_cmp(c.r, get_dis(c.o, p[k])) == -1) c = get_circle(p[i], p[j], p[k]);
}
}
}
}
}
c.r *= 1000;
c.r = round(c.r);
c.r /= 1000;
printf("%.3lf\n",c.r);
}