[关闭]
@Voleking 2017-03-13T10:57:35.000000Z 字数 24184 阅读 817

《挑战程序设计竞赛》Summary

ICPC Algorithms Unfinished

Tips

  1. pow 返回浮点数可能会导致误差，慎用！！！自己写。
  2. 样例数用 long long 可能会 TLE
  3. 多组数据记得初始化
  4. 小心过程量 overflow
  5. Type convertions
  6. 在取余的情况下，要避免减法运算的结果出现负数 (...+ MOD) % MOD; 多加点总不会错？
  7. fill 注意初始范围为从 1 到 n 的情况
  8. 提交换行
  9. c++11: long long int abs (long long int n);
  10. 字符串多输出 0
  11. 输入读取判断其中一个为零
  12. 搜索 flag 标记可能要加上方案数
  13. 数据范围特别大时考虑不可能情况
  14. SegmentTree 开大点以防 RE （1e5 -> 1 << 17)
  15. 滚动数组未全部清零

Basic Algorithm

Sort

  1. // Bubble Sort 
  2. for (int i = 0; i < N; i++)
  3.     for (int j = 0; j < N - i - 1; j++)
  4.         if (A[j] > A[j + 1]) swap(A[j], A[j + 1]);
  1. // Insertion Sort
  2. for (int i = 1; i < N; i++) {
  3.     int tmp = A[i], j;
  4.     for (j = i - 1; j >= 0 && A[j] > tmp; j--)
  5.         A[j + 1] = A[j];
  6.     A[++j] = tmp; 
  7. }
  1. // Selection Sort
  2. for (int i = 0; i < N; i++)
  3.     for (int j = i + 1; j < N; j++)
  4.         if (A[i] > A[j]) swap(A[i], A[j])
  1. // Merge Sort
  2. void merge_sort(int A[], int l, int r) {
  3.     if (l >= r) return ;
  4.     int mid = (l + r) / 2;
  5.     merge_sort(A, l, mid);
  6.     merge_sort(A, mid + 1, r); 
  7.     int i, j, k;
  8.     i = l; j = mid + 1; k = l;
  9.     while (i <= mid && j <= r) {
  10.         if (A[i] <= A[j]) {
  11.             B[k] = A[i++];
  12.         }
  13.         else {
  14.             B[k] = A[j++];
  15.         }
  16.         k++;
  17.     }
  18.     while (i <= mid) {
  19.         B[k] = A[i++];
  20.         k++;
  21.     }
  22.     while (j <= r) {
  23.         B[k] = A[j++];
  24.         k++;
  25.     }
  26.     memcpy(A, B, r - l + 1);
  27.     return ; 
  28. }
  1. // Quick Sort
  2. void quicksort(int A[], int l, int r) {
  3.     int i = l, j = r, mid = A[(r - l) / 2 + l];
  4.     while (i <= j) {
  5.         while (A[i] < mid) i++;
  6.         while (A[j] > mid) j--;
  7.         if (i <= j) {
  8.             swap(A[i], A[j]);
  9.             ++i; --j;
  10.         }
  11.     }
  12.     if (i < r) quicksort(A, i, r);
  13.     if (l < j) quicksort(A, l, j);        
  14.     return ;
  15. }

DP

LIS

  1. int A[MAX_N];
  2. long lis(int n) {
  3.     int dp[MAX_N];
  4.     fill(dp, dp + n, INF);
  5.     for (int i = 0; i < n; ++i)
  6.         *lower_bound(dp, dp + n, A[i]) = A[i];// lds: -A[i]; ln: upper_bound
  7.     return lower_bound(dp, dp + n, INF) - dp;
  8. }

Knapsack Problem

  1. for (int i = 0; i < N; ++i)
  2.     for (int j = W; j >= w[i]; --j)
  3.         f[j] = max(f[j - w[i]] + c[i], f[j]);

  1. for (int i = 0; i < N; ++i)
  2.     for (int j = w[i]; j <= W; ++j)
  3.         f[j] = max(f[j - w[i]] + c[i], f[v]);

注意循环顺序的不同背后思路。 

  1. for (int i = 0; i < N; ++i) {
  2.     int num = m[i];
  3.     for (int k = 1; num > 0; k <<= 1) {
  4.         int mul = min(k, num);
  5.         for (int j = W; j >= w[i] * mul; --j) {
  6.             f[j] = max(f[j - w[i] * mul] + v[i] * mul, f[j]);
  7.         }
  8.         num -= mul;
  9.     }
  10. }

  弄清楚上面三种背包后分情况就好 

  1. int n;
  2. ll w[MAX_N], v[MAX_N], W;
  3. Pll ps[1 << (MAX_N / 2)]; // (w, v);
  5. void solve() {
  6.     int n2 = n / 2;
  7.     for (int i = 0; i < 1 << n2; ++i) {
  8.         ll sw = 0, sv = 0;
  9.         for (int j = 0; j < n2; ++j)
  10.             if (i >> j & 1) {
  11.                 sw += w[j];
  12.                 sv += v[j];
  13.             }
  14.         ps[i] = Pll(sw, sv);
  15.     }
  16.     sort(ps, ps + (1 << n2));
  17.     int m = 1;
  18.     for (int i = 1; i < 1 << n2; ++i)
  19.         if (ps[m - 1].second < ps[i].second)
  20.             ps[m++] = ps[i];
  22.     ll res = 0;
  23.     for (int i = 0; i < 1 << (n - n2); ++i) {
  24.         ll sw = 0, sv = 0;
  25.         for (int j = 0; j < n - n2; ++j)
  26.             if (i >> j & 1) {
  27.                 sw += w[n2 + j];
  28.                 sv += v[n2 + j];
  29.             }
  30.         if (sw <= W) {
  31.             ll tv = (lower_bound(ps, ps + m, make_pair(W - sw, INF)) - 1)->second;
  32.             res = max(res, sv + tv);
  33.         }
  34.     }
  35.     printf("%lld\n", res);
  36. }

二维费用可由最多取 m 件等方式隐蔽给出。 

  1. for (int k = 0; k < K; ++k)
  2.     for (j = W; j >= 0; --j)
  3.         for (int i = 0; i <= m[k]; ++i)
  4.             f[j] = max(f[j - w[i]] + v[i], f[j]);

显然可以对每组中物品应用完全背包中“一个简单有效的优化”

由 NOIP2006 金明的预算方案引申，对每个附件先做一个 01 背包，再与组件得到一个 个物品组。
更一般问题，依赖关系由「森林」形式给出，涉及到树形 DP 以及泛化物品，这里不表。

更多内容详见「背包九讲」

Maximum Subarray Sum

  1. int max_subarray_sum(int A[], int n) {
  2.     int res, cur;
  3.     if (!A || n <= 0) return 0; 
  4.     res = cur = a[0];
  5.     for (int i = 0; i < n; ++i) {
  6.         if (cur < 0) cur = a[i];
  7.         else cur += a[i];
  8.         res = max(cur, res);
  9.     }
  10.     return res;
  11. }

