++w++ ICPCライブラリ

++w++(amylase, y3eadgbe, atetubou(is2012))で使用する ICPC ライブラリ置き場。 (last updated: 11/8 20:42)

TODO
- 線型計画を二段階単体法にする。
- 幾何(線分アレンジメント)
- 文字列(KMP method, Aho-Corasick algorithm)
- (ここには載せないが、蟻本から最大流・最小費用流・拡張ユークリッド・連立線型合同式を持っていく)
- 余裕があれば有名問題を増やす。
  - グラフ彩色

++w++ ICPCライブラリ

テンプレート

#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
#include <queue>
#include <cstring>
#include <complex>

using namespace std;

#define rep(i,n) for(int i=0;i<n;++i)
#define pb(a) push_back(a)
#define mp(a,b) make_pair(a,b)
#define F first
#define S second
#define SZ(a) (int)((a).size())

だいたいコードの先頭にあると考えておいてね。

幾何

基本要素

// primitives
typedef long double Real;
typedef complex<Real> Pt;
typedef pair<Pt, Pt> Line;
typedef pair<Pt, Real> Circle;
typedef vector<Pt> Poly;

Real eps = 1e-9;

inline istream& operator>>(istream& s, Pt& p) {return s >> p.real() >> p.imag();}
inline bool near(const Pt& p, const Pt& q){return abs(p - q) < eps;}

Real : long double, double の切り替え用に。
Pt : 点もしくはベクトル。abs(x) で長さ。arg(x) で偏角[-pi, +pi]。a + bi = pt(a, b)
Line : 線分・直線。通る2点（線分では両端）。
Circle : 円。中心点と半径。
Poly : 多角形。反時計回りに頂点を持つ。
eps : 許容誤差。問題によって値を修正する。
operator>> : (x座標) (y座標) と来た時に一発で読み取る。
near : 誤差を許容して同じ点かどうか見るよ。

内積・外積の大きさ

inline Real dot(const Pt& p, const Pt& q){return p.real() * q.real() + p.imag() * q.imag();}
inline Real cross(const Pt& p, const Pt& q){return p.real() * q.imag() - p.imag() * q.real();}

dot : a, b の内積。射影を取ったりするときに使う。
cross : a, b が張る平行四辺形の面積。三角形の面積とか。a, b を通る直線上の点 x が cross(x, b) = cross(a, b) を満たすことから、abの法線を軸としたパラメータとの解釈もある。

3点の位置関係

// ccw_b : old version of ccw.
inline Real ccw_b(const Pt& p, const Pt& q, const Pt& r) {return cross(q - p, r - p);}

/* ccw :
CD  : counter direction
CW  : clock wise
OS  : on segment
CCW : counter clock wise
D   : direction
 */
enum LPposit { P_CD = -2, P_CW = -1, P_OS = 0, P_CCW = 1, P_D = 2};
LPposit ccw(const Pt& p, const Pt& q, const Pt& r) {
  Real c = ccw_b(p, q, r);
  if (c < -eps) return P_CW;
  if (c >  eps) return P_CCW;
  if (dot(q - p, r - p) < -eps) return P_CD;
  if (dot(p - q, r - q) < -eps) return P_D;
  return P_OS;
}

ccw_b : 上で述べた cross のパラメータ解釈での値を計算。
ccw : コメント参照!

点と直線の距離

inline Real Sabs(const Line& l) {return abs(l.first - l.second); }
inline Real LPdist(const Line& l, const Pt& p) {return abs(ccw_b(l.first, l.second, p)) / Sabs(l); }
inline Real SPdist(const Line& l, const Pt& p) {
  if ((Labs(l) > abs(l.first - p) && Labs(l) > abs(l.second - p)) ||
      cos(arg(p - l.first) - arg(p - l.second)) > 0){

    return LPdist(l, p);
  }

  return min(abs(p - l.first), abs(p - l.second));
}

Sabs : 線分の長さ。
LPdist : 点と直線の距離。
SPdist : 点と線分の距離。
- L : line, S : segment, P : point

