dyck 列（正しい括弧列）の数え上げ問題に、深さ上限を与えた問題。dyck 列の総数はカタラン数として知られる。解法の本質はカタラン数の二項係数表現の導出と同じで、エラーが起きた地点から移動方向を反転させるとよい。カタラン数については分かりやすいサイトが見つからなかったため参考リンクは用意していないが、『数学ガール』の 1 巻に書いてあった説明が分かりやすかった覚えがある（読んだの高 1 くらいなので、具体的に何が書かれてたか覚えてないけど）。定番ネタなので探せばいくらでもあるだろう。

たとえば左のような移動経路は、右の移動経路に対応させる。E (エラー地点）から移動方向を反転させている。

        *                 
 *     * G       *        
S-*---*---  <=> S-*-*------     
   E *             E *    
    *                 *   
                       * G'
                        *

S→Gへの経路のうち、下側エラーが起きるパターン数は、S から G' への移動経路数と一致することが言える。ここまではカタラン数の計算方法として有名だろう。

今回は下側だけでなく上側も制約がつく。上側エラーも同様に計算できるが、下側と上側で二重カウントしてしまう。そのため包除原理的に重複を除去していく。

下側エラー→上側エラーとなるパターン数と上側エラー→下側エラーとなるパターン数もまた、先程と同様のテクニックにより計算できる（二回反転させる）。下側エラー→上側エラー→下側エラーというのも同様に計算できる（三階反転させる）。これらが計算できれば、包除原理的に解ける。下→上→下→上→…というのを永遠に求め続ける必要はなく、m の制約からして有限個である（m/n個程度であり、これは計算量に影響を与える）

計算量はわかりづらいが、O(n+m) になる。競技中には計算量解析が雑で、n≤500 のとき別の解法を取るような O( (n+m) sqrt(n) ) の解法を設計していたが、その必要はなかった。

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

using namespace std;

constexpr int mod = 1e9 + 7;

struct modint {
  int n;
  modint(int n = 0) : n(n) {}
};

modint operator+(modint a, modint b) { return modint((a.n += b.n) >= mod ? a.n - mod : a.n); }
modint operator-(modint a, modint b) { return modint((a.n -= b.n) < 0 ? a.n + mod : a.n); }
modint operator*(modint a, modint b) { return modint(1LL * a.n * b.n % mod); }
modint &operator+=(modint &a, modint b) { return a = a + b; }
modint &operator-=(modint &a, modint b) { return a = a - b; }
modint &operator*=(modint &a, modint b) { return a = a * b; }

modint fact[505050];
modint ifact[505050];
modint inv[505050];

modint C(int n, int r) {
  if (n < 0 || r < 0 || n < r) return 0;
  return fact[n] * ifact[r] * ifact[n - r];
}

modint g(int w, int y) {
  y = abs(y);
  if (w < y) return 0;
  if ((w - y) % 2 != 0) return 0;
  return C(w, (w - y) / 2);
}

int mirror(int y, int u) {
  return -(y - u) + u;
}

modint f(int w, int h, int y) {
  modint ret;
  ret += g(w, y);
  int y1 = y;
  int y2 = y;
  int d1 = -1;
  int d2 = h + 1;
  for (int ii = 0; min(-d1, d2) <= w; ii++) {
    y1 = mirror(y1, d1);
    y2 = mirror(y2, d2);
    if (ii & 1) {
      ret += g(w, y1);
      ret += g(w, y2);
    } else {
      ret -= g(w, y1);
      ret -= g(w, y2);
    }
    d1 -= h + 2;
    d2 += h + 2;
  }
  return ret;
}

int main() {
  fact[0] = 1;
  ifact[0] = 1;
  inv[1] = 1;
  for (int i = 2; i < 505050; i++) {
    inv[i] = inv[mod % i] * (mod - mod / i);
  }
  for (int i = 1; i < 505050; i++) {
    fact[i] = i * fact[i - 1];
    ifact[i] = inv[i] * ifact[i - 1];
  }
  int n, m;
  cin >> n >> m;

  modint ans = 0;
  for (int i = 0; i <= n - 1; i++) {
    ans += f(m, n - 1, i);
  }
  for (int i = 0; i <= n - 2; i++) {
    ans -= f(m, n - 2, i);
  }
  cout << ans.n << endl;
}

そういえば通り数って聞き馴染みがない（twitter 上でしか見ない）んだけど、よく使われる言葉なのだろうか？場合の数という表現は見覚えがある。

2017-09-16

yukicoder No.569 3 x N グリッドのパスの数

https://yukicoder.me/problems/no/569

解法そのものは自明だが、真面目にやるべきものではない。行列による解法の存在に気づければ、解は線形漸化式に従うことが分かるので、その性質を使って楽をする。

解の列に対する母関数 a(x) を考える。ここである多項式 f(x) が存在して、f(x)a(x) の次数が L で抑えられる（つまり次数が高い場所では係数が 0 になる）。このような多項式を線形回帰列に対する annihilator と呼ぶらしい。annihilator を求めるアルゴリズムとして Berlekamp-Massey アルゴリズムが知られている [1]。具体値が既知で漸化式が未知のときに漸化式を求めるアルゴリズムとも考えられる。

annihilator を計算する際に、数列の具体値を計算する必要がある。これは graphillion などの ZDD 関連のライブラリを用いて計算できるが、自分が書いたコードは何故か遅かった（何かテクニックがあるのだろうか？）。以前に simpath を書いたことがあったので、そちらを使うことにする [2]。フロンティアが極めて小さいため、各フェーズにおける状態数は10～20個程度にしかならない。そのため、n=30 程度までなら一瞬で計算できる。試しに n=10^5 で計算させてみたところ 1.58 秒で結果が得られた。

[ 解答コード ]

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

using namespace std;

constexpr int mod = 1e9 + 7;

struct modint {
  int n;
  modint(int n = 0) : n(n) {}
};

modint operator+(modint a, modint b) { return modint((a.n += b.n) >= mod ? a.n - mod : a.n); }
modint operator-(modint a, modint b) { return modint((a.n -= b.n) < 0 ? a.n + mod : a.n); }
modint operator*(modint a, modint b) { return modint(1LL * a.n * b.n % mod); }
bool operator<(modint a, modint b) { return a.n < b.n; }
modint &operator+=(modint &a, modint b) { return a = a + b; }
modint &operator-=(modint &a, modint b) { return a = a - b; }
modint &operator*=(modint &a, modint b) { return a = a * b; }

modint modinv(modint n) {
  if (n.n == 1) return 1;
  return modinv(mod % n.n) * (mod - mod / n.n);
}
modint operator/(modint a, modint b) { return a * modinv(b); }

modint simpath(int n) {
  constexpr int H = 4;
  vector<pair<int, int>> es;
  for (int x = 0; x < n; x++) {
    for (int y = 0; y < H; y++) {
      if (y + 1 < H) es.emplace_back(x * H + y, x * H + y + 1);
      if (x + 1 < n) es.emplace_back(x * H + y, x * H + y + H);
    }
  }
  vector<int> deg(n * H);
  for (auto e : es) {
    deg[e.first]++;
    deg[e.second]++;
  }
  map<map<int, int>, modint> dp0, dp1;
  dp0[{{0, 0}}] = 1;
  for (auto e : es) {
    const int u = e.first;
    const int v = e.second;
    deg[u]--;
    deg[v]--;
    dp1.clear();
    for (auto kv : dp0) {
      map<int, int> mate1 = kv.first;
      if (!kv.first.count(v)) mate1[v] = v;
      const int mu = mate1[u];
      const int mv = mate1[v];
      const bool del = deg[u] == 0 && u != 0;
      if (del) mate1.erase(u);
      if (deg[u] != 0 || (u != 0 && mu == u) || (u == 0 && mu != u) || mu == -1) dp1[mate1] += kv.second;
      if (mu == v || mu == -1 || mv == -1) continue;
      mate1[u] = -1;
      mate1[v] = -1;
      mate1[mu] = mv;
      mate1[mv] = mu;
      if (mate1[0] == -1) continue;
      if (del) mate1.erase(u);
      if (!del || mu != u) dp1[mate1] += kv.second;
    }
    swap(dp0, dp1);
  }
  return dp0.rbegin()->second;
}

// ref: https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm 
vector<modint> berlekamp_massey(vector<modint> s) {
  const int N = s.size();
  vector<modint> C(N);
  vector<modint> B(N);
  C[0] = 1;
  B[0] = 1;
  int L = 0;
  int m = 1;
  modint b = 1;
  for (int n = 0; n < N; n++) {
    modint d = s[n];
    for (int i = 1; i <= L; i++) d += C[i] * s[n - i];
    if (d.n == 0) {
      m++;
    } else if (2 * L <= n) {
      auto T = C;
      for (int i = 0; i + m < N; i++) C[i + m] -= B[i] * (d / b);
      L = n + 1 - L;
      B = T;
      b = d;
      m = 1;
    } else {
      for (int i = 0; i + m < N; i++) C[i + m] -= B[i] * (d / b);
      m++;
    }
  }
  C.resize(L + 1);
  reverse(C.begin(), C.end());
  assert(L < N - 1);
  return C;
}

vector<modint> poly_mod(vector<modint> a, const vector<modint> &m) {
  const int n = m.size();
  for (int i = a.size() - 1; i >= m.size(); i--) {
    for (int j = 0; j < m.size(); j++) {
      a[i - n + j] += a[i] * m[j];
    }
  }
  a.resize(m.size());
  return a;
}

// a*b mod m(x)
vector<modint> poly_mul(const vector<modint> &a, const vector<modint> &b, const vector<modint> &m) {
  vector<modint> ret(a.size() + b.size() - 1);
  for (int i = 0; i < a.size(); i++) {
    for (int j = 0; j < b.size(); j++) {
      ret[i + j] += a[i] * b[j];
    }
  }
  return poly_mod(ret, m);
}

// x^n mod m(x)
vector<modint> nth_power(long long n, const vector<modint> &m) {
  vector<modint> ret(1);
  vector<modint> x(2);
  ret[0] = x[1] = 1;
  while (n > 0) {
    if (n & 1) ret = poly_mul(ret, x, m);
    x = poly_mul(x, x, m);
    n /= 2;
  }
  return poly_mod(ret, m);
}

int main() {
  vector<modint> a(30);
  for (int i = 0; i < a.size(); i++) {
    a[i] = simpath(i + 1);
  }

  vector<modint> m = berlekamp_massey(a);
  m.pop_back();
  for (int i = 0; i < m.size(); i++) {
    m[i] *= mod - 1;
  }

  long long n;
  cin >> n;
  auto x = nth_power(n, m);

  modint ans;
  for (int i = 0; i < x.size(); i++) {
    ans += x[i] * a[i];
  }
  cout << ans.n << endl;
}

問題に関係のない simpath 関連のコードは以下の gist にまとめた。速さよりもコンパクトさを求めて書いた。まだ短くなる気がするので、劇的な改善策を見つけたい。コードもまだ綺麗とは言えない状態。『超高速グラフ列挙アルゴリズム』をまだ読んでないんだけど、もしかしたらコード、あるいは実装方法が載ってたりする？流石に載ってない？

https://gist.github.com/pekempey/80fb0978306cc077afa2800dd214cddd

[1] https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm
[2] http://pekempey.hatenablog.com/entry/2017/01/26/203424

2017-09-03

ARC 082 E: ConvexScore

Keywords: convex hull DP
Time Complexity: \(O(N^3)\)

汎用的な処理が多数含まれているのでメモ。

#include <iostream>
#include <algorithm>
#include <vector>
#include <complex>
#include <cmath>
#include <cassert>
#include <tuple>

using namespace std;

constexpr long long mod = 998244353;

using P = complex<long long>;

long long dot(P a, P b) {
  return real(conj(a) * b);
}

long long cross(P a, P b) {
  return imag(conj(a) * b);
}

long long gcd(long long a, long long b) {
  a = abs(a);
  b = abs(b);
  if (b == 0) return a;
  return gcd(b, a % b);
}

// Get minimum element pointing in the same direction.
P tomato(P a) {
  long long g = gcd(a.real(), a.imag());
  return P(a.real() / g, a.imag() / g);
}

// Lexical order (argument, norm)
bool argcomp(P p1, P p2) {
  bool a = p1.imag() > 0 || (p1.imag() == 0 && p1.real() >= 0);
  bool b = p2.imag() > 0 || (p2.imag() == 0 && p2.real() >= 0);
  if (a != b) return a;
  long long c = p2.real() * p1.imag();
  long long d = p1.real() * p2.imag();
  if (c != d) return c < d;
  return norm(p1) < norm(p2);
}

namespace std {
  bool operator<(const P &a, const P &b) {
    return argcomp(a, b);
  }
}

int main() {
  int n;
  cin >> n;

  vector<P> ps(n);
  for (int i = 0; i < n; i++) {
    int x, y;
    cin >> x >> y;
    ps[i] = P(x, y);
  }

  vector<pair<int, int>> es;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      if (i == j) continue;
      es.emplace_back(i, j);
    }
  }
  // If two vectors are pointing same direction,
  // these are arranged such that both aren't taken.
  sort(es.begin(), es.end(), [&](pair<int, int> a, pair<int, int> b) {
    P d = ps[a.second] - ps[a.first];
    P e = ps[b.second] - ps[b.first];
    if (cross(d, e) == 0 && dot(d, e) >= 0) {
      return dot(ps[a.second] - ps[a.first], ps[b.first] - ps[a.first]) < 0;
    } else {
      return d < e;
    }
  });

  // Count the points inside (NOT on) the triangle(i,j,k)
  static int in[200][200][200];
  for (int j = 0; j < n; j++) {
    vector<pair<P, int>> qs;
    for (int i = 0; i < n; i++) {
      if (i == j) continue;
      qs.emplace_back(ps[i] - ps[j], i);
    }
    sort(qs.begin(), qs.end());
    for (int i = 0; i < n; i++) {
      if (i == j) continue;
      P s = tomato(ps[j] - ps[i]);
      int f = lower_bound(qs.begin(), qs.end(), make_pair(s, 0)) - qs.begin();
      for (int l = 0; l < n - 1; l++) {
        int k = qs[(f + l) % (n - 1)].second;
        in[i][j][k] -= l + 1;
        in[j][k][i] -= l + 1;
        in[k][i][j] -= l + 1;
      }
    }
  }
  static int on[200][200];
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      if (i == j) continue;
      for (int k = 0; k < n; k++) {
        if (i == k || j == k) continue;
        on[i][j] += cross(ps[k] - ps[i], ps[j] - ps[i]) == 0 && dot(ps[k] - ps[i], ps[j] - ps[i]) > 0 && dot(ps[k] - ps[j], ps[i] - ps[j]) > 0;
        in[i][j][k] = max(0, in[i][j][k] + n);
      }
    }
  }

  vector<long long> two(n + 1);
  two[0] = 1;
  for (int i = 1; i <= n; i++) {
    two[i] = two[i - 1] * 2 % mod;
  }

  // Easy
  long long ans = 0;
  for (int i = 0; i < n; i++) {
    vector<long long> dp1(n);
    vector<long long> dp2(n);
    for (auto e : es) {
      int u = e.first;
      int v = e.second;

      if (u == i) {
        assert(dp1[v] == 0);
        dp1[v] = two[on[u][v]];
      } else if (v == i) {
        dp2[v] += dp2[u];
      } else {
        int c = on[i][v] + on[u][v] + in[i][u][v];
        dp2[v] += dp1[u] * two[c];
        dp2[v] += dp2[u] * two[c];
      }
      dp2[v] %= mod;
    }
    ans += dp2[i];
    ans %= mod;
  }

  cout << ans << endl;
}

2017-08-25

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

http://codeforces.com/contest/843/problem/D

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 というらしいですね。

pekempey's blog

ACPC2017Day1 F Steps