Set

  1. // 子集枚举
  2. int sub = sup;
  3. do {
  4.     sub = (sub - 1) & sup;
  5. } while (sub != sup); // -1 & sup = sup;
  7. // 势为 k 的集合枚举
  8. int comb = (1 << k) - 1;
  9. while (comb < 1 << n) {
  10.     int x = comb & -comb, y = comb + x;
  11.     comb = ((comb & ~y) / x >> 1) | y;
  12. }
  14. // 排列组合
  15. do {
  17. } while (next_permutation(A, A + N)); // prev_permutation

Bit operation

  1. int __builtin_ffs (unsigned int x)
  2. //返回x的最后一位1的是从后向前第几位，比如7368（1110011001000）返回4。
  3. int __builtin_clz (unsigned int x)
  4. // 返回前导的0的个数。
  5. int __builtin_ctz (unsigned int x)
  6. // 返回后面的0个个数，和__builtiin_clz相对。
  7. int __builtin_popcount (unsigned int x)
  8. // 返回二进制表示中1的个数。
  9. int __builtin_parity (unsigned int x)
  10. // 返回x的奇偶校验位，也就是x的1的个数模2的结果。

Data Structure

  1. // Heap
  2. int heap[MAX_N], sz = 0;
  3. void push(int x) {
  4.     int i = sz++;
  6.     while (i > 0) {
  7.         int p = (i - 1) / 2;
  8.         if (heap[p] <= x) break;
  9.         heap[p] = heap[i];
  10.         i = p;
  11.     }
  12.     heap[i] = x;
  13. }
  14. int pop() {
  15.     int ret = heap[0];
  16.     int x = heap[--sz];
  17.     int i = 0;
  18.     while (i * 2 + 1 < sz) {
  19.         int a = i * 2 + 1, b = i * 2 + 2;
  20.         if (b < sz && heap[b] < heap[a]) a = b;
  21.         if (heap[a] >= x) break;
  22.         heap[i] = heap[a];
  23.         i = a;
  24.     }
  25.     heap[i] = x;
  26.     return ret;
  27. }
  1. // Binary Search Tree
  2. struct node {
  3.     int val;
  4.     node *lch, rch;
  5. };
  7. node *insert(node *p, int x) {
  8.     if (p == NULL) {
  9.         node *q = new node;
  10.         q->val = x;
  11.         q->lch = q->rch = NULL;
  12.         return q;
  13.     } else {
  14.         if (x < p->val) p->lch = insert(p->lch, x);
  15.         else p->rch = insert(p->rch, x);
  16.         return p;
  17.     }
  18. }
  19. bool find(node *p, int x) {
  20.     if (p == NULL) return false;
  21.     else if (x == p->val) return true;
  22.     else if (x < p->val) return find(p->lch, x);
  23.     else return find(p->rch, x);
  24. }
  25. node *remove(node *p, int x) {
  26.     if (p == NULL) return NULL;
  27.     else if (x < p->val) p->lch = remove(p->lch, x);
  28.     else if (x > p->val) p->rch = remove(p->rch, x);
  29.     else if (p->lch == NULL) {
  30.         node *q = p->rch;
  31.         delete p;
  32.         return q;
  33.     } else if (p->lch->rch == NULL) {
  34.         node *q = p->lch;
  35.         q->rch = p->rch;
  36.         delete p;
  37.         return q;
  38.     }  else {
  39.         // 把左儿子子孙中最大的节点提到需要删除的节点上
  40.         node *q;
  41.         for (q = p->lch; q->rch->rch != NULL; q = q->rch);
  42.         node *r = q->rch;
  43.         q->rch = r->lch;
  44.         r->lch = p->lch;
  45.         r->rch = p->rch;
  46.         delete p;
  47.         return r;
  48.     }
  49.     return p;
  50. }
  1. // Union-find Set
  2. int par[MAX_N];
  3. int rnk[MAX_N];
  4. void init(int n) {
  5.     for (int i = 0; i < n; ++i) {
  6.         par[i] = i;
  7.         rnk[i] = 0;
  8.     }
  9. }
  10. int find(int x) {
  11.     return par[x] == x? x : par[x] = find(par[x]);
  12. }
  13. bool same(int x, int y) {
  14.     return find(x) == find(y);
  15. }
  16. void unite(int x, int y) {
  17.     x = find(x);
  18.     y = find(y);
  19.     if (x == y) return;
  20.     if (rnk[x] < rnk[y]) {
  21.         par[x] = y;
  22.     } else {
  23.         par[y] = x;
  24.         if (rnk[x] == rnk[y]) rnk[x]++;
  25.     }
  26. }

当然，更快捷简单的做法，是使用 C++ 的 container。

  1. // Segment Tree
  2. const int MAX_N = 1 << 17;
  3. int n, dat[2 * MAX_N - 1];
  4. void init(int _n) {
  5.     n = 1;
  6.     while (n < _n) n <<= 1;
  7.     for (int i = 0; i < 2 * n - 1; ++i)
  8.         dat[i] = INF;
  9. }
  10. void update(int k, int a) {
  11.     k += n - 1;
  12.     dat[k] = a;
  13.     while (k > 0) {
  14.         k = (k - 1) / 2;
  15.         dat[k] = min(dat[2 * k + 1], dat[2 * k + 2]);
  16.     }
  17. }
  18. // query [a, b), index k in [l, r)
  19. // query(a, b, 0, 0, n)
  20. int query(int a, int b, int k, int l, int r) {
  21.     if (r <= a || b <= l) return INF;
  22.     if (a <= l && r <= b) return dat[k];
  23.     else {
  24.         int v1 = query(a, b, k * 2 + 1, l, (l + r) / 2);
  25.         int v2 = query(a, b, k * 2 + 2, (l + r) / 2, r);
  26.         return min(v1, v2);
  27.     }
  28. }
  1. // RMQ
  2. int n, dat[2 * MAX_N - 1];
  3. void init(int _n) {
  4.     n = 1;
  5.     while (n < _n) n <<= 1;
  6.     for (int i = 0; i < 2 * n - 1; ++i)
  7.         dat[i] = INF;
  8. }
  9. void update(int k, int a) {
  10.     k += n - 1;
  11.     dat[k] = a;
  12.     while (k > 0) {
  13.         k = (k - 1) / 2;
  14.         dat[k] = min(dat[2 * k + 1], dat[2 * k + 2]);
  15.     }
  16. }
  17. // query [a, b), index k in [l, r)
  18. // query(a, b, 0, 0, n)
  19. int query(int a, int b, int k, int l, int r) {
  20.     if (r <= a || b <= l) return INF;
  21.     if (a <= l && r <= b) return dat[k];
  22.     else {
  23.         int v1 = query(a, b, k * 2 + 1, l, (l + r) / 2);
  24.         int v2 = query(a, b, k * 2 + 2, (l + r) / 2, r);
  25.         return min(v1, v2);
  26.     }
  27. }
  1. //Sparse Table
  2. const int MAX_N = 1e5 + 10;
  3. const int MAX_K = 31 - __builtin_clz(MAX_N);
  5. int n, ST[MAX_N][MAX_K + 1], A[MAX_N];
  6. void build(int N) {
  7.     for (int i = 0; i < N; ++i) 
  8.         ST[i][0] = A[i];
  9.     int k = 31 - __builtin_clz(N);
  10.     for (int j = 1; j <= k; ++j)
  11.         for (int i = 0; i <= N - (1 << j); ++i)
  12.             ST[i][j] = min(ST[i][j - 1i], ST[i + (1 << (j - 1))][j - 1]);
  13. }
  15. int query(int l, int r) {
  16.     if (l >= r) return 0;
  17.     int ans = INF, k = 31 - __builtin_clz(r - l);
  18.     for (int j = k; j >= 0; --j)
  19.         if (l + (1 << j) - 1 <= r) {
  20.             ans = min(ans, ST[l][j]);
  21.             l += 1 << j; 
  22.         }
  23.     return ans;
  24. }
  26. int RMQ(int l, int r) {
  27.     if (l >= r) return 0;
  28.     int k = 31 - __builtin_clz(r - l);
  29.     return min(ST[l][k], ST[r - (1 << k)][k]);
  30. }