線分の交差判定

bool crossS(const Line& p, const Line& q){
  return
    ccw(p.first, p.second, q.first) * ccw(p.first, p.second, q.second) <= 0 &&
    ccw(q.first, q.second, p.first) * ccw(q.first, q.second, p.second) <= 0;
}

線分 p と q が交差するかどうか。

直線の交点

Pt intersect(const Line& p, const Line& q) {
  Pt vp = p.second - p.first;
  Pt vq = q.second - q.first;
  Pt c(cross(vp, p.first), cross(vq, q.first));
  return Pt(cross(c, Pt(vp.real(), vq.real())), cross(c, Pt(vp.imag(), vq.imag()))) / cross(vp, vq);
}

多角形の面積

Real area(const Poly& p) {
  Real ret = 0.0;
  for(Poly::iterator it = p.begin(); it != p.end(); it++) ret += cross(*it, *(it + 1));
  return ret / 2;
}

多角形の面積。凸性は仮定しない。O(n)。
時計回りの時は面積が -1 倍される。abs 取ろう。

円と直線の交点

vector<pt> circle_line_intersect(line l,circle c){
  vector<pt> ret;
  LD di = dist(l,c.first);
  LD r=c.second;
  if(di+EPS > r) return ret;
  pt v=(l.second-l.first);
  v/=abs(v);  
  pt rv=v*pt(0,1);
  rv*=di;  
  if(dist(l,c.first+rv) > di+EPS) rv = -rv;
  v*=sqrt(r*r-di*di);
  ret.push_back(c.first+rv-v);
  ret.push_back(c.first+rv+v);
  return ret;
}

円と直線の交点を求める。
円が直線と接しない場合はemptyを返し、それ以外は取り敢えず二つの要素を持つvectorを返す。

円の共通接線

vector<line> contact(circle p, circle q){
  vector<line> ret;
  if(p.second < q.second) swap(p, q);
  long double d = abs(p.first - q.first);
  pt n = q.first - p.first;
  n /= abs(n);

  if(d + eps < abs(p.second - q.second)){
    ret.clear();
  } else if(eq(d, abs(p.second - q.second))){
    pt t, u;
    t = p.first + p.second * n;
    u = t + n * pt(0, 1);
    ret.push_back(make_pair(t, u));
  } else {
    if(!eq(p.second, q.second)){
      pt t = p.first + (p.second * d / (p.second - q.second)) * n;
      long double theta = asin((p.second - q.second) / d);
      pt u = n * pt(cos(theta), sin(theta));
      pt v = n * pt(cos(-theta), sin(-theta));
      u += t;
      v += t;
      ret.push_back(make_pair(t, u));
      ret.push_back(make_pair(t, v));
    } else {
      pt t = p.first + n * pt(0, 1) * p.second;
      pt u = p.first - n * pt(0, 1) * p.second;
      ret.push_back(make_pair(t, t+n));
      ret.push_back(make_pair(u, u+n));
    }

    if(eq(d, p.second + q.second)){
      pt t, u;
      t = p.first + p.second * n;
      u = t + n * pt(0, 1);
      ret.push_back(make_pair(t, u));
    } else if(d > p.second + q.second){
      pt t = p.first + (p.second * d / (p.second + q.second)) * n;
      long double theta = asin((p.second + q.second) / d);
      pt u = n * pt(cos(theta), sin(theta));
      pt v = n * pt(cos(-theta), sin(-theta));
      u += t;
      v += t;
      ret.push_back(make_pair(t, u));
      ret.push_back(make_pair(t, v));
    }
  }
  return ret;
}

vector<line> の各要素が接線。

凸包

// convex-hull                                                                  
// return minimum convex polygon that contains every point in vs.               
// depend on ccw                                                                
namespace std{
  bool operator<(const Pt& a, const Pt& b) {
    return a.real() != b.real() ? a.real() < b.real() : a.imag() < b.imag();
  }
}

