Template of Simplex Algorithm

任务

给定一个$n$维向量$a$、一个$m \times n$的矩阵$b$和一个$m$维向量$c$。要求构造一个$n$维向量$x$,使得在满足
$$
\forall i \in [1, m], \; \sum_{j=1}^n b_{i, j}x_j \le c_i
$$
的前提下
$$
\sum_{i=1}^n a_ix_i
$$
最大。

接口

int solve()

返回$-1$表示无解,$1$表示无界,$0$表示存在最大值。

输入

int n, m

如题。

double a[][]

$(m+1) \times (n+1)$的矩阵。第$0$行表示$a$向量,第$0$列表示$c$向量,第$1 \ldots m$行、$1 \ldots n$列表示$b$矩阵。

输出

double ans[]

solve()返回$0$时,表示$n$维向量$x$。此时的最大值为a[0][0]的相反数。

代码

/*
 * Life is like a diaper -- short and loaded.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>

static int const N = 45;
static long double const EPS = 1e-8;
int t;

int pm(long double const &n) {
  return (n < 0 ? -n : n) < EPS ? 0 : n < 0 ? -1 : 1;
}

class Array {
private:
  long double arr[N];

public:
  int size;

long double &operator[](int const &n) { return arr[n]; }

long double max() { return *std::max_element(arr + 1, arr + size + 1); }

void operator/=(long double n) {
    for (int i = 0; i <= size; ++i)
      arr[i] /= n;
  }

void add(Array a, long double n) {
    for (int i = 0; i <= size; ++i)
      arr[i] += a[i] * n;
  }
};

class Simplex {
private:
  int id[N];
  Array tmp;

void pivot(int e, int v) {
    a[e] /= a[e][v];
    for (int i = 0; i <= m; ++i)
      if (i != e)
        a[i].add(a[e], -a[i][v]);
    id[e] = v;
  }

int simplex(int n, int m) {
    for (; pm(a[0].max()) > 0;) {
      int col = -1, row = -1;
      for (int i = 1; i <= n; ++i)
        if (pm(a[0][i]) > 0) {
          col = i;
          break;
        }
      long double mn = 1e100, tmp;
      for (int i = 1; i <= m; ++i)
        if (pm(a[i][col]) > 0 && ((tmp = pm(a[i][0] / a[i][col] - mn)) < 0 ||
                                  tmp == 0 && id[i] < id[row]))
          mn = a[i][0] / a[i][col], row = i;
      if (row == -1)
        return 1;
      pivot(row, col);
    }
    return 0;
  }

public:
  int n, m;
  Array ans, a[N];

int solve() {
    for (int i = 1; i <= m; ++i)
      a[i][n + i] = 1, id[i] = n + i;
    int row = -1;
    long double mn = 1e100;
    for (int i = 1; i <= m; ++i)
      if (a[i][0] < mn)
        mn = a[i][0], row = i;
    if (pm(mn) < 0) {
      for (int i = 1; i <= n; ++i)
        tmp[i] = a[0][i], a[0][i] = 0;
      for (int i = 0; i <= m; ++i)
        a[i][n + m + 1] = -1;
      for (int i = 0; i <= m; ++i)
        a[i].size = n + m + 1;
      pivot(row, n + m + 1), simplex(n + m + 1, m);
      if (pm(a[0][0]) != 0)
        return -1;
      for (int i = 1; i <= m; ++i)
        if (id[i] == n + m + 1) {
          for (int j = 1; j <= n + m; ++j)
            if (pm(a[i][j]) != 0) {
              pivot(i, j);
              break;
            }
          break;
        }
      for (int i = 1; i <= n; ++i)
        a[0][i] = tmp[i];
      for (int i = 1; i <= m; ++i)
        a[0].add(a[i], -a[0][id[i]]);
    }
    for (int i = 0; i <= m; ++i)
      a[i].size = n + m;
    int ret = simplex(n + m, m);
    if (!ret)
      for (int i = 1; i <= m; ++i)
        if (id[i] <= n)
          ans[id[i]] = a[i][0];
    return ret;
  }
} solver;

int main() {
  scanf("%d%d%d", &solver.n, &solver.m, &t);
  for (int i = 1; i <= solver.n; ++i)
    scanf("%Lf", &solver.a[0][i]);
  for (int i = 1; i <= solver.m; scanf("%Lf", &solver.a[i++][0]))
    for (int j = 1; j <= solver.n; ++j)
      scanf("%Lf", &solver.a[i][j]);
  int res = solver.solve();
  if (res == -1)
    puts("Infeasible");
  else if (res == 1)
    puts("Unbounded");
  else {
    printf("%.10Lf\n", -solver.a[0][0]);
    if (t == 1) {
      for (int i = 1; i <= solver.n; ++i)
        printf("%.10Lf ", solver.ans[i]);
      puts("");
    }
  }
  return 0;
}

Hide

题意概述

有一棵$N$个节点的树,每个节点都是黑色或白色。有$Q$次操作,动态修改点的颜色,询问最远黑点对之间的距离。
数据范围:$1 \le N \le 10^5, \; 1 \le Q \le 5 \times 10^5$。

算法分析

考虑点分治,构造一棵重心树。为了方便,下面所有距离均表示在原树上的距离,堆均为大根堆。
设节点$i$在重心树上的父亲为$p_i$。对于每个节点$u$维护两个堆,第一个堆用来维护重心树上$u$子树中的所有黑点到$p_u$的距离,第二个堆用来维护重心树上$u$各个子树中距离$u$最远的黑点到$u$的距离。那么第二个堆可以借助第一个堆来维护。
接着需要一个全局的堆,用来维护所有节点第二个堆的前两个元素之和(经过该点的路径)以及所有黑点第二个堆的第一个元素(以该点为端点的路径)。
修改时,沿着重心树往上依次处理经过的节点,先在全局的堆中删除,再在第二个堆中删除,在第一个堆中修改后,维护第二个堆和全局的堆。

代码

/*
 * Computer programmers do it byte by byte.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <queue>

template <typename T> void read(T &n) {
  char c;
  for (; (c = getchar()) < '0' || c > '9';)
    ;
  for (n = c - '0'; (c = getchar()) >= '0' && c <= '9'; (n *= 10) += c - '0')
    ;
}

typedef int const ic;

static ic N = 100005;

class heap {
private:
  std::priority_queue<int> que, buf;

void clear() {
    for (; !que.empty() && !buf.empty() && que.top() == buf.top();
         que.pop(), buf.pop())
      ;
  }

public:
  void push(ic &n) {
    if (n)
      que.push(n);
  }

void erase(ic &n) {
    if (n)
      buf.push(n);
  }

void pop() {
    clear();
    if (!que.empty())
      que.pop();
  }

int top() {
    clear();
    return que.empty() ? 0 : que.top();
  }

int top2() {
    clear();
    if (que.empty())
      return 0;
    int tmp = que.top(), ret = tmp;
    que.pop(), clear();
    if (que.empty()) {
      que.push(tmp);
      return 0;
    }
    ret += que.top(), que.push(tmp);
    return ret;
  }
} global;

namespace tree {
int n, nume, h[N], col[N];
int tim, fi[N], dep[N], lca[N << 1][20];
int power[N << 1];
int root, up[N], size[N], mx[N], vis[N];
heap toup[N], me[N];
struct Edge {
  int v, nxt;
} e[N << 1];

void add_edge(ic &u, ic &v) {
  e[++nume] = (Edge){v, h[u]}, h[u] = nume;
  e[++nume] = (Edge){u, h[v]}, h[v] = nume;
}

void dfs(ic &t, ic &fa) {
  lca[++tim][0] = t, fi[t] = tim, dep[t] = dep[fa] + 1;
  for (int i = h[t]; i; i = e[i].nxt)
    if (e[i].v != fa)
      dfs(e[i].v, t), lca[++tim][0] = t;
}

void init() {
  for (int i = 2; i <= tim; i <<= 1)
    ++power[i];
  for (int i = 1; i <= tim; ++i)
    power[i] += power[i - 1];
  for (int i = 1; i < 20; ++i)
    for (int j = 1; j <= tim; ++j) {
      ic k = std::min(tim, j + (1 << i - 1));
      lca[j][i] = dep[lca[j][i - 1]] < dep[lca[k][i - 1]] ? lca[j][i - 1]
                                                          : lca[k][i - 1];
    }
}

int get_lca(ic &u, ic &v) {
  ic l = std::min(fi[u], fi[v]), r = std::max(fi[u], fi[v]);
  ic len = power[r - l + 1], k = r - (1 << len) + 1;
  return dep[lca[l][len]] < dep[lca[k][len]] ? lca[l][len] : lca[k][len];
}

int get_dist(ic &u, ic &v) {
  ic lca = get_lca(u, v);
  return dep[u] + dep[v] - 2 * dep[lca];
}

int get_size(ic &t, ic &fa) {
  int ret = 1;
  for (int i = h[t]; i; i = e[i].nxt)
    if (!vis[e[i].v] && e[i].v != fa)
      ret += get_size(e[i].v, t);
  return ret;
}

void get_root(ic &t, ic &fa, ic &tot) {
  size[t] = 1, mx[t] = 0;
  for (int i = h[t]; i; i = e[i].nxt)
    if (!vis[e[i].v] && e[i].v != fa)
      get_root(e[i].v, t, tot), size[t] += size[e[i].v],
          mx[t] = std::max(mx[t], size[e[i].v]);
  (mx[t] = std::max(mx[t], tot - size[t])) < mx[root] && (root = t);
}

void init_heap(ic &t, ic &fa, ic &dep, heap &hp) {
  hp.push(dep);
  for (int i = h[t]; i; i = e[i].nxt)
    if (!vis[e[i].v] && e[i].v != fa)
      init_heap(e[i].v, t, dep + 1, hp);
}

void build(int t) {
  vis[t] = true;
  for (int i = h[t]; i; i = e[i].nxt)
    if (!vis[e[i].v])
      root = 0,
      get_root(e[i].v, 0,
               size[e[i].v] < size[t] ? size[e[i].v] : get_size(e[i].v, t)),
      up[root] = t, init_heap(e[i].v, t, 1, toup[root]),
      me[t].push(toup[root].top()), build(root);
  global.push(me[t].top()), global.push(me[t].top2());
}

void modify(int t) {
  ic p = t;
  if (col[t])
    global.erase(me[t].top());
  else
    global.push(me[t].top());
  for (; up[t]; t = up[t]) {
    if (col[up[t]])
      global.erase(me[up[t]].top());
    global.erase(me[up[t]].top2());
    me[up[t]].erase(toup[t].top());
    if (col[p])
      toup[t].erase(get_dist(p, up[t]));
    else
      toup[t].push(get_dist(p, up[t]));
    me[up[t]].push(toup[t].top());
    if (col[up[t]])
      global.push(me[up[t]].top());
    global.push(me[up[t]].top2());
  }
  col[p] = !col[p];
}
} // namespace tree

int main() {
  read(tree::n), tree::mx[0] = tree::n;
  for (int i = 1, u, v; i < tree::n; ++i)
    read(u), read(v), tree::add_edge(u, v);
  tree::dfs(1, 0), tree::init(), tree::get_root(1, 0, tree::n),
      tree::build(tree::root);
  for (int i = 1; i <= tree::n; ++i)
    tree::col[i] = 1;
  int q;
  for (read(q); q--;) {
    char c;
    scanf(" %c", &c);
    if (c == 'G')
      printf("%d\n", global.top());
    else {
      int p;
      read(p), tree::modify(p);
    }
  }
  return 0;
}

Count on a Tree II

题目描述

You are given a tree with $N$ nodes. The tree nodes are numbered from $1$ to $N$. Each node has an integer weight.
We will ask you to perform the following operation:

  • $u$ $v$: ask for how many different integers that represent the weight of nodes there are on the path from $u$ to $v$.

题意概述

给定一棵$N$个节点的树,每个节点都有权值。有$M$次询问,每次询问给出$u, v$,求从$u$到$v$的路径上有多少种不同的权值。
数据范围:$1 \le N \le 40000, \; 1 \le M \le 10^5$。

算法分析

统计权值种类数,容易想到莫队,但这是在一棵树上,所以需要稍微做些转化。
求出树的括号序(DFS,访问和离开节点时均记录),那么每个询问都可以转化为对括号序上某个区间的询问。但是,如果区间内的某个节点出现了两次,那么这个节点不应该对答案产生贡献。因此,算法执行时可以用两个数组维护,一个表示节点的出现次数,另一个表示权值的出现次数。
具体来讲,令$s_i, e_i$分别表示访问、离开节点$i$的时间。假设$s_u \le s_v$。若$u$是$v$祖先,那么就相当于询问区间$[s_u, s_v]$;否则,相当于询问区间$[e_u, s_v]$,但此时$u, v$的LCA并没有被计算到答案中,因此需要把它加上。

代码

/*
 * Trap full -- please empty.
 */

