竞赛模板整理

数学

高斯消元 O(n^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
#include <bits/stdc++.h>
using namespace std;
const double EPS = 1e-9;
const int INF = 0x3f3f3f3f;

int gauss (vector < vector<double> > a, vector<double> & ans) {
int n = (int) a.size();
int m = (int) a[0].size() - 1;

vector<int> where (m, -1);
for (int col=0, row=0; col<m && row<n; ++col) {
int sel = row;
for (int i=row; i<n; ++i)
if (abs (a[i][col]) > abs (a[sel][col]))
sel = i;
if (abs (a[sel][col]) < EPS)
continue;
for (int i=col; i<=m; ++i)
swap (a[sel][i], a[row][i]);
where[col] = row;

for (int i=0; i<n; ++i)
if (i != row) {
double c = a[i][col] / a[row][col];
for (int j=col; j<=m; ++j)
a[i][j] -= a[row][j] * c;
}
++row;
}

ans.assign (m, 0);
for (int i=0; i<m; ++i)
if (where[i] != -1)
ans[i] = a[where[i]][m] / a[where[i]][i];
for (int i=0; i<n; ++i) {
double sum = 0;
for (int j=0; j<m; ++j)
sum += ans[j] * a[i][j];
if (abs (sum - a[i][m]) > EPS)
return 0;
}

for (int i=0; i<m; ++i)
if (where[i] == -1)
return INF;
return 1;
}

int main() {
int n, m;
cin >> n >> m;
vector<vector<double> > a(n, vector<double>(m + 1));
vector<double> ans(m);
for(int i = 0; i < n; ++i) {
for(int j = 0; j <= m; ++j) {
cin >> a[i][j];
}
}
int s = gauss(a, ans);
if(s == 0) cout << "NO" << endl;
else if(s == INF) cout << "MANY" << endl;
else {
for(double &i : ans) {
cout << i << ' ';
}
}
cout << endl;
return 0;
}

等价矩阵 O(n^3)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
const double EPS = 1E-9;

int rank = max(n,m);
vector<char> line_used (n);
for (int i=0; i<m; ++i) {
int j;
for (j=0; j<n; ++j)
if (!line_used[j] && abs(a[j][i]) > EPS)
break;
if (j == n)
--rank;
else {
line_used[j] = true;
for (int p=i+1; p<m; ++p)
a[j][p] /= a[j][i];
for (int k=0; k<n; ++k)
if (k != j && abs (a[k][i]) > EPS)
for (int p=i+1; p<m; ++p)
a[k][p] -= a[j][p] * a[k][i];
}
}

行列式的值 O(n^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
const double EPS = 1E-9;
int n;
vector < vector<double> > a (n, vector<double> (n));
... чтение n и a ...

double det = 1;
for (int i=0; i<n; ++i) {
int k = i;
for (int j=i+1; j<n; ++j)
if (abs (a[j][i]) > abs (a[k][i]))
k = j;
if (abs (a[k][i]) < EPS) {
det = 0;
break;
}
swap (a[i], a[k]);
if (i != k)
det = -det;
det *= a[i][i];
for (int j=i+1; j<n; ++j)
a[i][j] /= a[i][i];
for (int j=0; j<n; ++j)
if (j != i && abs (a[j][i]) > EPS)
for (int k=i+1; k<n; ++k)
a[j][k] -= a[i][k] * a[j][i];
}

cout << det;

离散数对

$a ^ x ≡ b\ \mod\ {m}$ 以知 $a,b,m$ 求解 $x$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int solve (int a, int b, int m) {
int n = (int) sqrt (m + .0) + 1;

int an = 1;
for (int i=0; i<n; ++i)
an = (an * a) % m;

map<int,int> vals;
for (int i=1, cur=an; i<=n; ++i) {
if (!vals.count(cur))
vals[cur] = i;
cur = (cur * an) % m;
}

for (int i=0, cur=b; i<=n; ++i) {
if (vals.count(cur)) {
int ans = vals[cur] * n - i;
if (ans < m)
return ans;
}
cur = (cur * a) % m;
}
return -1;
}

离散根提取

$x ^ k ≡ a \ \mod\ {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
67
68
69
70
71
72
73
int gcd (int a, int b) {
return a ? gcd (b%a, a) : b;
}

int powmod (int a, int b, int p) {
int res = 1;
while (b)
if (b & 1)
res = int (res * 1ll * a % p), --b;
else
a = int (a * 1ll * a % p), b >>= 1;
return res;
}

int generator (int p) {
vector<int> fact;
int phi = p-1, n = phi;
for (int i=2; i*i<=n; ++i)
if (n % i == 0) {
fact.push_back (i);
while (n % i == 0)
n /= i;
}
if (n > 1)
fact.push_back (n);

for (int res=2; res<=p; ++res) {
bool ok = true;
for (size_t i=0; i<fact.size() && ok; ++i)
ok &= powmod (res, phi / fact[i], p) != 1;
if (ok) return res;
}
return -1;
}

int main() {
int n, k, a;
cin >> n >> k >> a;
if (a == 0) {
puts ("1\n0");
return 0;
}

int g = generator (n);

int sq = (int) sqrt (n + .0) + 1;
vector < pair<int,int> > dec (sq);
for (int i=1; i<=sq; ++i)
dec[i-1] = make_pair (powmod (g, int (i * sq * 1ll * k % (n - 1)), n), i);
sort (dec.begin(), dec.end());
int any_ans = -1;
for (int i=0; i<sq; ++i) {
int my = int (powmod (g, int (i * 1ll * k % (n - 1)), n) * 1ll * a % n);
vector < pair<int,int> >::iterator it =
lower_bound (dec.begin(), dec.end(), make_pair (my, 0));
if (it != dec.end() && it->first == my) {
any_ans = it->second * sq - i;
break;
}
}
if (any_ans == -1) {
puts ("0");
return 0;
}

int delta = (n-1) / gcd (k, n-1);
vector<int> ans;
for (int cur=any_ans%delta; cur<n-1; cur+=delta)
ans.push_back (powmod (g, cur, n));
sort (ans.begin(), ans.end());
printf ("%d\n", ans.size());
for (size_t i=0; i<ans.size(); ++i)
printf ("%d ", ans[i]);

格雷码

1
2
3
4
5
6
7
8
int g (int n) {
return n ^ (n >> 1);
}
int rev_g (int g) {
int n = 0;
for (; g; g>>=1) n ^= g;
return n;
}

求阶乘 $O(p \log_p n).$

1
2
3
4
5
6
7
8
9
10
int factmod (int n, int p) {
int res = 1;
while (n > 1) {
res = (res * ((n/p) % 2 ? p-1 : 1)) % p;
for (int i=2; i<=n%p; ++i)
res = (res * i) % p;
n /= p;
}
return res % p;
}

原根

如果它 $g$ 是一个原始根模数 $n$ ,那么对于任何一这样$gcd(a,n)= 1$的整数,都有一个 $ķ$ 这样的整数 $g ^ k ≡ a \ \mod\ {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
int powmod (int a, int b, int p) {
int res = 1;
while (b)
if (b & 1)
res = int (res * 1ll * a % p), --b;
else
a = int (a * 1ll * a % p), b >>= 1;
return res;
}

int generator (int p) {
vector<int> fact;
int phi = p-1, n = phi;
for (int i=2; i*i<=n; ++i)
if (n % i == 0) {
fact.push_back (i);
while (n % i == 0)
n /= i;
}
if (n > 1)
fact.push_back (n);

for (int res=2; res<=p; ++res) {
bool ok = true;
for (size_t i=0; i<fact.size() && ok; ++i)
ok &= powmod (res, phi / fact[i], p) != 1;
if (ok) return res;
}
return -1;
}

卡特兰数应用

  • 由 $n$ 开 $n$ 括号和右括号组成的有效括号序列的数量。
  • 带 $n + 1$ 叶子的根二叉树的数量(顶点没有编号)。
  • 完全划分 $n + 1$ 乘数的方法的数量。
  • 凸 $n + 2$ 角的三角测量的数量(即,通过将对角线不相交成三角形的多边形的分区数)。
  • 连接 $2N$ 圆上的点与 $n$ 非交叉和弦的方法的数量。
  • 具有ñ内部顶点的非同构完整二叉树的数量(即,至少有一个子)。
  • 从点数量单调路径 $(0,0)$ 指向 $(n,n)$ 在一个方形网格尺寸 $n \times n$ 上方的主对角线没有上升。
  • $n$ 可以通过堆栈排序的长度排列的数量(当且仅当没有这样的索引时 $i<j <k$,可以显示排列按堆栈排序 $a_k <a_i <a_j$ )。
  • 一组 $n$ 元素的连续分区数(即分区为连续块)。
  • $1 \ldots n$ 使用 $n$ 矩形覆盖梯子的方式的数量(意味着由 $n$ 列组成的图,其中一个我具有高度我)。

康托展开(求任意排列的排名)

康托展开可以用来求一个 $1\sim n$ 的任意排列的排名。

什么是排列的排名?

把 $1\sim n$ 的所有排列按字典序排序,这个排列的位次就是它的排名。

时间复杂度?

康托展开可以在 $O(n^2)$ 的复杂度内求出一个排列的排名,在用到树状数组优化时可以做到 $O(n\log n)$ 。

怎么实现?

因为排列是按字典序排名的,因此越靠前的数字优先级越高。也就是说如果两个排列的某一位之前的数字都相同,那么如果这一位如果不相同,就按这一位排序。

比如 $4$ 的排列,$[2,3,1,4]<[2,3,4,1]$,因为在第 $3$ 位出现不同,则 $[2,3,1,4]$ 的排名在 $[2,3,4,1]$ 前面。

举个栗子

我们知道长为 $5$ 的排列 $[2,5,3,4,1]$ 大于以 $1$ 为第一位的任何排列,以 $1$ 为第一位的 $5$ 的排列有 $4!$ 种。这是非常好理解的。但是我们对第二位的 $5$ 而言,它大于第一位与这个排列相同的,而这一位比 $5$ 小的所有排列。不过我们要注意的是,这一位不仅要比 $5$ 小,还要满足没有在当前排列的前面出现过,不然统计就重复了。因此这一位为 $1,3$ 或 $4$ ,第一位为 $2$ 的所有排列都比它要小,数量为 $3\times 3!$。

按照这样统计下去,答案就是 $1+4!+3\times 3!+2!+1=46$。注意我们统计的是排名,因此最前面要 $+1$。

注意到我们每次要用到当前有多少个小于它的数还没有出现,这里用树状数组统计比它小的数出现过的次数就可以了。

逆康托展开

因为排列的排名和排列是一一对应的,所以康托展开满足双射关系,是可逆的。可以通过类似上面的过程倒推回来。

如果我们知道一个排列的排名,就可以推出这个排列。因为 $4!$ 是严格大于 $3\times 3!+2\times 2!+1\times 1!$ 的,所以可以认为对于长度为 $5$ 的排列,排名 $x$ 除以 $4!$ 向下取整就是有多少个数小于这个排列的第一位。

引用上面展开的例子

首先让 $46-1=45$,$45$ 代表着有多少个排列比这个排列小。$\lfloor\frac {45}{4!}\rfloor=1$,有一个数小于它,所以第一位是 $2$。

此时让排名减去 $1\times 4!$得到$21$,$\lfloor\frac {21}{3!}\rfloor=3$,有 $3$ 个数小于它,去掉已经存在的 $2$,这一位是 $5$。

$21-3\times 3!=3$,$\lfloor\frac {3}{2!}\rfloor=1$,有一个数小于它,那么这一位就是 $3$。

让 $3-1\times 2!=1$,有一个数小于它,这一位是剩下来的第二位,$4$,剩下一位就是 $1$。即 $[2,5,3,4,1]$。

实际上我们得到了形如有两个数小于它这一结论,就知道它是当前第 $3$ 个没有被选上的数,这里也可以用线段树维护,时间复杂度为 $O(n\log n)$。

约瑟夫

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int joseph (int n, int k) {
return n>1 ? (joseph (n-1, k) + k - 1) % n + 1 : 1;
}
int joseph (int n, int k) {
if (n == 1) return 0;
if (k == 1) return n-1;
if (k > n) return (joseph (n-1, k) + k) % n;
int cnt = n / k;
int res = joseph (n - cnt, k);
res -= n % k;
if (res < 0) res += n;
else res += res / (k - 1);
return res;
}

Stirling 数(子集划分)

根据例题来讲解:

(2007 普及)将$n$个数$(1,2,…,n)$分成$r$个部分。每个部分至少一个数。将不同划分方法的总数记为$S_n^r$。例如,$S_4^2=7$,这 7 种不同的划分方法依次为 ${\ (1) , (234) }\,{\ (2) , (134) }\,{\ (3) , (124) }\,{\ (4) , (123) }\,{\ (12) , (34) }\,{\ (13) , (24) }\,{\ (14) , (23) }$。当$n=6,r=3$时,$S_6^3$=( )

提示:先固定一个数,对于其余的 5 个数考虑$S_5^3$与$S_5^2$,再分这两种情况对原固定的数进行分析。

题解:在近几年算法竞赛中,递推算法越来越重要:

$$
S_6^3=3 \times S_5^3 + S_5^2
$$

$$
S_5^3=3 \times S_4^3 + S_4^2
$$

$$
S_5^2=2 \times S_4^2 + S_4^1
$$

第二类 stirling 数,显然拥有这样的性质:

$$
S_n^m = m \times S_{n-1}^{m} + S_{n-1}^{m-1}
$$

$$
S_n^1 = 1,S_n^0 = 0,S_n^n = 1
$$

而这些性质就可以总结成:

$$
S_n^3 = \frac{1}{2} \times (3^{n-1}+1) - 2^{n-1}
$$

整数拆分

正整数n的有序k分拆的个数等于

$$p_k(n)=\binom{n-1}{k-1}$$
$p_0(0)=1, p_k(0)=0, p_1(n)=1, p_n(n)=1$
$p_k(n)=p_1(n-k)+p_2(n-k)+\ldots +p_k(n-k) =p_{k-1}(n-1)+p_k(n-k)$

n个球 r 个盒子 是否允许有空盒 分配方案数
不同 不同 允许 $r^n$
不同 不同 不允许 $r!S(n,r)$
不同 相同 允许 $\sum_{k=1}^{r}S(n,k)$
不同 相同 不允许 $S(n,r)$
相同 不同 允许 $C(n+r-1, n)$
相同 不同 不允许 $C(n-1, r-1)$
相同 相同 允许 $\sum_{k=1}^{r}p_k(n)$
相同 相同 不允许 $p_r(n)$

$S(n,r)$ 为斯特林数

博弈

任意图表上的博弈

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
vector<int> g [100];
bool win [100];
bool loose [100];
bool used[100];
int degree[100];

void dfs (int v) {
used[v] = true;
for (vector<int>::iterator i = g[v].begin(); i != g[v].end(); ++i)
if (!used[*i]) {
if (loose[v])
win[*i] = true;
else if (--degree[*i] == 0)
loose[*i] = true;
else
continue;
dfs (*i);
}
}

FFT

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
#include <bits/stdc++.h>
using namespace std;

typedef complex<double> complex_t;
const double PI = acos(-1.0);
const int MAXN = 1 << 21;
complex_t A[MAXN], B[MAXN];

complex_t w[MAXN + 5];
int bitrev[MAXN + 5];
#define rep(i, l, r) for(int i = l; i < r; ++i)
int N, L;

void FFT(complex_t *a, int op = 1, int n = N) { // op 1 时是 DFT,on == -1 时是 IDFT
rep(i, 0, n) if (i < bitrev[i]) swap(a[i], a[bitrev[i]]);
for (int i = 2, lyc = n >> 1; i <= n; i <<= 1, lyc >>= 1)
for (int j = 0; j < n; j += i) {
complex_t *l = a + j, *r = a + j + (i >> 1), *p = w;
rep(k, 0, i >> 1) {
complex_t tmp = *r * *p;
*r = *l - tmp, *l = *l + tmp;
++l, ++r, p += lyc;
}
}
if(op == -1) {
for(int i = 0; i < N; ++i) {
a[i] = complex_t(a[i].real() / N, a[i].imag() / N);
}
}
}

inline void fft_prepare(int n) { //init
L = 0;
for ( ; (1 << L) < n - 1; ++L);
N = 1 << L;
rep(i, 0, N) bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (L - 1));
rep(i, 0, N) w[i] = complex_t(cos(2 * PI * i / N), sin(2 * PI * i / N));
}

int main() {
int n; cin >> n;
fft_prepare(n);
FFT(A, 1);
FFT(A, 1);
for(int i = 0; i < N; ++i) A[i] *= B[i];
FFT(A, -1);
// (int)(A[i].real() + 0.5)
return 0;
}

计算几何

O(N) 求简单多边形面积

1
2
3
4
5
6
7
8
9
10
double sq (const vector<point> & fig) {
double res = 0;
for (unsigned i=0; i<fig.size(); i++) {
point
p1 = i ? fig[i-1] : fig.back(),
p2 = fig[i];
res += (p1.x - p2.x) * (p1.y + p2.y);
}
return fabs (res) / 2;
}

凸多边形内接圆

先三分圆心的 $x$ 坐标,再三分圆心的 $y$ 坐标,$O(N \log ^2 C)$

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
const double EPS = 1E-9;
int steps = 60;

struct pt {
double x, y;
};

struct line {
double a, b, c;
};

double dist (double x, double y, line & l) {
return abs (x * l.a + y * l.b + l.c);
}

double radius (double x, double y, vector<line> & l) {
int n = (int) l.size();
double res = INF;
for (int i=0; i<n; ++i)
res = min (res, dist (x, y, l[i]));
return res;
}

double y_radius (double x, vector<pt> & a, vector<line> & l) {
int n = (int) a.size();
double ly = INF, ry = -INF;
for (int i=0; i<n; ++i) {
int x1 = a[i].x, x2 = a[(i+1)%n].x, y1 = a[i].y, y2 = a[(i+1)%n].y;
if (x1 == x2) continue;
if (x1 > x2) swap (x1, x2), swap (y1, y2);
if (x1 <= x+EPS && x-EPS <= x2) {
double y = y1 + (x - x1) * (y2 - y1) / (x2 - x1);
ly = min (ly, y);
ry = max (ry, y);
}
}
for (int sy=0; sy<steps; ++sy) {
double diff = (ry - ly) / 3;
double y1 = ly + diff, y2 = ry - diff;
double f1 = radius (x, y1, l), f2 = radius (x, y2, l);
if (f1 < f2)
ly = y1;
else
ry = y2;
}
return radius (x, ly, l);
}

int main() {

int n;
vector<pt> a (n);
... чтение a ...

vector<line> l (n);
for (int i=0; i<n; ++i) {
l[i].a = a[i].y - a[(i+1)%n].y;
l[i].b = a[(i+1)%n].x - a[i].x;
double sq = sqrt (l[i].a*l[i].a + l[i].b*l[i].b);
l[i].a /= sq, l[i].b /= sq;
l[i].c = - (l[i].a * a[i].x + l[i].b * a[i].y);
}

double lx = INF, rx = -INF;
for (int i=0; i<n; ++i) {
lx = min (lx, a[i].x);
rx = max (rx, a[i].x);
}

for (int sx=0; sx<stepsx; ++sx) {
double diff = (rx - lx) / 3;
double x1 = lx + diff, x2 = rx - diff;
double f1 = y_radius (x1, a, l), f2 = y_radius (x2, a, l);
if (f1 < f2)
lx = x1;
else
rx = x2;
}

double ans = y_radius (lx, a, l);
printf ("%.7lf", ans);

}

通过“缩小边”的方法在凸多边形中找到内切圆 O(n \ log 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
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

const double EPS = 1E-9;
const double PI = ...;

struct pt {
double x, y;
pt() { }
pt (double x, double y) : x(x), y(y) { }
pt operator- (const pt & p) const {
return pt (x-p.x, y-p.y);
}
};

double dist (const pt & a, const pt & b) {
return sqrt ((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

double get_ang (const pt & a, const pt & b) {
double ang = abs (atan2 (a.y, a.x) - atan2 (b.y, b.x));
return min (ang, 2*PI-ang);
}

struct line {
double a, b, c;
line (const pt & p, const pt & q) {
a = p.y - q.y;
b = q.x - p.x;
c = - a * p.x - b * p.y;
double z = sqrt (a*a + b*b);
a/=z, b/=z, c/=z;
}
};

double det (double a, double b, double c, double d) {
return a * d - b * c;
}

pt intersect (const line & n, const line & m) {
double zn = det (n.a, n.b, m.a, m.b);
return pt (
- det (n.c, n.b, m.c, m.b) / zn,
- det (n.a, n.c, m.a, m.c) / zn
);
}

bool parallel (const line & n, const line & m) {
return abs (det (n.a, n.b, m.a, m.b)) < EPS;
}

double get_h (const pt & p1, const pt & p2,
const pt & l1, const pt & l2, const pt & r1, const pt & r2)
{
pt q1 = intersect (line (p1, p2), line (l1, l2));
pt q2 = intersect (line (p1, p2), line (r1, r2));
double l = dist (q1, q2);
double alpha = get_ang (l2 - l1, p2 - p1) / 2;
double beta = get_ang (r2 - r1, p1 - p2) / 2;
return l * sin(alpha) * sin(beta) / sin(alpha+beta);
}

struct cmp {
bool operator() (const pair<double,int> & a, const pair<double,int> & b) const {
if (abs (a.first - b.first) > EPS)
return a.first < b.first;
return a.second < b.second;
}
};

int main() {
int n;
vector<pt> p;
... чтение n и p ...

vector<int> next (n), prev (n);
for (int i=0; i<n; ++i) {
next[i] = (i + 1) % n;
prev[i] = (i - 1 + n) % n;
}

set < pair<double,int>, cmp > q;
vector<double> h (n);
for (int i=0; i<n; ++i) {
h[i] = get_h (
p[i], p[next[i]],
p[i], p[prev[i]],
p[next[i]], p[next[next[i]]]
);
q.insert (make_pair (h[i], i));
}

double last_time;
while (q.size() > 2) {
last_time = q.begin()->first;
int i = q.begin()->second;
q.erase (q.begin());

next[prev[i]] = next[i];
prev[next[i]] = prev[i];
int nxt = next[i], nxt1 = (nxt+1)%n,
prv = prev[i], prv1 = (prv+1)%n;
if (parallel (line (p[nxt], p[nxt1]), line (p[prv], p[prv1])))
break;

q.erase (make_pair (h[nxt], nxt));
q.erase (make_pair (h[prv], prv));

h[nxt] = get_h (
p[nxt], p[nxt1],
p[prv1], p[prv],
p[next[nxt]], p[(next[nxt]+1)%n]
);
h[prv] = get_h (
p[prv], p[prv1],
p[(prev[prv]+1)%n], p[prev[prv]],
p[nxt], p[nxt1]
);

q.insert (make_pair (h[nxt], nxt));
q.insert (make_pair (h[prv], prv));
}

cout << last_time << endl;
}

找到三角形联合的区域。垂直分解法

给定N个三角形。需要找到他们关联的区域。

决定

这里我们考虑垂直分解的方法,这在几何问题中通常非常重要。

所以,我们有N个三角形可以以任何方式相互交叉。我们将在垂直分解的帮助下摆脱这些交叉点:我们将找到所有段的所有交叉点(形成三角形),并通过其横坐标对找到的点进行排序。假设我们得到了一些数组B.让我们沿着这个数组移动。在第i步,我们考虑元素B [i]和B [i + 1]。我们在直线X = B [i]和X = B [i + 1]之间具有垂直条带,并且根据阵列B的结构,在该条带内,这些区段彼此不重叠。因此,在该条带内部,三角形被切割成梯形,并且条带内的这些梯形的侧面根本不相交。我们将从下往上沿着这些梯形的两侧移动,并折叠梯形区域,确保每个部分只教一次。实际上,这个过程与嵌套括号的处理非常相似。通过在每个条带内添加梯形区域,并为所有条带添加结果,我们将找到答案 - 三角形的并集区域。

让我们再一次考虑从实现的角度来看添加梯形区域的过程。我们迭代所有三角形的所有边,如果任何边(不垂直,我们不需要垂直边,甚至反之亦然,它们将阻挡)进入这个垂直条(完全或部分),然后我们把这边放在一些矢量以这种形式进行此操作最方便:侧面与垂直条带边界的交点处的Y坐标,以及三角形的编号。在我们构建了包含边的部分的这个向量之后,我们按Y的值排序:首先是左边的Y,然后是右边的Y.结果,向量中的第一个元素将包含最低梯形的下边。现在我们来看看收到的矢量。让我成为当前的元素; 这意味着第i件是一些梯形的底面,一些块(可能包含几个空格),我们想要立即添加到答案的区域。因此,我们将某个三角形计数器设置为1,然后向上移动,如果我们第一次遇到三角形的边,则增加计数器,如果我们第二次遇到三角形,则减少计数器。如果在某个段j上计数器变为等于零,那么我们找到了块的上边界 - 我们停在此处,添加由段i和j限定的梯形区域,并且我分配j + 1并再次重复整个过程。如果我们第一次遇到三角形的一侧,并减少计数器,如果我们第二次遇到三角形。如果在某个段j上计数器变为等于零,那么我们找到了块的上边界 - 我们停在此处,添加由段i和j限定的梯形区域,并且我分配j + 1并再次重复整个过程。如果我们第一次遇到三角形的一侧,并减少计数器,如果我们第二次遇到三角形。如果在某个段j上计数器变为等于零,那么我们找到了块的上边界 - 我们停在此处,添加由段i和j限定的梯形区域,并且我分配j + 1并再次重复整个过程。

因此,由于垂直分解方法,我们使用仅使用两个段的交集的几何图元来解决这个问题。

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
struct segment {
int x1, y1, x2, y2;
};

struct point {
double x, y;
};

struct item {
double y1, y2;
int triangle_id;
};

struct line {
int a, b, c;
};

const double EPS = 1E-7;

void intersect (segment s1, segment s2, vector<point> & res) {
line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.a*s1.x1+l1.b*s1.y1 },
l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.a*s2.x1+l2.b*s2.y1 };
double det1 = l1.a * l2.b - l1.b * l2.a;
if (abs (det1) < EPS) return;
point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,
(l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };
if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x <= s2.x2+EPS)
res.push_back (p);
}

double segment_y (segment s, double x) {
return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);
}

bool eq (double a, double b) {
return abs (a-b) < EPS;
}

vector<item> c;

bool cmp_y1_y2 (int i, int j) {
const item & a = c[i];
const item & b = c[j];
return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;
}

int main() {
int n;
cin >> n;
vector<segment> a (n*3);
for (int i=0; i<n; ++i) {
int x1, y1, x2, y2, x3, y3;
scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3);
segment s1 = { x1,y1,x2,y2 };
segment s2 = { x1,y1,x3,y3 };
segment s3 = { x2,y2,x3,y3 };
a[i*3] = s1;
a[i*3+1] = s2;
a[i*3+2] = s3;
}

for (size_t i=0; i<a.size(); ++i)
if (a[i].x1 > a[i].x2)
swap (a[i].x1, a[i].x2), swap (a[i].y1, a[i].y2);

vector<point> b;
b.reserve (n*n*3);
for (size_t i=0; i<a.size(); ++i)
for (size_t j=i+1; j<a.size(); ++j)
intersect (a[i], a[j], b);

vector<double> xs (b.size());
for (size_t i=0; i<b.size(); ++i)
xs[i] = b[i].x;
sort (xs.begin(), xs.end());
xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());

double res = 0;
vector<char> used (n);
vector<int> cc (n*3);
c.resize (n*3);
for (size_t i=0; i+1<xs.size(); ++i) {
double x1 = xs[i], x2 = xs[i+1];
size_t csz = 0;
for (size_t j=0; j<a.size(); ++j)
if (a[j].x1 != a[j].x2)
if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) {
item it = { segment_y (a[j], x1), segment_y (a[j], x2), (int)j/3 };
cc[csz] = (int)csz;
c[csz++] = it;
}
sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);
double add_res = 0;
for (size_t j=0; j<csz; ) {
item lower = c[cc[j++]];
used[lower.triangle_id] = true;
int cnt = 1;
while (cnt && j<csz) {
char & cur = used[c[cc[j++]].triangle_id];
cur = !cur;
if (cur) ++cnt; else --cnt;
}
item upper = c[cc[j-1]];
add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;
}
res += add_res * (x2 - x1) / 2;
}

cout.precision (8);
cout << fixed << res;
}

