[关闭]
@oujunsong 2022-04-26T21:12:45.000000Z 字数 48386 阅读 100

各种代码模板

代码模板

卖价： 元

最大公约数 (gcd) [#nt01]:

  1. int gcd(int a, int b) {return b ? gcd(b, a % b) : a;}

快速幂 (递归版本) [#nt02]:

  1. ll PowerMod(ll a, ll n, ll m = mod) {
  2.     if (!n || a == 1) return 1ll;
  3.     ll s = PowerMod(a, n >> 1, m);
  4.     s = s * s % m;
  5.     return n & 1 ? s * a % m : s;
  6. }

快速幂 (非递归版本) [#nt03]:

  1. ll PowerMod(ll a, int n, ll c = 1) {for (; n; n >>= 1, a = a * a % mod) if (n & 1) c = c * a % mod; return c;}

扩展 Euclid 算法 (exgcd) [#nt04]:

  1. ll exgcd(ll a, ll b, ll &x, ll &y) {
  2.     if (b) {ll d = exgcd(b, a % b, y, x); return y -= a / b * x, d;}
  3.     else return x = 1, y = 0, a;
  4. }

基本预处理 [#nt05]:

  1. // 组合数
  2. // factorials, fact-inverses, unbounded
  3. void init() {
  4.     int i;
  5.     for (*fact = i = 1; i < N; ++i) fact[i] = (ll)fact[i - 1] * i % mod;
  6.     --i, finv[i] = PowerMod(fact[i], mod - 2);
  7.     for (; i; --i) finv[i - 1] = (ll)finv[i] * i % mod;
  8. }
  10. // factorials, fact-inverses, bounded
  11. void init(int n) {
  12.     int i;
  13.     for (*fact = i = 1; i <= n; ++i) fact[i] = (ll)fact[i - 1] * i % mod;
  14.     finv[n] = PowerMod(fact[n], mod - 2);
  15.     for (i = n; i; --i) finv[i - 1] = (ll)finv[i] * i % mod;
  16. }
  18. // inverses, factorials, fact-inverses, unbounded
  19. void init() {
  20.     int i;
  21.     for (inv[1] = 1, i = 2; i < N; ++i) inv[i] = ll(mod - mod / i) * inv[mod % i] % mod;
  22.     for (*finv = *fact = i = 1; i < N; ++i) fact[i] = (ll)fact[i - 1] * i % mod, finv[i] = (ll)finv[i - 1] * inv[i] % mod;
  23. }
  25. // inverses, factorials, fact-inverses, bounded
  26. void init(int n) {
  27.     int i;
  28.     for (inv[1] = 1, i = 2; i <= n; ++i) inv[i] = ll(mod - mod / i) * inv[mod % i] % mod;
  29.     for (*finv = *fact = i = 1; i <= n; ++i) fact[i] = (ll)fact[i - 1] * i % mod, finv[i] = (ll)finv[i - 1] * inv[i] % mod;
  30. }

线性筛素数 (筛最小素因子) [#nt06]:

  1. int pn = 0, c[N], p[78540];
  3. void sieve(int n) {
  4.     int i, j, v;
  5.     memset(c, -1, sizeof c);
  6.     for (i = 2; i <= n; ++i) {
  7.         if (!~c[i]) p[pn] = i, c[i] = pn++;
  8.         for (j = 0; (v = i * p[j]) <= n && j <= c[i]; ++j) c[v] = j;
  9.     }
  10. }

Miller-Rabin 素数测试 [#nt07]:

  1. // No pseudoprime below 2^64 to base 2, 3, 5, 7, 11, 13, 82, 373.
  2. const int test[8] = {2, 3, 5, 7, 11, 13, 82, 373};
  3. bool Miller_Rabin(ll n) {
  4.     if (n < MAX) return !c[n] && n > 1;
  5.     int c, i, j; ll s = n - 1, u, v;
  6.     for (c = 0; !(s & 1); ++c, s >>= 1);
  7.     for (i = 0; i < 8; ++i) {
  8.         if (!(n % test[i])) return false;
  9.         u = PowerMod(test[i], s, n);
  10.         for (j = 0; j < c; ++j, u = v)
  11.             if (v = MulMod(u, u, n), u != 1 && u != n - 1 && v == 1) return false;
  12.         if (u != 1) return false;
  13.     }
  14.     return true;
  15. }

分解质因数 [#nt08]:

  1. // version 1
  2. inline void push(ll prime, int alpha) {pr[cnt] = prime, al[cnt++] = alpha;}
  4. // version 2
  5. inline void push(ll prime, int alpha) {
  6.     map::iterator it = result.find(prime);
  7.     it == result.end() ? (void)result.emplace(prime, alpha) : (void)(it->second += alpha);
  8. }
  10. void factor(ll n) {
  11.     int i, j; cnt = 0; // result.clear();
  12.     for (i = 0; n > 1; ) {
  13.         if (n >= N) {
  14.             for (; (ll)p[i] * p[i] <= n && n % p[i]; ++i);
  15.             if (n % p[i]) return push(n, 1);
  16.         } else i = c[n];
  17.         for (j = 0; !(n % p[i]); n /= p[i], ++j); push(p[i], j);
  18.     }
  19. }

线性筛 Euler φ 函数 [#nt09]:

  1. void sieve(int n) {
  2.     int i, j, v; phi[1] = 1;
  3.     memset(c, -1, sizeof c);
  4.     for (i = 2; i <= n; ++i) {
  5.         if (!~c[i]) p[pn] = i, c[i] = pn++, phi[i] = i - 1;
  6.         for (j = 0; (v = i * p[j]) <= n && j < c[i]; ++j) c[v] = j, phi[v] = phi[i] * (p[j] - 1);
  7.         if (v <= n) c[v] = j, phi[v] = phi[i] * p[j];
  8.     }
  9. }

线性筛 Möbius μ 函数 [#nt10]:

  1. void sieve(int n) {
  2.     int i, j, v; mu[1] = 1;
  3.     memset(c, -1, sizeof c);
  4.     for (i = 2; i <= n; ++i) {
  5.         if (!~c[i]) p[pn] = i, c[i] = pn++, mu[i] = -1;
  6.         for (j = 0; (v = i * p[j]) <= n && j < c[i]; ++j) c[v] = j, mu[v] = -mu[i];
  7.         if (v <= n) c[v] = j;
  8.     }
  9. }

遍历 ⌊n / i⌋ 的所有可能值 (旧版整除分块) [#nt11]:

  1. #define next(n, i) (n / (n / i + 1))
  2. for(i = n; i >= 1; i = j){
  3.     j = next(n, i);
  4.     // do something
  5. }
  8. #define _next(n, i) (n / (n / (i + 1)))
  9. for(i = 0; i < n; i = j){
  10.     j = _next(n, i);
  11.     // do something
  12. }
  14. // 0 1 100 1 2 50 2 3 33 3 4 25 4 5 20 5 6 16 16 6 7 14 7 8 12 8 9 11 9 10 10 10 11 9 11 12 8 12 13 7 14 15 6 16 17 5 20 21 4 25 26 3 33 34 2 50 51 1 100 101 0 ?

新版整除分块 [#nt12]:

  1. // int32 ver.
  2. int n, cnt, srn;
  3. int v[G];
  5. inline int ID(int x) {return x <= srn ? x : cnt + 1 - n / x;}
  7. {
  8.     int i;
  9.     srn = sqrt(n), cnt = srn * 2 - (srn * (srn + 1) > n);
  10.     for (i = 1; i <= srn; ++i) v[i] = i, v[cnt - i + 1] = n / i;
  11. }
  13. // int64 ver.
  14. ll n, v[G];
  15. int cnt, srn;
  17. inline int ID(ll x) {return x <= srn ? x : cnt + 1 - n / x;}
  19. {
  20.     int i;
  21.     srn = sqrt(n), cnt = srn * 2 - (srn * (srn + 1ll) > n);
  22.     for (i = 1; i <= srn; ++i) v[i] = i, v[cnt - i + 1] = n / i;
  23. }

杜氏求和法 / 杜教筛 [#nt13]:

  1. void Dsum() {
  2.     int i, j, k, l, n, sq, _;
  3.     for (i = 1; i <= cnt && v[i] < N; ++i) Mu[i] = smu[v[i]], Phi[i] = sphi[v[i]];
  4.     for (; i <= cnt; ++i) {
  5.         int &m = Mu[i]; ll &p = Phi[i];
  6.         n = v[i], sq = sqrt(n), m = 1, p = n * (n + 1ll) / 2, l = n / (sq + 1);
  7.         for (j = sq; j; --j) {
  8.             _ = (k = n / j) - l, m -= _ * smu[j], p -= _ * sphi[j], l = k;
  9.             if (j != 1 && k > sq) m -= Mu[_ = ID(k)], p -= Phi[_];
  10.         }
  11.     }
  12. }

大小步离散对数 [#nt14]:

  1. int Boy_Step_Girl_Step(){
  2.     mii :: iterator it;
  3.     ll boy = 1, girl = 1, i, t;
  4.     for(i = 0; i < B; i++){
  5.         m.insert(pr((int)boy, i)); // if there exists in map, then ignore it
  6.         (boy *= a) %= p;
  7.     }
  8.     for(i = 0; i * B < p; i++){
  9.         t = b * inv(girl) % p;
  10.         it = m.find((int)t);
  11.         if(it != m.end())
  12.             return i * B + it->second;
  13.         (girl *= boy) %= p;
  14.     }
  15.     return -1;
  16. }
  18. // input p and then B = (int)sqrt(p - 1);
  19. // inv(x) is the multiplicative inverse of x in Z[p].
  20. // mii: map <int, int> and m is a variable of which.

Pollard-Rho 素因数分解 [#nt15]:

  1. ll Pollard_Rho(ll n, int c) {
  2.     ll i = 1, k = 2, x = rand() % (n - 1) + 1, y = x, p;
  3.     for (; k <= 131072; ) {
  4.         x = (MulMod(x, x, n) + c) % n;
  5.         p = gcd<ll>((ull)(y - x + n) % n, n);
  6.         if (p != 1 && p != n) return p;
  7.         if (++i == k) y = x, k <<= 1;
  8.     }
  9.     return 0;
  10. }
  12. inline void push(ll prime, int alpha) {
  13.     map::iterator it = result.find(prime);
  14.     it == result.end() ? (void)result.insert(pr(prime, alpha)) : (void)(it->second += alpha);
  15. }
  17. void factor(ll n) {
  18.     int i, j, k; ll p;
  19.     for (i = 0; n != 1; )
  20.         if (n >= MAX) {
  21.             if (Miller_Rabin(n)) {push(n, 1); return;}
  22.             for (k = 97; !(p = Pollard_Rho(n, k)); --k);
  23.             factor(p); n /= p;
  24.         } else {
  25.             i = (c[n] ? c[n] : n);
  26.             for (j = 0; !(n % i); n /= i, ++j);
  27.             push(i, j);
  28.         }
  29. }

字符串 Hash [#str01]:

  1. const ll mod[3] = {900000011, 998244353, 1000000007};
  2. const ll bas[3] = {4493, 8111, 8527};
  3. // you can choose your bases and modulos
  5. char s[S];
  6. ll pw[3][S], Hash[3][S];
  8. inline ll getHash(int id, int L, int R){ // str[L..R-1]
  9.     ll J = (Hash[id][R] - Hash[id][L] * pw[id][R - L]) % mod[id];
  10.     return J < 0 ? J + mod[id] : J;
  11. }
  13. // following is the pretreatment
  14. ll *f, *g;
  15. for(j = 0; j < 3; j++){
  16.     f = Hash[j]; f[0] = 0;
  17.     g = pw[j]; g[0] = 1;
  18.     for(i = 0; i < n; i++){
  19.         f[i + 1] = (f[i] * bas[j] + (s[i] - 'a')) % mod[j];
  20.         g[i + 1] = g[i] * bas[j] % mod[j];
  21.     }
  22. }

KMP 模式匹配 [#str02]:

  1. // pretreatment (border) 
  2. for (j = *f = -1, i = 1; i < n; f[++i] = ++j)
  3.     for (; ~j && s[j] != s[i]; j = f[j]);
  5. // KMP
  6. for (j = i = 0; i < _n; ++i) {
  7.     for (; ~j && s[j] != m[i]; j = f[j]);
  8.     if (++j == n) printf("%d ", i - s_n + 2);
  9. }

字典树 [#str03]:

  1. void append(char *s){
  2.     char *p = s;
  3.     int t = 1, id; // t = 0 is also OK
  4.     for(; *p; ++p){
  5.         id = *p - 97;
  6.         t = d[t][id] ? d[t][id] : (d[t][id] = ++V);
  7.     }
  8.     ++val[t];
  9. }
  11. // the process of matching is just like going on a DFA, so omit it.

Aho-Corasick 自动机 [#str04]:

  1. void build(){
  2.     int h, ta = 1, i, t, id;
  3.     que[0] = 1;
  4.     f[1] = 0;
  5.     for(h = 0; h < ta; ++h)
  6.         for(i = que[h], id = 0; id < 26; ++id){ // 26 is the size of char-set
  7.             t = (f[i] ? d[f[i]][id] : 1); // 1 or 0 depend on the root of trie
  8.             int &u = d[i][id];
  9.             if(!u) {u = t; continue;}
  10.             f[u] = t; val[u] |= val[t]; que[ta++] = u;
  11.             la[u] = (v[t] ? t : la[t]);
  12.         }
  13. }

后缀数组 (倍增构造) 与最长公共前缀 [#str05]:

  1. struct LCP {
  2.     int n, *sa, *rnk, *st[LN];
  3.     LCP () : n(0), sa(NULL), rnk(NULL) {memset(st, 0, sizeof st);}
  4.     ~LCP () {
  5.         if (sa) delete [] sa;
  6.         if (rnk) delete [] rnk;
  7.         for (int i = 0; i < LN; ++i) if (st[i]) delete [] st[i];
  8.     }
  9.     void construct(const char *s) {
  10.         int i, j, k, m = 256, p, limit; n = strlen(s);
  11.         int *x = new int[n + 1], *y = new int[n + 1], *buf = new int[max(n, m)], *f, *g = new int[n + 1];
  12.         auto cmp = [this] (const int *a, const int u, const int v, const int l) {return a[u] == a[v] && (u + l >= n ? 0 : a[u + l]) == (v + l >= n ? 0 : a[v + l]);};
  13.         if (sa) delete [] sa; sa = new int[n];
  14.         if (rnk) delete [] rnk;
  15.         for (i = 0; i < LN; ++i) if (st[i]) delete [] st[i], st[i] = 0;
  16.         for (i = 0; i < n; ++i) sa[i] = i, x[i] = (unsigned char)s[i];
  17.         std::sort(sa, sa + n, [s] (const int u, const int v) {return (unsigned char)s[u] < (unsigned char)s[v];});
  18.         for (j = 1; j < n; j <<= 1, m = ++p) {
  19.             std::iota(y, y + j, n - j), p = j;
  20.             for (i = 0; i < n; ++i) if (sa[i] >= j) y[p++] = sa[i] - j;
  21.             memset(buf, 0, m << 2);
  22.             for (i = 0; i < n; ++i) ++buf[ x[y[i]] ];
  23.             for (i = 1; i < m; ++i) buf[i] += buf[i - 1];
  24.             for (i = n - 1; i >= 0; --i) sa[ --buf[ x[y[i]] ] ] = y[i];
  25.             std::swap(x, y);
  26.             x[*sa] = p = 1, x[n] = 0;
  27.             for (i = 1; i < n; ++i) x[sa[i]] = (cmp(y, sa[i - 1], sa[i], j) ? p : ++p);
  28.             if (p >= n) break;
  29.         }
  30.         if (rnk = x, n == 1) *x = 0;
  31.         else for (i = 0; i < n; ++i) --x[i];
  32.         delete [] buf, delete [] y;
  33.         for (p = i = 0; i < n; ++i) {
  34.             if (p && --p, !x[i]) continue;
  35.             for (j = sa[x[i] - 1], limit = n - max(i, j); p < limit && s[i + p] == s[j + p]; ++p);
  36.             g[ x[i] - 1 ] = p;
  37.         }
  38.         *st = g, k = n - 1;
  39.         for (j = 0; 1 << (j + 1) < n; ++j) {
  40.             k -= 1 << j, f = g, g = st[j + 1] = new int[k + 1];
  41.             for (i = 0; i < k; ++i)
  42.                 g[i] = min(f[i], f[i + (1 << j)]);
  43.         }
  44.     }
  45.     inline int operator () (const int u, const int v) {
  46.         assert((unsigned)u < (unsigned)n && (unsigned)v < (unsigned)n);
  47.         if (u == v) return n - u;
  48.         int L, R, c; std::tie(L, R) = std::minmax(rnk[u], rnk[v]), c = lg2(R - L);
  49.         return min(st[c][L], st[c][R - (1 << c)]);
  50.     }
  51. };

后缀自动机 [#str06]:

  1. #define q d[p][x]
  2. int extend(int x) {
  3.     for (p = np, val[np = ++cnt] = val[p] + 1; p && !q; q = np, p = pa[p]);
  4.     if (!p) pa[np] = 1;
  5.     else if (val[p] + 1 == val[q]) pa[np] = q;
  6.     else {
  7.         int nq = ++cnt;
  8.         val[nq] = val[p] + 1, memcpy(d[nq], d[q], 104);
  9.         pa[nq] = pa[q], pa[np] = pa[q] = nq;
  10.         for (int Q = q; p && q == Q; q = nq, p = pa[p]);
  11.     }
  12.     return f[np] = 1, np;
  13. }
  14. #undef q

根据后缀自动机构造后缀树 [#str07]:

  1. void build() {
  2.     int i, j, c; used[1] = true;
  3.     for (i = cnt; i; --i) if (~pos[i])
  4.         for (j = i; !used[j]; j = pa[j])
  5.             c = pos[j] - val[pa[j]], pos[pa[j]] = pos[j],
  6.             child[pa[j]][int(s[c])] = j, used[j] = true;
  7. //  dfs(1);
  8. }

Z 算法 [#str08]:

  1. void Z() {
  2.     int i, Max = 0, M = 0;
  3.     for (i = 1; i < n; ++i) {
  4.         z[i] = (i < Max ? std::min(z[i - M], Max - i) : 0);
  5.         for (; s[z[i]] == s[i + z[i]]; ++z[i]);
  6.         if (i + z[i] > Max) Max = i + z[i], M = i;
  7.     }
  8. }

Manacher 回文串 [#str09]:

  1. void Manacher() {
  2.     int n = 2, i, Max = 0, M = 0;
  3.     t[0] = 2, t[1] = 1;
  4.     for (i = 0; i < S; ++i) t[n++] = s[i], t[n++] = 1; t[n++] = 3;
  5.     for (i = 0; i < n; ++i){
  6.         f[i] = (i < Max ? std::min(f[M * 2 - i], Max - i) : 1);
  7.         for(; t[i + f[i]] == t[i - f[i]]; ++f[i]);
  8.         if (i + f[i] > Max) Max = i + f[i], M = i;
  9.     }
  10. }

回文自动机 [#str10]:

  1. int get_fail(int x) {for (; ptr[~val[x]] != *ptr; x = fail[x]); return x;}
  3. int extend(int x) {
  4.     int &q = d[p = get_fail(p)][x];
  5.     if (!q) fail[++cnt] = d[get_fail(fail[p])][x], val[q = cnt] = val[p] + 2;
  6.     return p = q;
  7. }
  9. val[1] = -1, p = 0, *fail = cnt = 1;

多串后缀自动机 / 广义后缀自动机 [#str11]:

  1. #define q d[p][x]
  2. #define try_split(v) { \
  3.     if (val[p] + 1 == val[q]) v = q; \
  4.     else { \
  5.         int nq = ++cnt; \
  6.         val[nq] = val[p] + 1, memcpy(d[nq], d[q], 104); \
  7.         pa[nq] = pa[q], v = pa[q] = nq; \
  8.         for (int Q = q; p && q == Q; q = nq, p = pa[p]); \
  9.     } \
  10. }
  12. int extend(int x) {
  13.     if (p = np, q) try_split(np)
  14.     else {
  15.         for (val[np = ++cnt] = val[p] + 1; p && !q; q = np, p = pa[p]);
  16.         if (p) try_split(pa[np]) else pa[np] = 1;
  17.     }
  18.     return np;
  19. }
  20. #undef q

深度优先搜索，dfs 序，树上基本操作 [#tr01]:

  1. void dfs(int x) {
  2.     int i, y; o[++cnt] = x, id[x] = cnt, size[x] = 1;
  3.     for (i = first[x]; i; i = next[i])
  4.         if ((y = to[i]) != p[x])
  5.             p[y] = x, dep[y] = dep[x] + 1, dfs(y), size[x] += size[y];
  6.     eid[x] = cnt;
  7. }

最近公共祖先 (LCA) 的倍增算法 (树上倍增的基本模板) [#tr02]:

  1. void dfs(int x) {
  2.     int i, y;
  3.     for (i = 0; P[i][x]; ++i) P[i + 1][x] = P[i][P[i][x]];
  4.     for (i = first[x]; i; i = next[i])
  5.         if ((y = to[i]) != p[x])
  6.             p[y] = x, dep[y] = dep[x] + 1, dfs(y);
  7. }
  9. int jump_until(int x, int d) {
  10.     for (int S = dep[x] - d; S; S &= S - 1) x = P[ctz(S)][x];
  11.     return x;
  12. }
  14. int LCA(int x, int y) {
  15.     if (dep[x] < dep[y]) std::swap(x, y);
  16.     if (x = jump_until(x, dep[y]), x == y) return x;
  17.     for (int i = LN - 1; i >= 0; --i)
  18.         if (P[i][x] != P[i][y])
  19.             x = P[i][x], y = P[i][y];
  20.     return p[x];
  21. }

最近公共祖先 (LCA) 的 ST 表算法 [#tr03]:

  1. int cnt = 0, id[N], o[N], st[LN][N];
  3. void dfs(int x, int px = 0) {
  4.     int i, y; st[0][cnt] = id[px], o[++cnt] = x, id[x] = cnt;
  5.     for (i = first[x]; i; i = next[i])
  6.         if ((y = to[i]) != px) dfs(y, x);
  7. }
  9. void build_sparse_table() {
  10.     int *f, *g = *st, i, j, k = n - 1;
  11.     for (j = 0; 1 << (j + 1) < n; ++j) {
  12.         f = g, g = st[j + 1], k -= 1 << j;
  13.         for (i = 1; i <= k; ++i)
  14.             g[i] = min(f[i], f[i + (1 << j)]);
  15.     }
  16. }
  18. inline int LCA_relabel(int x, int y) {
  19.     if (x == y) return x;
  20.     if (x > y) std::swap(x, y);
  21.     int c = lg2(y - x);
  22.     return min(st[c][x], st[c][y - (1 << c)]);
  23. }
  25. inline int LCA(int x, int y) {
  26.     if (x == y) return x;
  27.     int L, R, c; std::tie(L, R) = std::minmax(id[x], id[y]), c = lg2(R - L);
  28.     return o[min(st[c][L], st[c][R - (1 << c)])];
  29. }

轻重链剖分 (Heavy-Light Decomposition) [#tr04]:

  1. int p[N], dep[N], size[N];
  2. int cnt = 0, id[N], prf[N], top[N];
  4. void dfs_wt(int x) {
  5.     int i, y, &z = prf[x]; size[x] = !next[first[x]];
  6.     for (i = first[x]; i; i = next[i])
  7.         if ((y = to[i]) != p[x]) {
  8.             p[y] = x, dep[y] = dep[x] + 1;
  9.             dfs_wt(y), size[x] += size[y];
  10.             size[y] > size[z] ? z = y : 0;
  11.         }
  12. }
  14. void dfs_hld(int x, int r) {
  15.     int i, y; id[x] = ++cnt, top[x] = r;
  16.     if (!prf[x]) return;
  17.     dfs_hld(prf[x], r);
  18.     for (i = first[x]; i; i = next[i])
  19.         if (!top[y = to[i]]) dfs_hld(y, y);
  20. }
  22. int solve(int u, int v) {
  23.     int x = top[u], y = top[v];
  24.     for (; x != y; u = p[x], x = top[u]) {
  25.         if (dep[x] < dep[y]) {swap(u, v); swap(x, y);}
  26.         // do something on [x, u]
  27.     }
  28.     if (dep[u] > dep[v]) {swap(u, v); swap(lx, ly);}
  29.     // do something on [u, v]
  30.     // return something
  31. }

树分治 (静态点分治) [#tr05]:

  1. namespace Centroid {
  2.     int V, Gm, G, size[N];
  4.     void init(int _V) {V = _V, Gm = INT_MAX;}
  6.     int get(int x, int px = 0) {
  7.         int i, y, M = 0; size[x] = 1;
  8.         for (i = first[x]; i; i = next[i])
  9.             if ((y = to[i]) != px && !banned[y])
  10.                 get(y, x), up(M, size[y]), size[x] += size[y];
  11.         return up(M, V - size[x]), M <= Gm ? (Gm = M, G = x) : G;
  12.     }
  13. }
  15. #define get_centroid(x, y) (Centroid::init(y), Centroid::get(x))
  17. void dfs(int x, int px = 0, int dep = 1) {
  18.     int i, y; size[x] = 1;
  19.     for (i = first[x]; i; i = next[i])
  20.         if ((y = to[i]) != px && !banned[y])
  21.             dfs(y, x, dep + 1), size[x] += size[y];
  22. }
  24. void solve(int x) {
  25.     int i, y, G;
  26.     banned[x] = true, dfs(x);
  27.     // do something
  28.     for (i = first[x]; i; i = next[i])
  29.         if (!fy[y = to[i]])
  30.             G = get_centroid(y, size[y]), solve(G);
  31. }
  33. // main
  34. G = get_centroid(1, n), solve(G);

虚树 (深度栈算法) [#tr06]:

  1. int cnt_vir, vir[N], virp[N];
  3. inline bool idcmp(const int x, const int y) {return id[x] < id[y];}
  5. void DSA(int n, int *v) {
  6.     #define ins(x) (virp[x] = stack[top], stack[++top] = vir[cnt_vir++] = x)
  7.     int i, x, y, top = 0; cnt_vir = 0;
  8.     for (i = 0; i < n; ++i)
  9.         if (x = v[i], !top) {
  10.             for (; top; --top) stack[top] = 0; ins(x);
  11.         } else {
  12.             stack[top + 1] = 0;
  13.             for (y = LCA(x, stack[top]); dep[ stack[top] ] > dep[y]; --top);
  14.             virp[ stack[top + 1] ] = y;
  15.             if (stack[top] != y) ins(y); ins(x);
  16.         }
  17.     for (; top; --top) stack[top] = 0;
  18.     std::sort(vir, vir + cnt_vir, idcmp);
  19. }

长链剖分 (Long-Short Decomposition) [#tr07]:

  1. int p[N], dep[N], f[N];
  2. int cnt = 0, id[N], prf[N], top[N];
  3. int len[N], buf[N * 3], *alloc = buf, *near[N];
  5. void dfs_len(int x) {
  6.     int i, y, &z = prf[x];
  7.     for (i = 0; i < LN && P[i][x]; ++i) P[i + 1][x] = P[i][P[i][x]];
  8.     for (i = first[x]; i; i = next[i])
  9.         if ((y = to[i]) != p[x]) {
  10.             p[y] = x, dep[y] = dep[x] + 1;
  11.             dfs_len(y), up(f[x], f[y] + 1);
  12.             f[y] > f[z] ? z = y : 0;
  13.         }
  14. }
  16. void dfs_lsd(int x, int r, int l = 1) { // long-short decomposition
  17.     int i, y; id[x] = ++cnt, top[x] = r;
  18.     if (!prf[x]) {len[x] = l; return;}
  19.     dfs_lsd(prf[x], r, l + 1), len[x] = len[prf[x]];
  20.     for (i = first[x]; i; i = next[i])
  21.         if (!top[y = to[i]]) dfs_lsd(y, y);
  22. }
  24. int ancestor(int x, int k) {
  25.     if (dep[x] < k) return 0;
  26.     if (!k) return x;
  27.     x = P[i][x], near[top[x]][dep[x] - dep[top[x]] - (k - (1 << lg2(k)))];
  28. }
  30. void init() {
  31.     int i, j, x, *pt;
  32.     *f = *dep = -1, dfs_len(1), dfs_lsd(1, 1);
  33.     for (i = 1; i < N; ++i)
  34.         if (top[i] == i) {
  35.             pt = alloc + len[i], alloc += len[i] * 2 + 1, *pt = i;
  36.             for (j = 1, x = p[i]; j <= len[i] && x; ++j, x = p[x]) pt[-j] = x;
  37.             for (j = 1, x = prf[i]; j <= len[i] && x; ++j, x = prf[x]) pt[j] = x;
  38.         }
  39. }
  40. Trie Tree [#tr08]:
  41. int son[N][26], cnt[N], idx;
  43. void insert(char *str)
  44. {
  45.     int p = 0;
  46.     for (int i = 0; str[i]; i ++ )
  47.     {
  48.         int u = str[i] - 'a';
  49.         if (!son[p][u]) son[p][u] = ++ idx;
  50.         p = son[p][u];
  51.     }
  52.     cnt[p] ++ ;
  53. }
  55. int query(char *str)
  56. {
  57.     int p = 0;
  58.     for (int i = 0; str[i]; i ++ )
  59.     {
  60.         int u = str[i] - 'a';
  61.         if (!son[p][u]) return 0;
  62.         p = son[p][u];
  63.     }
  64.     return cnt[p];
  65. }

双向邻接表 [#gr01]:

  1. inline void addedge(int u, int v) {
  2.     to[++E] = v, next[E] = first[u], first[u] = E;
  3.     to[++E] = u, next[E] = first[v], first[v] = E;
  4. }

有向图的强连通分量 [#gr02]:

  1. void dfs(int x) {
  2.     int i, y;
  3.     id[x] = low[x] = ++cnt, instack[x] = true, sta[stc++] = x;
  4.     for(i = first[x]; i; i = next[i])
  5.         if (!id[y = to[i]])
  6.             dfs(y), down(low[x], low[y]);
  7.         else if (instack[y])
  8.             down(low[x], id[y]);
  9.     if (id[x] == low[x])
  10.         for (y = 0; y != x; )
  11.             y = sta[--stc], instack[y] = false, top[y] = x;
  12. }

有向无环图的拓扑排序 [#gr03]:

  1. void toposort() {
  2.     int i, h, t = 0, x, y;
  3.     for (i = 1; i <= V; ++i) if (!deg[i]) que[t++] = i;
  4.     for (h = 0; h < t; ++h)
  5.         for (i = first[x = que[h]]; i; i = next[i])
  6.             if (!--deg[y = to[i]]) que[t++] = y;
  7. }

无向图的桥边 [#gr04]:

  1. void dfs(int x, int px = 0) {
  2.     int i, y;
  3.     id[x] = low[x] = ++cnt;
  4.     for (i = first[x]; i; i = next[i])
  5.         if (!id[y = to[i]]) {
  6.             dfs(y, i), down(low[x], low[y]);
  7.             if (id[x] < low[y]) // edge i is a bridge
  8.         } else if (px - 1 ^ i - 1 ^ 1)
  9.             down(low[x], id[y]);
  10. }

无向图的割点 [#gr05]:

  1. void dfs(int x, int px = 0) {
  2.     int i, y, c = 0; bool is_cut = false;
  3.     id[x] = low[x] = ++cnt;
  4.     for (i = first[x]; i; i = next[i])
  5.         if (!id[y = to[i]]) {
  6.             dfs(y, i), down(low[x], low[y]), ++c;
  7.             if (id[x] <= low[y]) is_cut = true;
  8.         } else if (px - 1 ^ i - 1 ^ 1)
  9.             down(low[x], id[y]);
  10.     if (!px) is_cut = c > 1;
  11. }

无向图的边双缩点 [#gr06]:

  1. void dfs(int x, int px = 0) {
  2.     int i, y;
  3.     id[x] = low[x] = ++cnt, stack[stc++] = x;
  4.     for (i = first[x]; i; i = next[i])
  5.         if (!id[y = to[i]]) {
  6.             dfs(y, i), down(low[x], low[y]);
  7.             if (id[x] < low[y]) bridge[ad(i)] = bridge[i] = true;
  8.         } else if (px - 1 ^ i - 1 ^ 1)
  9.             down(low[x], id[y]);
  10.     if (id[x] == low[x])
  11.         for (y = 0; y != x; )
  12.             y = stack[--stc], top[y] = x;
  13. }

无向图的点双缩点 / 圆方树 [#gr07]:

  1. void dfs(int x, int px = 0) {
  2.     int i, y, z;
  3.     id[x] = low[x] = ++cnt, stack[top++] = x;
  4.     for (i = first[x]; i; i = next[i])
  5.         if (!id[y = to[i]]) {
  6.             dfs(y, i), down(low[x], low[y]);
  7.             if (id[x] <= low[y])
  8.                 for (link(++bcc_cnt + V, x), z = 0; z != y; )
  9.                     link(z = stack[--top], bcc_cnt + V);
  10.         } else if (px - 1 ^ i - 1 ^ 1)
  11.             down(low[x], id[y]);
  12. }

带权 edge 结构体 [#gr08]:

  1. struct edge {
  2.     int u, v, w;
  3.     edge (int u0 = 0, int v0 = 0, int w0 = 0) : u(u0), v(v0), w(w0) {}
  4. };

单源最短路 (Dijkstra) [#gr09]:

  1. typedef std::pair <int, int> pr;
  3. int d[N];
  4. std::priority_queue <pr, std::vector <pr>, std::greater <pr>> pq;
  6. void Dijkstra(int s) {
  7.     int i, x, y, D;
  8.     memset(d, 63, (V + 1) << 2);
  9.     for (pq.emplace(d[s] = 0, s); !pq.empty(); ) {
  10.         std::tie(D, x) = pq.top(), pq.pop();
  11.         if (d[x] != D) continue;
  12.         for (i = first[x]; i; i = next[i])
  13.             if (D + e[i].w < d[y = e[i].v])
  14.                 pq.emplace(d[y] = D + e[i].w, y);
  15.     }
  16. }

所有点对的最短路 (Floyd) [#gr10]:

  1. void Floyd() {
  2.     int i, j, k;
  3.     for (k = 1; k <= V; ++k)
  4.         for (i = 1; i <= V; ++i)
  5.             for (j = 1; j <= V; ++j)
  6.                 down(d[i][j], d[i][k] + d[k][j]);
  7. }

最小生成树 (Kruskal) [#gr11]:

  1. struct edge {
  2.     int u, v, w;
  3.     edge (int u0 = 0, int v0 = 0, int w0 = 0): u(u0), v(v0), w(w0) {}
  4.     bool operator < (const edge &b) const {return w < b.w;}
  5. };
  7. void Kruskal(){
  8.     uf.resize(V);
  9.     sort(e + 1, e + (E + 1));
  10.     for(i = 1; i <= E; i++)
  11.         if(!uf.test(e[i].u, e[i].v, true)){
  12.             // e is an edge of minimum spanning tree
  13.             if(++Es >= V - 1) return;
  14.         }
  15. }

网络流 edge 结构体 [#gr12]:

  1. struct edge {
  2.     int u, v, f; // f is remaining flow
  3.     edge(int u0 = 0, int v0 = 0, int f0 = 0): u(u0), v(v0), f(f0) {}
  4. };

最大流 (Dinic) [#gr13]:

  1. namespace Flow {
  2.     #define ad(x) ((x - 1 ^ 1) + 1)
  4.     const int N = 2000, M = 100000;
  6.     struct edge {
  7.         int u, v, f;
  8.         edge (int u0 = 0, int v0 = 0, int f0 = 0) : u(u0), v(v0), f(f0) {}
  9.     } e[M];
  11.     int V = 2, E = 0, si = 1, ti = 2, flow;
  12.     int first[N], next[M];
  13.     int dep[N], cur[N], que[N];
  15.     inline void addedge(int u, int v, int f) {
  16.         e[++E] = edge(u, v, f), next[E] = first[u], first[u] = E;
  17.         e[++E] = edge(v, u), next[E] = first[v], first[v] = E;
  18.     }
  20.     bool bfs() {
  21.         int h, t = 1, i, x, y;
  22.         memset(dep, -1, sizeof dep);
  23.         que[0] = si, dep[si] = 0;
  24.         for (h = 0; h < t; ++h) {
  25.             if ((x = que[h]) == ti) return true;
  26.             for (i = first[x]; i; i = next[i])
  27.                 if (dep[y = e[i].v] == -1 && e[i].f)
  28.                     que[t++] = y, dep[y] = dep[x] + 1;
  29.         }
  30.         return false;
  31.     }
  33.     int dfs(int x, int lim) {
  34.         int a, c, f = 0;
  35.         if (x == ti || !lim) return lim;
  36.         for (int &i = cur[x]; i; i = next[i])
  37.             if (dep[e[i].v] == dep[x] + 1 && e[i].f) {
  38.                 a = std::min(lim - f, e[i].f);
  39.                 c = dfs(e[i].v, a);
  40.                 e[i].f -= c, e[ad(i)].f += c;
  41.                 if ((f += c) == lim) return f;
  42.             }
  43.         return f;
  44.     }
  46.     int Dinic() {
  47.         for (flow = 0; bfs(); flow += dfs(si, INT_MAX))
  48.             memcpy(cur, first, sizeof cur);
  49.         return flow;
  50.     }
  51. }

最小费用最大流 (Dinic) [#gr14]:

  1. namespace CF {
  2.     #define ad(x) ((x - 1 ^ 1) + 1)
  4.     const int N = 2000, M = 100000, INF = 0x7f7f7f7f;
  6.     struct edge {
  7.         int u, v, c, f;
  8.         edge (int u0 = 0, int v0 = 0, int c0 = 0, int f0 = 0) : u(u0), v(v0), c(c0), f(f0) {}
  9.     } e[M];
  11.     int V = 2, E = 0, si = 1, ti = 2, flow, cost;
  12.     int first[N], next[M];
  13.     int dep[N], cur[N], que[M << 1];
  14.     bool in_que[N], used[N];
  16.     inline void addedge(int u, int v, int c, int f) {
  17.         e[++E] = edge(u, v, c, f), next[E] = first[u], first[u] = E;
  18.         e[++E] = edge(v, u, -c), next[E] = first[v], first[v] = E;
  19.     }
  21.     bool bfs() {
  22.         int h = M, t = h + 1, i, x, y;
  23.         memset(dep, 127, sizeof dep);
  24.         que[h] = ti, dep[ti] = 0, in_que[ti] = true;
  25.         for (; h < t; ) {
  26.             x = que[h++], in_que[x] = false;
  27.             for (i = first[x]; i; i = next[i])
  28.                 if (dep[y = e[i].v] > dep[x] - e[i].c && e[ad(i)].f) {
  29.                     dep[y] = dep[x] - e[i].c;
  30.                     if (!in_que[y])
  31.                         in_que[y] = true, (dep[y] >= dep[que[h]] ? que[t++] : que[--h]) = y;
  32.                 }
  33.         }
  34.         return dep[si] < INF;
  35.     }
  37.     int dfs(int x, int lim) {
  38.         int a, c, f = 0;
  39.         if (x == ti || !lim) return lim;
  40.         used[x] = true;
  41.         for (int &i = cur[x]; i; i = next[i])
  42.             if (dep[e[i].v] == dep[x] - e[i].c && e[i].f && !used[e[i].v]) {
  43.                 a = std::min(lim - f, e[i].f);
  44.                 c = dfs(e[i].v, a);
  45.                 e[i].f -= c, e[ad(i)].f += c;
  46.                 if ((f += c) == lim) return f;
  47.             }
  48.         return f;
  49.     }
  51.     int Dinic() {
  52.         int f;
  53.         for (cost = flow = 0; bfs(); ) {
  54.             memcpy(cur, first, sizeof cur);
  55.             memset(used, 0, sizeof used);
  56.             flow += f = dfs(si, INF);
  57.             cost += dep[si] * f;
  58.         }
  59.         return cost;
  60.     }
  61. }

二分图最大匹配 (增广路算法) [#gr15]:

  1. bool dfs(int x){
  2.     for(int i = 1; i <= n_g; i++)
  3.         if(e[x][i] && !used[i]){
  4.             used[i] = 1;
  5.             if(!boy[i] || dfs(boy[i])){
  6.                 boy[i] = x; girl[x] = i; return true;
  7.             }
  8.         }
  9.     return false;
  10. }

一般图最大匹配 (带花树算法) [#gr16]:

  1. #define unknown -1
  2. #define boy 0
  3. #define girl 1
  5. int LCA(int x, int y) {
  6.     for (++hash_cnt; x; y ? swap(x, y) : (void)0) {
  7.         x = ancestor(x);
  8.         if (hash[x] == hash_cnt) return x; // visited
  9.         hash[x] = hash_cnt;
  10.         x = prev[match[x]];
  11.     }
  12.     return 0x131a371;
  13. }
  15. void blossom(int x, int y, int root, int &t) { // vertices in blossoms are all boys !
  16.     for (int z; ancestor(x) != root; y = z, x = prev[y]) {
  17.         prev[x] = y; z = match[x];
  18.         if (col[z] == girl) que[t++] = z, col[z] = boy;
  19.         if (ancestor(x) == x) p[x] = root;
  20.         if (ancestor(z) == z) p[z] = root;
  21.     }
  22. }
  24. bool bfs(int st) {
  25.     int h, t = 1, i, x, y, b, g;
  26.     que[0] = st; col[st] = boy;
  27.     for (h = 0; h < t; ++h)
  28.         for (i = first[x = que[h]]; i; i = next[i])
  29.             if (col[y = to[i]] == unknown) { // common step
  30.                 prev[y] = x; col[y] = girl;
  31.                 if (!match[y]) { // augment (change mates) !!!
  32.                     for (g = y; g; swap(match[b], g))
  33.                         match[g] = b = prev[g];
  34.                     return true;
  35.                 }
  36.                 col[que[t++] = match[y]] = boy;
  37.             } else if(col[y] == boy && ancestor(x) != ancestor(y)) { // flower !!!
  38.                 b = LCA(x, y); blossom(x, y, b, t); blossom(y, x, b, t);
  39.             }
  40.     return false;
  41. }
  43. // main process
  44. for (i = 1; i <= V; ++i) {
  45.     for (v = 1; v <= V; ++v) col[v] = unknown, p[v] = v;
  46.     if (!match[i] && bfs(i)) ++ans;
  47. }

最小树形图 (朱-刘-Edmonds) [#gr17]:

  1. namespace DMST {
  2.     const int N = 108, M = 10054, INF = 0x3f3f3f3f;
  4.     int V, E, r;
  5.     int p[N], best[N];
  6.     int bel[N], used[N], cyc[N];
  7.     int min[N], delta[M];
  8.     bool ans[M];
  9.     edge e[M];
  10.     pr elim[N * M];
  12.     inline void reset(int n, int root) {V = n, E = 0, r = root;}
  14.     inline void link(int u, int v, int w) {e[E++] = edge(u, v, w);}
  16.     int CLE() {
  17.         int i, u, v, w, n, r, cnt, tot = 0, w0 = 0, w1 = 0;
  18.         n = V, r = DMST::r, std::iota(bel, bel + (V + 1), 0);
  19.         memset(ans, false, E), memset(delta, 0, E << 2);
  20.         for (; ; n = cnt, r = cyc[r]) {
  21.             memset(min, 63, (V + 1) << 2),
  22.             memset(used, 0, (V + 1) << 2), used[r] = -1;
  23.             memset(cyc, cnt = 0, (V + 1) << 2);
  24.             for (i = 0; i < E; ++i) {
  25.                 u = bel[ e[i].u ], v = bel[ e[i].v ], w = e[i].w - delta[i];
  26.                 if (u != v && w < min[v]) p[v] = u, min[v] = w, best[v] = i;
  27.             }
  28.             for (i = 1; i <= n; ++i) if (i != r) {
  29.                 if (min[i] >= INF) return -1;
  30.                 for (u = i; !(u == r || used[u]); u = p[u]) used[u] = i;
  31.                 if (used[u] == i) {
  32.                     ++cnt, v = u;
  33.                     do cyc[v] = cnt, v = p[v]; while (v != u);
  34.                 }
  35.             }
  36.             if (!(w = cnt)) break;
  37.             for (i = 1; i <= n; ++i)
  38.                 if (cyc[i]) ans[ best[i] ] = true;
  39.                 else cyc[i] = ++cnt;
  40.             for (i = 0; i < E; ++i) {
  41.                 u = bel[ e[i].u ], v = bel[ e[i].v ];
  42.                 if (cyc[u] != cyc[v])
  43.                     if (delta[i] += min[v], cyc[v] <= w)
  44.                         elim[tot++] = pr(best[v], i);
  45.             }
  46.             for (i = 1; i <= V; ++i) bel[i] = cyc[bel[i]];
  47.         }
  48.         for (i = 1; i <= n; ++i) if (i != r) ans[ best[i] ] = true;
  49.         for (i = tot - 1; i >= 0; --i) if (ans[elim[i].second]) ans[elim[i].first] = false;
  50.         for (i = 0; i < E; ++i) if (ans[i]) ++w0, w1 += e[i].w;
  51.         return assert(w0 == V - 1), w1;
  52.     }
  53. }

并查集 (不封装，按秩合并) [#ds01]:

  1. int ancestor(int x) {return p[x] == x ? x : (p[x] = ancestor(p[x]));}
  3. bool test(int x, int y, bool un = false) {
  4.     if ((x = ancestor(x)) == (y = ancestor(y))) return true;
  5.     if (un) size[x] > size[y] ? std::swap(x, y) : (void)0, p[x] = y, size[y] += size[x];
  6.     return false;
  7. }

并查集 (封装) [#ds02]:

  1. struct UFind{
  2.     int sz, *p;
  3.     UFind (): sz(0) {p = NULL;}
  4.     ~UFind () {if(p) delete [] (p);}
  5.     void resize(int size){
  6.         if(p) delete [] (p); p = new int[(sz = size) + 1];
  7.         for(int i = 0; i <= sz; i++) p[i] = i;
  8.     }
  9.     int ancestor(int x) {return x == p[x] ? x : p[x] = ancestor(p[x]);}
  10.     bool test(int x, int y, bool un = false){
  11.         if((x = ancestor(x)) == (y = ancestor(y))) return true;
  12.         if(un) p[x] = y; return false;
  13.     }
  14. };

树状数组 (不封装) [#ds03]:

  1. #define lowbit(x) (x & -x)
  3. int sum(int h){
  4.     int s = 0;
  5.     while(h){
  6.         s += a[h];
  7.         h -= lowbit(h);
  8.     }
  9.     return s;
  10. }
  12. int add(int h, int v){
  13.     while(h <= n){
  14.         a[h] += v;
  15.         h += lowbit(h);
  16.     }
  17.     return v;
  18. }

树状数组 (封装，可清空) [#ds04]:

  1. struct BIT {
  2.     #define lowbit(x) (x & -x)
  3.     int n, ti, tag[N]; ll x[N];
  4.     BIT () : ti(0) {}
  5.     inline void resize(int size) {n = size;}
  6.     inline void clear() {++ti;}
  7.     inline ll & recover(int id) {return tag[id] == ti ? x[id] : (tag[id] = ti, x[id] = 0);}
  8.     ll sum(int h) {ll s = 0; for (; h > 0; h -= lowbit(h)) s += recover(h); return s;}
  9.     void add(int h, ll v) {for (; h <= n; h += lowbit(h)) recover(h) += v;}
  10. };

线段树 (不封装) [#ds05]:

  1. #define segc int M = L + R - 1 >> 1, lc = id << 1, rc = lc | 1
  3. void add(int id, int L, int R, int h, int v){
  4.     if(L == R) return void(st[id] += v);
  5.     segc; h <= M ? add(lc, L, M, h, v) : add(rc, M + 1, R, h, v);
  6.     x[id].v = x[lc].v + x[rc].v;
  7. }
  9. int range(int id, int L, int R, int ql, int qr){
  10.     if(ql <= L && qr >= R) return st[id];
  11.     segc, s = 0;
  12.     if(ql <= M) s += range(lc, L, M, ql, min(qr, M));
  13.     if(qr > M) s += range(rc, M + 1, R, max(ql, M + 1), qr);
  14.     return s;
  15. }

线段树 (封装) [#ds06]:

  1. struct ST{
  2.     #define segc int M = L + R - 1 >> 1, lc = id << 1, rc = lc | 1
  3.     int sz;
  4.     struct node {int v, f; bool zero;} *x;
  5.     ST (): sz(0) {x = NULL;}
  6.     ~ST () {if(x) delete [] (x);}
  7.     void resize(int size) {sz = size; int sz0 = sz << 3; if(x) delete [] (x);
  8.         x = new node[sz0]; memset(x, 0, sz0 * sizeof(node));}
  10.     void add(int h, int v) {add(1, 1, sz, h, v);}
  11.     int range(int l, int r) {return query(1, 1, sz, l, r);}
  13.     void add(int id, int L, int R, int h, int v){
  14.         if(L == R) return void(x[id].v += v);
  15.         segc; h <= M ? add(lc, L, M, h, v) : add(rc, M + 1, R, h, v);
  16.         x[id].v = x[lc].v + x[rc].v;
  17.     }
  19.     int query(int id, int L, int R, int ql, int qr){
  20.         if(ql <= L && R <= qr) return x[id].v;
  21.         segc, s = 0;
  22.         if(ql <= M) s += query(lc, L, M, ql, min(qr, M));
  23.         if(qr > M) s += query(rc, M + 1, R, max(ql, M + 1), qr);
  24.         return s;
  25.     }
  26. };

可持久化线段树 (add 操作) [#ds07]:

  1. int add(int _id, int L, int R, int h, int v){
  2.     int id = ++cnt; x[id] = x[_id]; x[id].v += v;
  3.     if(L == R) return id;
  4.     int M = L + R - 1 >> 1;
  5.     if(h <= M) x[id].lc = add(x[id].lc, L, M, h, v);
  6.     else x[id].rc = add(x[id].rc, M + 1, R, h, v);
  7.     return id;
  8. }

伸展树 (Splay) [#ds08]:

  1. #define pa p[nd]
  2. #define root nd[0].c[0]
  3. struct node {int v, sz, p, c[2]; ll sum;} nd[N];
  5. inline int dir(int x) {return x == x[nd].pa.c[1];}
  7. inline void update(int x){
  8.     x[nd].sz = x[nd].c[0][nd].sz + x[nd].c[1][nd].sz + 1;
  9.     x[nd].sum = x[nd].c[0][nd].sum + x[nd].c[1][nd].sum + x[nd].v;
  10. }
  12. void rotate(int x){
  13.     int y = x[nd].p, d = !dir(x);
  14.     nd[y[nd].c[!d] = x[nd].c[d]].p = y;
  15.     x[nd].p = y[nd].p;
  16.     y[nd].pa.c[dir(y)] = x;
  17.     nd[x[nd].c[d] = y].p = x;
  18.     update(y);
  19. }
  21. void splay(int x, int g = 0){
  22.     for(; x[nd].p != g; rotate(x))
  23.         if(x[nd].pa.p != g) rotate(dir(x) ^ dir(x[nd].p) ? x : x[nd].p);
  24.     update(x);
  25. }
  27. void insert(int x){
  28.     int y = 0, d = 0;
  29.     if(root) for(y = root; d = (val[x] < y[nd].v), y[nd].c[d]; y = y[nd].c[d]);
  30.     y[nd].c[d] = x;
  31.     nd[x].v = nd[x].sum = val[x]; nd[x].sz = 1; nd[x].p = y; nd[x].c[0] = nd[x].c[1] = 0;
  32.     splay(x);
  33. }
  35. void erase(int x){
  36.     if(x[nd].p) splay(x);
  37.     if(!(x[nd].c[0] && x[nd].c[1])){
  38.         int d = (x[nd].c[1] > 0);
  39.         x[nd].c[d][nd].p = 0;
  40.         root = x[nd].c[d];
  41.     }else{
  42.         int y;
  43.         for(y = x[nd].c[1]; y[nd].c[0]; y = y[nd].c[0]);
  44.         splay(y, x);
  45.         nd[y[nd].c[0] = x[nd].c[0]].p = y;
  46.         y[nd].p = 0; root = y;
  47.         update(y);
  48.     }
  49. }
  51. int kth(int x, int v){
  52.     if(x[nd].sz < v) return -1;
  53.     for(int j; ; ){
  54.         j = x[nd].c[0][nd].sz;
  55.         if(v == j + 1) return x;
  56.         x = x[nd].c[v > j]; v > j ? (v -= j + 1) : v;
  57.     }
  58. }

动态树 (Link-Cut Tree) [#ds09]:

  1. namespace LCT {
  2.     #define pa p[nd]
  3.     struct node {bool rev; int v, p, c[2];} nd[N];
  4.     inline int dir(int x) {return !nd[x].p ? -1 : x == nd[x].pa.c[0] ? 0 : x == nd[x].pa.c[1] ? 1 : -1;}
  5.     inline void set(int x, int px, int c) {if (nd[x].p = px, ~c) nd[px].c[c] = x;}
  6.     inline void reverse(int x) {x && (std::swap(nd[x].c[0], nd[x].c[1]), nd[x].rev = !nd[x].rev);}
  7.     void push_down(int x) {if (nd[x].rev) reverse(nd[x].c[0]), reverse(nd[x].c[1]), nd[x].rev = false;}
  8.     void pull_down(int x) {if (~dir(x)) pull_down(nd[x].p); push_down(x);}
  9.     inline void update(int x) {const int l = nd[x].c[0], r = nd[x].c[1]; nd[x].v = (l ? nd[l].v : 0) + (r ? nd[r].v : 0);}
  10.     void rotate(int x) {int y = nd[x].p, d = !dir(x); set(nd[x].c[d], y, !d), set(x, nd[y].p, dir(y)), set(y, x, d);}
  11.     void splay(int x) {for (pull_down(x); ~dir(x); rotate(x)) if (~dir(nd[x].p)) rotate(dir(x) ^ dir(nd[x].p) ? x : nd[x].p); update(x);}
  12.     void access(int x) {for (int y = 0; x; y = x, x = nd[x].p) splay(x), nd[x].c[1] = y, update(x);}
  13.     void make_root(int x) {access(x), splay(x), reverse(x);}
  14.     int find_root(int x) {for (access(x), splay(x); push_down(x), nd[x].c[0]; x = nd[x].c[0]); return splay(x), x;}
  15.     int split(int x, int y) {return make_root(x), access(y), splay(y), y;}
  16.     void link(int x, int y) {make_root(x), nd[x].p = y;}
  17.     void cut(int x, int y) {split(x, y), nd[x].p = nd[y].c[0] = 0, update(y);}
  18.     void trylink(int x, int y) {x == y || (split(x, y), ~dir(x)) || (nd[x].p = y);}
  19.     void trycut(int x, int y) {split(x, y), nd[y].c[0] == x && !nd[x].c[1] && (nd[x].p = nd[y].c[0] = 0, update(y), 0);}
  20. }

树堆 (Treap) [#ds10]:

  1. struct random_t{
  2.     ull seed;
  3.     static const ull multiplier = 0x5deece66dll;;
  4.     static const ull addend = 0xbll;
  5.     static const ull mask = 0xffffffffffffll;
  7.     random_t () {char *x = new char; seed = (ull)x; delete x;}
  9.     unsigned int next(){
  10.         seed = (seed * multiplier + addend) & mask;
  11.         return seed >> 16;
  12.     }
  14.     unsigned int next(unsigned int n){
  15.         return n * (ull)next() >> 32;
  16.     }
  17. }rnd;
  19. #define pa p[nd]
  20. #define C(x) c[x][nd]
  21. #define root nd[0].c[0]
  23. struct node {int v, sz, pt; int c[2], p; unsigned priority;} nd[N];
  24. int cnt = 0;
  26. inline int dir(int x) {return x == nd[x].pa.c[1];}
  28. inline void update(int x){
  29.     nd[x].sz = nd[x].C(0).sz + nd[x].C(1).sz + nd[x].pt;
  30. }
  32. void rotate(int x){
  33.     int y = nd[x].p, d = !dir(x);
  34.     nd[nd[y].c[!d] = nd[x].c[d]].p = y;
  35.     nd[x].p = nd[y].p;
  36.     nd[y].pa.c[dir(y)] = x;
  37.     nd[nd[x].c[d] = y].p = x;
  38.     update(y); update(x);
  39. }
  41. void push_up(int x) {for(; x; x = nd[x].p) update(x);}
  43. void rotation(int x) {for(; nd[x].p && nd[x].priority < nd[x].pa.priority; ) rotate(x);}
  45. void insert(int v){
  46.     int x = 0, y = 0, d = 0;
  47.     if(root) for(y = root; nd[y].v != v && nd[y].c[d = v > nd[y].v]; y = nd[y].c[d]);
  48.     if(nd[y].v == v) {++nd[y].pt; push_up(y); return;}
  49.     nd[y].c[d] = x = ++cnt;
  50.     nd[x].v = v; nd[x].sz = nd[x].pt = 1; nd[x].p = y; nd[x].c[0] = nd[x].c[1] = 0;
  51.     nd[x].priority = rnd.next(); push_up(y);
  52.     rotation(x);
  53. }
  55. void erase(int v){
  56.     int x = 0, y = 0, d = 0;
  57.     if(root) for(x = root; nd[x].v != v && nd[x].c[d = v > nd[x].v]; x = nd[x].c[d]);
  58.     if(nd[x].v != v) return;
  59.     if(nd[x].pt > 1) {--nd[x].pt; push_up(x); return;}
  60.     for(; nd[x].c[0] && nd[x].c[1]; rotate(nd[x].c[d]))
  61.         d = nd[x].C(0).priority > nd[x].C(1).priority;
  62.     d = nd[x].c[1] > 0; y = nd[x].c[d];
  63.     nd[y].p = nd[x].p; nd[x].pa.c[dir(x)] = y; push_up(nd[y].p);
  64. }
  66. int rank(int v){
  67.     int x = 0, d = 0, k = 0;
  68.     if(root) for(x = root; nd[x].v != v && nd[x].c[d = v > nd[x].v]; x = nd[x].c[d])
  69.         if(d) k += nd[x].C(0).sz + nd[x].pt;
  70.     if(nd[x].v != v) return -1;
  71.     return k + nd[x].C(0).sz;
  72. }
  74. int kth(int k){
  75.     int x = root;
  76.     if(nd[x].sz < k) return -1;
  77.     for(int j, z; ; ){
  78.         j = nd[x].C(0).sz; z = nd[x].pt;
  79.         if(j <= k && k < j + z) return x;
  80.         x = nd[x].c[k > j]; k > j ? (k -= j + z) : k;
  81.     }
  82. }
  84. int prev(int v){
  85.     int x = 0, r = -1, d = 0;
  86.     if(root) for(x = root; x; x = nd[x].c[d]) if(d = v > x[nd].v) up(r, nd[x].v); // too little
  87.     return r;
  88. }
  90. int succ(int v){
  91.     int x = 0, r = INT_MAX, d = 0;
  92.     if(root) for(x = root; x; x = nd[x].c[d]) if(!(d = v >= x[nd].v)) down(r, nd[x].v); // too large
  93.     return r;
  94. }

析合树 (构造) [#ds11]:

  1. int n, p[N];
  3. namespace DCTree {
  4.     typedef std::pair <int, int> pr;
  5.     const int N = ::N * 2;
  7.     enum type {leaf, disjunct, conjunct} I[N];
  9.     pr st[20][N];
  10.     int stack1[N], stack2[N], stack[N];
  11.     int cnt, root, left[N], mid[N], right[N];
  13.     inline void up(pr &x, const pr &y) {x.first > y.first ? x.first = y.first : 0, x.second < y.second ? x.second = y.second : 0;}
  15.     void build_sparse_table() {
  16.         int i, j, k = n; pr *f, *g = *st;
  17.         for (j = 0; 1 << (j + 1) <= n; ++j) {
  18.             f = g, g = st[j + 1], k -= 1 << j;
  19.             for (i = 1; i <= k; ++i)
  20.                 up(g[i] = f[i], f[i + (1 << j)]);
  21.         }
  22.     }
  24.     inline bool is_consecutive(int L, int R) {
  25.         int c = lg2(R - L); pr ans = st[c][L]; up(ans, st[c][R - (1 << c)]);
  26.         return ans.second - ans.first == R - L - 1;
  27.     }
  29.     namespace ST {
  30.         #define segc int M = (L + R - 1) >> 1, lc = id << 1, rc = lc | 1
  31.         struct node {int v, tag;} x[N * 4];
  33.         void build(int id, int L, int R) {
  34.             x[id].v = L, x[id].tag = 0;
  35.             if (L == R) return;
  36.             segc; build(lc, L, M), build(rc, M + 1, R);
  37.         }
  39.         void add(int id, int L, int R, int ql, int qr, int v) {
  40.             if (ql <= L && R <= qr) {x[id].v += v, x[id].tag += v; return;}
  41.             segc;
  42.             if (ql <= M) add(lc, L, M, ql, qr, v);
  43.             if (qr > M) add(rc, M + 1, R, ql, qr, v);
  44.             x[id].v = std::min(x[lc].v, x[rc].v) + x[id].tag;
  45.         }
  47.         int find_suf(int id, int L, int R, int v, int cv = 0) {
  48.             if (cv + x[id].v > v) return -1;
  49.             if (L == R) return L;
  50.             segc, p = find_suf(lc, L, M, v, cv += x[id].tag);
  51.             return ~p ? p : find_suf(rc, M + 1, R, v, cv);
  52.         }
  53.     }
  55.     int pa[N], fc[N], nc[N], deg[N];
  57.     inline void link(int x, int px) {pa[x] = px, nc[x] = fc[px], fc[px] = x, ++deg[px];}
  59.     void build() {
  60.         int i, l, top1 = 0, top2 = 0, top = 0, &v = root, u; cnt = n;
  61.         for (i = 1; i <= n; ++i) st[0][i] = pr(p[i], p[i]), left[i] = right[i] = i, I[i] = leaf;
  62.         build_sparse_table(), ST::build(1, 1, n);
  63.         for (i = 1; i <= n; ++i) {
  64.             for (; top1 && p[i] > p[ stack1[top1] ]; --top1)
  65.                 ST::add(1, 1, n, stack1[top1 - 1] + 1, stack1[top1], p[i] - p[ stack1[top1] ]);
  66.             for (; top2 && p[i] < p[ stack2[top2] ]; --top2)
  67.                 ST::add(1, 1, n, stack2[top2 - 1] + 1, stack2[top2], p[ stack2[top2] ] - p[i]);
  68.             stack1[++top1] = stack2[++top2] = i;
  69.             l = ST::find_suf(1, 1, n, i);
  70.             for (v = i; top && left[ u = stack[top] ] >= l; --top)
  71.                 if (I[u] == conjunct && is_consecutive(mid[u], i + 1))
  72.                     right[u] = i, link(v, u), v = u;
  73.                 else if (is_consecutive(left[u], i + 1)) {
  74.                     I[++cnt] = conjunct, link(u, cnt), link(v, cnt);
  75.                     left[cnt] = left[u], right[cnt] = i, mid[cnt] = left[v], v = cnt;
  76.                 } else {
  77.                     I[++cnt] = disjunct, link(v, cnt);
  78.                     for (; top && !is_consecutive(left[ u = stack[top] ], i + 1); --top)
  79.                         link(u, cnt); link(u, cnt);
  80.                     left[cnt] = left[u], right[cnt] = i, v = cnt;
  81.                 }
  82.             stack[++top] = v;
  83.         }
  84.     }
  85. }

Chtholly 树 (ODT) / set 维护连续段 [#ds12]:

  1. namespace CTree {
  2.     map C;
  4.     map::iterator split(int pos) {
  5.         map::iterator it = C.lower_bound(pos), jt = it;
  6.         return it->first == pos ? it : C.emplace_hint(it, pos, (--jt)->second);
  7.     }
  9.     void modify(int l, int r, int v) {
  10.         map::iterator it = split(l), jt = split(r), i, j;
  11.         for (j = it; (i = j++) != jt; ) // do something on [i, j]
  12.         C.erase(it, jt), C.emplace(l, v); // do something
  13.     }
  14. }

可删除的堆 [#ds13]:

  1. struct heap {
  2.     std::priority_queue <int, std::vector <int>, std::greater <int>> A, B;
  4.     inline void emplace(int x) {A.emplace(x);}
  5.     inline void erase(int x) {B.emplace(x);}
  6.     inline void normalize() {for (; !B.empty() && A.top() == B.top(); A.pop(), B.pop());}
  7.     inline int top() {return normalize(), A.empty() ? INF : A.top();}
  8.     inline bool pop() {return normalize(), !A.empty() && (A.pop(), true);}
  9.     inline void clear() {for (; !A.empty(); A.pop()); for (; !B.empty(); B.pop());}
  10. };

IO 优化 [#ed01]:

  1. #define ID isdigit(c = *next++)
  2. #define IS isspace(c = *next++)
  4. struct Istream {
  5.     char *next, buf[20030731];
  6.     Istream & init(FILE *f = stdin) {fread(buf, 1, sizeof buf, f); next = buf; return *this;}
  7.     Istream & operator >> (int &x) {
  8.         int c; x = 0;
  9.         for (; !ID; ) if (!~c) return *this;
  10.         for (x = c & 15; ID; x = x * 10 + (c & 15));
  11.         return *this;
  12.     }
  13.     Istream & operator >> (char *x) {
  14.         int c;
  15.         for (; IS; ) if (!~c) return *this;
  16.         for (*x++ = c; !IS; *x++ = c) if (!~c) break;
  17.         return *x = 0, *this;
  18.     }
  19.     char get() {return *next++;}
  20. } cin;
  22. struct Ostream {
  23.     char *next, buf[20030731], _buf[34];
  24.     Ostream () {next = buf;}
  25.     void flush(FILE *f = stdout) {fwrite(buf, 1, next - buf, f); next = buf;}
  26.     Ostream & operator << (int x) {
  27.         if (!x) return put(48), *this;
  28.         int i;
  29.         for (i = 0; x; x /= 10) _buf[++i] = x % 10 | 48;
  30.         for (; i; --i) put(_buf[i]);
  31.         return *this;
  32.     }
  33.     Ostream & operator << (char c) {return put(c), *this;}
  34.     Ostream & operator << (const char *s) {for (char *p = (char*)s; *p; ++p) put(*p); return *this;}
  35.     void put(char c) {*next++ = c;}
  36. } cout;

动态规划的转移 [#ed02]:

  1. inline void up(int &x, const int y) {x < y ? x = y : 0;}
  2. inline void down(int &x, const int y) {x > y ? x = y : 0;}
  3. inline int min(const int x, const int y) {return x < y ? x : y;}
  4. inline int max(const int x, const int y) {return x < y ? y : x;}
  5. inline int & reduce(int &x) {return x += x >> 31 & mod;}
  6. inline void add(int &x, const int y) {x += y - mod, x += x >> 31 & mod;}
  7. inline void sub(int &x, const int y) {x -= y, x += x >> 31 & mod;}
  8. inline int & half(int &x) {return x = (x >> 1) + (-(x & 1) & iv2);}
  9. inline int & neg(int &x) {return x = (!x - 1) & (mod - x);}
  11. // another ver.
  12. inline bool up(int &x, const int y) {return x < y ? x = y, 1 : 0;}
  13. inline bool down(int &x, const int y) {return x > y ? x = y, 1 : 0;}
  14. inline int & add(int &x, const int y) {return x += y - mod, x += x >> 31 & mod;}
  15. inline int & sub(int &x, const int y) {return x -= y, x += x >> 31 & mod;}
  17. // templated ver.
  18. #define templated template <typename T>
  19. inline bool up(T &x, const T y) {return x < y ? x = y, 1 : 0;}
  20. inline bool down(T &x, const T y) {return x > y ? x = y, 1 : 0;}
  21. inline T min(const T x, const T y) {return x < y ? x : y;}
  22. inline T max(const T x, const T y) {return x < y ? y : x;}

乘模函数 (64 bit 大整数相乘取模) [#ed03]:

  1. inline ll MulMod(ll a, ll b, ll m) {
  2.     ll t = (a * b - (ll)((ld)a * b / m) * m) % m;
  3.     return t + (t >> 63 & m);
  4. }

ST 表 [#ed04]:

  1. namespace RMQ_simple {
  2.     void build_sparse_table() {
  3.         int *f, *g = *st, i, j, k = n;
  4.         for (j = 0; 1 << (j + 1) <= n; ++j) {
  5.             f = g, g = st[j + 1], k -= 1 << j;
  6.             for (i = 0; i < k; ++i)
  7.                 g[i] = min(f[i], f[i + (1 << j)]);
  8.         }
  9.     }
  11.     inline int range(int L, int R) {
  12.         int c = lg2(R - L);
  13.         return min(st[c][L], st[c][R - (1 << c)]);
  14.     }
  15. }
  17. namespace RMQ_fake {
  18.     const int LN = 18, TH = 18;
  20.     int n, *a, pre[N], suf[N], bel[N];
  21.     int L, st[LN][N / TH], *lay = *st;
  23.     void build(int n_, int *a_) {
  24.         int i, j, k, I, *f, *g = lay; n = n_, a = a_;
  25.         for (I = i = 0; I < n; ++i, I += TH) {
  26.             pre[I] = a[I], bel[I] = i;
  27.             for (j = 1; j < TH && I + j < n; ++j) pre[I + j] = min(pre[I + j - 1], a[I + j]), bel[I + j] = i;
  28.             for (--j, suf[I + j] = a[I + j]; --j >= 0; suf[I + j] = min(suf[I + j + 1], a[I + j]));
  29.             lay[i] = suf[I];
  30.         }
  31.         for (k = L = i, j = 0; 1 << (j + 1) <= L; ++j)
  32.             for (f = g, g = st[j + 1], k -= 1 << j, i = 0; i < k; ++i)
  33.                 g[i] = min(f[i], f[i + (1 << j)]);
  34.     }
  36.     inline int range(int L, int R) {
  37.         if (R <= L + TH) return *std::min_element(a + L, a + R);
  38.         int x = bel[L] + 1, y = bel[--R], c = lg2(y - x);
  39.         return min(min(suf[L], pre[R]), x < y ? min(st[c][x], st[c][y - (1 << c)]) : INT_MAX);
  40.     }
  41. }

离散化 [#ed05]:

  1. namespace DC {
  2.     int F[N]; pr D[N];
  4.     int Discretize(int n) {
  5.         int i, cnt = 0; std::sort(D, D + n);
  6.         for (i = 0; i < n; ++i)
  7.             F[D[i].second] = (i && D[i].first == D[i - 1].first ? cnt - 1 : (D[cnt] = D[i], cnt++));
  8.         return cnt;
  9.     }
  10. }
  11. Hash Map [#ed06]:
  12. class hash_map{
  13. public:
  14.     static const int HASH_MAX = 0xffffff, N = 8000000;
  15.     int cnt, first[HASH_MAX + 2], next[N]; data z[N];
  16.     inline int getHash(int key) {return (key ^ key << 3 ^ key >> 2) & HASH_MAX;}
  18.     void clear() {for(; cnt > 0; --cnt) first[z[cnt].hash] = 0;}
  20.     data * find(int key, bool inserted){
  21.         int x = getHash(key), i;
  22.         for(i = first[x]; i; i = next[i]) if(z[i].key == key) return z + i;
  23.         if(!inserted) return NULL;
  24.         z[++cnt] = data(key, 0, x); next[cnt] = first[x]; first[x] = cnt;
  25.         return z + cnt;
  26.     }
  27. };

快速 Fourier 变换 (Fast Fourier Transform) [#ed07]:

  1. // 'Fast Number Theory Transform' is in memos/12.html.
  2. typedef std::complex <double> C;
  4. namespace rpoly_base {
  5.     int l, n; double iv; C w2[N];
  7.     void init(int n = N) {
  8.         int i, t = min(n > 1 ? lg2(n - 1) : 0, 21); double angle = pi;
  9.         for (*w2 = C(1), i = 0; i <= t; ++i) w2[1 << i] = std::polar(1., angle *= .5);
  10.         for (i = 3; i < n; ++i) if (i & (i - 1)) w2[i] = w2[i & (i - 1)] * w2[i & -i];
  11.     }
  13.     inline void FFT_init(int len) {n = 1 << (l = len), iv = 1. / n;}
  15.     void DIF(C *a) {
  16.         int i, len = n >> 1; C *j, *k, *o, R;
  17.         for (i = 0; i < l; ++i, len >>= 1)
  18.             for (j = a, o = w2; j != a + n; j += len << 1, ++o)
  19.                 for (k = j; k != j + len; ++k)
  20.                     R = *o * k[len], k[len] = *k - R, *k += R;
  21.     }
  23.     void DIT(C *a) {
  24.         int i, len = 1; C *j, *k, *o, R;
  25.         for (i = 0; i < l; ++i, len <<= 1)
  26.             for (j = a, o = w2; j != a + n; j += len << 1, ++o)
  27.                 for (k = j; k != j + len; ++k)
  28.                     R = *k + k[len], k[len] = (*k - k[len]) * *o, *k = R;
  29.     }
  30. }
  32. namespace rpoly {
  33.     using namespace rpoly_base;
  35.     C B1[N], B2[N], B3[N], B4[N];
  37.     void mulf(int A, int B, double *a, double *b, double *c) {
  38.         if (!(A || B)) {*c = *a * *b; return;}
  39.         int i; FFT_init(lg2(A + B) + 1);
  40.         for (i = 0; i <= A; ++i) B1[i] = C(a[i]); memset(B1 + i, 0, (n - i) << 4);
  41.         for (i = 0; i <= B; ++i) B2[i] = C(b[i]); memset(B2 + i, 0, (n - i) << 4);
  42.         DIF(B1), DIF(B2);
  43.         for (i = 0; i < n; ++i) B1[i] *= B2[i];
  44.         DIT(B1), std::reverse(B1 + 1, B1 + n);
  45.         for (i = 0; i <= A + B; ++i) c[i] = B1[i].real() * iv;
  46.     }
  48.     void muli(int A, int B, int *a, int *b, int *c) {
  49.         if (!(A || B)) {*c = (ll)*a * *b % mod; return;}
  50.         int i, j, k; FFT_init(lg2(A + B) + 1);
  51.         for (i = 0; i <= A; ++i) B1[i] = C(a[i] & 32767, a[i] >> 15); memset(B1 + i, 0, (n - i) << 4);
  52.         for (i = 0; i <= B; ++i) B2[i] = C(b[i] & 32767, b[i] >> 15); memset(B2 + i, 0, (n - i) << 4);
  53.         DIF(B1), DIF(B2);
  54.         *B3 = B1->real() * *B2, *B4 = B1->imag() * *B2;
  55.         for (k = 0; k < l; ++k)
  56.             for (i = 1 << k, j = (2 << k) - 1; i < 2 << k; ++i, --j)
  57.                 B3[j] = (B1[i] + std::conj(B1[j])) * B2[i] * C(.5, 0.),
  58.                 B4[j] = (B1[i] - std::conj(B1[j])) * B2[i] * C(0., -.5);
  59.         DIT(B3), DIT(B4);
  60.         for (i = 0; i <= A + B; ++i)
  61.             c[i] = (ll(B3[i].real() * iv + .5) + (ll((B3[i].imag() + B4[i].real()) * iv + .5) % mod << 15) + (ll(B4[i].imag() * iv + .5) % mod << 30) + mod) % mod;
  62.     }
  63. }

Gauss 消元法 [#ed08]:

  1. // Gauss elimination with type 'double'
  2. struct LnEqn{
  3.     int sz;
  4.     double **m, *b;
  5.     LnEqn (): sz(0) {m = NULL; b = NULL;}
  6.     void resize(int size){
  7.         sz = size; m = new double *[sz];
  8.         for(int i = 0; i < sz; i++){
  9.             m[i] = new double[sz];
  10.             memset(m[i], 0, sz << 3);
  11.         }
  12.         b = new double[sz];
  13.         memset(b, 0, sz << 3);
  14.     }
  15.     ~LnEqn (){
  16.         if(m) {for(int i = 0; i < sz; i++) delete [] (m[i]); delete [] (m);}
  17.         if(b) delete [] (b);
  18.     }
  19.     bool solve(){
  20.         int i, j, k, maxi; double coe;
  21.         for(k = 0; k < sz; k++){
  22.             maxi = k;
  23.             for(i = k + 1; i < sz; i++)
  24.                 if(fabs(m[i][k]) > fabs(m[maxi][k]))
  25.                     maxi = i;
  26.             if(fabs(m[maxi][k]) < 1e-8) return false;
  27.             if(maxi != k){
  28.                 swap(m[maxi], m[k]);
  29.                 swap(b[maxi], b[k]);
  30.             }
  31.             coe = 1.0 / m[k][k];
  32.             for(j = 0; j < sz; j++)
  33.                 m[k][j] *= coe;
  34.             m[k][k] = 1.0;
  35.             b[k] *= coe;
  36.             for(i = 0; i < sz; i++){
  37.                 if(i == k) continue;
  38.                 coe = m[i][k];
  39.                 for(j = 0; j < sz; j++)
  40.                     m[i][j] -= coe * m[k][j];
  41.                 m[i][k] = 0.0;
  42.                 b[i] -= coe * b[k];
  43.             }
  44.         }
  45.         return true;
  46.     }
  47. };

线性规划 (xₚ ≥ 0, bₚ ≥ 0) [#ed09]:

  1. const double eps = 1e-8;
  3. int id[N + E];
  5. double m[E][N], b[E], *c = m[0], ans[N + E];
  7. void pivot(int toFree, int toBasic); // basic(1,m) -> free, free(1,n) -> basic
  9. void pivot(int r, int c){
  10.     int i, j; double coe = 1.0 / m[r][c];
  11.     swap(id[n + r], id[c]);
  12.     m[r][c] = 1.0;
  13.     for(j = 1; j <= n; ++j)
  14.         m[r][j] *= coe;
  15.     b[r] *= coe;
  16.     for(i = 0; i <= e; ++i){
  17.         if(i == r) continue;
  18.         coe = m[i][c]; m[i][c] = 0.0;
  19.         for(j = 1; j <= n; ++j)
  20.             m[i][j] -= coe * m[r][j];
  21.         b[i] -= coe * b[r];
  22.     }
  23. }
  25. bool simplex(){
  26.     int i, j, basic, free; double G;
  27.     for(; ; ){
  28.         basic = free = 0; G = INFINITY;
  29.         for(i = 1; i <= n; ++i) // free (nonbasic) variable
  30.             if(c[i] > c[free]) free = i;
  31.         if(!free) return true;
  32.         for(j = 1; j <= e; ++j) // basic variable
  33.             if(m[j][free] > eps && b[j] < G * m[j][free]){
  34.                 G = b[j] / m[j][free]; basic = j;
  35.             }
  36.         if(!basic) return puts("Unbounded"), false;
  37.         pivot(basic, free);
  38.     }
  39. }
  41. // initialize :
  42. for(j = 1; j <= n + e; ++j) id[j] = j;
  43. c[0] = eps;
  45. // output:
  46. if(simplex()){
  47.     printf("%.8lg\n", -*b);
  48.     for(i = 1; i <= e; ++i) ans[id[n + i]] = b[i];
  49.     for(j = 1; j <= n; ++j) printf("%.8lg%c", ans[j], j == n ? 10 : 32);
  50. }

Gauss 消元求行列式 (模意义) [#ed10]:

  1. int gauss(int n) {
  2.     int i, j, k, det = 1; ll coe;
  3.     static int *m[N];
  4.     for (i = 0; i < n; ++i) m[i] = mat[i];
  5.     for (i = 0; i < n; ++i) {
  6.         for (j = i; j < n && !m[j][i]; ++j);
  7.         if (j == n) return 0;
  8.         if (i != j) det = mod - det, std::swap(m[i], m[j]);
  9.         det = (ll)det * m[i][i] % mod;
  10.         coe = PowerMod(m[i][i], mod - 2);
  11.         for (j = 0; j < n; ++j) m[i][j] = m[i][j] * coe % mod;
  12.         for (k = i + 1; k < n; ++k)
  13.             for (coe = mod - m[k][i], j = i; j < n; ++j) m[k][j] = (m[k][j] + coe * m[i][j]) % mod;
  14.     }
  15.     return det;
  16. }

Nim 积 (32/64 bit) [#ed11]:

  1. namespace nimbers {
  2.     constexpr u32 n2f[16] = {0x0001u, 0x8ff5u, 0x6cbfu, 0xa5beu, 0x761au, 0x8238u, 0x4f08u, 0x95acu, 0xf340u, 0x1336u, 0x7d5eu, 0x86e7u, 0x3a47u, 0xe796u, 0xb7c3u, 0x0008u},
  3.                   f2n[16] = {0x0001u, 0x2827u, 0x392bu, 0x8000u, 0x20fdu, 0x4d1du, 0xde4au, 0x0a17u, 0x3464u, 0xe3a9u, 0x6d8du, 0x34bcu, 0xa921u, 0xa173u, 0x0ebcu, 0x0e69u};
  4.     inline u32 nimber2field(u32 x) {u32 y = 0; for (; x; x &= x - 1) y ^= n2f[ctz(x)]; return y;}
  5.     inline u32 field2nimber(u32 x) {u32 y = 0; for (; x; x &= x - 1) y ^= f2n[ctz(x)]; return y;}
  6.     inline u32 __builtin_double(u32 x) {return x << 1 ^ (x < 0x8000u ? 0 : 0x1681fu);}
  8.     u16 ln[65536], exp[131075], *Hexp = exp + 3, *H2exp = exp + 6;
  10.     inline void init() {
  11.         int i;
  12.         for (*exp = i = 1; i < 65535; ++i) exp[i] = __builtin_double(exp[i - 1]);
  13.         for (i = 1; i < 65535; ++i) exp[i] = field2nimber(exp[i]), ln[exp[i]] = i;
  14.         memcpy(exp + 65535, exp, 131070);
  15.         memcpy(exp + 131070, exp, 10);
  16.     }
  18.     inline u16 product(u16 A, u16 B) {return A && B ? exp[ln[A] + ln[B]] : 0;}
  19.     inline u16 H(u16 A) {return A ? Hexp[ln[A]] : 0;}
  20.     inline u16 H2(u16 A) {return A ? H2exp[ln[A]] : 0;}
  21.     inline u16 Hproduct(u16 A, u16 B) {return A && B ? Hexp[ln[A] + ln[B]] : 0;}
  23.     inline u32 product(u32 A, u32 B) {
  24.         u16 a = A & 65535, b = B & 65535, c = A >> 16, d = B >> 16, e = product(a, b);
  25.         return u32(product(u16(a ^ c), u16(b ^ d)) ^ e) << 16 | (Hproduct(c, d) ^ e);
  26.     }
  28.     inline u32 H(u32 A) {
  29.         u16 a = A & 65535, b = A >> 16;
  30.         return H(u16(a ^ b)) << 16 | H2(b);
  31.     }
  33.     inline u64 product(u64 A, u64 B) {
  34.         u32 a = A & UINT_MAX, b = B & UINT_MAX, c = A >> 32, d = B >> 32, e = product(a, b);
  35.         return u64(product(a ^ c, b ^ d) ^ e) << 32 | (H(product(c, d)) ^ e);
  36.     }
  37. }

01 背包问题 [#ed12]:

  1. // 01 背包
  2. // 有 N 件物品和一个容量为 V 的背包
  3. // 第 i 件物品费用是 c[i]，价值是 w[i]
  4. // 求解将哪些物品装入背包，可使这些物品的总体积不超过背包容量，且总价值最大。
  6. // 无优化
  7. for (int i = 1; i <= N; i++)
  8. {
  9.     for (int j = 1; j <= V; j++)
  10.     {
  11.         if (j < c[i]) dp[i][j] = dp[i - 1][j];
  12.         else dp[i][j] = std::max(dp[i - 1][j], dp[i - 1][j - c[i]] + w[i]);
  13.     }
  14. }
  16. // 空间优化
  17. for (int i = 1; i <= N; i++)
  18. {
  19.     for (int j = V; j >= c[i]; j--)
  20.     {
  21.         dp[j] = std::max(dp[j], dp[j - c[i]] + w[i]);
  22.     }
  23. }

多重背包问题 [#ed13]:

  2. // 多重背包问题
  3. // 有 N 件物品和一个容量为 V 的背包，每种物品都有无限件可用
  4. // 第 i 件物品费用是 c[i]，价值是 w[i]
  5. // 求解将哪些物品装入背包，可使这些物品的总体积不超过背包容量，且总价值最大。
  7. // 无优化
  8. for (int i = 1; i <= N; i++)
  9. {
  10.     for (int j = 1; j <= V; j++)
  11.     {
  12.         dp[i][j] = dp[i - 1][j];
  13.         if (j >= c[i]) dp[i][j] = std::max(dp[i - 1][j], dp[i][j - c[i]] + w[i]);
  14.     }
  15. }
  17. // 空间优化
  18. for (int i = 1; i <= N; i++)
  19. {
  20.     for (int j = c[i]; j <= V; j++)
  21.     {
  22.         dp[j] = std::max(dp[j], dp[j - c[i]] + w[i]);
  23.     }
  24. }

整数二分 [#ed14]:

  2. bool check(int x) {/* ... */} // 检查x是否满足某种性质
  4. // 区间[l, r]被划分成[l, mid]和[mid + 1, r]时使用：
  5. int bsearch(int l, int r)
  6. {
  7.     while (l < r)
  8.     {
  9.         int mid = l + r << 1;
  10.         if (check(mid)) r = mid;    // check()判断mid是否满足性质
  11.         else l = mid + 1;
  12.     }
  13.     return l;
  14. }
  16. // 区间[l, r]被划分成[l, mid - 1]和[mid, r]时使用：
  17. int bsearch(int l, int r)
  18. {
  19.     while (l < r)
  20.     {
  21.         int mid = l + r + 1 << 1;
  22.         if (check(mid)) l = mid;
  23.         else r = mid - 1;
  24.     }
  25.     return l;
  26. }

高精度加法 [#ed15]:

  1. std::vector <int> add(std::vector <int> &A, std::vector <int> &B)
  2. {
  3.     std::vector <int> C;
  5.     int t = 0;
  6.     for (int i = 0; i < A.size() || i < B.size(); i++)
  7.     {
  8.         if (i < A.size()) t += A[i];
  9.         if (i < B.size()) t += B[i];
  10.         C.push_back(t % 10);
  11.         t /= 10;
  12.     }
  14.     if (t) C.push_back(1);
  15.     return C;
  16. }

二维向量/点、计算几何基础 [#cg01]:

  1. const double eps = 1e-8;
  3. #define lt(x, y) ((x) < (y) - eps)
  4. #define gt(x, y) ((x) > (y) + eps)
  5. #define le(x, y) ((x) <= (y) + eps)
  6. #define ge(x, y) ((x) >= (y) - eps)
  7. #define eq(x, y) (le(x, y) && ge(x, y))
  8. #define dot(x, y, z) (((y) - (x)) * ((z) - (x)))
  9. #define cross(x, y, z) (((y) - (x)) ^ ((z) - (x)))
  11. struct vec2 {
  12.     double x, y;
  13.     vec2 (double x0 = 0.0, double y0 = 0.0) : x(x0), y(y0) {}
  14.     vec2 * read() {scanf("%lf%lf", &x, &y); return this;}
  15.     inline vec2 operator - () const {return vec2(-x, -y);}
  16.     inline vec2 operator + (const vec2 &B) const {return vec2(x + B.x, y + B.y);}
  17.     inline vec2 operator - (const vec2 &B) const {return vec2(x - B.x, y - B.y);}
  18.     inline vec2 operator * (double k) const {return vec2(x * k, y * k);}
  19.     inline vec2 operator / (double k) const {return *this * (1.0 / k);}
  20.     inline double operator * (const vec2 &B) const {return x * B.x + y * B.y;}
  21.     inline double operator ^ (const vec2 &B) const {return x * B.y - y * B.x;}
  22.     inline double norm2() const {return x * x + y * y;}
  23.     inline double norm() const {return sqrt(x * x + y * y);}
  24.     inline bool operator < (const vec2 &B) const {return lt(x, B.x) || le(x, B.x) && lt(y, B.y);}
  25.     inline bool operator == (const vec2 &B) const {return eq(x, B.x) && eq(y, B.y);}
  26.     inline bool operator << (const vec2 &B) const {return lt(y, 0) ^ lt(B.y, 0) ? lt(B.y, 0) : gt(*this ^ B, 0) || ge(*this ^ B, 0) && ge(x, 0) && lt(B.x, 0);}
  27.     inline vec2 trans(double a11, double a12, double a21, double a22) const {return vec2(x * a11 + y * a12, x * a21 + y * a22);}
  28. };
  30. /*
  31. operator * : Dot product
  32. operator ^ : Cross product
  33. norm2() : |v|^2 = v.v
  34. norm() : |v| = sqrt(v.v)
  35. operator < : Two-key compare
  36. operator << : Polar angle compare
  37. trans : Transition with a 2x2 matrix
  38. */

直线及其函数 [#cg02]:

  1. struct line {
  2.     double A, B, C; // Ax + By + C = 0, > 0 if it represents halfplane.
  3.     line (double A0 = 0.0, double B0 = 0.0, double C0 = 0.0) : A(A0), B(B0), C(C0) {}
  4.     line (const vec2 &u, const vec2 &v) : A(u.y - v.y), B(v.x - u.x), C(u ^ v) {} // left > 0
  5.     inline vec2 normVec() const {return vec2(A, B);}
  6.     inline double norm2() const {return A * A + B * B;}
  7.     inline double operator () (const vec2 &P) const {return A * P.x + B * P.y + C;}
  8. };
  10. inline vec2 intersection(const line u, const line v) {return vec2(u.B * v.C - u.C * v.B, u.C * v.A - u.A * v.C) / (u.A * v.B - u.B * v.A);}
  12. inline bool parallel(const line u, const line v) {double p = u.normVec() ^ v.normVec(); return eq(p, 0);}
  14. inline bool perpendicular(const line u, const line v) {double p = u.normVec() * v.normVec(); return eq(p, 0);}
  16. inline bool sameDir(const line u, const line v) {return parallel(u, v) && gt(u.normVec() * v.normVec(), 0);}
  18. inline line bisector(const vec2 u, const vec2 v) {return line(v.x - u.x, v.y - u.y, 0.5 * (u.norm2() - v.norm2()));}
  20. inline double dis2(const vec2 P, const line l) {return l(P) * l(P) / l.norm2();}
  22. inline vec2 projection(const vec2 P, const line l) {return P - l.normVec() * (l(P) / l.norm2());}
  24. inline vec2 symmetry(const vec2 P, const line l) {return P - l.normVec() * (2 * l(P) / l.norm2());}

多边形操作 [#cg03]:

  1. // Relation of 3 points. (2 inside, 1 outside, 0 not collinear)
  2. inline int collinear(const vec2 u, const vec2 v, const vec2 P) {double p = cross(P, u, v); return eq(p, 0) ? 1 + le(dot(P, u, v), 0) : 0;}
  4. // Perimeter of a polygon
  5. double perimeter(int n, vec2 *poly) {
  6.     double ret = (poly[n - 1] - *poly).norm();
  7.     for (int i = 1; i < n; ++i) ret += (poly[i - 1] - poly[i]).norm();
  8.     return ret;
  9. }
  11. // Directed area of a polygon (positive if CCW)
  12. double area(int n, vec2 *poly) {
  13.     double ret = poly[n - 1] ^ *poly;
  14.     for (int i = 1; i < n; ++i) ret += poly[i - 1] ^ poly[i];
  15.     return ret * 0.5;
  16. }
  18. // Point in polygon (2 on boundary, 1 inside, 0 outside)
  19. int contain(int n, vec2 *poly, const vec2 P) {
  20.     int in = 0; vec2 *r = poly + (n - 1), *u, *v;
  21.     for (int i = 0; i < n; r = poly, ++poly, ++i) {
  22.         if (collinear(*r, *poly, P) == 2) return 2;
  23.         gt(r->y, poly->y) ? (u = poly, v = r) : (u = r, v = poly);
  24.         if (ge(P.y, v->y) || lt(P.y, u->y)) continue;
  25.         in ^= gt(cross(P, *u, *v), 0);
  26.     }
  27.     return in;
  28. }

平面凸包 (Graham Scan) [#cg04]:

  1. int graham(int n, vec2 *p, vec2 *dest) {
  2.     int i; vec2 *ret = dest;
  3.     std::iter_swap(p, std::min_element(p, p + n));
  4.     std::sort(p + 1, p + n, [p] (const vec2 A, const vec2 B) {double r = cross(*p, A, B); return gt(r, 0) || (ge(r, 0) && lt((A - *p).norm2(), (B - *p).norm2()));});
  5.     for (i = 0; i < 2 && i < n; ++i) *ret++ = p[i];
  6.     for (; i < n; *ret++ = p[i++])
  7.         for (; ret != dest + 1 && ge(cross(ret[-2], p[i], ret[-1]), 0); --ret);
  8.     return *ret = *p, ret - dest;
  9. }

旋转卡壳求凸集直径 [#cg05]:

  1. double convDiameter(int n, vec2 *poly) {
  2.     int l = std::min_element(poly, poly + n) - poly, r = std::max_element(poly, poly + n) - poly, i = l, j = r;
  3.     double diam = (poly[l] - poly[r]).norm2();
  4.     do {
  5.         (ge(poly[(i + 1) % n] - poly[i] ^ poly[(j + 1) % n] - poly[j], 0) ? ++j : ++i) %= n;
  6.         up(diam, (poly[i] - poly[j]).norm2());
  7.     } while (i != l || j != r);
  8.     return diam;
  9. }

三角形外心 & 最小圆覆盖 [#cg06]:

  1. inline vec2 circumCenter(const vec2 A, const vec2 B, const vec2 C) {
  2.     vec2 a = B - A, b = C - A, AO;
  3.     double det = 0.5 / (a ^ b), na = a.norm2(), nb = b.norm2();
  4.     AO = vec2((na * b.y - nb * a.y) * det, (nb * a.x - na * b.x) * det);
  5.     return A + AO;
  6. }
  8. double minCircleCover(int n, vec2 *p, vec2 *ret = NULL) {
  9.     int i, j, k; double r2 = 0.0;
  10.     std::random_shuffle(p + 1, p + (n + 1));
  11.     vec2 C = p[1];
  12.     for (i = 2; i <= n; ++i)
  13.         if (gt((p[i] - C).norm2(), r2))
  14.             for (C = p[i], r2 = 0, j = 1; j < i; ++j)
  15.                 if (gt((p[j] - C).norm2(), r2))
  16.                     for (C = (p[i] + p[j]) * 0.5, r2 = (p[j] - C).norm2(), k = 1; k < j; ++k)
  17.                         if (gt((p[k] - C).norm2(), r2))
  18.                             C = circumCenter(p[i], p[j], p[k]), r2 = (p[k] - C).norm2();
  19.     return ret ? *ret = C : 0, r2;
  20. }

半平面交 (平行处理版) [#cg07]:

  1. inline bool HPIcmp(const line u, const line v) {return sameDir(u, v) ? gt((fabs(u.A) + fabs(u.B)) * v.C, (fabs(v.A) + fabs(v.B)) * u.C) : u.normVec() << v.normVec();}
  3. inline bool geStraight(const vec2 A, const vec2 B) {return lt(A ^ B, 0) || le(A ^ B, 0) && lt(A * B, 0);}
  5. inline bool para_nega_test(const line u, const line v) {
  6.     return parallel(u, v) && lt(u.normVec() * v.normVec(), 0) && (fabs(u.A) + fabs(u.B)) * v.C + (fabs(v.A) + fabs(v.B)) * u.C < -7e-14;
  7. }
  9. int HPI(int n, line *l, vec2 *poly, vec2 *&ret) { // -1 : Unbounded, -2 : Infeasible
  10.     int h = 0, t = -1, i, j, que[n + 5];
  11.     std::sort(l, l + n, HPIcmp);
  12.     n = std::unique(l, l + n, sameDir) - l;
  13.     for (j = i = 0; i < n && j < n; ++i) {
  14.         for (up(j, i + 1); j < n && !geStraight(l[i].normVec(), l[j].normVec()); ++j);
  15.         if (para_nega_test(l[i], l[j])) return -2;
  16.     }
  17.     if (n <= 2 || geStraight(l[n - 1].normVec(), l->normVec())) return -1;
  18.     for (i = 0; i < n; ++i) {
  19.         if (geStraight(l[que[t]].normVec(), l[i].normVec())) return -1;
  20.         for (; h < t && le(l[i](poly[t - 1]), 0); --t);
  21.         for (; h < t && le(l[i](poly[h]), 0); ++h);
  22.         que[++t] = i; h < t ? poly[t - 1] = intersection(l[ que[t - 1] ], l[ que[t] ]) : 0;
  23.     }
  24.     for (; h < t && le(l[ que[h] ](poly[t - 1]), 0); --t);
  25.     return h == t ? -2 : (poly[t] = intersection(l[ que[t] ], l[ que[h] ]), ret = poly + h, t - h + 1);
  26. }

平面上的圆几何 [#cg08]:

  1. const double pi = M_PI, pi2 = pi * 2., pi_2 = M_PI_2;
  3. inline double angle(const vec2 u, const vec2 v) {return atan2(u ^ v, u * v);}
  5. // intersection of circle and line
  6. int intersection(double r2, const vec2 O, const line l, vec2 *ret) {
  7.     double d2 = dis2(O, l); vec2 j = l.normVec();
  8.     if (gt(d2, r2)) return ret[0] = ret[1] = vec2(INFINITY, INFINITY), 0;
  9.     if (ge(d2, r2)) return ret[0] = ret[1] = projection(O, l), 1;
  10.     if (le(d2, 0)) {
  11.         j = j * sqrt(r2 / j.norm2());
  12.         ret[0] = O + j.trans(0, -1, 1, 0);
  13.         ret[1] = O + j.trans(0, 1, -1, 0);
  14.     } else {
  15.         double T = copysign(sqrt((r2 - d2) / d2), l(O)); j = j * (-l(O) / j.norm2());
  16.         ret[0] = O + j.trans(1, T, -T, 1);
  17.         ret[1] = O + j.trans(1, -T, T, 1);
  18.     }
  19.     return 2;
  20. }
  22. // area of intersection(x^2 + y^2 = r^2, triangle OBC)
  23. double interArea(double r2, const vec2 B, const vec2 C) {
  24.     if (eq(B ^ C, 0)) return 0;
  25.     vec2 is[2];
  26.     if (!intersection(r2, vec2(), line(B, C), is)) return 0.5 * r2 * angle(B, C);
  27.     bool b = gt(B.norm2(), r2), c = gt(C.norm2(), r2);
  28.     if (b && c) return 0.5 * (lt(dot(*is, B, C), 0) ? r2 * (angle(B, *is) + angle(is[1], C)) + (is[0] ^ is[1]) : r2 * angle(B, C));
  29.     else if (b) return 0.5 * (r2 * angle(B, *is) + (*is ^ C));
  30.     else if (c) return 0.5 * ((B ^ is[1]) + r2 * angle(is[1], C));
  31.     else return 0.5 * (B ^ C);
  32. }
  34. // tangents to circle((0, 0), r) through P
  35. int tangent(double r, const vec2 P, line *ret) {
  36.     double Q = P.norm2() - r * r;
  37.     if (lt(Q, 0)) return ret[0] = ret[1] = line(INFINITY, INFINITY, INFINITY), 0;
  38.     if (le(Q, 0)) return ret[0] = ret[1] = line(P.x, P.y, -P.norm2()), 1;
  39.     Q = sqrt(Q) / r;
  40.     ret[0] = line(P.x + Q * P.y, P.y - Q * P.x, -P.norm2());
  41.     ret[1] = line(P.x - Q * P.y, P.y + Q * P.x, -P.norm2());
  42.     return 2;
  43. }
  45. // tangets to circle(O, r) through P
  46. int tangent(double r, const vec2 O, const vec2 P, line *ret) {
  47.     int R = tangent(r, P - O, ret);
  48.     if (R)
  49.         ret[0].C -= ret[0].A * O.x + ret[0].B * O.y,
  50.         ret[1].C -= ret[1].A * O.x + ret[1].B * O.y;
  51.     return R;
  52. }

三维计算几何基础 [#cg09]:

  1. #define triple(x, y, z) ((x) * ((y) ^ (z)))
  3. struct vec3 {
  4.     double x, y, z;
  5.     vec3 (double x0 = 0, double y0 = 0, double z0 = 0) : x(x0), y(y0), z(z0) {}
  6.     vec3 * read() {scanf("%lf%lf%lf", &x, &y, &z); return this;}
  7.     inline vec3 operator - () const {return vec3(-x, -y, -z);}
  8.     inline vec3 operator + (const vec3 &B) const {return vec3(x + B.x, y + B.y, z + B.z);}
  9.     inline vec3 operator - (const vec3 &B) const {return vec3(x - B.x, y - B.y, z - B.z);}
  10.     inline vec3 operator * (double k) const {return vec3(x * k, y * k);}
  11.     inline vec3 operator / (double k) const {return *this * (1.0 / k);}
  12.     inline double operator * (const vec3 &B) const {return x * B.x + y * B.y + z * B.z;}
  13.     inline vec3 operator ^ (const vec3 &B) const {return vec3(y * B.z - z * B.y, z * B.x - x * B.z, x * B.y - y * B.x);}
  14.     inline double norm2() const {return x * x + y * y + z * z;}
  15.     inline double norm() const {return sqrt(x * x + y * y + z * z);}
  16.     inline bool operator < (const vec3 &B) const {return lt(x, B.x) || le(x, B.x) && (lt(y, B.y) || le(y, B.y) && lt(z, B.z));}
  17.     inline bool operator == (const vec3 &B) const {return eq(x, B.x) && eq(y, B.y) && eq(z, B.z);}
  18. };
  20. // Positive if Right-hand rule
  21. inline double volume(const vec3 A, const vec3 B, const vec3 C, const vec3 D) {return triple(B - A, C - A, D - A);}
  23. struct line3 {
  24.     vec3 P, t;
  25.     line3 (vec3 _P = vec3(), vec3 _t = vec3()) : P(_P), t(_t) {}
  26. };
  28. inline double dis2(const vec3 P, const line3 l) {return ((P - l.P) ^ l.t).norm2() / l.t.norm2();}
  30. struct plane {
  31.     double A, B, C, D; // Ax + By + Cz + D = 0
  32.     plane (double A0 = 0.0, double B0 = 0.0, double C0 = 0.0, double D0 = 0.0) : A(A0), B(B0), C(C0), D(D0) {}
  33.     plane (const vec3 &u, const vec3 &v, const vec3 &w) {vec3 t = (v - u) ^ (w - u); A = t.x, B = t.y, C = t.z, D = -triple(u, v, w);} // > 0 if it follows Right-hand rule
  34.     inline vec3 normVec() const {return vec3(A, B, C);}
  35.     inline double norm2() const {return A * A + B * B + C * C;}
  36.     inline double operator () (const vec3 &P) const {return A * P.x + B * P.y + C * P.z + D;}
  37. };
  39. inline double dis2(const vec3 P, const plane F) {return F(P) * F(P) / F.norm2();}

三维凸包 (卷包裹法) [#cg10]:

  1. namespace CH3D {
  2.     typedef std::pair <int, int> pr;
  3.     typedef std::set <pr> set;
  5.     struct triangle {vec3 A, B, C;} tr[N];
  7.     vec3 p[N], q[N], r[N];
  8.     set E;
  9.     int n, cnt;
  11.     inline vec3 randvec3() {return vec3((double)rand() / RAND_MAX, (double)rand() / RAND_MAX, (double)rand() / RAND_MAX);}
  13.     void wrap(int u, int v) {
  14.         if (E.find({u, v}) == E.end()) {
  15.             int i, w = -1;
  16.             for (i = 0; i < n; ++i)
  17.                 if (i != u && i != v && (!~w || ge(volume(q[w], q[u], q[v], q[i]), 0))) w = i;
  18.             if (~w) {
  19.                 tr[cnt++] = (triangle){p[u], p[v], p[w]};
  20.                 E.emplace(u, v); E.emplace(v, w); E.emplace(w, u);
  21.                 wrap(w, v); wrap(u, w);
  22.             }
  23.         }
  24.     }
  26.     void main(int _n, vec3 *_p) {
  27.         int i;
  28.         n = _n; cnt = 0; E.clear();
  29.         memcpy(p, _p, n * sizeof(vec3));
  30.         std::iter_swap(p, std::min_element(p, p + n));
  31.         for (i = 0; i < n; ++i) q[i] = p[i] + randvec3() * 1e-6;
  32.         for (i = 2; i < n; ++i)
  33.             if (ge(cross(q[0], q[i], q[1]).z, 0)) std::swap(p[1], p[i]), std::swap(q[1], q[i]);
  34.         wrap(0, 1);
  35.     }
  36. }

自适应 Simpson 积分 [#cg11]:

  1. double Simpson(double L, double M, double R, double fL, double fM, double fR, double eps) {
  2.     double LM = (L + M) * 0.5, fLM = f(LM), MR = (M + R) * 0.5, fMR = f(MR);
  3.     double A = (fL + fM * 4.0 + fR) * (R - L) * sixth,
  4.            AL = (fL + fLM * 4.0 + fM) * (M - L) * sixth,
  5.            AR = (fM + fMR * 4.0 + fR) * (R - M) * sixth;
  6.     if (fabs(AL + AR - A) < eps) return (2.0 * (AL + AR) + A) * third;
  7.     return Simpson(L, LM, M, fL, fLM, fM, eps * 0.6) + Simpson(M, MR, R, fM, fMR, fR, eps * 0.6);
  8. }
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注