#include <map>
#include <cstdio>
#include <cstring>
#include <algorithm>

template <typename T>
void read(T &n) {
  char c; for (; (c = getchar()) < '0' || c > '9'; ) ;
  for (n = c - '0'; (c = getchar()) >= '0' && c <= '9'; (n *= 10) += c - '0') ;
}

typedef int const ic;
typedef long long ll;

static ic N = 40005;
static ic M = 100005;
static ic K = 300;
std::map <int, int> num;
int nume, h[N], w[N], base[N], seq[N << 1];
int tim, st[N], ed[N], dep[N], up[N][16];
int mans, ml, mr, mvis[N], mrec[N], ans[M];
struct Edge {
  int v, nxt;
} e[N << 1];
struct Query {
  int l, r, lca, id;
  bool operator < (const Query &q) const {
    return l / K == q.l / K ? r < q.r : l / K < q.l / K;
  }
} q[M];

void add_edge(ic &u, ic &v) {
  e[++ nume] = (Edge) { v, h[u] }, h[u] = nume;
  e[++ nume] = (Edge) { u, h[v] }, h[v] = nume;
}

void dfs(ic &t, ic &fa) {
  seq[st[t] = ++ tim] = t, dep[t] = dep[fa] + 1, up[t][0] = fa;
  for (int i = h[t]; i; i = e[i].nxt)
    if (e[i].v != fa) dfs(e[i].v, t);
  seq[ed[t] = ++ tim] = t;
}