半平面交

  • 求凸多边形相交的面积
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
#include <stdio.h>
#include <math.h>
#include <algorithm>
using namespace std;
const double eps = 1e-8;
const int maxn = 20010;

int tail, head, it;
struct Point {
double x, y;
} P[maxn];
struct Line {
Point a, b;
double angle;
void getAngle() {angle = atan2(b.y-a.y, b.x-a.x);}
} L[maxn], deq[maxn];

int dcmp(double x) {
return x < -eps ? -1 : x > eps;
}
double xmult(Point a, Point b, Point c) {
return (a.x-c.x)*(b.y-c.y)-(a.y-c.y)*(b.x-c.x);
}
bool cmp(Line u, Line v) {
int d = dcmp(u.angle-v.angle);
if(d) return d > 0;
return dcmp(xmult(u.a, v.a, v.b)) > 0;
///Clockwise:大于0取向量左半部分为半平面,小于0,取右半部分
}
Point intersection(Line u, Line v)
{
Point ret = u.a;
double t = ((u.a.x-v.a.x)*(v.a.y-v.b.y)
-(u.a.y-v.a.y)*(v.a.x-v.b.x))
/ ((u.a.x-u.b.x)*(v.a.y-v.b.y)
-(u.a.y-u.b.y)*(v.a.x-v.b.x));
ret.x += (u.b.x-u.a.x)*t, ret.y += (u.b.y-u.a.y)*t;
return ret;
}
bool judge(Line l1, Line l2, Line l3) {
Point p = intersection(l2, l3);
return dcmp(xmult(p, l1.a, l1.b)) < 0;
///Clockwise:大于小于符号与上面cmp()中注释处相反
}
void HPI(Line L[], int n){
sort(L, L+n, cmp);
int tmp = 1;
for(int i = 1; i < n; ++i) {
if(dcmp(L[i].angle-L[tmp-1].angle) != 0) {
L[tmp++] = L[i];
}
}
n = tmp;
deq[0] = L[0], deq[1] = L[1];
head = 0, tail = 1;
for(int i = 2; i < n; ++i) {
while(head < tail && judge(L[i], deq[tail-1], deq[tail]))
tail--;
while(head < tail && judge(L[i], deq[head+1], deq[head]))
head++;
deq[++tail] = L[i];
}
while(head < tail && judge(deq[head], deq[tail-1], deq[tail]))
tail--;
while(head < tail && judge(deq[tail], deq[head+1], deq[head]))
head++;
if(head == tail) return ;
it = 0;
for(int i = head; i < tail; ++i) {
P[it++] = intersection(deq[i], deq[i+1]);
}
if(tail > head+1) {
P[it++] = intersection(deq[head], deq[tail]);
}
}
double getArea(Point p[], int n) {
double area = 0;
for(int i = 1; i < n-1; ++i) {
area += xmult(P[0], P[i], P[i+1]);
}
return fabs(area)/2.0;
}