Graph

  1. struct edge {
  2.     int from;
  3.     int to, dis;
  4. };
  5. vector<edge> G[MAX_V];
  6. vector<edge> es;
  7. bool vis[MAX_V];
  8. int V, E, pre[MAX_V], dist[MAX_V];
  9. // int cost[MAX_V][MAX_V];
  1. // Shortest Way
  2. void dijkstra(int s) {
  3.     priority_queue<Pii, vector<Pii>, greater<Pii> > que;// fisrt 是最短距离，second 是顶点编号
  4.     fill(dist, dist + V, INF);
  5.     dist[s] = 0; que.push(Pii(0, s));
  6.     while (!que.empty()) {
  7.         Pii p = que.top(); que.pop();
  8.         int v = p.second;
  9.         if (dist[v] < p.first) continue;
  10.         for (int i = 0; i < G[v].size(); i++) {
  11.             edge e = G[v][i];
  12.             if (dist[e.to] > dist[v] + e.dis) {
  13.                 dist[e.to] = dist[v] + e.dis;
  14.                 que.push(Pii(dist[e.to], e.to));
  15.             }
  16.         }
  17.     } 
  18. }
  19. void bellman_ford(int s) {
  20.     fill(dist, dist + V, INF);
  21.     dist[s] = 0;
  22.     while (true) {
  23.         bool update = false;
  24.         for (int i = 0; i < E; ++i) {
  25.             edge e = es[i];
  26.             if (dist[e.from] != INF && dist[e.from] + e.dis < dist[e.to]) {
  27.                 update = true;
  28.                 dist[e.to] = dist[e.from] + e.dis;
  29.             }
  30.         }
  31.         if (!update) break;
  32.     }
  33. }
  34. bool find_negative_loop() {
  35.     memset(dist, 0, sizeof dist);
  36.     for (int i = 0; i < V; ++i) 
  37.         for (int j = 0; j < E; ++j) {
  38.             edge e = es[j];
  39.             if (d[e.to] > d[e.from] + e.dis) {
  40.                 d[e.to] = d[e.from] + e.dis;
  41.                 if (i == V - 1) return true;
  42.             }
  43.         }
  44.     return false;
  45. }
  46. void spfa(int s) {
  47.     queue<int> que;
  48.     fill(dist, dist + V, INF);
  49.     fill(vis, vis + V, false);
  50.     dist[s] = 0; que.push(s); vis[s] = true;
  51.     while (!que.empty()) {
  52.         int v = que.front(); que.pop();
  53.         vis[v] = false;
  54.         for (int i = 0; i < G[v].size(); ++i) {
  55.             int u = G[v][i].to;
  56.             if (dist[u] > dist[v] + G[v][i].dis) {
  57.                 dist[u] = dist[v] + G[v][i].dis;
  58.                 if (!vis[u]) {
  59.                     que.push(u);
  60.                     vis[u] = true;
  61.                 }
  62.             }
  63.         }
  64.     }
  65. }
  1. // Spanning Tree
  2. int prime() {
  3.     /*
  4.     fill(dist, dist + V, INF);
  5.     fill(vis, vis + V, false);
  6.     dist[0] = 0;
  7.     int res = 0;
  8.     while (true) {
  9.         int v = -1;
  10.         for (int u = 0; u < V; ++u) {
  11.             if(!vis[u] && (v == -1 || dist[u] < dist[v])) v = u;
  12.         }
  13.         if (v == -1) break;
  14.         vis[v] = true;
  15.         res += dist[v];
  16.         for (int u = 0; u < V; u++)
  17.             dist[u] = min(dist[u], cost[v][u]);
  18.     }
  19.     //*/
  20.     priority_queue<Pii, vector<Pii>, greater<Pii> > que;
  21.     int res = 0;
  22.     fill(dist, dist + V, INF);
  23.     fill(vis, vis + V, false);
  24.     dist[0] = 0;
  25.     que.push(Pii(0, 0));
  26.     while (!que.empty()) {
  27.         Pii p = que.top(); que.pop();
  28.         int v = p.second;
  29.         if (vis[v] || dist[v] < p.first) continue;
  30.         res += dist[v]; vis[v] = true;
  31.         for (int i = 0; i < G[v].size(); ++i) {
  32.             edge e = G[v][i];
  33.             if (dist[e.to] > e.dis) {
  34.                 dist[e.to] = e.dis;
  35.                 que.push(Pii(dist[e.to], e.to));
  36.             }
  37.         }
  38.     }
  39.     return res;
  40. }
  42. bool cmp(const edge e1, const edge e2) {
  43.     return e1.dis < e2.dis;
  44. }
  45. int kruskal() {
  46.     sort(es.begin(), es.end(), cmp);
  47.     init(V);
  48.     int res = 0;
  49.     for (int i = 0; i < E; ++i) {
  50.         edge e = es[i];
  51.         if (!same(e.from, e.to)) {
  52.             unite(e.from, e.to);
  53.             res += e.dis;
  54.         }
  55.     }
  56.     return res;
  57. }
  1. // SCC
  2. int V, cmp[MAX_V];
  3. vector<int> G[MAX_V], rG[MAX_V], vs;
  4. bool used[MAX_V];
  6. void add_edge(int from, int to) {
  7.     G[from].push_back(to); rG[to].push_back(from);
  8. }
  9. void dfs(int v) {
  10.     used[v] = true;
  11.     for (int i = 0; i < G[v].size(); ++i)
  12.         if (!used[G[v][i]]) dfs(G[v][i]);
  13.     vs.push_back(v);
  14. }
  15. void rdfs(int v, int k) {
  16.     used[v] = true; 
  17.     cmp[v] = k;
  18.     for (int i = 0; i < rG[v].size(); ++i) 
  19.         if (!used[rG[v][i]]) rdfs(rG[v][i], k);
  20. }
  21. int scc() {
  22.     memset(used, 0, sizeof used);
  23.     vs.clear();
  24.     for (int v = 0; v < V; ++v)
  25.         if (!used[v]) dfs(v);
  26.     memset(used, 0, sizeof used);
  27.     int k = 0;
  28.     for (int i = vs.size() - 1; i >= 0; --i)
  29.         if (!used[vs[i]]) rdfs(vs[i], k++);
  30.     return k;
  31. }
  1. // tarjan
  2. // flag: in stack or not;
  3. void dfs(int v) {
  4.     int u;
  5.     dfn[v] = low[v] = ++it;
  6.     flag[v] = true;
  7.     st.push(v);
  8.     rep(i, 0, G[v].size()) {
  9.         u = G[v][i];
  10.         if (!dfn[u]) {
  11.             dfs(u);
  12.             low[v] = min(low[v], low[u]);
  13.             if (dfn[v] <= low[u]) { // 割点
  14.             }
  15.             if (dfn[v] < low[u]) { // 割边
  16.             }
  17.         } else if (flag[u])
  18.             low[v] = min(low[v], dfn[u]);
  19.     }
  20.     if (dfn[v] == low[v]) {
  21.         ++cnt;
  22.         do {
  23.             tp[u = st.top()] = cnt;
  24.             flag[u] = false;
  25.             st.pop();
  26.         } while (u != v);
  27.     }
  28. }
  29. void tarjan() {
  30.     cnt = it = 0;
  31.     memset(dfn, 0, sizeof dfn);
  32.     rep(i, 0, V) if (!dfn[i])
  33.         dfs(i);
  34. }
  1. // Bipartite Matching
  2. void add_edge(int u, int v) {
  3.     G[u].push_back(v); G[v].push_back(u);
  4. }
  5. bool dfs(int v) {
  6.     used[v] = true;
  7.     rep(i, 0, G[v].size()) {
  8.         int u = G[v][i], w = match[u];
  9.         if (w < 0 || (!used[w] && dfs(w))) {
  10.             match[v] = u; match[u] = v;
  11.             return true;
  12.         }
  13.     }
  14.     return false;
  15. }
  16. int bipartite_matching() {
  17.     int res = 0;
  18.     memset(match, -1, sizeof match);
  19.     rep (v, 0, V) 
  20.         if (match[v] < 0) {
  21.             memset(used, false, sizeof used);
  22.             if (dfs(v)) ++res;
  23.         }
  24.     return res;
  25. }
  1. // Network Flow
  2. struct edge{
  3.     int to, cap, rev;
  4. };
  5. vector<edge> G[MAX_V];
  6. int level[MAX_V], iter[MAX_V];
  7. void add_edge(int from, int to, int cap) {
  8.     G[from].push_back((edge){to, cap, static_cast<int>(G[to].size())});
  9.     G[to].push_back((edge){from, 0, static_cast<int>(G[from].size() - 1)});
  10. }
  11. // Ford-Fulkerson
  12. int dfs(int v, int t, int f) {
  13.     if (v == t) return f;
  14.     flag[v] = true;
  15.     rep(i, 0, G[v].size()) {
  16.         edge &e = G[v][i];
  17.         if (!flag[e.to] && e.cap > 0) {
  18.             int d = dfs(e.to, t, min(f, e.cap));
  19.             if (d > 0) {
  20.                 e.cap -= d;
  21.                 G[e.to][e.rev].cap += d;
  22.                 return d;
  23.             }
  24.         }
  25.     }
  26.     return 0;
  27. }
  28. int max_flow(int s, int t) {
  29.     int flow = 0;
  30.     for(;;) {
  31.         memset(flag, false, sizeof flag);
  32.         int f = dfs(s, t, INF);
  33.         if (!f) return flow;
  34.         flow += f;
  35.     }
  36. }
  37. // Dinic
  38. void bfs(int s) {
  39.     memset(level, -1, sizeof(level));
  40.     queue<int> que;
  41.     level[s] = 0; que.push(s);
  42.     while (!que.empty()) {
  43.         int v = que.front(); que.pop();
  44.         for (int i = 0; i < G[v].size(); ++i) {
  45.             edge &e = G[v][i];
  46.             if (e.cap > 0 && level[e.to] < 0) {
  47.                 level[e.to] = level[v] + 1;
  48.                 que.push(e.to);
  49.             }
  50.         }
  51.     }
  52. }
  53. int dfs(int v, int t, int f) {
  54.     if (v == t) return f;
  55.     for (int &i = iter[v]; i < G[v].size(); ++i) {
  56.         edge &e = G[v][i];
  57.         if (e.cap > 0 && level[v] < level[e.to]) {
  58.             int d = dfs(e.to, t, min(f, e.cap));
  59.             if (d > 0) {
  60.                 e.cap -= d;
  61.                 G[e.to][e.rev].cap += d;
  62.                 return d;
  63.             }
  64.         }
  65.     }
  66.     return 0;
  67. }
  68. int max_flow(int s, int t) {
  69.     int flow = 0;
  70.     for (;;) {
  71.         bfs(s);
  72.         if (level[t] < 0) return flow;
  73.         memset(iter, 0, sizeof iter);
  74.         int f;
  75.         while ((f = dfs(s, t, INF)) > 0) {
  76.             flow += f;
  77.         }
  78.     }
  79. }
  1. // min_cost_flow
  2. void add_edge(int from, int to, int cap, int cost) {
  3.     G[from].push_back((edge){to, cap, cost, (int)G[to].size()});
  4.     G[to].push_back((edge){from, 0, -cost, (int)G[from].size() - 1});
  5. }
  6. int min_cost_flow(int s, int t, int f) {
  7.     int res = 0;
  8.     fill(h, h + V, 0);
  9.     while (f > 0) {
  10.         priority_queue<Pii, vector<Pii>, greater<Pii> > que;
  11.         fill(dist, dist + V, INF);
  12.         dist[s] = 0; que.push(Pii(0, s));
  13.         while (!que.empty()) {
  14.             Pii p = que.top(); que.pop();
  15.             int v = p.second;
  16.             if (dist[v] < p.first) continue;
  17.             rep(i, 0, G[v].size()) {
  18.                 edge &e = G[v][i];
  19.                 if (e.cap > 0 && dist[e.to] > dist[v] + e.cost + h[v] - h[e.to]) {
  20.                     dist[e.to] = dist[v] + e.cost + h[v] - h[e.to];
  21.                     prevv[e.to] = v;
  22.                     preve[e.to] = i;
  23.                     que.push(Pii(dist[e.to], e.to));
  24.                 }
  25.             }
  26.         }
  27.         if (dist[t] == INF) return -1;
  28.         rep(v, 0, V) h[v] += dist[v];
  29.         int d = f;
  30.         for (int v = t; v != s; v = prevv[v]) 
  31.             d = min(d, G[prevv[v]][preve[v]].cap);
  32.         f -= d;
  33.         res += d * h[t];
  34.         for (int v = t; v != s; v = prevv[v]) {
  35.             edge &e = G[prevv[v]][preve[v]];
  36.             e.cap -= d;
  37.             G[v][e.rev].cap += d;
  38.         }
  39.     }
  40.     return res;
  41. }
  1. // stoer_wagner 全局最小割
  2. void search() {
  3.     memset(vis, false, sizeof vis);
  4.     memset(wet, 0, sizeof wet);
  5.     S = T = -1;
  6.     int imax, tmp;
  7.     rep(i, 0, V) {
  8.         imax = -INF;
  9.         rep(j, 0, V)
  10.             if (!cmb[j] && !vis[j] && wet[j] > imax) {
  11.                 imax = wet[j];
  12.                 tmp = j;
  13.             }
  14.         if (T == tmp) return;
  15.         S = T; T = tmp;
  16.         mc = imax;
  17.         vis[tmp] = true;
  18.         rep(j, 0, V) 
  19.             if (!cmb[j] && !vis[j])
  20.                 wet[j] += G[tmp][j];
  21.     }
  22. }
  23. int stoer_wagner() {
  24.     memset(cmb, false, sizeof cmb);
  25.     int ans = INF;
  26.     rep(i, 0, V - 1) {
  27.         search();
  28.         ans = min(ans, mc);
  29.         if (ans == 0) return 0;
  30.         cmb[T] = true;
  31.         rep(j, 0, V)
  32.             if (!cmb[j]) {
  33.                 G[S][j] += G[T][j];
  34.                 G[j][S] += G[j][T];
  35.             }
  36.     }
  37.     return ans;
  38. }
  1. // 
  2. const int MAX_LOG_V = 32 - __builtin_clz(MAX_V);
  4. vector<int> G[MAX_V];
  5. int root, parent[MAX_LOG_V][MAX_V], depth[MAX_V];
  7. void dfs(int v, int p, int d) {
  8.     parent[0][v] = p;
  9.     depth[v] = d;
  10.     for (int i = 0; i < G[v].size(); i++)
  11.         if (G[v][i] != p) dfs(G[v][i], v, d + 1);
  12. }
  13. void init(int V) {
  14.     dfs(root, -1, 0);
  15.     for (int k = 0; k + 1 < MAX_LOG_V; k++)
  16.         for (int v = 0; v < V; v++)
  17.             if (parent[k][v] < 0) parent[k + 1][v] = -1;
  18.             else parent[k + 1][v] = parent[k][parent[k][v]];
  19. }
  20. int lca(int u, int v) {
  21.     if (depth[u] > depth[v]) swap(u, v);
  22.     for (int k = 0; k < MAX_LOG_V; k++)
  23.         if ((depth[v] - depth[u]) >> k & 1)
  24.             v = parent[k][v];
  25.     if (u == v) return u;
  26.     for (int k = MAX_LOG_V - 1; k >= 0; k--)
  27.         if (parent[k][u] != parent[k][v])
  28.             u = parent[k][u], v = parent[k][v];
  29.     return parent[0][u];
  30. }
  31. // LCA--RMQ
  32. vector<int> G[MAX_V];
  33. int root, vs[MAX_V * 2 - 1], depth[MAX_V * 2 - 1], id[MAX_V];
  35. void dfs(int v, int p, int d, int &k) {
  36.     id[v] = k;
  37.     vs[k] = v;
  38.     depth[k++] = d;
  39.     for (int i = 0; i < G[v].size(); i++) {
  40.         if (G[v][i] != p) {
  41.             dfs(G[v][i], v, d + 1, k);
  42.             vs[k] = v;
  43.             depth[k++] = d;
  44.         }
  45.     }
  46. }
  47. void init(int V) {
  48.     int k = 0;
  49.     dfs(root, -1, 0, k);
  50.     rmq_init(depth, 2 * V - 1);
  51. }
  52. int lca(int u, int v) {
  53.     return vs[query(min(id[u], id[v]), max(id[u], id[v]) + 1)];
  54. }

