icpcgeolib - arosh/arosh.github.com GitHub Wiki

#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
using namespace std;

typedef double D;
const D EPS = 1e-8;
const D PI = M_PI;
struct P {
  D x, y;
  P() { }
  P(D x_, D y_) : x(x_), y(y_) { }
};
struct L {
  P a, b;
  L() { }
  L(P a_, P b_) : a(a_), b(b_) { }
};
struct C {
  P p; D r;
  C() { }
  C(P p_, D r_) : p(p_), r(r_) { }
};
P operator +(P a, P b) { return P(a.x + b.x, a.y + b.y); }
P operator -(P a, P b) { return P(a.x - b.x, a.y - b.y); }
P operator *(P a, D b) { return P(a.x * b, a.y * b); }
P operator /(P a, D b) { return P(a.x / b, a.y / b); }
P vec(P st, P en) { return en - st; }
D inp(P a, P b) { return a.x*b.x + a.y*b.y; }
D outp(P a, P b) { return a.x*b.y - a.y*b.x; }
D abs(P p) { return sqrt(inp(p,p)); }
int sgn(D a, D b) {
  if(a < b - EPS) return -1;
  if(a > b + EPS) return +1;
  return 0;
}
// !!! とりあえずここまで写してください !!!
D arg(P p) { return atan2(p.y,p.x); }
P rot90(P p) { return P(-p.y,p.x); }
P rot(P p, D rad) {
  return P(cos(rad)*p.x - sin(rad)*p.y,
           sin(rad)*p.x + cos(rad)*p.y);
}
int ccw(P a, P b, P c) { // 重なっている点があるとうまく動かないと思われる
  b = vec(a,b); c = vec(a,c);

  // a - b - c が折れ曲がるとき
  if(sgn(outp(b,c),0.0) > 0) return +1; // 反時計回り
  if(sgn(outp(b,c),0.0) < 0) return -1; // 時計回り

  // a - b - c が直線上に並ぶとき
  // sgnを使わない実装が主流っぽい？
  if(sgn(inp(b,c),0.0) < 0)  return +2; // c - a - b
  if(sgn(abs(b),abs(c)) < 0) return -2; // a - b - c
  return 0;                             // a - c - b
}
P proj(L l, P p) { // 直線lに対する点pの写像
  P a = vec(l.a,l.b);
  P b = vec(l.a,p);
  D t = inp(a,b)/abs(a);
  return l.a + (a/abs(a))*t;
}
P refl(L l, P p) { // 直線lに対する点pの反射
  return p + vec(p, proj(l, p)) * 2;
}
// 線分と点
bool iSP(L s, P p) {
  return ccw(s.a, s.b, p) == 0;
}
D dSP(L s, P p) {
  P r = proj(s, p);
  if(iSP(s, r)) return abs(p - r); // 写像がs上にある
  return min(abs(p - s.a), abs(p - s.b)); // 写像がs上にない
}
// 直線と点
bool iLP(L l, P p) {
  return abs(ccw(l.a, l.b, p)) != 1;
}
D dLP(L l, P p) {
  return abs(p - proj(l, p));
}
// 線分と線分
bool iSS(L s, L t) { // 端を含む
  return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
    ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}