int main()
{
int n;
while(~scanf("%d", &n)) {
n += 4;
L[0] = (Line){(Point){0, 10000}, (Point){0, 0}};
L[1] = (Line){(Point){10000, 10000}, (Point){0, 10000}};
L[2] = (Line){(Point){10000, 0}, (Point){10000, 10000}};
L[3] = (Line){(Point){0, 0}, (Point){10000, 0}};
L[0].getAngle(), L[1].getAngle(), L[2].getAngle(), L[3].getAngle();
for(int i = 4; i < n; ++i) {
scanf("%lf%lf%lf%lf", &L[i].a.x, &L[i].a.y, &L[i].b.x, &L[i].b.y);
L[i].getAngle();
}
HPI(L, n);
printf("%.1f\n", getArea(P, it));
}
return 0;
}

数据结构

笛卡尔树

笛卡尔树是一种同时满足二叉搜索树和堆的性质的数据结构。 可在一个数组上构造出来(时间复杂度可以达到O(n))。树中节点有几个属性, key(节点元素的大小)、index(节点在原数组中的索引)、left(左子节点)、right(右子节点)、parent(父节点)。

性质

  • 树中的元素满足二叉搜索树性质,要求按照中序遍历得到的序列为原数组序列
  • 树中节点满足堆性质,节点的key值要大于其左右子节点的key值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
