/*
      [幾何ライブラリ]
      ・円周率
      ・誤差チェッカー
      ・点(ベクトル)を表わす構造体
      ・線分(直線)を表わす構造体
      ・ベクトルのノルムの2乗
      ・ベクトルの大きさ
      ・ベクトルの内積
      ・ベクトルの外積
      ・ベクトルの直交判定
      ・ベクトルの平行判定
      ・直線sに対して点pから垂線を引いたときの交点x(直線sに対する点pの射影)
      ・直線sに対して点pと線対称の位置にある点p'
      ・点と点の距離
      ・直線と点の距離
      ・線分と点の距離
      ・線分p0p1と点p2の位置関係を求める
      ・2本の線分の交差判定
      ・線分と線分の距離
      ・交差する2本の直線の交点
      ・交差する2本の線分の交点
      ・座標cを中心とする半径rの円を表わす構造体
      ・円と直線の交点の組(接点なら同じものを2つ)を求める(小さい順に入れる)
      ・ベクトルの偏角
      ・極座標から直交座標に変換
      ・円と円の交点の組(接点なら同じものを2つ)を求める(小さい順に入れる)
      ・点の内包判定
      ・凸包
      */
//円周率
static const double pi = 3.141592653589793;
//誤差チェッカー
#define EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS)
//点を表わす構造体
struct Point
{
public:
    double x, y;
    Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
    /*
          以下ベクトルに対する演算子の定義
          +:ベクトルの足し算
          -:ベクトルの引き算
          *:スカラー倍
          /:スカラーの逆数倍
          */
    Point operator+(const Point p) { return Point(x + p.x, y + p.y); }
    Point operator-(const Point p) { return Point(x - p.x, y - p.y); }
    Point operator*(const double k) { return Point(x * k, y * k); }
    Point operator/(const double k) { return Point(x / k, y / k); }
    double abs() { return sqrt(norm()); }
    double norm() { return x * x + y * y; }
    //ベクトルの比較(x軸方向のベクトルから比較)
    bool operator<(const Point &p) const
    {
        return x != p.x ? x < p.x : y < p.y;
    }
    bool operator>(const Point &p) const
    {
        return x != p.x ? x > p.x : y > p.y;
    }
    bool operator==(const Point &p) const
    {
        return equals(x, p.x) && equals(y, p.y);
    }
    bool operator!=(const Point &p) const
    {
        return (!equals(x, p.x)) || (!equals(y, p.y));
    }
};
//ベクトルとして扱う場合
typedef Point Vector;
//線分を表わす構造体
struct Segment
{
    Point p1, p2;
    Segment() = default;
    Segment(Point _p1, Point _p2) : p1(_p1), p2(_p2) {}
};
//直線として扱う場合
typedef Segment Line;
//ベクトルのノルムの2乗
double norm(Vector a)
{
    return a.x * a.x + a.y * a.y;
}
//ベクトルの大きさ
double abs(Vector a)
{
    return sqrt(norm(a));
}
//ベクトルの内積
double dot(Vector a, Vector b)
{
    return a.x * b.x + a.y * b.y;
}
//ベクトルの外積
double cross(Vector a, Vector b)
{
    return a.x * b.y - a.y * b.x;
}
//内積≓0かどうかによるベクトルの直交判定(様々な引数で答えられるように複数用意)
bool Orthogonal(Vector a, Vector b)
{
    return equals(dot(a, b), 0.0);
}
bool Orthogonal(Point a1, Point a2, Point b1, Point b2)
{
    return Orthogonal(a1 - a2, b1 - b2);
}
bool Orthogonal(Segment s1, Segment s2)
{
    return equals(dot(s1.p2 - s1.p1, s2.p2 - s2.p1), 0.0);
}
//外積≓0かどうかによるベクトルの平行判定(様々な引数で答えられるように複数用意)
bool Parallel(Vector a, Vector b)
{
    return equals(cross(a, b), 0.0);
}
bool Parallel(Point a1, Point a2, Point b1, Point b2)
{
    return Parallel(a1 - a2, b1 - b2);
}
bool Parallel(Segment s1, Segment s2)
{
    return equals(cross(s1.p2 - s1.p1, s2.p2 - s2.p1), 0.0);
}
//直線sに対して点pから垂線を引いたときの交点x(直線sに対する点pの射影)を求める
Point Project(Line s, Point p)
{
    Vector base = s.p2 - s.p1, hypo = p - s.p1;
    double scalar = (dot(base, hypo) / norm(base));
    return s.p1 + base * scalar;
}
//直線sに対して点pと線対称の位置にある点p'を求める
Point Reflect(Line s, Point p)
{
    return p + (Project(s, p) - p) * 2.0;
}
//点と点の距離
double GetDist(Point a, Point b)
{
    return abs(a - b);
}
//直線と点の距離
double GetDistLP(Line l, Point p)
{
    return abs(Project(l, p) - p);
}
//線分と点の距離
double GetDistSP(Segment s, Point p)
{
    if (dot(s.p2 - s.p1, p - s.p1) < 0.0)
        return abs(p - s.p1);
    if (dot(s.p1 - s.p2, p - s.p2) < 0.0)
        return abs(p - s.p2);
    return GetDistLP(s, p);
}
//Counter-Clockwise(AOJ:CGL_1_Cより)
static const int COUNTER_CLOCKWISE = 1; //反時計回り
static const int CLOCKWISE = -1;        //時計回り
static const int ONLINE_BACK = 2;       //線分の後ろにある
static const int ONLINE_FRONT = -2;     //線分の前にある
static const int ON_SEGMENT = 0;        //線分上にある
//線分p0p1と点p2の位置関係を求める
int CCW(Point p0, Point p1, Point p2)
{
    Vector a = p1 - p0, b = p2 - p0; //a:p0->p1,b:p0->p2
    if (cross(a, b) > EPS)
        return COUNTER_CLOCKWISE; //外積が正なら反時計回り(sinθ>0)
    if (cross(a, b) < -EPS)
        return CLOCKWISE; //外積が負なら時計回り(sinθ<0)
    if (dot(a, b) < -EPS)
        return ONLINE_BACK; //外積が0かつ内積が負(cosθ<0)
    if (a.norm() < b.norm())
        return ONLINE_FRONT; //p2がp0p1に含まれないときp0p2=p0p1+p1p2となり|a|^2<|b|^2が成立
    return ON_SEGMENT;       //最後にp0p2p1が残る
}
//2本の線分の交差判定(点と線分で2種類用意)
bool Intersect(Point p1, Point p2, Point p3, Point p4)
{
    return (CCW(p1, p2, p3) * CCW(p1, p2, p4) <= 0 && CCW(p3, p4, p1) * CCW(p3, p4, p2) <= 0);
}
bool Intersect(Segment s1, Segment s2)
{
    return Intersect(s1.p1, s1.p2, s2.p1, s2.p2);
}
//線分と線分の距離
double GetDist(Segment s1, Segment s2)
{
    //線分が交差していた場合、距離は0.0
    if (Intersect(s1, s2))
        return 0.0;
    return min(min(GetDistSP(s1, s2.p1), GetDistSP(s1, s2.p2)), min(GetDistSP(s2, s1.p1), GetDistSP(s2, s1.p2)));
}
//交差する2本の直線の交点
Point CrossPointL(Line s1, Line s2)
{
    double a = s1.p1.x, b = s1.p1.y, c = s1.p2.x, d = s1.p2.y;
    double e = s2.p1.x, f = s2.p1.y, g = s2.p2.x, h = s2.p2.y;
    Point ret(((f * g - e * h) * (c - a) - (b * c - a * d) * (g - e)) / ((d - b) * (g - e) - (c - a) * (h - f)), ((f * g - e * h) * (d - b) - (b * c - a * d) * (h - f)) / ((d - b) * (g - e) - (c - a) * (h - f)));
    return ret;
}
//交差する2本の線分の交点
Point CrossPoint(Segment s1, Segment s2)
{
    if (!Intersect(s1, s2))
        cout << "2本の線分は交点を持たない" << endl; //デバッグ用
    Vector base = s2.p2 - s2.p1, hypo1 = s1.p1 - s2.p1, hypo2 = s1.p2 - s2.p1;
    double d1 = fabs(cross(base, hypo1)) / abs(base);
    double d2 = fabs(cross(base, hypo2)) / abs(base);
    Point x = s1.p1 + (s1.p2 - s1.p1) * (d1 / (d1 + d2));
    return x;
}
//座標cを中心とする半径rの円を表わす構造体
struct Circle
{
public:
    Point c;
    double r;
    Circle(Point _c, double _r) : c(_c), r(_r) {}
};
//円と直線の交点の組(接点なら同じものを2つ)を求める(小さい順に入れる)
pair<Point, Point> CrossPoints(Circle c, Line l)
{
    if (GetDistLP(l, c.c) - c.r > EPS)
        cout << "円と直線は交点を持たない" << endl;      //デバッグ用
    Point mid = Project(l, c.c);                         //交点の中間
    Vector e = (l.p2 - l.p1) / abs(l.p2 - l.p1);         //直線lの単位ベクトル
    double halfdist = sqrt(c.r * c.r - norm(c.c - mid)); //三平方の定理で2交点の距離の半分を求める
    Point s = mid + e * halfdist, t = mid - e * halfdist;
    if (s > t)
        swap(s, t);
    return make_pair(s, t); //s<tの順
}
//ベクトルの角度
double arg(Vector p)
{
    return atan2(p.y, p.x);
}
//極座標から直交座標に変換
Vector polar(double r, double theta)
{
    return Point(cos(theta) * r, sin(theta) * r);
}
//円と円の交点の組(接点なら同じものを2つ)を求める(小さい順に入れる)
pair<Point, Point> CrossPoints(Circle c1, Circle c2)
{
    if (GetDist(c1.c, c2.c) - (c1.r + c2.r) > EPS)
        cout << "2円は交点を持たない" << endl;                             //デバッグ用
    double d = abs(c2.c - c1.c);                                           //中心間距離
    double a = acos((d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d * c1.r)); //余弦定理+arccos
    double t = arg(c2.c - c1.c);                                           //偏角
    Point p1 = c1.c + polar(c1.r, t + a), p2 = c1.c + polar(c1.r, t - a);
    if (p1 > p2)
        swap(p1, p2);
    return make_pair(p1, p2);
}
//多角形は点の列として扱う
typedef vector<Point> Polygon;
//点の内包判定
/*
          IN 2
          ON 1
          OUT 0
      */
