Keywords: SSSP | reweighting technique
Time Complexity: \( O(Q(N+M)) \)

Explanation

最小費用流や、負辺を含む全点対最短経路[1]などで用いられるポテンシャルを用いたreweightテクニックを使う。

Reweight

グラフ \(G\) における単一始点最短経路長を \(d(u)\) とする。このとき \(w'( (u,v))=w( (u,v))-d(u)+d(v)\) と置き直したグラフ \(G'\) における最短経路と \(G\) における最短経路は一致する。

元グラフ\(G\)の重みが\(w_2\)に変わった（非減少に変化）とき、新たな最短経路長を計算したいとする。このグラフの最短経路長を \(d_2(u)\) としたとき、先程と同様に \(w_2'( (u,v))=w_2( (u,v))-d(u)+d(v)\) と変形したとしても、やはり最短経路は変わらない。

最短経路が変わらないというのは、\(d_2'(u)=d_2(u)+d(u) \) だからである。\(-d(u)+d(v)\)の補正が最短経路長に与える影響が経路に依存していない（これがポテンシャルと呼ばれる由来だろう）ことから、最短経路が変化しないことが言える。

この問題を解く。reweighted graphにおける最短経路長はすべて 0 である。ここで重み増加クエリが来たとしても reweighted graph の最短経路は \(n\) でバウンドできる。つまり \(n\) 個 queue を持っておけば線形 dijkstra が可能である。

TLEしたため、type 2のクエリをやや溜めている。

#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>

using namespace std;

constexpr int N = 1e5;
int n, m, q;
int us[N], vs[N], ww[N];
vector<int> g[N];
long long dist[N];
int dist2[N];
vector<int> que[N*2];

void dijkstra() {
  for (int i = 0; i < n; i++) {
    dist[i] = 1e16;
  }
  priority_queue<pair<long long, int>> q;
  q.emplace(0, 0);
  dist[0] = 0;
  while (!q.empty()) {
    long long d = -q.top().first;
    int u = q.top().second;
    q.pop();
    if (dist[u] < d) continue;
    for (int id : g[u]) {
      int v = vs[id];
      if (dist[v] > dist[u] + ww[id]) {
        dist[v] = dist[u] + ww[id];
        q.emplace(-dist[v], v);
      }
    }
  }
}

int cap;

void dijkstra2() {
  if (cap == 0) return;
  que[0].push_back(0);
  for (int i = 0; i < n; i++) {
    dist2[i] = cap;
  }
  dist2[0] = 0;
  for (int d = 0; d < cap; d++) {
    while (!que[d].empty()) {
      int u = que[d].back(); que[d].pop_back();
      if (dist2[u] < d) continue;
      for (int id : g[u]) {
        int v = vs[id];
        int w = min<long long>(cap, ww[id] + dist[u] - dist[v]);
        if (dist2[v] > dist2[u] + w && dist2[u] + w < cap) {
          dist2[v] = dist2[u] + w;
          que[dist2[v]].push_back(v);
        }
      }
    }
  }
  for (int i = 0; i < n; i++) {
    dist[i] += dist2[i];
  }
  cap = 0;
}

int main() {
  cin >> n >> m >> q;
  for (int i = 0; i < m; i++) {
    scanf("%d %d %d", &us[i], &vs[i], &ww[i]);
    us[i]--;
    vs[i]--;
    g[us[i]].push_back(i);
  }
  dijkstra();
  while (q--) {
    int type;
    scanf("%d", &type);
    if (type == 1) {
      dijkstra2();
      int v;
      scanf("%d", &v);
      if (dist[v - 1] >= 1e16) {
        puts("-1");
      } else {
        printf("%lld\n", dist[v - 1]);
      }
    } else {
      int c;
      scanf("%d", &c);
      while (c--) {
        int l;
        scanf("%d", &l);
        ww[l - 1]++;
      }
      cap += n;
      if (cap == n*2) dijkstra2();
    }
  }
}

ポテンシャルは何でも良いが、ポテンシャルに最短経路長を用いることで三角不等式により全辺の重みを非負にすることができる。APSP では SSSP を n 回行うという方法が考えられるが、負辺が存在する場合に Dijkstra 法が使えない。そこで Bellman-Ford 法を用いて負辺を除去し、Dijkstra 法が使える状態に落ちし込む。これが Johnson's algorithm である。

References

[1]
https://en.wikipedia.org/wiki/Johnson%27s_algorithm. Retrieved 2017-08-29

2017-08-21

Goldman Sachs CodeSprint: Transaction Certificates

keywords: rolling hash | birthday attack

解法

ハッシュ値が重複するまでひたすら生成し続けるだけで高速に発見できることが知られている。これは誕生日攻撃と呼ばれている。nが小さい場合は重複値が存在しない可能性があるので、100以上になるまで大きくしておく。

#include <iostream>
#include <algorithm>
#include <map>
#include <vector>

using namespace std;

int main() {
  int n, k, p;
  unsigned long long m;
  cin >> n >> k >> p >> m;

  int u = 1;
  while (n * u < 100) {
    u++;
  }
  n *= u;

  auto gen = [&](int seed) {
    unsigned long long x = seed;
    vector<int> val;
    for (int i = 0; i < n; i++) {
      x ^= x << 13;
      x ^= x >> 7;
      x ^= x << 17;
      val.push_back(x % k + 1);
    }
    return val;
  };

  auto calc_hash = [&](vector<int> a) {
    long long hash = 0;
    for (long long x : a) {
      hash = (hash * p + x) & (m - 1);
    }
    return hash;
  };

  map<unsigned, int> mp;
  
  for (int i = 1;; i++) {
    long long hash = calc_hash(gen(i));
    if (mp.count(hash)) {
      for (long long x : gen(i)) {
        cout << x << ' ';
      }
      cout << endl;
      for (long long x : gen(mp[hash])) {
        cout << x << ' ';
      }
      return 0;
    } else {
      mp[hash] = i;
    }
  }
}

二冪なので他にも攻略法があって、それはhos氏のブログが参考になる。

文字列 \(110000\) と文字列 \(000011\) の rolling hash はそれぞれ

\begin{equation}
1x^5+1x^4+0x^3+0x^2+0x^1+0 \\
0x^5+0x^4+0x^3+0x^2+1x^1+1
\end{equation}

であり、その差は \(x^5+x^4-x-1\) である。この多項式は \(x\) が奇数のとき常に \(32\) の倍数になる。つまり法を \(32\)、基数を奇数としたとき必ず衝突する。法が二冪であればこのような多項式が容易に構成できる。

\(v(n)\) を \(2\) が \(n\) を割り切る回数とし、多項式 \(f\) に対して

\begin{equation}
V(f)=\min_{n \in 2\mathbb{Z}+1}v(f(n))
\end{equation}

とする。\(V(f) \ge 32\) となる多項式 \(f\) を見つければ良い。\(V(x^{2k} - 1) = V( (x^k - 1) (x^k+1) ) \ge 1+V(x^k - 1)\) なので、

\begin{equation}
V( (x-1)(x^2-1)(x^4-1)\cdots(x^{128}-1) ) \ge 32
\end{equation}

である。これは展開時の係数がすべて \(\pm 1\) であり、さらにその係数は組合せ的に求まる。

n, k, p, m = map(int, input().split())
# (x-1)(x^2-1)(x^4-1)(x^8-1)(x^16-1)(x^32-1)(x^64-1)(x^128-1)
#  1    2      3      4      5       6       7      8        >32
n = (256+n-1)//n*n
ans = [0]*(n-256) + [(-1)**bin(i).count('1') for i in range(256)]

mp1 = {0: 1, 1: 2, -1: 1}
mp2 = {0: 1, 1: 1, -1: 2}
print(*[mp1[v] for v in ans])
print(*[mp2[v] for v in ans])

補足

よく考えると基数が偶数 \((p=2)\) のときバグる。
\(x\) が奇数のとき \(x^2+1\) は \(4k+2\) 型の整数になるので、\(v(x^2+1)=1\) である。これにより厳密に \( v(x^{2^k} - 1)=v(x^2-1)+k \) である。この結果は LTE の補題に似ているが、LTE の補題は偶素数と奇素数で振る舞いがわずかに異なることに注意。
対二冪の多項式の係数列を Thue-Morse sequence というらしいですね。

2017-08-18

August Long Challenge 2017: Flower Pots

CHTのO(BN^2)解を最初投げたんだけど通らなかった。解法は難しくないにも関わらず通している人が少ないのは集団心理か何かだろうか。

計算量: \(O(BN^2)\)
キーワード: DP

解説

まずDPで解くことを考える。状態は着火区間と経過秒数とする。\( [L,R] \rightarrow [L',R'] \)に遷移するパターンは3つある。

\(L\)から\(L'\)に伝播、\(R\)から\(R'\)に伝播
\(L\)から\(L'\)と\(R'\)に伝播
\(R\)から\(L'\)と\(R'\)に伝播

この問題を解く鍵は遷移2or3は複数回起きないということである。最適構造の雰囲気を図1に示す。

f:id:pekempey:20170817154309p:plain — 図1. 最適構造

最適解にならないケースを図2に示す。\( C\to y\)に直接繋げた方がコストは小さくなる。

f:id:pekempey:20170817161026p:plain — 図2. 最適解にならないケース

ということでDPしなくても分岐地点を全探索すれば良いことがわかる。\(x\)を固定したとき、そのコストを知るには \( dist_b(1,\ast),dist_b(n,\ast),dist_b(c,\ast) \) が必要である――何手掛かったかを知る必要があるため下付きの b として表現した――。これは多少手を抜いたO(BN^2)のDPで計算したとしても間に合う。

起点 \(x\) を決めたとき、伝播先(l,r)はO(N)通りしかないことに注意したい。

#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;

void upd(long long &x, long long y) {
  x = min(x, y);
}

constexpr int N = 3000;
int n, b, c;
long long a[N];
long long ll[N][N], rr[N][N], cc[N][N];

void solve() {
  cin >> n >> b >> c;
  for (int i = 0; i < n; i++) {
    cin >> a[i];
  }
  c--;

  auto dist = [&](int i, int j) {
    return (a[i] - a[j]) * (a[i] - a[j]);
  };

  constexpr long long inf = 1.1e16;
  fill_n(*ll, N*N, inf);
  fill_n(*rr, N*N, inf);
  fill_n(*cc, N*N, inf);
  ll[0][0] = 0;
  rr[0][n - 1] = 0;
  cc[0][c] = 0;
  for (int i = 0; i + 1 < b; i++) {
    for (int j = 0; j < n; j++) {
      for (int k = 0; k < n; k++) {
        upd(ll[i + 1][k], ll[i][j] + dist(k, j));
        upd(rr[i + 1][k], rr[i][j] + dist(k, j));
        upd(cc[i + 1][k], cc[i][j] + dist(k, j));
      }
    }
  }

  long long ans = inf;
  for (int i = 0; i < n; i++) {
    int l = i;
    int r = i;
    while (l != 0 || r != n - 1) {
      if (l == 0) {
        r++;
      } else if (r == n - 1) {
        l--;
      } else if (a[i] - a[l - 1] == a[r + 1] - a[i]) {
        l--;
        r++;
      } else if (a[i] - a[l - 1] < a[r + 1] - a[i]) {
        l--;
      } else {
        r++;
      }
      for (int j = 0; j < b; j++) {
        upd(ans, cc[j][i] + max(dist(i, l), dist(i, r)) + ll[b - j - 1][l] + rr[b - j - 1][r]);
      }
    }
  }
  cout << ans << endl;
}

int main() {
  int T;
  cin >> T;
  while (T--) solve();
}

2017-08-18

August Long Challenge 2017: Walks on the binary tree

計算量: \(O(Q \log^2 N)\)
キーワード: persistent segment tree | zobrist hash | binary counter | binary search on segment tree | lcp

考察

これはそもそもbinary stringの問題である。
文字列数が2であればlen(s)+len(t)-lcp(s,t)で計算できる。複数文字列があるならソートして隣接要素間のLCPをすべて引けば良い。
sortはLCPがあればいいので、LCPの方法を考える。
ハッシュを用いてLCPを求めることにする。zobrist hash化したものをsegtreeに載せておけば、LCPは二分探索でいける。
二進カウンタと同じ考え方で、更新クエリで書き換わる桁数はO(1) amortizedである。文字列をsegtreeで持つことにしていたので、永続により空間O(Q log N)で保持できるということになる。
segtree上の[0,k)型クエリの二分探索は工夫することでO(log^2 n)からO(log N)に落とせる。fenwick tree上の二分探索(lower_bound)として有名？
先読みソート+linked listを用いて動的更新に対応させた。

#include <iostream>
#include <algorithm>
#include <vector>
#include <cassert>
#include <string>
 
using namespace std;
 
struct node {
  node *l = nullptr;
  node *r = nullptr;
  uint64_t sum = 0;
};
 
constexpr int N = 1 << 17;
 
uint64_t rnd[N];
 
uint64_t sum(node *x) {
  return x ? x->sum : 0ULL;
}
 
node *build(int l, int r) {
  node *ret = new node();
  if (r - l == 1) return ret;
  ret->l = build(l, l + r >> 1);
  ret->r = build(l + r >> 1, r);
  return ret;
}
 
int getval(int k, node *x, int l = 0, int r = N) {
  if (r - l == 1) {
    return x->sum != 0;
  }
  int m = l + r >> 1;
  if (k < m) {
    return getval(k, x->l, l, m);
  } else {
    return getval(k, x->r, m, r);
  }
}
 
node *setval(int k, int v, node *x, int l = 0, int r = N) {
  x = x ? new node(*x) : new node();
  if (r - l == 1) {
    x->sum = v ? rnd[l] : 0;
    return x;
  }
  int m = l + r >> 1;
  if (k < m) {
    x->l = setval(k, v, x->l, l, m);
  } else {
    x->r = setval(k, v, x->r, m, r);
  }
  x->sum = sum(x->l) ^ sum(x->r);
  return x;
}
 
int lcp(node *x, node *y, int k = N) {
  if (k == 1) return 0;
  if (sum(x->l) == sum(y->l)) {
    return lcp(x->r, y->r, k / 2) + k / 2;
  } else {
    return lcp(x->l, y->l, k / 2);
  }
}
 
bool compare(node *x, node *y) {
  int l = lcp(x, y);
  if (l == N - 1) return false;
  return getval(l, x) < getval(l, y);
}
 
void solve() {
  int n, q;
  cin >> n >> q;
 
  node *curr = build(0, N);
 
  vector<node *> strs(q);
 
  vector<bool> query(q);
  for (int i = 0; i < q; i++) {
    char type;
    scanf(" %c", &type);
    if (type == '!') {
      int k;
      scanf("%d", &k);
      k = n - 1 - k;
      for (; k >= 0; k--) {
        bool b = getval(k, curr);
        curr = setval(k, !b, curr);
        if (!b) break;
      }
    } else {
      query[i] = true;
    }
    strs[i] = curr;
  }
 
  vector<int> perm(q);
  for (int i = 0; i < q; i++) {
    perm[i] = i;
  }
 
  sort(perm.begin(), perm.end(), [&](int i, int j) {
    return compare(strs[i], strs[j]);
  });
 
  vector<int> inv(q);
  for (int i = 0; i < q; i++) {
    inv[perm[i]] = i;
  }
 
  long long ans = 1LL * q * n + 1;
  for (int ii = 0; ii + 1 < q; ii++) {
    int i = perm[ii];
    int j = perm[ii + 1];
    ans -= min(n, lcp(strs[i], strs[j]));
  }
 
  vector<int> left(q, -1);
  vector<int> right(q, -1);
 
  for (int i = 0; i + 1 < q; i++) {
    left[i + 1] = i;
    right[i] = i + 1;
  }
 
  auto del = [&](int i) {
    int y = inv[i];
    int x = left[y];
    int z = right[y];
    if (x != -1) ans += min(n, lcp(strs[perm[x]], strs[perm[y]]));
    if (z != -1) ans += min(n, lcp(strs[perm[y]], strs[perm[z]]));
    if (x != -1 && z != -1) ans -= min(n, lcp(strs[perm[x]], strs[perm[z]]));
    if (x != -1) right[x] = z;
    if (z != -1) left[z] = x;
    ans -= n;
  };
  for (int i = 0; i < q; i++) {
    if (query[i]) {
      del(i);
    }
  }
 
  vector<long long> anss;
  for (int i = q - 1; i >= 0; i--) {
    if (query[i]) {
      anss.push_back(ans);
    } else {
      del(i);
    }
  }
 
  for (int i = (int)anss.size() - 1; i >= 0; i--) {
    printf("%lld\n", anss[i]);
  }
}
 
uint64_t xorshift() {
  static uint64_t x = 1234567;
  x ^= x << 13;
  x ^= x >> 7;
  x ^= x << 17;
  return x;
}
 
int main() {
  for (int i = 0; i < N; i++) {
    rnd[i] = xorshift();
  }
 
  int T;
  cin >> T;
  while (T--) {
    solve();
  }
}

2017-08-18

August Long Challenge 2017: Hill Jumping

正答数を見る限りだと正攻法ではないんだろなとは思う。

計算量: \(O(100 Q \log N)\)
キーワード: link-cut tree | envelope

考察

各点におけるジャンプ先をsucc[i]とする。変更クエリにより書き換わるsuccの要素数はたかだか200である。
succ[i]で繋いでいくと木になる。この問題は木上のkth-ancestorを求める問題と見做せる。
link-cut treeによりlink, cut, kth-ancestorはいずれもO(log N) amortizedで処理できる。
kth-ancestorは二分探索木のkth-minと同様

#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <deque>
#include <cassert>
 
using namespace std;
 
constexpr int N = 1e5 + 10;
 
struct node {
  node *l = nullptr;
  node *r = nullptr;
  node *p = nullptr;
  int size = 1;
  int key;
};
 
bool is_root(node *x) {
  return !x->p || (x != x->p->l && x != x->p->r);
}
 
int size(node *x) {
  return x ? x->size : 0;
}
 
void fix(node *x) {
  x->size = 1 + size(x->l) + size(x->r);
}
 
void rot(node *x) {
  node *y = x->p;
  node *z = y->p;
  if (z) {
    if (y == z->l) z->l = x;
    if (y == z->r) z->r = x;
  }
  x->p = z;
  y->p = x;
  if (x == y->l) {
    y->l = x->r;
    x->r = y;
    if (y->l) y->l->p = y;
  } else {
    y->r = x->l;
    x->l = y;
    if (y->r) y->r->p = y;
  }
  fix(y);
}
 
void splay(node *x) {
  while (!is_root(x)) {
    node *y = x->p;
    if (!is_root(y)) {
      node *z = y->p;
      rot((x == y->l) == (y == z->l) ? y : x);
    }
    rot(x);
  }
  fix(x);
}
 
void expose(node *x) {
  node *tmp = x;
  node *r = nullptr;
  while (x) {
    splay(x);
    x->r = r;
    fix(x);
    r = x;
    x = x->p;
  }
  splay(tmp);
}
 
void link(node *c, node *p) {
  expose(c);
  expose(p);
  p->r = c;
  c->p = p;
  fix(p);
}
 
void cut(node *x) {
  expose(x);
  x->l->p = nullptr;
  x->l = nullptr;
  fix(x);
}
 
int kth_element(node *x, int k) {
  expose(x);
  while (true) {
    if (size(x->r) == k) {
      break;
    } else if (k < size(x->r)) {
      x = x->r;
    } else {
      if (!x->l) {
        break;
      }
      k -= 1 + size(x->r);
      x = x->l;
    }
  }
  splay(x);
  return x->key;
}
 
long long ft[N];
 
void add(int k, int v) {
  for (k++; k < N; k += k & -k) {
    ft[k] += v;
  }
}
 
long long getval(int k) {
  long long ret = 0;
  for (k++; k > 0; k &= k - 1) {
    ret += ft[k];
  }
  return ret;
}
 
int main() {
  int n, q;
  cin >> n >> q;
  for (int i = 0; i < n; i++) {
    int a;
    scanf("%d", &a);
    add(i, a);
    add(i + 1, -a);
  }
 
  vector<int> succ(n, -1);
  for (int i = 0; i < n; i++) {
    succ[i] = i;
  }
  vector<node *> nodes(n);
  for (int i = 0; i < n; i++) {
    nodes[i] = new node();
    nodes[i]->key = i;
  }
 
  auto update = [&](int l, int r) {
    l = max(0, l);
    deque<int> mono;
    for (int i = min(n - 1, r + 100); i >= r + 1; i--) {
      while (!mono.empty() && getval(i) >= getval(mono.front())) {
        mono.pop_front();
      }
      mono.push_front(i);
    }
    for (int i = r; i >= l; i--) {
      while (!mono.empty() && mono.back() > i + 100) {
        mono.pop_back();
      }
      while (!mono.empty() && getval(i) >= getval(mono.front())) {
        mono.pop_front();
      }
      mono.push_front(i);
      if (mono.size() == 1) {
        if (succ[i] != i) {
          cut(nodes[i]);
        }
        succ[i] = i;
      } else {
        if (succ[i] != i) {
          cut(nodes[i]);
        }
        succ[i] = mono[1];
        link(nodes[i], nodes[mono[1]]);
      }
    }
  };
  update(0, n - 1);
 
  while (q--) {
    int type;
    scanf("%d", &type);
    if (type == 1) {
      int u, cnt;
      scanf("%d %d", &u, &cnt);
      u--;
      printf("%d\n", kth_element(nodes[u], cnt) + 1);
    } else {
      int l, r, x;
      scanf("%d %d %d", &l, &r, &x);
      l--;
      r--;
      add(l, x);
      add(r + 1, -x);
      update(r - 99, r);
      update(l - 100, l - 1);
    }
  }
}

pekempey's blog

AIM Tech Round 4 (Div. 1) D. Dynamic Shortest Path

Explanation

Reweight

References

Goldman Sachs CodeSprint: Transaction Certificates

解法

補足

August Long Challenge 2017: Flower Pots

解説

August Long Challenge 2017: Walks on the binary tree

考察

August Long Challenge 2017: Hill Jumping

考察