void build(int n) {
int top = 0;
for(int i = 1; i <= n; ++i) {
l[i] = 0, r[i] = 0, p[i] = 0;
}
for(int i = 1; i <= n; ++i) {
int k = top;
while(k > 0 && a[stk[k-1]] < a[i]) --k;
if(k) r[stk[k-1]] = i;
if(k < top) l[i] = stk[k];
stk[k++] = i;
top = k;
}
for(int i = 1; i <= n; ++i) {
p[l[i]] = p[r[i]] = 1;
}
int rt = 0;
for(int i = 1; i <= n; ++i) {
if(p[i] == 0) rt = i;
}
}

解析表达式

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
bool delim (char c) {
return c == ' ';
}

bool is_op (char c) {
return c=='+' || c=='-' || c=='*' || c=='/' || c=='%';
}

int priority (char op) {
return
op == '+' || op == '-' ? 1 :
op == '*' || op == '/' || op == '%' ? 2 :
-1;
}

void process_op (vector<int> & st, char op) {
int r = st.back(); st.pop_back();
int l = st.back(); st.pop_back();
switch (op) {
case '+': st.push_back (l + r); break;
case '-': st.push_back (l - r); break;
case '*': st.push_back (l * r); break;
case '/': st.push_back (l / r); break;
case '%': st.push_back (l % r); break;
}
}