int Contain(Polygon g, Point p)
{
    int n = g.size();
    bool x = false;
    for (int i = 0; i < n; i++)
    {
        Vector a = g[i] - p;           //p->g[i]
        Vector b = g[(i + 1) % n] - p; //p->g[i+1] 但しn-1の次は0
        if (abs(cross(a, b)) < EPS && dot(a, b) < EPS)
            return 1; //外積の大きさが0で内積0以下なら線分上
        if (a.y > b.y)
            swap(a, b);
        if (a.y < EPS && EPS < b.y && cross(a, b) > EPS)
            x = (!x); //xの真偽値を反転
    }
    return (x ? 2 : 0); //xが真なら2(内包)を、偽なら0(外)を返す
}
//凸包を求める
Polygon AndrewScan(Polygon s)
{
    Polygon u, l;
    if (s.size() < 3)
        return s;
    sort(s.begin(), s.end()); //x座標(y座標)の順でソート
    //xが小さい順に2つ追加
    u.push_back(s[0]);
    u.push_back(s[1]);
    //xが大きい順に2つ追加
    l.push_back(s[s.size() - 1]);
    l.push_back(s[s.size() - 2]);
    //凸包の上部生成
    for (int i = 2; i < s.size(); i++)
    {
        for (int n = u.size(); n >= 2 && CCW(u[n - 2], u[n - 1], s[i]) == COUNTER_CLOCKWISE; n--)
        {
            u.pop_back(); //凸にならなかったら頂点を削除
        }
        u.push_back(s[i]);
    }
    //凸包の下部生成
    for (int i = s.size() - 3; i >= 0; i--)
    {
        for (int n = l.size(); n >= 2 && CCW(l[n - 2], l[n - 1], s[i]) == COUNTER_CLOCKWISE; n--)
        {
            l.pop_back();
        }
        l.push_back(s[i]);
    }
    //反時計回りになるように凸包の点列を生成
    reverse(l.begin(), l.end());
    for (int i = u.size() - 2; i >= 1; i--)
        l.push_back(u[i]);
    return l;
}
    
    
© 2020 kacho65535