Computational Geometry

  1. const double EPS = 1e-8;
  2. int sgn(double x) { return x < -EPS ? -1 : x > EPS ? 1 : 0;}
  3. struct Point {
  4.     double x, y;
  5.     Point(double x = 0, double y = 0) : x(x), y(y) {}
  6.     Point operator + (Point p) { return Point(x + p.x, y + p.y); }
  7.     Point operator - (Point p) { return Point(x - p.x, y - p.y); }
  8.     Point operator * (double d) { return Point(x * d, y * d); }
  9.     bool operator < (Point p) {
  10.         return x != p.x? x < p.x : y < p.y;
  11.     }
  12.     double dot(Point p) { return add(x * p.x, y * p.y); }// 内积
  13.     double det(Point p) { return add(x * p.y, -y * p.x);}// 外积
  14. };
  15. bool on_seg(Point p1, Point p2, Point q) {
  16.     return (p1 - q).det(p2 - q) == 0 && (p1 - q).dot(p2 - q) <= 0;
  17. }
  18. P intersection(Point p1, Point p2, Point q1, Point q2) {
  19.     return p1 + (p2 - p1) * ((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1));
  20. }
  21. // 凸包
  22. int convex_hull(Point *ps, int n, Point *ch) {
  23.     sort(ps, ps + n);
  24.     int k = 0;
  25.     for (int i = 0; i < n; ++i) {
  26.         while (k > 1 && (ch[k - 1] - ch[k - 2].det(ps[i] - ch[k - 1])) <= 0) k--;
  27.         ch[k++] = ps[i];
  28.     }
  29.     for (int i = n - 2, t = k; i >= 0; --i) {
  30.         while (k > t && (ch[k - 1] - ch[k - 2].det(ps[i] - ch[k - 1])) <= 0) k--;
  31.         ch[k++] = ps[i];
  32.     }
  33.     return k - 1;
  34. }

