++w++ ICPCライブラリ

++w++(amylase, y3eadgbe, atetubou(is2012))で使用する ICPC ライブラリ置き場。 (last updated: 7/6 12:21)

テンプレート

#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())

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

幾何

基本要素

typedef long double LD;
typedef complex<long double> pt;
typedef pair<pt, pt> line;
typedef pair<pt, long double> circle;
typedef vector<pt> poly;

• LD : long doubleをこれで書く。定数高速化が必要なときとかはtypedef double LDとかに書き換える。
• pt : 点もしくはベクトル。abs(x) で長さ。arg(x) で偏角[-pi, +pi]。a + bi = pt(a, b)
• line : 線分・直線。通る2点（線分では両端）。
• circle : 円。中心点と半径。
• poly : 多角形。反時計回りに頂点を持つ。

内積・外積の大きさ

long double inprod(pt a, pt b){
  return a.real() * b.real() + a.imag() * b.imag();
}

long double det(pt p, pt q){
  return p.real() * q.imag() - p.imag() * q.real();
}

• inprod : a, b の内積。射影を取ったりするときに使う。
• det : a, b が張る平行四辺形の面積。三角形の面積とか。

3点の位置関係

long double ccw(pt p, pt q, pt r){
  pt n = (q-p) * pt(0, 1);
  long double c = inprod(n, p);
  return inprod(n, r) - c;
}

• ccw > 0 → 反時計回り。

• ccw < 0 → 時計回り。

• ccw = 0 → 一直線上にある。
• 直線と点の位置関係にも使えますよ。

点と直線の距離

long double dist(line l, pt c){
  pt a = l.first, b = l.second;
  return abs(det(b-a, c-a)) / abs(b-a);
}

• 直線 l と点 c の距離。

線分の交差判定

bool is_cross(line p, 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 が交差するかどうか。

線分と点の距離

ld stick_point_dist(line l, pt c){
  pt n = (l.second - l.first) * pt(0, 1);
  if(ccw(l.first, l.first + n, c) * ccw(l.second, l.second + n, c) <= 0){
    return lpdist(l, c);
  } else {
    return min(abs(l.first - c), abs(l.second - c));
  }
}

直線の交点

pt intersect(line l, line m){
  pt x = l.second - l.first;
  pt y = m.second - m.first;
  pt p = l.first;
  pt c = m.first - l.first;

  double d = det(x, y);
  double k = det(c, y) / d;

  return p + k * x;
}

多角形の面積

long double area(poly p){
  long double ret = 0.0;
  pt bp = p.back();
  rep(i, p.size()){
    ret += p[i].imag() * bp.real() - p[i].real() * bp.imag();
    bp = p[i];
  }
  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> の各要素が接線。

数論

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

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;
}

あとでけす

蟻本の最小費用流

const li MAX_V = 1200;

struct edge { li to, cap, cost, rev; };

li V;
vector<edge> G[MAX_V];
li dist[MAX_V];
li prevv[MAX_V], preve[MAX_V];

void add_edge(li from, li to, li cap, li cost){
  G[from].push_back((edge){to, cap, cost, G[to].size()});
  G[to].push_back((edge){from, 0, -cost, G[from].size() - 1});
}

li min_cost_flow(li s, li t, li f) {
  li res = 0;

  while (f > 0) {
    fill(dist, dist + V, inf);
    dist[s] = 0;
    bool update = true;
    while (update) {
      update = false;
      for (li v = 0; v < V; v++) {
        if (dist[v] == inf) continue;
        for (li i = 0; i < G[v].size(); i++) {
          edge &e = G[v][i];
          if (e.cap > 0 && dist[e.to] > dist[v] + e.cost) {
            dist[e.to] = dist[v] + e.cost;
            prevv[e.to] = v;
            preve[e.to] = i;
            update = true;
          }
        }
      }
    }

    if (dist[t] == inf) {
      return -1;
    }

    li d = f;
    for (li v = t; v != s; v = prevv[v]) {
      d = min(d, G[prevv[v]][preve[v]].cap);
    }
    f -= d;
    res += d * dist[t];
    for (li v = t; v != s; v = prevv[v]) {
      edge &e = G[prevv[v]][preve[v]];
      e.cap -= d;
      G[v][e.rev].cap += d;
    }
  }

  return res;
}

蟻本の最大二部マッチング

vector<vi> g;
vector<int> used, match;
int V;

void add_edge(int u, int v){
  g[v].push_back(u);
  g[u].push_back(v);
  return;
}

bool dfs(int v){
  used[v] = 1;
  for(int i=0; i<g[v].size(); i++){
    int u = g[v][i], w = match[u];
    if(w < 0 || !used[w] && dfs(w)){
      match[v] = u;
      match[u] = v;
      return true;
    }
  }
  return false;
}

int bipartite_matching(){
  int res = 0;
  match = vector<int>(V, -1);
  for(int v=0; v<V; v++){
    if(match[v] < 0){
      used = vector<int>(V, 0);
      if(dfs(v)){
        res++;
      }
    }
  }
  return res;
}

蟻本の最大流

struct edge{int to,cap,rev;};
const int MAX_V=600;
vector<edge> G[MAX_V];
int level[MAX_V];
int iter[MAX_V];

void add_edge(int from,int to,int cap){
  G[from].pb(((edge){to,cap,G[to].size()}));
  G[to].pb(((edge){from,0,G[from].size()-1}));
}

void bfs(int s){
  memset(level,-1,sizeof(level));
  queue<int> que;
  level[s]=0;
  que.push(s);
  while(!que.empty()){
    int v=que.front();que.pop();
    for(int i=0;i<SZ(G[v]);++i){
      edge &e=G[v][i];
      if(e.cap>0 && level[e.to]<0){
        level[e.to]=level[v]+1;
        que.push(e.to);
      }
    }
  }
}

int dfs(int v,int t,int f){
  if(v==t) return f;
  for(int &i=iter[v];i<SZ(G[v]);++i){
    edge &e=G[v][i];
    if(e.cap>0 && level[v]<level[e.to]){
      int d=dfs(e.to,t,min(f,e.cap));
      if(d>0){
        e.cap -=d;
        G[e.to][e.rev].cap += d;
        return d;
      }
    }
  }
  return 0;
}

int max_flow(int s,int t){
  int flow = 0;
  for(;;){
    bfs(s);
    if(level[t]<0) return flow;
    memset(iter,0,sizeof(iter));
    int f;
    while((f=dfs(s,t,1<<30))>0)
      flow += f;
  }
}