Poly convexHull(vector<Pt> ps) {
  int n = ps.size();
  sort(ps.begin(), ps.end());

  Poly ret(2 * n);
  int m = 0;
  for (int i = 0; i < n; i++) {
    while (m >= 2 && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
    ret[m++] = ps[i];
  }

  int t = m;
  for (int i = n-2; i >= 0; i--) {
    while (m >= t && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
    ret[m++] = ps[i];
  }

  ret.resize(m - 1);
  return ret;
}

Pt の座標軸ソートもついてる。
convexHull : ps に含まれる点をすべて含む最小の凸多角形を返す。spaghetthi source を参考にしたよ。

数論

ミラー・ラビン素数判定法

long long powmod(long long x, long long p, long long m){
  if(p == 0) return 1;
  long long rt = powmod(x, p/2, m);
  if(p % 2 == 0){
    return rt * rt % m;
  } else {
    return (rt * rt % m) * x % m;
  }
}

int miller_test[9] = {2,3,5,7,11,13,17,19,23};
bool isprime(long long n){
  if(n <= ISPR_MAX) return !isnpr[n];

  long long d = n-1, s = 0;
  while(d%2 == 0){
    d /= 2;
    s++;
  }

  for(int i=0; i<9; i++){
    bool iscomp = true;
    li x = powmod(miller_test[i], d, n);
    iscomp = iscomp && (x % n != 1);

    for(int r=0; r<s; r++){
      iscomp = iscomp && (x % n != n-1);
      x = x * x % n;
    }
    if(iscomp) return false;
  }

  return true;
}

ポラード・ロー素因数分解法

long long myrand(long long c, long long n, long long x){
  return (x * x + c) % n;
}
 
long long pollard(long long n){
  long long x = 2, y = 2, d = 1, c = rand();
  while(d == 1){
    x = myrand(c, n, x);
    y = myrand(c, n, myrand(c, n, y));
    d = __gcd(abs(x-y), n);
  }
  return d;
}

数学

行列乗算

typedef vector<vector<li> > matrix;
const li mod = 1000000007;

matrix ident(const int& n) {
    matrix ret(n, vi(n, 0));
    rep(i, n) ret[i][i] = 1;
    return ret;
}

matrix matadd(const matrix& p, const matrix& q) {
    int n = p.size(), m = p[0].size();
    matrix ret = p;
    rep(i, n) rep(j, m) ret[i][j] = (p[i][j] + q[i][j]) % mod;
    return ret;
}

matrix matmul(const matrix& p, const matrix& q) {
    int n = p.size(), m = q[0].size(), l = p[0].size();
    matrix ret = matrix(n, vi(m, 0));
    rep(i, n) rep(j, m) rep(k, l) ret[i][j] = (ret[i][j] + p[i][k] * q[k][j] % mod) % mod;
    return ret;
}

matrix matpow(const matrix& p, const li& x) {
    if (x == 0) return ident(p.size());
    if (x == 1) return p;
    matrix ret = matpow(p, x / 2);
    ret = matmul(ret, ret);
    if (x % 2) ret = matmul(ret, p);
    return ret;
}

2次元配列で行列を表現。エラーチェックなし。mod とってます。

畳み込み

typedef double Real;
typedef complex<Real> Complex;

vector<Complex> fastFourierTransform(const vector<Complex>& x) {
    if (__builtin_popcount(x.size()) > 1) {
        cerr << "fastFourierTransform: x.size() should be power of 2." << endl;
        return x;
    }

    const int n = x.size();
    if (n == 1) {
        return x;
    }

    vector<Complex> even(n / 2), odd(n / 2), ret(n);
    for (int i = 0; i < n; ++i) {
        if (i % 2 == 0) {
            even[i / 2] = x[i];
        } else {
            odd[i / 2] = x[i];
        }
    }

    even = fastFourierTransform(even);
    odd = fastFourierTransform(odd);

    for (int i = 0; i < n; ++i) {
        int j = i % (n / 2);
        ret[i] = even[j] + odd[j] * exp(Complex(0, -2 * M_PI * i / n));
    }

    return ret;
}

vector<Complex> inverseFastFourierTransform(vector<Complex> x) {
    const int n = x.size();
    for (int i = 0; i < n; ++i) {
        x[i] = conj(x[i]);
    }
    x = fastFourierTransform(x);
    for (int i = 0; i < n; ++i) {
        x[i] = conj(x[i]) / Complex(n);
    }
    return x;
}

int conv_minpow2(int x) {
    int ret = 1;
    while (ret < x) ret <<= 1;
    return ret;
}

vector<Complex> convolution(vector<Complex> x,
                            vector<Complex> y) {
    int l = 2 * max(conv_minpow2(x.size()), conv_minpow2(y.size()));
    x.resize(l);
    y.resize(l);

    x = fastFourierTransform(x);
    y = fastFourierTransform(y);

    vector<Complex> ret(l);
    for (int i = 0; i < l; i++) {
        ret[i] = x[i] * y[i];
    }

    ret = inverseFastFourierTransform(ret);
    while (not ret.empty() && ret.back() == Complex(0)) ret.pop_back();
    return ret;
}

高速フーリエ変換による畳み込み演算。O(n log n)。
多項式の掛け算がナイーブにやるとO(n^2)だが畳込みを使うと計算量が落ちたり。

有名問題

線型計画法

// linear programming by simplex method
// solve maximize(cx) s.t. Ax <= a
// if not bounded, return -1.
double simplex(vector<double>& c, 
               vector<vector<double> >& A,
               vector<double>& a) {
  // introduce slack variable
  int n = c.size(), m = a.size(), dim = n + m;
  vector<int> N(n), B(m);
  for (int i = 0; i < n; ++i) N[i] = i;
  for (int i = 0; i < m; ++i) {
    B[i] = i + n;
    c.push_back(0);
    for (int j = 0; j < n; j++) A[i][j] *= -1;
    for (int j = 0; j < m; j++) A[i].push_back(0);
  }

  double ans = 0;

  while (true) {
    // check optimized or not
    int s = -1;
    for (int i = 0; i < n; i++) if (c[N[i]] > 0) s = N[i];
    if (s < 0) break;

    // check bounded or not
    double bound = 1e300;
    int r = -1;
    for (int i = 0; i < m; i++) {
      if (A[i][N[s]] < 0) {
        double nbound = -a[i] / A[i][N[s]];
        if (nbound < bound) {
          bound = nbound;
          r = i;
        }
      }
    }

    if (r < 0) return -1;

    // pivotting
    for (int i = 0; i < dim; i++) if (i != N[s]) A[r][i] /= -A[r][N[s]];
    a[r] /= -A[r][N[s]];
    A[r][B[r]] = 1.0 / A[r][N[s]];
    A[r][N[s]] = 0;

    for (int i = 0; i < m; i++) {
      if (i == r) continue;

      for (int j = 0; j < dim; j++) if (j != N[s]) A[i][j] += A[i][N[s]] * A[r][j];
      a[i] += A[i][N[s]] * a[r];
      A[i][N[s]] = 0;
    }

    ans += c[N[s]] * a[r];
    for (int i = 0; i < dim; i++) if (i != N[s]) c[i] += c[N[s]] * A[r][i];
    c[N[s]] = 0;

    swap(N[s], B[r]);
  }

  return ans;
}

simplex(c, A, a) : maximize cx, where Ax <= a, x >= 0 なる線型計画問題を単体法で解く。
x = 0 が実行可能解であることを仮定している（つまり、TODOとして二段階単体法に改善することがあげられる）。
発散する場合には -1 を返す。有限な最適解は必ず正になっているはずなので、正負を見れば発散したかどうかがわかる。

amylase/icpc

++w++ ICPCライブラリ

テンプレート

幾何

基本要素

内積・外積の大きさ

3点の位置関係

点と直線の距離

線分の交差判定

直線の交点

多角形の面積

円と直線の交点

円の共通接線

凸包

数論

ミラー・ラビン素数判定法

ポラード・ロー素因数分解法

数学

行列乗算

畳み込み

有名問題

線型計画法