int get_lca(int u, int v) {
  if (dep[u] > dep[v]) std::swap(u, v);
  for (int i = 15; ~ i; -- i)
    if (dep[up[v][i]] >= dep[u]) v = up[v][i];
  if (u == v) return u;
  for (int i = 15; ~ i; -- i)
    if (up[u][i] != up[v][i]) u = up[u][i], v = up[v][i];
  return up[u][0];
}

void add(ic &u) {
  if (mvis[u]) {
    -- mrec[w[u]], mvis[u] = 0;
    if (! mrec[w[u]]) -- mans;
  } else {
    if (! mrec[w[u]]) ++ mans;
    ++ mrec[w[u]], mvis[u] = 1;
  }
}

void get(ic &l, ic &r) {
  for (; ml < l;) add(seq[ml ++]);
  for (; ml > l;) add(seq[-- ml]);
  for (; mr < r;) add(seq[++ mr]);
  for (; mr > r;) add(seq[mr --]);
}

int main() {
  int n, m; read(n), read(m);
  for (int i = 1; i <= n; ++ i) {
    read(w[i]);
    if (! num.count(w[i])) num[w[i]] = num.size();
    w[i] = num[w[i]];
  }
  for (int i = 1, u, v; i < n; ++ i) read(u), read(v), add_edge(u, v);
  dfs(1, 0);
  for (int i = 1; i < 16; ++ i)
    for (int j = 1; j <= n; ++ j) up[j][i] = up[up[j][i - 1]][i - 1];
  for (int i = 1, u, v; i <= m; ++ i) {
    read(u), read(v), q[i].id = i;
    if (st[u] > st[v]) std::swap(u, v);
    if (ed[u] > ed[v]) q[i].l = st[u] + 1, q[i].r = st[v], q[i].lca = u;
    else q[i].l = ed[u], q[i].r = st[v], q[i].lca = get_lca(u, v);
  }
  std::sort(q + 1, q + m + 1);
  mans = ml = mr = mvis[1] = mrec[w[1]] = 1;
  for (int i = 1; i <= m; ++ i) get(q[i].l, q[i].r), add(q[i].lca), ans[q[i].id] = mans, add(q[i].lca);
  for (int i = 1; i <= m; ++ i) printf("%d\n", ans[i]);
  return 0;
}