int calc (string & s) {
vector<int> st;
vector<char> op;
for (size_t i=0; i<s.length(); ++i)
if (!delim (s[i]))
if (s[i] == '(')
op.push_back ('(');
else if (s[i] == ')') {
while (op.back() != '(')
process_op (st, op.back()), op.pop_back();
op.pop_back();
}
else if (is_op (s[i])) {
char curop = s[i];
while (!op.empty() && priority(op.back()) >= priority(s[i]))
process_op (st, op.back()), op.pop_back();
op.push_back (curop);
}
else {
string operand;
while (i < s.length() && isalnum (s[i])))
operand += s[i++];
--i;
if (isdigit (operand[0]))
st.push_back (atoi (operand.c_str()));
else
st.push_back (get_variable_val (operand));
}
while (!op.empty())
process_op (st, op.back()), op.pop_back();
return st.back();
}

含一元操作码

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
bool delim (char c) {
return c == ' ';
}

bool is_op (char c) {
return c=='+' || c=='-' || c=='*' || c=='/' || c=='%';
}

int priority (char op) {
if (op < 0)
return 4; // op == -'+' || op == -'-'
return
op == '+' || op == '-' ? 1 :
op == '*' || op == '/' || op == '%' ? 2 :
-1;
}

void process_op (vector<int> & st, char op) {
if (op < 0) {
int l = st.back(); st.pop_back();
switch (-op) {
case '+': st.push_back (l); break;
case '-': st.push_back (-l); break;
}
}
else {
int r = st.back(); st.pop_back();
int l = st.back(); st.pop_back();
switch (op) {
case '+': st.push_back (l + r); break;
case '-': st.push_back (l - r); break;
case '*': st.push_back (l * r); break;
case '/': st.push_back (l / r); break;
case '%': st.push_back (l % r); break;
}
}
}