D dSS(L s, L t) {
  if(iSS(s, t)) return 0;
  return min(min(dSP(s, t.a), dSP(s, t.b)),
      min(dSP(t, s.a), dSP(t, s.b)));
}
// P cSS(L s, L t) はiSS(s,t)を確認してからcLLを使う
// 直線と直線
bool iLL(L l, L m) {
  return sgn(outp(vec(l.a, l.b), vec(m.a, m.b)), 0.0) != 0 || // 平行でない
    sgn(outp(vec(l.a, l.b), vec(l.a, m.a)), 0.0) == 0;   // 同じ直線
}
D dLL(L l, L m) {
  return iLL(l, m) ? 0.0 : dLP(l, m.a);
}
P cLL(L l, L m) {
  D d = outp(vec(m.a, m.b), vec(l.a, l.b));
  return l.a + vec(l.a, l.b) * outp(vec(m.a, m.b), vec(l.a, m.b)) / d;
}
// 直線と線分
bool iLS(L l, L s) {
  return sgn(outp(vec(l.a, l.b), vec(l.a, s.a)), 0.0) *
    sgn(outp(vec(l.a, l.b), vec(l.a, s.b)), 0.0) <= 0;
}
D dLS(L l, L s) {
  if(iSS(l, s)) return 0.0;
  return min(dLP(l, s.a), dLP(l, s.b));
}
// P cLS(L s, L t) はiSS(s,t)を確認してからcLLを使う
// 円と直線
int iCL(C c, L l) {
  D d = dLP(l, c.p);
  int s = sgn(d, c.r);
  if(s < 0) return 2; // 2点で交差する
  if(s == 0) return 1; // 1点で接する
  return 0; // 触れない
}
vector<P> cCL(C c, L l) {
  P a = proj(l, c.p);
  D s = abs(c.p - a);
  D t = sqrt(c.r*c.r-s*s);
  P v = vec(l.a,l.b)/abs(l.a-l.b);
  vector<P> res;
  res.push_back(a+v*t);
  res.push_back(a-v*t);
  return res;
}
// 円と点
vector<P> tCP(C c, P p) {
  vector<P> res;
  D d = abs(c.p - p);
  if(sgn(d,c.r) < 0) return res;
  D rc = c.r*c.r/d;
  D rs = sqrt(max(0.0, c.r*c.r - rc*rc));
  P v = (p - c.p) / abs(p - c.p);
  res.emplace_back(c.p + v*rc + rot90(v)*rs);
  res.emplace_back(c.p + v*rc - rot90(v)*rs);
  return res;
}
// 円と線分
int iCS(C c, L s) {
  if(sgn(abs(c.p - s.a), c.r) <= 0 && sgn(abs(c.p - s.b), c.r) <= 0)
    return -2; // 線分の両端が円の内側(境界含む)にある
  if(sgn(abs(c.p - s.a), c.r) <= 0 || sgn(abs(c.p - s.b), c.r) <= 0)
    return -1; // 線分の端の一方が円の内側(境界含む)、他方が円の外側にある
  if(sgn(dSP(s, c.p), c.r) < 0)
    return +2; // 線分と円周が2点で交わる
  if(sgn(dSP(s, c.p), c.r) == 0)
    return +1; // 線分と円周が1点で交わる
  return 0;
}
// vector<P> cCS(C c, L s) はiCSの状態によって定義が難しいが、
// p <- cCL(c,s)からiSP(s,p)==trueのものだけを抜き出せば良い気がする
//
// 円と円
// int iCC(C a, C b) はaがbの内側にある場合などもあるので定義が難しい
// 2点で接していることを確認すること
vector<P> cCC(C a, C b) {
  D d = abs(b.p - a.p);
  D x = (d*d + a.r*a.r - b.r*b.r) / (2*d);
  D y = sqrt(a.r*a.r - x*x);
  P v = (b.p - a.p) / d;
  vector<P> res;
  res.emplace_back(a.p + v*x + rot90(v)*y);
  res.emplace_back(a.p + v*x - rot90(v)*y);
  return res;
}
// tCCinに同じ円を食わせたり、
// tCCoutにどちらかが他方に囲まれている円を食わせたりすると破綻することは分かっているが、
// そのあたりが厳密にverifyできていない (AOJ 2201でACすることは確認した)
vector<L> tCCout(C a, C b) {
  vector<L> res;
  if(sgn(abs(a.r - b.r), abs(a.p - b.p)) >= 0) return res;
  if(sgn(a.r, b.r) == 0) {
    P v = (b.p - a.p) / abs(b.p - a.p);
    v = rot90(v);
    res.emplace_back(a.p + v*a.r, b.p + v*b.r);
    res.emplace_back(a.p - v*a.r, b.p - v*b.r);
    return res;
  }
  P p = (a.p*b.r - b.p*a.r) / (b.r - a.r);
  vector<P> at = tCP(a, p);
  vector<P> bt = tCP(b, p);
  for(int i = 0; i < (int)min(at.size(), bt.size()); ++i) {
    res.emplace_back(at[i], bt[i]);
  }
  return res;
}
vector<L> tCCin(C a, C b) {
  vector<L> res;
  if(sgn(abs(a.r + b.r), abs(a.p - b.p)) >= 0) return res;
  P p = (a.p*b.r + b.p*a.r) / (a.r + b.r);
  vector<P> at = tCP(a, p);
  vector<P> bt = tCP(b, p);
  for(int i = 0; i < (int)min(at.size(), bt.size()); ++i) {
    res.emplace_back(at[i], bt[i]);
  }
  return res;
}
vector<L> tCC(C a, C b) {
  vector<L> res;
  for(L l : tCCout(a, b)) res.push_back(l);
  for(L l : tCCin(a, b)) res.push_back(l);
  return res;
}
// 半時計回りを仮定している (時計回りならabsを取る)
D area(vector<P> g) {
  D res = 0.0;
  for(int i = 0; i < (int)g.size(); i++) {
    res += outp(g[i], g[(i+1)%g.size()]);
  }
  return res / 2.0;
}
// ON = 0, IN = 1, OUT = -1
int containsGP(vector<P> g, P p) {
  int side = -1;
  for(int i = 0; i < (int)g.size(); i++) {
    if(ccw(g[i], g[(i+1)%g.size()], p) == 0) return 0;
    P a = vec(p, g[i]);
    P b = vec(p, g[(i+1)%g.size()]);
    if(a.y > b.y) swap(a, b);
    if(sgn(a.y, 0.0) <= 0 && sgn(b.y, 0.0) > 0 && sgn(outp(a, b), 0.0) > 0) side *= -1;
  }
  return side;
}
bool operator <(P a, P b) {
  if(sgn(a.x, b.x) != 0) return a.x < b.x;
  return a.y < b.y;
}
// 凸包を構成する点を得る。半時計回り
// 辺上の点も含めるときは ccw(..) == -1 とすること
vector<P> convex_hull(vector<P> ps) {
  int N = ps.size();
  int k = 0; // 凸包を構成する点の数
  sort(ps.begin(), ps.end());
  vector<P> res(N * 2);
  for(int i = 0; i < N; i++) {
    // 時計回りの折れ曲がりがあったら、折れ点を削除
    while(k >= 2 && ccw(res[k - 2], res[k - 1], ps[i]) <= 0) k--;
    res[k++] = ps[i];
  }
  int t = k + 1;
  // 右端は取らない (重複するから)
  // 左端は取る (重複するけど、あとから取り除いたほうが楽)
  for(int i = N - 2; i >= 0; i--) {
    while(k >= t && ccw(res[k - 2], res[k - 1], ps[i]) <= 0) k--;
    res[k++] = ps[i];
  }
  res.resize(k - 1);
  return res;
}
// l.a -> l.b の厳密に右側の領域を切り落とす
vector<P> convex_cut(vector<P> g, L l) {
  vector<P> res;
  for(int i = 0; i < (int)g.size(); ++i) {
    P a = g[i];
    P b = g[(i+1)%g.size()];
    if(ccw(l.a, l.b, a) != -1) res.push_back(a);
    // 端の点を含まないiLS
    if(ccw(l.a, l.b, a)*ccw(l.a, l.b, b) < 0) {
      res.push_back(cLL(L(a, b), l));
    }
  }
  return res;
}
// aとbの垂直二等分線をaが左側に来るように計算する
L bisector(P a, P b) {
  P p = (a + b) / 2;
  return L(p, p + rot90(vec(a,b)));
}
// 外枠をgとして点集合vのk番目の点のボロノイ領域を返す
vector<P> voronoi_cell(vector<P> g, vector<P> v, int k) {
  for(int i = 0; i < (int)v.size(); ++i) {
    if(i == k) continue;
    g = convex_cut(g, bisector(v[k], v[i]));
  }
  return g;
}
struct Edge {
  int from, to;
  double cost;
  Edge(int from_, int to_, double cost_) : from(from_), to(to_), cost(cost_) { }
};
// uniqueするのに必要 (悪用厳禁)
// もしかしたらEPSを考慮したほうが良いかもしれない
bool operator ==(P a, P b) {
  return a.x == b.x && a.y == b.y;
}
vector<vector<Edge>> segment_arrangement(vector<L> ss, vector<P> &ps) {
  for(int i = 0; i < (int)ss.size(); ++i) {
    ps.push_back(ss[i].a);
    ps.push_back(ss[i].b);
    for(int j = i + 1; j < (int)ss.size(); ++j) {
      // iSSは端を含むが、後でuniqueするので気にしない
      if(iSS(ss[i], ss[j])) {
        ps.push_back(cLL(ss[i], ss[j]));
      }
    }
  }
  sort(ps.begin(), ps.end());
  ps.erase(unique(ps.begin(), ps.end()), ps.end());

  vector<vector<Edge>> G(ps.size());
  for(int i = 0; i < (int)ss.size(); ++i) {
    vector<pair<double,int>> li;
    // ps中にssの端点も含まれていることに注意
    for(int j = 0; j  < (int)ps.size(); ++j) {
      if(iSP(ss[i], ps[j])) {
        li.emplace_back(abs(vec(ss[i].a,ps[j])), j);
      }
    }
    sort(li.begin(), li.end());
    for(int j = 0; j < (int)li.size() - 1; ++j) {
      int a = li[j].second;
      int b = li[j+1].second;
      D len = abs(vec(ps[a], ps[b]));
      G[a].emplace_back(a, b, len);
      G[b].emplace_back(b, a, len);
    }
  }
  return G;
}

cerr << "line(" << a.x << "," << a.y << "," << b.x << "," << b.y << ");" << endl;
cerr << "circle(" << c.p.x << "," << c.p.y << "," << c.r << ");" << endl;