Little Pony and Lord Tirek

题目描述

Lord Tirek is a centaur and the main antagonist in the season four finale episodes in the series “My Little Pony: Friendship Is Magic”. In “Twilight’s Kingdom” (Part 1), Tirek escapes from Tartarus and drains magic from ponies to grow stronger.
The core skill of Tirek is called Absorb Mana. It takes all mana from a magic creature and gives them to the caster.
Now to simplify the problem, assume you have $n$ ponies (numbered from $1$ to $n$). Each pony has three attributes:

  • $s_i$: amount of mana that the pony has at time $0$;
  • $m_i$: maximum mana that the pony can have;
  • $r_i$: mana regeneration per unit time.

Lord Tirek will do $m$ instructions, each of them can be described with three integers: $t_i, l_i, r_i$. The instruction means that at time $t_i$, Tirek will use Absorb Mana on ponies with numbers from $l_i$ to $r_i$ (both borders inclusive). We’ll give you all the m instructions in order, count how much mana Tirek absorbs for each instruction.

题意概述

有$n$匹马,第$i$匹马初始时有$s_i$点魔法,魔法上限为$m_i$,每秒增加$r_i$点魔法。有$m$次操作,第$i$次操作在第$t_i$秒进行,将区间$[l_i, r_i]$中马的魔法清零。求每次操作清除了多少点魔法。
数据范围:$1 \le n, m \le 10^5, \; 0 \le s_i \le m_i \le 10^5, \; 0 \le r_i \le 10^5, \; 0 \le t_i \le 10^9$。