Simpson 公式——二次函数近似原函数积分:

Math Problem

  1. // returning count of nk in range [l, r], from Infinity
  2. template<typename T> T mps(T l, T r, T k) { 
  3.     return ((r - (r % k + k) % k) - (l + (k - l % k) % k)) / k + 1;
  4. }
  5. template<typename T> T gcd(T a, T b) {
  6.     //return (b)? gcd(b, a  % b) : a;
  7.     while (b) { T t = a % b; a = b; b = t; } return a;
  8. }
  9. template<typename T> T lcm(T a, T b) { 
  10.     return a / gcd(a, b) * b;
  11. }
  12. // find (x, y) s.t. a x + b y = gcd(a, b) = d
  13. template<typename T> T exgcd(T a, T b, T &x, T &y) { 
  14.     T d = a; 
  15.     if (b) { 
  16.         d = exgcd(b, a % b, y, x); 
  17.         y -= a / b * x; 
  18.     } else {
  19.         x = 1; y = 0; 
  20.     } 
  21.     return d;
  22. }
  23. template<typename T> T modular_linear(T a, T b, T n) {
  24.     T d, e, x, y;
  25.     d = exgcd(a, n, x, y);
  26.     if (b % d)
  27.         return -1;
  28.     e = x * (b / d) % n + n;
  29.     return e % (n / d);
  30. }
  31. template<typename T> T mod_mult(T a, T b, T mod) {
  32.   T res = 0;
  33.   while (b) {
  34.     if (b & 1) {
  35.         res = (res + a) % mod;
  36.         // res += a;
  37.         // if (res >= mod) res -= mod;
  38.     }
  39.     a = (a + a) % mod;
  40.     // a <<= 1;
  41.     // if (a >= mod) a -= mod;
  42.     b >>= 1;
  43.   }
  44.   return res;
  45. }
  46. template<typename T> T mod_pow(T x, T n, T mod) {
  47.     T res = 1;
  48.     while (n) {
  49.         if (n & 1) res = mod_mult(res, x, mod);
  50.         x = mod_mult(x, x, mod);
  51.         n >>= 1;
  52.     }
  53.     return res;
  54.     // return b ? mod_pow(a * a % mod, b >> 1, mod) * (b & 1 ? a : 1) % mod : 1;
  55. }
  56. template<typename T> T mod_inverse(T a, T m) {
  57.     T x, y;
  58.     exgcd(a, m, x, y);
  59.     return (m + x % m) % m;
  60. }
  61. // 线性求逆元 小心乘法溢出？
  62. void init_inverse() {
  63.     inv[1] = 1;
  64.     for (int i = 2; i < MAX_N; i++)
  65.         inv[i] = (MOD - (MOD / i) * inv[MOD % i] % MOD) % MOD;
  66. }
  67. //A[i] * x % M[i] = B[i];
  68. std::pair<int, int> linear_congruence(const std::vector<int> &A, const std::vector<int> &B, const std::vector<int> &M) {
  69.         int x = 0, m = 1;
  70.         for(int i = 0; i < A.size(); i++) {
  71.                 int a = A[i] * m, b = B[i] - A[i] * x, d = gcd(M[i], a);
  72.                 if(b % d != 0) return std::make_pair(0, -1); // no solutioin
  73.                 int t = b / d * mod_inverse(a / d, M[i] / d) % (M[i] / d);
  74.                 x = x + m * t;
  75.                 m *= M[i] / d;
  76.         }
  77.         x = (x + m) % m;
  78.         return std::make_pair(x % m, m);
  79. }u
  80. ll CRT(vector<ll> &a, vector<ll> &m) {
  81.     ll M = 1LL, res = 0;
  82.     for (int i = 0; i < m.size(); ++i)
  83.         M *= m[i];
  84.     for (int i = 0; i < m.size(); ++i) {
  85.         ll Mi, Ti;
  86.         Mi = M / m[i]; Ti = mod_inverse(Mi, mi);
  87.         res = (res + a[i] * (Mi * Ti % M) % M) % M;
  88.     }
  89.     return res;
  90. }
  91. // only for MOD < 1e6;
  92. ll fact[MOD + 10];
  93. void init() {
  94.     fact[0] = 1;
  95.     for (int i = 1; i <= MOD; ++i)
  96.         fact[i] = fact[i - 1] * i % MOD;
  97. }
  98. int mod_fact(int n, int p, int &e) {
  99.     e = 0;
  100.     if (n == 0) return 1;
  101.     int res = mod_fact(n / p, p, e);
  102.     e += n / p;
  103.     if (n / p % 2 != 0) return res * (p - fact[n % p]) % p;
  104.     return res * fact[n % p] % p;
  105. }
  106. int mod_comb(int n, int k, int p) {
  107.     if (n < 0 || k < 0 || n < k) return 0;
  108.     if (n == 0) return 1;
  109.     int e1, e2, e3;
  110.     int a1 = mod_fact(n, p, e1), a2 = mod_fact(k, p, e2), a3 = mod_fact(n - k, p, e3);
  111.     if (e1 > e2 + e3) return 0;
  112.     return a1 * mod_inverse(a2 * a3 % p, p) % p;
  113. }
  114. ll lucas(ll n, ll k, const ll &p) {
  115.     if (n < 0 || k < 0 || n < k) return 0;
  116.     if (n == 0) return 1;
  117.     return lucas(n / p, k / p, p) * mod_comb(n % p, k % p, p) % p;
  118. }
  1. // 矩阵快速幂
  2. typedef vector<int> vec;
  3. typedef vector<vec> mat;
  4. mat G(MAX_N);
  6. mat mat_mul(mat &A, mat &B) {
  7.     mat C(A.size(), vec(B[0].size()));
  8.     for (int i = 0; i < A.size(); ++i)
  9.         for (int k = 0; k < B.size(); ++k)
  10.             for (int j = 0; j < B[0].size(); ++j)
  11.                 C[i][j] = (C[i][j] + A[i][k] % MOD * B[k][j] % MOD + MOD) % MOD;
  12.     return C;
  13. }
  14. mat mat_pow(mat A, ll n) {
  15.     mat B(A.size(), vec(A.size()));
  16.     for (int i = 0; i < A.size(); ++i)
  17.         B[i][i] = 1;
  18.     while (n > 0) {
  19.         if (n & 1) B = mat_mul(B, A);
  20.         A = mat_mul(A, A);
  21.         n >>= 1;
  22.     }
  23.     return B;
  24. }
  1. // prime number
  2. bool is_prime(int n) {
  3.     for (int i = 2; i * i <= n; ++i)
  4.         if (n % i == 0) return false;
  5.     return n != 1;
  6. }
  7. vector<int> divisor(int n) {
  8.     vector<int> res;
  9.     for (int i = 1; i * i <= n; ++i) {
  10.         if (n % i == 0) {
  11.             res.push_back(i);
  12.             if (i != n / i) res.push_back(n / i);
  13.         }
  14.     }
  15.     return res;
  16. }
  17. map<int, int> prime_factor(int n) {
  18.     map<int, int> res;
  19.     for (int i = 2; i * i <= n; ++i) {
  20.         while (n % i == 0) {
  21.             ++res[i];
  22.             n /= i;
  23.         }
  24.     }
  25.     if (n != 1) res[n] = 1;
  26.     return res;
  27. }
  28. int prime[MAX_N];
  29. bool isPrime[MAX_N + 1];
  30. int seive(int n) {
  31.     int p = 0;
  32.     fill(isPrime, isPrime + n + 1, true);
  33.     isPrime[0] = isPrime[1] = false;
  34.     for (int i = 2; i <= n; ++i) 
  35.         if (isPrime[i]) {
  36.             prime[p++] = i;
  37.             for (int j = 2 * i; j <= n; j += i) isPrime[j] = false;
  38.         }
  39.     return p;
  40. }
  41. // the number of prime in [L, r)
  42. // 对区间 [l, r）内的整数执行筛法，prime[i - l] = true <=> i 是素数
  43. bool segPrimeSmall[MAX_L];
  44. bool segPrime[MAX_SQRT_R];
  45. void segment_sieve(ll l, ll r) {
  46.     for (int i = 0; (ll)i * i < r; ++i) segPrimeSmall[i] = true;
  47.     for (int i = 0; i < r - l; ++i) segPrime[i] = true;
  48.     for (int i = 2; (ll)i * i < r; ++i) {
  49.         if (segPrimeSmall[i]) {
  50.             for (int j = 2 * i; (ll)j * j <= r; j += i) segPrimeSmall[j] = false;
  51.             for (ll j = max(2ll, (l + i - 1) / i) * i; j < r; j += i) segPrime[j - l] = false;
  52.         }
  53.     }
  54. }
  55. // Miller_Rabin
  56. bool check(ll a, ll n, ll x, ll t) {
  57.     ll res = mod_pow(a, x, n);
  58.     ll last = res;
  59.     for (int i = 1; i <= t; ++i) {
  60.         res = mod_mult(res, res, n);
  61.         if (res == 1 && last != 1 && last != n - 1) return true;
  62.         last = res;
  63.     }
  64.     if (res != 1) return true;
  65.     return false;
  66. }
  67. bool Miller_Rabin(ll n) {
  68.     if (n < MAX_N) return isPrime[n]; // small number may get wrong answer?!
  69.     if (n < 2) return false;
  70.     if (n == 2) return true;
  71.     if ((n & 1) == 0) return false;
  72.     ll x = n - 1, t = 0;
  73.     while ((x & 1) == 0) {
  74.         x >>= 1;
  75.         ++t;
  76.     }
  77.     for (int i = 0; i < S; ++i) {
  78.         ll a = rand() % (n - 1) + 1;
  79.         if (check(a, n, x, t))
  80.             return false;
  81.     }
  82.     return true;
  83. }
  84. // find factors
  85. vector<ll> factor;
  86. ll Pollard_rho(ll x, ll c) {
  87.     ll i = 1, k = 2;
  88.     ll x0 = rand() % x;
  89.     ll y = x0;
  90.     while (true) {
  91.         ++i;
  92.         x0 = (mod_mult(x0, x0, x) + c) % x;
  93.         ll d;
  94.         if (y == x0) d = 1;
  95.         else 
  96.             if (y > x0)
  97.                 d = gcd(y - x0, x);
  98.             else d = gcd(x0 - y, x);
  99.         if (d != 1 && d != x) return d;
  100.         if (y == x0) return x;
  101.         if (i == k) {
  102.             y = x0;
  103.             k += k;
  104.         }
  105.     }
  106. }
  107. void find_factor(ll n) {
  108.     if (n == 1) return ;
  109.     if (Miller_Rabin(n)) {
  110.         factor.push_back(n);
  111.         return ;
  112.     }
  113.     ll p = n;
  114.     while (p >= n) p = Pollard_rho(p, rand() % (n - 1) + 1);
  115.     find_factor(p);
  116.     find_factor(n / p);
  117. }
  1. #include<bits/stdc++>
  2. //Meisell-Lehmer
  3. const int MAX_N = 5e6 + 2;
  4. bool np[MAX_N];
  5. int prime[MAX_N], pi[MAX_N];
  6. int getprime()
  7. {
  8.     int cnt = 0;
  9.     np[0] = np[1] = true;
  10.     pi[0] = pi[1] = 0;
  11.     for(int i = 2; i < MAX_N; ++i)
  12.     {
  13.         if(!np[i]) prime[++cnt] = i;
  14.         pi[i] = cnt;
  15.         for(int j = 1; j <= cnt && i * prime[j] < MAX_N; ++j)
  16.         {
  17.             np[i * prime[j]] = true;
  18.             if(i % prime[j] == 0)   break;
  19.         }
  20.     }
  21.     return cnt;
  22. }
  23. const int M = 7;
  24. const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
  25. int phi[PM + 1][M + 1], sz[M + 1];
  26. void init() {
  27.     getprime();
  28.     sz[0] = 1;
  29.     for(int i = 0; i <= PM; ++i)  phi[i][0] = i;
  30.     for(int i = 1; i <= M; ++i) {
  31.         sz[i] = prime[i] * sz[i - 1];
  32.         for(int j = 1; j <= PM; ++j) phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
  33.     }
  34. }
  35. int sqrt2(ll x) {
  36.     ll r = (ll)sqrt(x - 0.1);
  37.     while(r * r <= x)   ++r;
  38.     return int(r - 1);
  39. }
  40. int sqrt3(ll x) {
  41.     ll r = (ll)cbrt(x - 0.1);
  42.     while(r * r * r <= x)   ++r;
  43.     return int(r - 1);
  44. }
  45. ll getphi(ll x, int s)
  46. {
  47.     if(s == 0)  return x;
  48.     if(s <= M)  return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
  49.     if(x <= prime[s]*prime[s])   return pi[x] - s + 1;
  50.     if(x <= prime[s]*prime[s]*prime[s] && x < MAX_N) {
  51.         int s2x = pi[sqrt2(x)];
  52.         ll ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
  53.         for(int i = s + 1; i <= s2x; ++i) ans += pi[x / prime[i]];
  54.         return ans;
  55.     }
  56.     return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
  57. }
  58. ll getpi(ll x) {
  59.     if(x < MAX_N)   return pi[x];
  60.     ll ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
  61.     for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) ans -= getpi(x / prime[i]) - i + 1;
  62.     return ans;
  63. }
  64. ll lehmer_pi(ll x) {
  65.     if(x < MAX_N)   return pi[x];
  66.     int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
  67.     int b = (int)lehmer_pi(sqrt2(x));
  68.     int c = (int)lehmer_pi(sqrt3(x));
  69.     ll sum = getphi(x, a) +(ll)(b + a - 2) * (b - a + 1) / 2;
  70.     for (int i = a + 1; i <= b; i++) {
  71.         ll w = x / prime[i];
  72.         sum -= lehmer_pi(w);
  73.         if (i > c) continue;
  74.         ll lim = lehmer_pi(sqrt2(w));
  75.         for (int j = i; j <= lim; j++) sum -= lehmer_pi(w / prime[j]) - (j - 1);
  76.     }
  77.     return sum;
  78. }
  79. int main() {
  80.     init();
  81.     ll n;
  82.     while(~scanf("%lld",&n))
  83.     {
  84.         printf("%lld\n",lehmer_pi(n));
  85.     }
  86.     return 0;
  87. }
  1. // 欧拉函数
  2. int euler_phi(int n) {
  3.     int res = n;
  4.     for (int i = 2; i * i <= n; ++i) {
  5.         if (n % i == 0) {
  6.             res = res / i * (i - 1);
  7.             for (; n % i == 0; n /= i);
  8.         }
  9.     }
  10.     if (n != 1) res = res / n * (n - 1);
  11.     return res;
  12. }
  13. int euler[MAX_N];
  14. void euler_phi_sieve() {
  15.     for (int i = 0; i < MAX_N; ++i) euler[i] = i;
  16.     for (int i = 2; i < MAX_N; ++i)
  17.         if (euler[i] == i)
  18.             for (int j = i; j < MAX_N; j += i) euler[j] = euler[j] / i * (i - 1);
  19. }