int calc (string & s) {
bool may_unary = true;
vector<int> st;
vector<char> op;
for (size_t i=0; i<s.length(); ++i)
if (!delim (s[i]))
if (s[i] == '(') {
op.push_back ('(');
may_unary = true;
}
else if (s[i] == ')') {
while (op.back() != '(')
process_op (st, op.back()), op.pop_back();
op.pop_back();
may_unary = false;
}
else if (is_op (s[i])) {
char curop = s[i];
if (may_unary && isunary (curop)) curop = -curop;
while (!op.empty() && (
curop >= 0 && priority(op.back()) >= priority(curop)
|| curop < 0 && priority(op.back()) > priority(curop))
)
process_op (st, op.back()), op.pop_back();
op.push_back (curop);
may_unary = true;
}
else {
string operand;
while (i < s.length() && isalnum (s[i])))
operand += s[i++];
--i;
if (isdigit (operand[0]))
st.push_back (atoi (operand.c_str()));
else
st.push_back (get_variable_val (operand));
may_unary = false;
}
while (!op.empty())
process_op (st, op.back()), op.pop_back();
return st.back();
}
//值得注意的是,在最简单的情况下,例如,当一元操作只被允许+和- ,
//右结合的戏剧没有作用,因此,在这种情况下,在方案无并发症无法进入。即 周期:

while (!op.empty() && (
curop >= 0 && priority(op.back()) >= priority(curop)
|| curop < 0 && priority(op.back()) > priority(curop))
)
process_op (st, op.back()), op.pop_back();

//可以替换为:

while (!op.empty() && priority(op.back()) >= priority(curop))
process_op (st, op.back()), op.pop_back();