算法分析

显然,答案只与两次询问的时间差有关。
用set维护一些不相交的区间(它们的并是整个序列),每段区间内所有马的上次操作时间相同。每次操作就相当于切割一些区间然后合并一些区间。由于每次至多切割两个区间,因此总共操作的区间个数是$O(n)$的。
考虑区间询问。每匹马都可以表示为一个分段函数,我们用主席树来维护,第$i$棵线段树表示经过$i$秒后每匹马的函数。因为有些马有初始魔法,所以需要建两棵主席树,一棵表示第$i$匹马的初始魔法为$s_i$,另一棵表示所有马的初始魔法为$0$。
那么只要预处理出主席树,对于每次操作,在set中找到其对应的一些区间,在主席树上询问答案,并更新区间。

代码

/*
 * So, is the glass half empty, half full, or just twice as large as it needs to
 * be?
 */

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <set>

template <typename T> void read(T &n) {
  char c;
  for (; (c = getchar()) < '0' || c > '9';)
    ;
  for (n = c - '0'; (c = getchar()) >= '0' && c <= '9'; (n *= 10) += c - '0')
    ;
}

typedef long long ll;
typedef int const ic;

static ic N = 100005;
int numn, root[N << 1];
struct Pony {
  int s, m, r;
} p[N];
struct Event {
  int t, p, b;
} e[N];
struct Query {
  int t, l, r;
  bool operator<(Query const &a) const { return l < a.l; }
} q[N];
struct SegmentTree {
  int child[2];
  ll a, b;
} node[N << 6];
std::set<Query> st;

void update(ic &root) {
  if (!root)
    return;
  node[root].a = node[node[root].child[0]].a + node[node[root].child[1]].a;
  node[root].b = node[node[root].child[0]].b + node[node[root].child[1]].b;
}

