跳转至

隨機增量法

引入

隨機增量算法是計算幾何的一個重要算法,它對理論知識要求不高,算法時間複雜度低,應用範圍廣大。

增量法 (Incremental Algorithm) 的思想與第一數學歸納法類似,它的本質是將一個問題化為規模剛好小一層的子問題。解決子問題後加入當前的對象。寫成遞歸式是:

\[ T(n)=T(n-1)+g(n) \]

增量法形式簡潔,可以應用於許多的幾何題目中。

增量法往往結合隨機化,可以避免最壞情況的出現。

最小圓覆蓋問題

題意描述

在一個平面上有 \(n\) 個點,求一個半徑最小的圓,能覆蓋所有的點。

過程

假設圓 \(O\) 是前 \(i-1\) 個點的最小覆蓋圓,加入第 \(i\) 個點,如果在圓內或邊上則什麼也不做。否則,新得到的最小覆蓋圓肯定經過第 \(i\) 個點。

然後以第 \(i\) 個點為基礎(半徑為 \(0\)),重複以上過程依次加入第 \(j\) 個點,若第 \(j\) 個點在圓外,則最小覆蓋圓必經過第 \(j\) 個點。

重複以上步驟。(因為最多需要三個點來確定這個最小覆蓋圓,所以重複三次)

遍歷完所有點之後,所得到的圓就是覆蓋所有點得最小圓。

性質

時間複雜度 \(O(n)\),證明詳見參考資料。

空間複雜度 \(O(n)\)

實現

代碼實現
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>

using namespace std;

int n;
double r;

struct point {
  double x, y;
} p[100005], o;

double sqr(double x) { return x * x; }

double dis(point a, point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)); }

bool cmp(double a, double b) { return fabs(a - b) < 1e-8; }

point geto(point a, point b, point c) {
  double a1, a2, b1, b2, c1, c2;
  point ans;
  a1 = 2 * (b.x - a.x), b1 = 2 * (b.y - a.y),
  c1 = sqr(b.x) - sqr(a.x) + sqr(b.y) - sqr(a.y);
  a2 = 2 * (c.x - a.x), b2 = 2 * (c.y - a.y),
  c2 = sqr(c.x) - sqr(a.x) + sqr(c.y) - sqr(a.y);
  if (cmp(a1, 0)) {
    ans.y = c1 / b1;
    ans.x = (c2 - ans.y * b2) / a2;
  } else if (cmp(b1, 0)) {
    ans.x = c1 / a1;
    ans.y = (c2 - ans.x * a2) / b2;
  } else {
    ans.x = (c2 * b1 - c1 * b2) / (a2 * b1 - a1 * b2);
    ans.y = (c2 * a1 - c1 * a2) / (b2 * a1 - b1 * a2);
  }
  return ans;
}

int main() {
  scanf("%d", &n);
  for (int i = 1; i <= n; i++) scanf("%lf%lf", &p[i].x, &p[i].y);
  for (int i = 1; i <= n; i++) swap(p[rand() % n + 1], p[rand() % n + 1]);
  o = p[1];
  for (int i = 1; i <= n; i++) {
    if (dis(o, p[i]) < r || cmp(dis(o, p[i]), r)) continue;
    o.x = (p[i].x + p[1].x) / 2;
    o.y = (p[i].y + p[1].y) / 2;
    r = dis(p[i], p[1]) / 2;
    for (int j = 2; j < i; j++) {
      if (dis(o, p[j]) < r || cmp(dis(o, p[j]), r)) continue;
      o.x = (p[i].x + p[j].x) / 2;
      o.y = (p[i].y + p[j].y) / 2;
      r = dis(p[i], p[j]) / 2;
      for (int k = 1; k < j; k++) {
        if (dis(o, p[k]) < r || cmp(dis(o, p[k]), r)) continue;
        o = geto(p[i], p[j], p[k]);
        r = dis(o, p[i]);
      }
    }
  }
  printf("%.10lf\n%.10lf %.10lf", r, o.x, o.y);
  return 0;
}

練習

最小圓覆蓋

「HNOI2012」射箭

CodeForces 442E

參考資料與擴展閲讀

http://www.doc88.com/p-007257893177.html

https://www.cnblogs.com/aininot260/p/9635757.html

https://wenku.baidu.com/view/162699d63186bceb19e8bbe6.html

https://blog.csdn.net/u014609452/article/details/62039612