Loj 6360 复燃「恋之埋火」

Posted on | Views:

高维空间的最小圆覆盖.

题意就是要求 $m$ 维空间内 $n$ 个点的最小圆覆盖.

把二维的随机增量算法直接拓展到高维,不难发现复杂度的分析是类似的.

当我们已经加入 $k$ 个点时,需要在这 $k$ 个点所在的 $k-1$ 维平面上找一个圆心,使得它到这 $k$ 个点距离相等.

从这 $k$ 个点中选一个点作为原点,它到另外 $k-1$ 个点的 $k-1$ 个向量显然就是这个平面的一组基底.

只需要求出这些向量各自的系数即可确定圆心.

根据圆心到原点和到其他任意一个点的距离相同,可以列出 $k-1$ 个方程,高斯消元求出系数,进而确定圆心位置.

时间复杂度 $O(nm^3)$ .

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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
//%std
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out = 0, fh = 1;
	char jp = getchar();
	while ((jp > '9' || jp < '0') && jp != '-')
		jp = getchar();
	if (jp == '-')
		fh = -1, jp = getchar();
	while (jp >= '0' && jp <= '9')
		out = out * 10 + jp - '0', jp = getchar();
	return out * fh;
}
const double eps = 1e-8;
double sq(double x)
{
	return x * x;
}
const int N = 2e4 + 10, M = 10;
int n, m;
typedef vector<double> vm;
vm operator + (vm A, vm B)
{
	vm C(m);
	for (int i = 0; i < m; ++i)
		C[i] = A[i] + B[i];
	return C;
}
vm operator - (vm A, vm B)
{
	vm C(m);
	for (int i = 0; i < m; ++i)
		C[i] = A[i] - B[i];
	return C;
}
vm operator * (vm A, double lambda)
{
	vm B(m);
	for (int i = 0; i < m; ++i)
		B[i] = A[i] * lambda;
	return B;
}
double dist(vm A, vm B)
{
	double s = 0;
	for (int i = 0; i < m; ++i)
		s += sq(A[i] - B[i]);
	return s;
}
int id[M], k = 0;
vm p[N], centre, vec[M];
double radius, a[M][M];
void Gauss()
{
	for (int i = 1; i < k; ++i)
	{
		int pos = i;
		for (int j = i + 1; j < k; ++j)
			if (fabs(a[j][i]) > fabs(a[pos][i]))
				pos = j;
		if (pos != i)
			swap(a[i], a[pos]);
		for (int j = 1; j < k; ++j)
			if (i != j)
			{
				double t = a[j][i] / a[i][i];
				if (fabs(t) <= eps)
					continue;
				for (int l = 1; l <= k; ++l)
					a[j][l] -= a[i][l] * t;					
			}
	}
}
void pr(vm x)
{
	for (int i = 0; i < m; ++i)
		printf("%.6lf,", x[i]);
	puts("");
}
void find_centre()
{
	for (int i = 0; i <= k; ++i)
		for (int j = 0; j <= k; ++j)
			a[i][j] = 0;
	for (int i = 1; i < k; ++i) 
		vec[i] = p[id[i]] - p[id[0]];
	for (int i = 1; i < k; ++i)
		for (int j = 0; j < m; ++j)
		{
			a[i][k] += sq(vec[i][j]);			
			for (int l = 1; l < k; ++l)
				a[i][l] += vec[i][j] * 2 * vec[l][j];  				
		}
	Gauss();
	for (int i = 0; i < m; ++i)
		centre[i] = 0;
	for (int i = 1; i < k; ++i)
		centre = centre + vec[i] * (a[i][k] / a[i][i]);
	centre = centre + p[id[0]];
}
void dfs(int x, int bound)
{
	if (x > m)
		return;
	for (int i = 0; i < bound; ++i)
	{
		if (dist(centre, p[i]) - radius > eps)
		{
			id[k++] = i;
			find_centre();
			radius = dist(centre, p[i]);
			dfs(x + 1, i);
			k--;
		}
	}
}
int main()
{
	srand(1919810);
	n = read(), m = read();
	for (int i = 0; i < n; ++i)
	{
		p[i].resize(m);
		for (int j = 0; j < m; ++j)
			scanf("%lf", &p[i][j]);
	}
	random_shuffle(p, p + n);
	centre = p[0], radius = -1;
	dfs(0, n);
	for (int i = 0; i < m; ++i)
		printf("%.6lf ", centre[i]);
	puts("");
	return 0;
}