对于 函数，有如下性质：

  1. int mu[MAX_N];
  2. void moebius() {
  3.     int cnt = 0; mu[1] = 1;
  4.     memset(vis, 0, sizeof vis);
  5.     for (int i = 2; i < MAX_N; ++i) {
  6.         if (!vis[i]) {
  7.             prime[cnt++] = i;
  8.             mu[i] = -1;
  9.         }
  10.         for (int j = 0; j < cnt && i * prime[j] < MAX_N; ++j) {
  11.             vis[i * prime[j]] = true;
  12.             if (i % prime[j])
  13.                 mu[i * prime[j]] = -mu[i];
  14.             else 
  15.                 mu[i * prime[j]] = 0, break;
  16.         }
  17.     }
  18. }
  20. map<int, int> moebius(int n) {
  21.     map<int, int> res;
  22.     vector<int> primes;
  23.     for (int i = 2; i * i <= n; ++i) {
  24.         if (n % i == 0) {
  25.             primes.push_back(i);
  26.             while (n % i == 0) n /= i;
  27.         }
  28.     }
  29.     if (n != 1) primes.push_back(n);
  31.     int m = primes.size();
  32.     for (int i = 0; i < (1 << m); ++i) {
  33.         int mu = 1, d = 1;
  34.         for (int j = 0; j < m; ++j) {
  35.             if (i >> j & 1) {
  36.                 mu *= -1;
  37.                 d *= primes[j];
  38.             }
  39.         }
  40.         res[d] = mu;
  41.     }
  42.     return res;
  43. }
  1. // Guass_jordan
  2. const double eps = 1e-8;
  3. typedef vector<double> vec;
  4. typedef vector<vec> mat;
  6. vec gauss_joedan(const mat &A, const vec& b) {
  7.     int n = A.size();
  8.     mat B(n, vec(n + 1));
  9.     for (int i = 0; i < n; ++i)
  10.         for (int j = 0; j < n; ++j) B[i][j] = A[i][j];
  11.     for (int i = 0; i < n; ++i) B[i][n] = b[i];
  13.     for (int i = 0; i < n; ++i) {
  14.         int pivot = i;
  15.         for (int j = i; j < n; ++j)
  16.             if (abs(B[j][i]) > abs(B[pivot][i])) pivot = j;
  17.         if (i != pivot) swap(B[i], B[pivot]);
  19.         if (abs(B[i][i]) < eps) return vec();
  21.         for (int j = i + 1; j <= n; ++j) B[i][j] /= B[i][i];
  22.         for (int j = 0; j < n; ++j) if (i != j)
  23.             for (int k = i + 1; k <= n; ++k) B[j][k] -= B[j][i] * B[i][k];
  24.     }
  26.     vec x(n);
  27.     for (int i = 0; i < n; ++i) x[i] = B[i][n];
  28.     return x;
  29. }
  31. vec gauss_joedan_xor(const mat& A, const vec& b) {
  32.     int n = A.size();
  33.     mat B(n, vec(n + 1));
  34.     for (int i = 0; i < n; ++i)
  35.         for (int j = 0; j < n; ++j) B[i][j] = A[i][j];
  36.     for (int i = 0; i < n; ++i) B[i][n] = b[i];
  38.     for (int i = 0; i < n; ++i) {
  39.         int pivot = i;
  40.         for (int j = i; j < n; ++j)
  41.             if (B[j][i]) { 
  42.                 pivot = j;
  43.                 break;
  44.             }
  45.         if (pivot != i) swap(B[i], B[pivot]);
  47.         for (int j = 0; j < n; ++j) if (i != j && B[j][i])
  48.                 for (int k = i + 1; k <= n; ++k) B[j][k] ^= B[i][k];
  49.             }
  50.     }
  52.     vec x(n);
  53.     for (int i = 0; i < n; ++i) x[i] = B[i][n];
  54.     return x;
  55. }