void modify(int &root, ic &nl, ic &nr, ic &p, ic &a, ic &b, ic &flag) {
  if (p > nr || nl > p)
    return;
  if (flag)
    node[++numn] = node[root], root = numn;
  else if (!root)
    root = ++numn;
  if (nl == nr) {
    node[root].a = a, node[root].b = b;
    return;
  }
  ic mid = (nl + nr) >> 1;
  modify(node[root].child[0], nl, mid, p, a, b, flag);
  modify(node[root].child[1], mid + 1, nr, p, a, b, flag);
  update(root);
}

ll query(ic &root, ic &nl, ic &nr, ic &t, ic &l, ic &r) {
  if (!root || l > nr || nl > r)
    return 0;
  if (l <= nl && nr <= r)
    return node[root].a * t + node[root].b;
  ic mid = (nl + nr) >> 1;
  return query(node[root].child[0], nl, mid, t, l, r) +
         query(node[root].child[1], mid + 1, nr, t, l, r);
}

int main() {
  int n, m;
  read(n);
  for (int i = 0; i < n; ++i)
    read(p[i].s), read(p[i].m), read(p[i].r);
  read(m);
  for (int i = 0; i < m; ++i)
    read(q[i].t), read(q[i].l), read(q[i].r), --q[i].l, --q[i].r;
  for (int i = 0; i < n; ++i)
    modify(root[0], 0, n - 1, i, p[i].r, 0, 0),
        e[i] = (Event){p[i].r ? p[i].m / p[i].r + 1 : N, i, p[i].m};
  std::sort(e, e + n, [](Event const &a, Event const &b) { return a.t < b.t; });
  for (int i = 1, p = 0; i < N; ++i) {
    root[i] = root[i - 1];
    for (; p < n && e[p].t == i; ++p)
      modify(root[i], 0, n - 1, e[p].p, 0, e[p].b, 1);
  }
  for (int i = 0; i < n; ++i)
    modify(root[N], 0, n - 1, i, p[i].r, p[i].s, 0),
        e[i] = (Event){p[i].r ? (p[i].m - p[i].s) / p[i].r + 1 : N, i, p[i].m};
  std::sort(e, e + n, [](Event const &a, Event const &b) { return a.t < b.t; });
  for (int i = 1, p = 0; i < N; ++i) {
    root[N + i] = root[N + i - 1];
    for (; p < n && e[p].t == i; ++p)
      modify(root[N + i], 0, n - 1, e[p].p, 0, e[p].b, 1);
  }
  st.insert((Query){0, 0, n - 1});
  for (int i = 0; i < m; ++i) {
    auto rec = st.lower_bound((Query){q[i].t, q[i].l, q[i].r});
    if (rec == st.end())
      --rec;
    for (; rec != st.begin() && rec->r >= q[i].l; --rec)
      ;
    if (rec->r < q[i].l)
      ++rec;
    auto tmp = rec;
    Query l = (Query){rec->t, rec->l, q[i].l - 1};
    ll ans = 0;
    for (; rec != st.end() && rec->l <= q[i].r; ++rec)
      ans += rec->t ? query(root[std::min(q[i].t - rec->t, N - 1)], 0, n - 1,
                            q[i].t - rec->t, std::max(rec->l, q[i].l),
                            std::min(rec->r, q[i].r))
                    : query(root[std::min(q[i].t - rec->t, N - 1) + N], 0,
                            n - 1, q[i].t - rec->t, std::max(rec->l, q[i].l),
                            std::min(rec->r, q[i].r));
    printf("%lld\n", ans), --rec;
    Query r = (Query){rec->t, q[i].r + 1, rec->r};
    for (; tmp != st.end() && tmp->l <= q[i].r;)
      st.erase(tmp++);
    st.insert((Query){q[i].t, q[i].l, q[i].r});
    if (l.l <= l.r)
      st.insert(l);
    if (r.l <= r.r)
      st.insert(r);
  }
  return 0;
}