刷題的時候發現了這么一個新的東西：Voronoi圖和Delaunay三角剖分
發現這個東西可以$O(nlogn)$解決平面圖最小生成樹問題感覺非常棒
然后就去學了..
看的n+e的blog，感謝n+e的耐心教導..
Voronoi圖是個啥
Delaunay三角剖分
最優三角剖分就是使每一個三角形的外接圓都不包含其他的點的三角剖分
這個算法就是求最優三角剖分的
簡單來說就是分治合並
對於點數小於等於$3$的可以直接連邊
合並的時候
1）先找到兩邊最下面的點，這個可以用凸包求，然后連邊
2）對於現在得到的兩個點$p_1$、$p_2$，找到一個點連接着$p_1$且由這三個點的外接圓不包含別的任何點，並刪除這個外接圓經過的邊，$p_2$也是如此
3）看現在找出來的兩個點$y_1$、$y_2$，找其中一個點使得它與$p_1$、$p_2$的外接圓不包含另外一個點，使其與對應的點連邊
4）重復（2）（3）直到無邊可連
對n+e代碼的修改
一開始壓根不知道怎么實現這個算法的時候去看了看n+e的代碼..
發現他的代碼中每次都要遍歷$p_1$和$p_2$的所有邊，這樣的做法在特殊的圖里面是$O(n^2)$的
可以說他自己出的數據偏水？？（霧
然后看了上面那篇詳細的文章，發現它的邊是按照順時針或者逆時針的方向進行取邊的，這樣遇到的邊要不刪掉要不就是最優的
所以我開了一個雙向鏈表來儲存邊，使得邊是按照$(-\dfrac{\pi}{2},\dfrac{3\pi}{2}]$順時針排列
然后就會發現細節一大堆..
調了我3天的東西就直接放上來好了
判斷一個點是否在某三個點的外接圓內
在這篇blog里面有講，但是貌似圖都爆掉了？？
就是把點映射到$z=x^2+y^2$的拋物面上，這三個點就會形成一個新的平面，然后再判斷剩下的那個點和這個平面的關系即可
下面放幾張n+e給我的圖片方便理解？
拋物面$z=x^2+y^2$
其實第三張圖應該很清晰了吧..
在圓內的點都會在平面下方，而在圓外的都會在平面上方
這個用三維叉積點積判一下就好了
Code
下面就貼一下我3天的成果吧..
還有什么問題請指教..
bzoj4219
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <vector>
#define eps 1e-10
using namespace std;
const int Maxn = 100010;
double _max(double x, double y) { return x > y ? x : y; }
double _min(double x, double y) { return x < y ? x : y; }
struct node {
int y, nxt, frt, opp;
}a[Maxn<<3]; int first[Maxn], last[Maxn], num[Maxn<<3], now;
int sta[Maxn], tp;
int pb(int x, int k, int y) {
int u = num[now++];
a[u].y = y; a[u].frt = k;
if(k){
a[u].nxt = a[k].nxt;
if(a[k].nxt) a[a[k].nxt].frt = u;
else last[x] = u;
a[k].nxt = u;
} else {
if(first[x]) a[first[x]].frt = u;
else last[x] = u;
a[u].nxt = first[x]; first[x] = u;
}
return u;
}
int pf(int x, int k, int y) {
int u = num[now++];
a[u].y = y; a[u].nxt = k;
if(k){
a[u].frt = a[k].frt;
if(a[k].frt) a[a[k].frt].nxt = u;
else first[x] = u;
a[k].frt = u;
} else {
if(last[x]) a[last[x]].nxt = u;
else first[x] = u;
a[u].frt = last[x]; last[x] = u;
}
return u;
}
void del(int x, int k) {
num[--now] = k;
if(a[k].nxt) a[a[k].nxt].frt = a[k].frt;
else last[x] = a[k].frt;
if(a[k].frt) a[a[k].frt].nxt = a[k].nxt;
else first[x] = a[k].nxt;
}
double _abs(double x) { return x < 0 ? -x : x; }
int zero(double x) { return _abs(x) < eps ? 1 : 0; }
struct Point {
double x, y;
Point(double x = 0, double y = 0) : x(x), y(y) {}
bool operator<(const Point &A) const { return zero(x-A.x) ? y < A.y : x < A.x; }
Point operator-(const Point &A) const { return Point(x-A.x, y-A.y); }
}list[Maxn]; int n; double X, Y;
double Cross(Point A, Point B) { return A.x*B.y-B.x*A.y; }
double Dot(Point A, Point B) { return A.x*B.x+A.y*B.y; }
double dis(Point A) { return sqrt(Dot(A, A)); }
double dis(int x, int y) { return dis(list[y]-list[x]); }
struct Point3 {
double x, y, z;
Point3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
Point3 operator-(const Point3 &A) const { return Point3(x-A.x, y-A.y, z-A.z); }
};
double Dot(Point3 A, Point3 B) { return A.x*B.x+A.y*B.y+A.z*B.z; }
Point3 Cross(Point3 A, Point3 B) { return Point3(A.y*B.z-A.z*B.y, A.z*B.x-A.x*B.z, A.x*B.y-A.y*B.x); }
Point3 t(Point A) { return Point3(A.x, A.y, A.x*A.x+A.y*A.y); }
bool incir(Point D, Point A, Point B, Point C) {
if(Cross(B-A, C-A) < -eps) swap(B, C);
Point3 aa = t(A), bb = t(B), cc = t(C), dd = t(D);
return Dot(Cross(bb-aa, cc-aa), dd-aa) < -eps;
}
bool incir(int D, int A, int B, int C) { return incir(list[D], list[A], list[B], list[C]); }
void divi(int L, int R) {
if(L == R) return;
if(L+1 == R){
int k1 = pb(L, 0, R); int k2 = pf(R, 0, L);
a[k1].opp = k2; a[k2].opp = k1;
return;
}
int mid = L+R>>1, i, j, k;
divi(L, mid); divi(mid+1, R);
int p1 = 0, p2 = 0; tp = 0;
for(i = L; i <= R; i++){
while(tp > 1 && Cross(list[i]-list[sta[tp-1]], list[sta[tp]]-list[sta[tp-1]]) > eps) tp--;
sta[++tp] = i;
}
for(i = 1; i < tp; i++) if(sta[i] <= mid && sta[i+1] > mid){ p1 = sta[i]; p2 = sta[i+1]; break; }
int kp1, kp2;
for(kp1 = last[p1]; kp1; kp1 = a[kp1].frt){
if(Cross(list[a[kp1].y]-list[p1], list[p2]-list[p1]) < eps || Cross(list[a[kp1].y]-list[p1], Point(0,1)) < -eps) break;
}
int k1 = pb(p1, kp1, p2);
for(kp2 = first[p2]; kp2; kp2 = a[kp2].nxt){
if(Cross(list[a[kp2].y]-list[p2], list[p1]-list[p2]) > -eps || Cross(list[a[kp2].y]-list[p2], Point(0,1)) > eps) break;
}
int k2 = pf(p2, kp2, p1);
a[k1].opp = k2; a[k2].opp = k1;
while(1){
int np1 = 0, np2 = 0;
for(; kp1; kp1 = a[kp1].frt){
if(Cross(list[a[kp1].y]-list[p1], list[p2]-list[p1]) > -eps){
if(Cross(list[a[kp1].y]-list[p1], Point(0,1)) > -eps) continue;
else break;
}
if(a[kp1].frt && incir(a[a[kp1].frt].y, p1, p2, a[kp1].y)) del(a[kp1].y, a[kp1].opp), del(p1, kp1);
else { np1 = kp1; break; }
}
for(; kp2; kp2 = a[kp2].nxt){
if(Cross(list[a[kp2].y]-list[p2], list[p1]-list[p2]) < eps){
if(Cross(list[a[kp2].y]-list[p2], Point(0,1)) < -eps) continue;
else break;
}
if(a[kp2].nxt && incir(a[a[kp2].nxt].y, p1, p2, a[kp2].y)) del(a[kp2].y, a[kp2].opp), del(p2, kp2);
else { np2 = kp2; break; }
}
if(!np1 && !np2) break;
if(!np2 || (np1 && !incir(a[kp2].y, p1, p2, a[kp1].y))){
p1 = a[kp1].y;
k2 = pf(p2, kp2, p1);
for(kp1 = last[p1]; kp1; kp1 = a[kp1].frt){
if(Cross(list[a[kp1].y]-list[p1], list[p2]-list[p1]) < eps || Cross(list[a[kp1].y]-list[p1], Point(0,1)) < -eps) break;
}
k1 = pb(p1, kp1, p2);
a[k1].opp = k2; a[k2].opp = k1;
} else {
p2 = a[kp2].y;
k1 = pb(p1, kp1, p2);
for(kp2 = first[p2]; kp2; kp2 = a[kp2].nxt){
if(Cross(list[a[kp2].y]-list[p2], list[p1]-list[p2]) > -eps || Cross(list[a[kp2].y]-list[p2], Point(0,1)) > eps) break;
}
k2 = pf(p2, kp2, p1);
a[k1].opp = k2; a[k2].opp = k1;
}
}
}
struct enode {
int x, y; double d;
enode(int x = 0, int y = 0, double d = 0) : x(x), y(y), d(d) {}
bool operator<(const enode &A) const { return d < A.d; }
}e[Maxn<<4]; int el;
int fa[Maxn];
int ff(int x) { return fa[x] == x ? x : fa[x] = ff(fa[x]); }
int main() {
int i, j, k;
scanf("%d%lf%lf", &n, &X, &Y);
for(i = 1; i <= n; i++) scanf("%lf%lf", &list[i].x, &list[i].y);
now = 1;
for(i = 1; i <= n<<3; i++) num[i] = i;
sort(list+1, list+n+1);
divi(1, n);
for(i = 1; i <= n; i++){
e[++el] = enode(i, n+1, _min(list[i].x, Y-list[i].y));
e[++el] = enode(i, n+2, _min(list[i].y, X-list[i].x));
for(k = first[i]; k; k = a[k].nxt) e[++el] = enode(i, a[k].y, dis(i, a[k].y)/2.0);
}
sort(e+1, e+el+1);
for(i = 1; i <= n+2; i++) fa[i] = i;
for(i = 1; i <= el; i++){
int fx = ff(e[i].x), fy = ff(e[i].y);
if(fx != fy){
fa[fx] = fy;
if(ff(n+1) == ff(n+2)){ printf("%lf\n", e[i].d); return 0; }
}
}
return 0;
}
小總結？？
雖然這次搞這個東西用的時間很長.. 但是收獲還是很大的..
最后差那么幾個點wa還是很想放棄的..
但是還是堅持下了來了嘛..
加油啊..