String

  1. Hash
  2. KMP
  3. Extend KMP
  4. trie树 poj2001 2503 3630 1056 hdu 1075 1251 1247 1298 1671
  5. Manacher 算法
  6. AC自动机
  7. 后缀数组
  8. 后缀树
  9. 后缀自动机
  10. 回文自动机
  1. // 最小最大表示法：
  2. int getMinString(const string &s) {  
  3.     int len = (int)s.length();
  4.     int i = 0, j = 1, k = 0;  
  5.     while(i < len && j < len && k < len) {
  6.         int t = s[(i+k)%len] - s[(j+k)%len];
  7.         if(t == 0) k++;
  8.         else {
  9.             if(t > 0) i += k + 1;//getMaxString: t < 0
  10.             else j += k + 1;
  11.             if(i==j) j++;
  12.             k = 0;
  13.         }
  14.     }
  15.     return min(i,j);
  16. }
  1. // KMP
  2. int nxt[MAX_N];
  3. void getNext(const string &str) {
  4.     int len = str.length();
  5.     int j = 0, k;
  6.     k = nxt[0] = -1;
  7.     while (j < len) {
  8.         if (k == -1 || str[j] == str[k]) 
  9.             nxt[++j] = ++k;
  10.         else k = nxt[k];
  11.     }
  12. }
  13. int kmp(const string &tar, const string &pat) {
  14.     getNext(pat);
  15.     int num, j, k;
  16.     int lenT = tar.length(), lenP = pat.length();
  17.     num = j = k = 0;
  18.     while (j < lenT) {
  19.         if(k == -1 || tar[j] == pat[k])
  20.             j++, k++;
  21.         else k = nxt[k];
  22.         if(k == lenP) {
  23.             // res = max(res, j - lenP);
  24.             k = nxt[k];
  25.             ++num;
  26.         }
  27.     }
  28.     return num;//lenP - res - 1;
  29. }
  1. // AC 自动机
  2. int ans[MAX_N], d[MAX_N];
  4. struct Trie {
  5.     int nxt[MAX_N][26], fail[MAX_N], end[MAX_N];
  6.     int root, L;
  7.     int newnode() {
  8.         for(int i = 0; i < 26; i++)
  9.             nxt[L][i] = -1;
  10.         end[L++] = 0;
  11.         return L-1;
  12.     }
  13.     void init() {
  14.         L = 0;
  15.         root = newnode();
  16.     }
  17.     void insert(char buf[]) {
  18.         int len = strlen(buf);
  19.         int now = root;
  20.         for(int i = 0; i < len; i++) {
  21.             if(nxt[now][buf[i]-'a'] == -1)
  22.                 nxt[now][buf[i]-'a'] = newnode();
  23.             now = nxt[now][buf[i]-'a'];
  24.         }
  25.         end[now] = 1;
  26.         d[now] = len;
  27.     }
  28.     void build() {
  29.         queue<int> Q;
  30.         fail[root] = root;
  31.         for(int i = 0; i < 26; i++)
  32.             if(nxt[root][i] == -1)
  33.                 nxt[root][i] = root;
  34.             else {
  35.                 fail[nxt[root][i]] = root;
  36.                 Q.push(nxt[root][i]);
  37.             }
  38.         while( !Q.empty() ) {
  39.             int now = Q.front(); Q.pop();
  40.             for(int i = 0; i < 26; i++)
  41.                 if(nxt[now][i] == -1)
  42.                     nxt[now][i] = nxt[fail[now]][i];
  43.                 else {
  44.                     fail[nxt[now][i]] = nxt[fail[now]][i];
  45.                     Q.push(nxt[now][i]);
  46.                 }
  47.         }
  48.     }
  49.     void solve(char buf[]) {
  50.         int cur = root;
  51.         int len = strlen(buf);
  52.         int index;
  53.         for(int i = 0; i < len; ++i) {
  54.             if(buf[i] >= 'A' && buf[i] <= 'Z')
  55.                 index = buf[i] - 'A';
  56.             else if(buf[i] >= 'a' && buf[i] <= 'z')
  57.                 index = buf[i] - 'a';
  58.             else continue;
  59.             cur = nxt[cur][index];
  60.             int x = cur;
  61.             while(x != root) {
  62.                 if(end[x]) {
  63.                     ans[i + 1] -= 1;
  64.                     ans[i - d[x] + 1] += 1;
  65.                     break;
  66.                 }
  67.                 x = fail[x];
  68.             }
  69.         }
  70.     }
  71. };
  73. Trie ac;
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注