数学 高斯消元 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 70 #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 74 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 * 1l l * a % p), --b; else a = int (a * 1l l * 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 * 1l l * 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 * 1l l * k % (n - 1 )), n) * 1l l * 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 * 1l l * a % p), --b; else a = int (a * 1l l * 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$,再分这两种情况对原固定的数进行分析。
题解:在近几年算法竞赛中,递推算法越来越重要:
第二类 stirling 数,显然拥有这样的性质:
而这些性质就可以总结成:
整数拆分 正整数n的有序k分拆的个数等于
$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) { 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) { 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 ); 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 124 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 103 #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 ; } 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 ; } 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 94 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 ; 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();