API tools faq
paste
Matrix_code

geo - Template

Matrix_code
Oct 9th, 2016 (edited)
184
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 10.01 KB | None | 0 0
raw download clone embed print report
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define PI          acos(-1.0)
  4. #define EPS         1e-6
  5. #define Z(x)        fabs(x)<EPS
  6. #define E(x,y)      Z(x-y)
  7. #define RAD(x)      (x)*PI/180
  8. #define PRE(x)      ((x)==0?n-1:x-1)
  9. #define NEX(x)      ((x)==n-1?0:x+1)
  10. int dcmp(double x) { return x > EPS ? 1: x < -EPS ? -1 : 0; }
  11. struct point {
  12.     double  x,y;
  13.     point(double x=0, double y=0) : x(x), y(y) {}
  14.     point(const point &p) : x(p.x), y(p.y)    {}
  15.     bool operator <  (const point &p) const { return (x<p.x)||(E(x,p.x) && y < p.y );}
  16.     bool operator == (const point &p) const { return E(x,p.x)&&E(y,p.y);    }
  17.     point operator + (const point &p) const { return point(x+p.x , y + p.y);}
  18.     point operator - (const point &p) const { return point(x-p.x , y - p.y);}
  19.     point operator * (double c) const { return point(x*c,y*c);}
  20.     point operator / (double c) const { return point(x/c,y/c);}
  21.     void input(){scanf("%lf %lf",&x,&y);}
  22. };
  23.  
  24. double dot(point p,point q){ return p.x * q.x + p.y * q.y;}
  25. double cross(point p,point q) {return p.x*q.y - p.y*q.x;}
  26. double sqdist(point p,point q) {return dot(p-q , p-q);}
  27. double dist(point p,point q)  { return sqrt(sqdist(p,q));}
  28. double mag(point p) {return sqrt(p.x*p.x+p.y*p.y);}
  29. double angle(point a,point b) { return acos(dot(a,b)/mag(a)/mag(b));}
  30.  
  31. // rotate a point cross or CW around the origin
  32. point rotccw90(point p){return point(-p.y,p.x);}
  33. point rotcw90(point p)   {return point(p.y,-p.x);}
  34. point rotccw(point p,double ang) {return point(p.x*cos(ang) - p.y*sin(ang) , p.x*sin(ang) + p.y * cos(ang));}
  35.  
  36. // project point c onto line through a and b
  37. point PPL (point a,point b,point c){return a + (b-a)* dot(c-a,b-a) / dot(b-a,b-a);}
  38.  
  39. // project point c onto line segment through a and b
  40. point PPS(point a, point b, point c) {
  41.     double r = dot(b-a,b-a);
  42.     if (fabs(r) < EPS) return a;
  43.     r = dot(c-a, b-a)/r;
  44.     if (r < 0) return a;
  45.     if (r > 1) return b;
  46.     return a + (b-a)*r;
  47. }
  48. // determine if lines from a to b and c to d are parallel or collinear
  49. bool IsLP(point a, point b, point c, point d) {
  50.     return fabs(cross(b-a, c-d)) < EPS;
  51. }
  52.  
  53. bool IsLC(point a, point b, point c, point d) {
  54.     return IsLP(a, b, c, d)
  55.     && fabs(cross(a-b, a-c)) < EPS
  56.     && fabs(cross(c-d, c-a)) < EPS;
  57. }
  58. // Checks if p lies on the segment connection ab
  59. bool onSegment(point a, point b, point p)
  60. {
  61.     return Z(cross(a-p,b-p)) && dot(a-p,b-p) < 0;
  62. }
  63. // determine if line segment from a to b intersects with
  64. // line segment from c to d
  65. bool SGIN(point a, point b, point c, point d) {
  66.     if (IsLC(a, b, c, d)) {
  67.         if (sqdist(a, c) < EPS || sqdist(a, d) < EPS ||
  68.             sqdist(b, c) < EPS || sqdist(b, d) < EPS) return true;
  69.         if (dot(c-a, c-b) > 0 && dot(d-a, d-b) > 0 && dot(c-b, d-b) > 0)
  70.             return false;
  71.         return true;
  72.     }
  73.     if (cross(d-a, b-a) * cross(c-a, b-a) > 0) return false;
  74.     if (cross(a-c, d-c) * cross(b-c, d-c) > 0) return false;
  75.     return true;
  76. }
  77.  
  78.  
  79. point LLIN(point a, point b, point c, point d) {
  80.     b=b-a; d=c-d; c=c-a;
  81.     return a + b*cross(c, d)/cross(b, d);
  82. }
  83.  
  84. // compute center of circle given three points
  85. point CCC(point a, point b, point c) {
  86.     b=(a+b)/2;
  87.     c=(a+c)/2;
  88.     return LLIN(b, b+rotcw90(a-b), c, c+rotcw90(a-c));
  89. }
  90.  
  91. //point in polygon(1-Strictly Interior,0-Strictly Exterior
  92. bool PIPoly(const vector<point> &p, point q) {
  93.     bool c = 0;
  94.     for (int i = 0; i < p.size(); i++){
  95.         int j = (i+1)%p.size();
  96.         if (((p[i].y <= q.y && q.y < p[j].y) ||
  97.              (p[j].y <= q.y && q.y < p[i].y)) &&
  98.             q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (p[j].y - p[i].y))
  99.             c = !c;
  100.     }
  101.     return c;
  102. }
  103.  
  104. // determine if point is on the boundary of a polygon
  105. bool POPoly(const vector<point> &p, point q) {
  106.     for (int i = 0; i < p.size(); i++)
  107.         if (sqdist(PPS(p[i], p[(i+1)%p.size()], q), q) < EPS)
  108.             return true;
  109.     return false;
  110. }
  111.  
  112. // compute intersection of line through points a and b with
  113. // circle centered at c with radius r > 0
  114. vector<point> CLIN(point a, point b, point c, double r) {
  115.     vector<point> ret;
  116.     b = b-a;
  117.     a = a-c;
  118.     double A = dot(b, b);
  119.     double B = dot(a, b);
  120.     double C = dot(a, a) - r*r;
  121.     double D = B*B - A*C;
  122.     if (D < -EPS) return ret;
  123.     ret.push_back(c+a+b*(-B+sqrt(D+EPS))/A);
  124.     if (D > EPS)
  125.         ret.push_back(c+a+b*(-B-sqrt(D))/A);
  126.     return ret;
  127. }
  128.  
  129. // compute intersection of circle centered at a with radius r
  130. // with circle centered at b with radius R
  131. vector<point> CCIN(point a, point b, double r, double R) {
  132.     vector<point> ret;
  133.     double d = sqrt(sqdist(a, b));
  134.     if (d > r+R || d+min(r, R) < max(r, R)) return ret;
  135.     double x = (d*d-R*R+r*r)/(2*d);
  136.     double y = sqrt(r*r-x*x);
  137.     point v = (b-a)/d;
  138.     ret.push_back(a+v*x + rotccw90(v)*y);
  139.     if (y > 0)
  140.         ret.push_back(a+v*x - rotccw90(v)*y);
  141.     return ret;
  142. }
  143. // p = a + t * (b-a). find t.
  144. double getT(point a, point dir, point p)
  145. {
  146.       if(dcmp(dir.x) == 0) return (p.y - a.y) / dir.y;
  147.       return (p.x - a.x)/ dir.x;
  148. }
  149. // Gives area of a circle sector passing through a,b
  150. double SectorArea(point r, point a, point b, double R)
  151. {
  152.       double ang = angle(a-r,b-r);
  153.       return R*R*ang/2;
  154. }
  155. // Common area of a circle and and a segment(a,b)
  156. double TRICA(point r, point a, point b, double R)
  157. {
  158.       double ra = dist(r, a) , rb = dist(r, b);
  159.       if(ra < R + EPS && rb < R + EPS) return cross(a - r, b - r) / 2;
  160.       if(dcmp(cross(a-r,b-r)) == 0) return 0;
  161.       double rtos = dist(r, PPS(a,b,r));
  162.       if(rtos > R - EPS) return SectorArea(r, a, b, R);
  163.       vector< point > ins = CLIN(a,b,r,R);
  164.       if(ins.size() < 2) return SectorArea(r, a, b, R);
  165.  
  166.       point ta = ins[0], tb = ins[1];
  167.       double t1 = getT(a, b - a, ta) , t2 = getT(a, b-a, tb);
  168.       if(t1 > t2) swap(ta, tb);
  169.  
  170.       if(ra < R + EPS) return cross(a - r, tb - r) / 2 + SectorArea(r, tb , b,R);
  171.       if(rb < R + EPS) return cross(ta - r, b - r) / 2 + SectorArea(r , a, ta ,R);
  172.       return cross(ta - r, tb - r) / 2 + SectorArea(r, a, ta, R) + SectorArea(r, tb, b,R);
  173. }
  174. // Simple polygon intersection area with a circle.
  175. double SPICA(vector< point > &p, point r, double R )
  176. {
  177.       double ret = 0;
  178.       int n = p.size();
  179.       for(int i= 0; i < n; i ++) {
  180.             int turn = dcmp(cross(p[i] - r, p[NEX(i)] -  r));
  181.             if(turn > 0 ) ret += TRICA(r, p[i] , p[NEX(i)],R);
  182.             else ret -= TRICA(r,p[NEX(i)], p[i],R);
  183.  
  184.       }
  185.       return fabs(ret);
  186.  
  187. }
  188. // This code computes the area or centroid of a (possibly nonconvex)
  189. // polygon, assuming that the coordinates are listed in a clockwise or
  190. // counterclockwise fashion.  Note that the centroid is often known as
  191. // the "center of gravity" or "center of mass".
  192.  
  193. double ComputeSignedArea(const vector<point> &p) {
  194.     double area = 0;
  195.     for(int i = 0; i < p.size(); i++) {
  196.         int j = (i+1) % p.size();
  197.         area += p[i].x*p[j].y - p[j].x*p[i].y;
  198.     }
  199.     return area / 2.0;
  200. }
  201.  
  202. double CAR(const vector<point> &p) {
  203.     return fabs(ComputeSignedArea(p));
  204. }
  205. point ComputeCentroid(const vector<point> &p) {
  206.   point c(0,0);
  207.   double scale = 6.0 * ComputeSignedArea(p);
  208.   for (int i = 0; i < p.size(); i++){
  209.     int j = (i+1) % p.size();
  210.     c = c + (p[i]+p[j])*(p[i].x*p[j].y - p[j].x*p[i].y);
  211.   }
  212.   return c / scale;
  213. }
  214.  
  215. // centroid
  216. point CCN(const vector<point> &p) {
  217.     point c(0,0);
  218.     double scale = 6.0 * ComputeSignedArea(p);
  219.     for (int i = 0; i < p.size(); i++){
  220.         int j = (i+1) % p.size();
  221.         c = c + (p[i]+p[j])*(p[i].x*p[j].y - p[j].x*p[i].y);
  222.     }
  223.     return c / scale;
  224. }
  225.  
  226. double myAngle(point a,point b)
  227. {
  228.     double ang = angle(a,b);
  229.     if(dcmp(cross(a,b)) >= 0) return ang;
  230.     return 2*PI - ang;
  231. }
  232. // tests whether or not a given polygon (in CW or cross order) is simple
  233. bool IsSimple(const vector<point> &p) {
  234.     for (int i = 0; i < p.size(); i++) {
  235.         for (int k = i+1; k < p.size(); k++) {
  236.             int j = (i+1) % p.size();
  237.             int l = (k+1) % p.size();
  238.             if (i == l || j == k) continue;
  239.             if (SGIN(p[i], p[j], p[k], p[l]))
  240.                 return false;
  241.         }
  242.     }
  243.     return true;
  244. }
  245.  
  246. bool cw(const point &a, const point &b, const point &c) {
  247.     return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x) < 0;
  248. }
  249. vector<point> convexHull(vector<point>& p) {
  250.     int n = (int)p.size();
  251.     if (n <= 1)
  252.         return p;
  253.     sort(p.begin(), p.end());
  254.     int cnt = 0;
  255.     vector<point> q(n * 2);
  256.     for (int i = 0; i < n; q[cnt++] = p[i++])
  257.         for (; cnt >= 2 && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt)
  258.             ;
  259.     for (int i = n - 2, t = cnt; i >= 0; q[cnt++] = p[i--])
  260.         for (; cnt > t && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt)
  261.             ;
  262.     q.resize(cnt - 1 - (q[0] == q[1]));
  263.     return q;
  264. }
  265. /*==============================================================================*/
  266.  
  267.  
  268.  
  269. // For integer Calculation
  270. struct Pint{
  271.     long long x,y;
  272.     Pint(){}
  273.     Pint(long long x,long long y):x(x),y(y){}
  274.     bool operator<(const Pint &p) const { return x<p.x||(x==p.x&& y<p.y);}
  275.     bool operator==(const Pint &p) const {return x==p.x && y==p.y; }
  276.     Pint operator-(const Pint &p) const { return {x-p.x,y-p.y};           }
  277. };
  278. long long dot(Pint a,Pint b) { return a.x*b.x + a.y*b.y; }
  279. long long cross(Pint a,Pint b) { return a.x*b.y - a.y*b.x;}
  280. bool cw(const Pint &a, const Pint &b, const Pint &c) {
  281.     return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x) < 0;
  282. }
  283. vector<Pint> convexHull(vector<Pint> &p)
  284. {
  285.  
  286.     int n = (int)p.size();
  287.     if (n <= 1)
  288.         return p;
  289.     sort(p.begin(), p.end());
  290.     int cnt = 0;
  291.     vector<Pint> q(n * 2);
  292.     for (int i = 0; i < n; q[cnt++] = p[i++])
  293.         for (; cnt >= 2 && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt)
  294.             ;
  295.     for (int i = n - 2, t = cnt; i >= 0; q[cnt++] = p[i--])
  296.         for (; cnt > t && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt)
  297.             ;
  298.     q.resize(cnt - 1 - (q[0] == q[1]));
  299.     return q;
  300. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!