API tools faq
paste
Advertisement
apiad

Untitled

apiad
May 20th, 2022
969
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 10.20 KB | None | 0 0
raw download clone embed print report
  1. typedef double db;
  2. const db EPS = 1e-9;
  3.  
  4. inline int sign(db a) { return a < -EPS ? -1 : a > EPS; }
  5.  
  6. inline int cmp(db a, db b){ return sign(a-b); }
  7.  
  8. struct P {
  9.     db x, y;
  10.     P() {}
  11.     P(db _x, db _y) : x(_x), y(_y) {}
  12.     P operator+(P p) { return {x + p.x, y + p.y}; }
  13.     P operator-(P p) { return {x - p.x, y - p.y}; }
  14.     P operator*(db d) { return {x * d, y * d}; }
  15.     P operator/(db d) { return {x / d, y / d}; }
  16.  
  17.     bool operator<(P p) const {
  18.         int c = cmp(x, p.x);
  19.         if (c) return c == -1;
  20.         return cmp(y, p.y) == -1;
  21.     }
  22.  
  23.     bool operator==(P o) const{
  24.         return cmp(x,o.x) == 0 && cmp(y,o.y) == 0;
  25.     }
  26.  
  27.     db dot(P p) { return x * p.x + y * p.y; }
  28.     db det(P p) { return x * p.y - y * p.x; }
  29.      
  30.     db distTo(P p) { return (*this-p).abs(); }
  31.     db alpha() { return atan2(y, x); }
  32.     void read() { cin>>x>>y; }
  33.     void write() {cout<<"("<<x<<","<<y<<")"<<endl;}
  34.     db abs() { return sqrt(abs2());}
  35.     db abs2() { return x * x + y * y; }
  36.     P rot90() { return P(-y,x);}
  37.     P unit() { return *this/abs(); }
  38.     int quad() const { return sign(y) == 1 || (sign(y) == 0 && sign(x) >= 0); }
  39.     P rot(db an){ return {x*cos(an)-y*sin(an),x*sin(an) + y*cos(an)}; }
  40. };
  41.  
  42. struct L{ //ps[0] -> ps[1]
  43.     P ps[2];
  44.     P& operator[](int i) { return ps[i]; }
  45.     P dir() { return ps[1] - ps[0]; }
  46.     L (P a,P b) {
  47.         ps[0]=a;
  48.         ps[1]=b;
  49.     }
  50.     bool include(P p) { return sign((ps[1] - ps[0]).det(p - ps[0])) > 0; }
  51.     L push(){ // push eps outward
  52.         const double eps = 1e-8;
  53.         P delta = (ps[1] - ps[0]).rot90().unit() * eps;
  54.         return {ps[0] + delta, ps[1] + delta};
  55.     }
  56. };
  57.  
  58. #define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
  59. #define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))
  60.  
  61. bool chkLL(P p1, P p2, P q1, P q2) {
  62.     db a1 = cross(q1, q2, p1), a2 = -cross(q1, q2, p2);
  63.     return sign(a1+a2) != 0;
  64. }
  65.  
  66. P isLL(P p1, P p2, P q1, P q2) {
  67.     db a1 = cross(q1, q2, p1), a2 = -cross(q1, q2, p2);
  68.     return (p1 * a2 + p2 * a1) / (a1 + a2);
  69. }
  70.  
  71. P isLL(L l1,L l2){ return isLL(l1[0],l1[1],l2[0],l2[1]); }
  72.  
  73. bool intersect(db l1,db r1,db l2,db r2){
  74.     if(l1>r1) swap(l1,r1); if(l2>r2) swap(l2,r2);
  75.     return !( cmp(r1,l2) == -1 || cmp(r2,l1) == -1 );
  76. }
  77.  
  78. bool isSS(P p1, P p2, P q1, P q2){
  79.     return intersect(p1.x,p2.x,q1.x,q2.x) && intersect(p1.y,p2.y,q1.y,q2.y) &&
  80.     crossOp(p1,p2,q1) * crossOp(p1,p2,q2) <= 0 && crossOp(q1,q2,p1)
  81.             * crossOp(q1,q2,p2) <= 0;
  82. }
  83.  
  84. bool isSS_strict(P p1, P p2, P q1, P q2){
  85.     return crossOp(p1,p2,q1) * crossOp(p1,p2,q2) < 0 && crossOp(q1,q2,p1)
  86.             * crossOp(q1,q2,p2) < 0;
  87. }
  88.  
  89. bool isMiddle(db a, db m, db b) {
  90.     return sign(a - m) == 0 || sign(b - m) == 0 || (a < m != b < m);
  91. }
  92.  
  93. bool isMiddle(P a, P m, P b) {
  94.     return isMiddle(a.x, m.x, b.x) && isMiddle(a.y, m.y, b.y);
  95. }
  96.  
  97. bool onSeg(P p1, P p2, P q){
  98.     return crossOp(p1,p2,q) == 0 && isMiddle(p1, q, p2);
  99. }
  100.  
  101. bool onSeg_strict(P p1, P p2, P q){
  102.     return crossOp(p1,p2,q) == 0 && sign((q-p1).dot(p1-p2)) * sign((q-p2).dot(p1-p2)) < 0;
  103. }
  104.  
  105. P proj(P p1, P p2, P q) {
  106.     P dir = p2 - p1;
  107.     return p1 + dir * (dir.dot(q - p1) / dir.abs2());
  108. }
  109.  
  110. P reflect(P p1, P p2, P q){
  111.     return proj(p1,p2,q) * 2 - q;
  112. }
  113.  
  114. db nearest(P p1,P p2,P q){
  115.     if (p1==p2) return p1.distTo(q);
  116.     P h = proj(p1,p2,q);
  117.     if(isMiddle(p1,h,p2))
  118.         return q.distTo(h);
  119.     return min(p1.distTo(q),p2.distTo(q));
  120. }
  121.  
  122. db disSS(P p1, P p2, P q1, P q2){
  123.     if(isSS(p1,p2,q1,q2)) return 0;
  124.     return min(min(nearest(p1,p2,q1),nearest(p1,p2,q2)), min(nearest(q1,q2,p1),nearest(q1,q2,p2)));
  125. }
  126.  
  127. db rad(P p1,P p2){
  128.     return atan2l(p1.det(p2),p1.dot(p2));
  129. }
  130.  
  131. db incircle(P p1, P p2, P p3){
  132.     db A = p1.distTo(p2);
  133.     db B = p2.distTo(p3);
  134.     db C = p3.distTo(p1);
  135.     return sqrtl(A*B*C/(A+B+C));
  136. }
  137.  
  138. //polygon
  139.  
  140. db area(vector<P> ps){
  141.     db ret = 0; rep(i,0,ps.size()) ret += ps[i].det(ps[(i+1)%ps.size()]);
  142.     return ret/2;
  143. }
  144.  
  145. int contain(vector<P> ps, P p){ //2:inside,1:on_seg,0:outside
  146.     int n = ps.size(), ret = 0;
  147.     rep(i,0,n){
  148.         P u=ps[i],v=ps[(i+1)%n];
  149.         if(onSeg(u,v,p)) return 1;
  150.         if(cmp(u.y,v.y)<=0) swap(u,v);
  151.         if(cmp(p.y,u.y) >0 || cmp(p.y,v.y) <= 0) continue;
  152.         ret ^= crossOp(p,u,v) > 0;
  153.     }
  154.     return ret*2;
  155. }
  156.  
  157. vector<P> convexHull(vector<P> ps) {
  158.     int n = ps.size(); if(n <= 1) return ps;
  159.     sort(ps.begin(), ps.end());
  160.     vector<P> qs(n * 2); int k = 0;
  161.     for (int i = 0; i < n; qs[k++] = ps[i++])
  162.         while (k > 1 && crossOp(qs[k - 2], qs[k - 1], ps[i]) <= 0) --k;
  163.     for (int i = n - 2, t = k; i >= 0; qs[k++] = ps[i--])
  164.         while (k > t && crossOp(qs[k - 2], qs[k - 1], ps[i]) <= 0) --k;
  165.     qs.resize(k - 1);
  166.     return qs;
  167. }
  168.  
  169. vector<P> convexHullNonStrict(vector<P> ps) {
  170.     //caution: need to unique the Ps first
  171.     int n = ps.size(); if(n <= 1) return ps;
  172.     sort(ps.begin(), ps.end());
  173.     vector<P> qs(n * 2); int k = 0;
  174.     for (int i = 0; i < n; qs[k++] = ps[i++])
  175.         while (k > 1 && crossOp(qs[k - 2], qs[k - 1], ps[i]) < 0) --k;
  176.     for (int i = n - 2, t = k; i >= 0; qs[k++] = ps[i--])
  177.         while (k > t && crossOp(qs[k - 2], qs[k - 1], ps[i]) < 0) --k;
  178.     qs.resize(k - 1);
  179.     return qs;
  180. }
  181.  
  182. db convexDiameter(vector<P> ps){
  183.     int n = ps.size(); if(n <= 1) return 0;
  184.     int is = 0, js = 0; rep(k,1,n) is = ps[k]<ps[is]?k:is, js = ps[js] < ps[k]?k:js;
  185.     int i = is, j = js;
  186.     db ret = ps[i].distTo(ps[j]);
  187.     do{
  188.         if((ps[(i+1)%n]-ps[i]).det(ps[(j+1)%n]-ps[j]) >= 0)
  189.             (++j)%=n;
  190.         else
  191.             (++i)%=n;
  192.         ret = max(ret,ps[i].distTo(ps[j]));
  193.     }while(i!=is || j!=js);
  194.     return ret;
  195. }
  196.  
  197. vector<P> convexCut(const vector<P>&ps, P q1, P q2) {
  198.     vector<P> qs;
  199.     int n = ps.size();
  200.     rep(i,0,n){
  201.         P p1 = ps[i], p2 = ps[(i+1)%n];
  202.         int d1 = crossOp(q1,q2,p1), d2 = crossOp(q1,q2,p2);
  203.         if(d1 >= 0) qs.pb(p1);
  204.         if(d1 * d2 < 0) qs.pb(isLL(p1,p2,q1,q2));
  205.     }
  206.     return qs;
  207. }
  208.  
  209. //min_dist
  210.  
  211. db min_dist(vector<P>&ps,int l,int r){
  212.     if(r-l<=5){
  213.         db ret = 1e100;
  214.         rep(i,l,r) rep(j,l,i) ret = min(ret,ps[i].distTo(ps[j]));
  215.         return ret;
  216.     }
  217.     int m = (l+r)>>1;
  218.     db ret = min(min_dist(ps,l,m),min_dist(ps,m,r));
  219.     vector<P> qs; rep(i,l,r) if(abs(ps[i].x-ps[m].x)<= ret) qs.pb(ps[i]);
  220.     sort(qs.begin(), qs.end(),[](P a,P b) -> bool {return a.y<b.y; });
  221.     rep(i,1,qs.size()) for(int j=i-1;j>=0&&qs[j].y>=qs[i].y-ret;--j)
  222.         ret = min(ret,qs[i].distTo(qs[j]));
  223.     return ret;
  224. }
  225.  
  226. int type(P o1,db r1,P o2,db r2){
  227.     db d = o1.distTo(o2);
  228.     if(cmp(d,r1+r2) == 1) return 4;
  229.     if(cmp(d,r1+r2) == 0) return 3;
  230.     if(cmp(d,abs(r1-r2)) == 1) return 2;
  231.     if(cmp(d,abs(r1-r2)) == 0) return 1;
  232.     return 0;
  233. }
  234.  
  235. vector<P> isCL(P o,db r,P p1,P p2){
  236.     if (cmp(abs((o-p1).det(p2-p1)/p1.distTo(p2)),r)>0) return {};
  237.     db x = (p1-o).dot(p2-p1), y = (p2-p1).abs2(), d = x * x - y * ((p1-o).abs2() - r*r);
  238.     d = max(d,(db)0.0); P m = p1 - (p2-p1)*(x/y), dr = (p2-p1)*(sqrt(d)/y);
  239.     return {m-dr,m+dr}; //along dir: p1->p2
  240. }
  241.  
  242. vector<P> isCC(P o1, db r1, P o2, db r2) { //need to check whether two circles are the same
  243.     db d = o1.distTo(o2);
  244.     if (cmp(d, r1 + r2) == 1) return {};
  245.     if (cmp(d,abs(r1-r2))==-1) return {};
  246.     d = min(d, r1 + r2);
  247.     db y = (r1 * r1 + d * d - r2 * r2) / (2 * d), x = sqrt(r1 * r1 - y * y);
  248.     P dr = (o2 - o1).unit();
  249.     P q1 = o1 + dr * y, q2 = dr.rot90() * x;
  250.     return {q1-q2,q1+q2};//along circle 1
  251. }
  252.  
  253. vector<P> tanCP(P o, db r, P p) {
  254.     db x = (p - o).abs2(), d = x - r * r;
  255.     if (sign(d) <= 0) return {}; // on circle => no tangent
  256.     P q1 = o + (p - o) * (r * r / x);
  257.     P q2 = (p - o).rot90() * (r * sqrt(d) / x);
  258.     return {q1-q2,q1+q2}; //counter clock-wise
  259. }
  260.  
  261.  
  262. vector<L> extanCC(P o1, db r1, P o2, db r2) {
  263.     vector<L> ret;
  264.     if (cmp(r1, r2) == 0) {
  265.         P dr = (o2 - o1).unit().rot90() * r1;
  266.         ret.pb(L(o1 + dr, o2 + dr)), ret.pb(L(o1 - dr, o2 - dr));
  267.     } else {
  268.         P p = (o2 * r1 - o1 * r2) / (r1 - r2);
  269.         vector<P> ps = tanCP(o1, r1, p), qs = tanCP(o2, r2, p);
  270.         rep(i,0,min(ps.size(),qs.size())) ret.pb(L(ps[i], qs[i])); //c1 counter-clock wise
  271.     }
  272.     return ret;
  273. }
  274.  
  275. vector<L> intanCC(P o1, db r1, P o2, db r2) {
  276.     vector<L> ret;
  277.     P p = (o1 * r2 + o2 * r1) / (r1 + r2);
  278.     vector<P> ps = tanCP(o1,r1,p), qs = tanCP(o2,r2,p);
  279.     rep(i,0,min(ps.size(),qs.size())) ret.pb(L(ps[i], qs[i])); //c1 counter-clock wise
  280.     return ret;
  281. }
  282.  
  283. db areaCT(db r, P p1, P p2){
  284.     vector<P> is = isCL(P(0,0),r,p1,p2);
  285.     if(is.empty()) return r*r*rad(p1,p2)/2;
  286.     bool b1 = cmp(p1.abs2(),r*r) == 1, b2 = cmp(p2.abs2(), r*r) == 1;
  287.     if(b1 && b2){
  288.         if(sign((p1-is[0]).dot(p2-is[0])) <= 0 &&
  289.             sign((p1-is[0]).dot(p2-is[0])) <= 0)
  290.         return r*r*(rad(p1,is[0]) + rad(is[1],p2))/2 + is[0].det(is[1])/2;
  291.         else return r*r*rad(p1,p2)/2;
  292.     }
  293.     if(b1) return (r*r*rad(p1,is[0]) + is[0].det(p2))/2;
  294.     if(b2) return (p1.det(is[1]) + r*r*rad(is[1],p2))/2;
  295.     return p1.det(p2)/2;
  296. }
  297.  
  298. bool parallel(L l0, L l1) { return sign( l0.dir().det( l1.dir() ) ) == 0; }
  299.  
  300. bool sameDir(L l0, L l1) { return parallel(l0, l1) && sign(l0.dir().dot(l1.dir()) ) == 1; }
  301.  
  302. bool cmp (P a,  P b) {
  303.     if (a.quad() != b.quad()) {
  304.         return a.quad() < b.quad();
  305.     } else {
  306.         return sign( a.det(b) ) > 0;
  307.     }
  308. }
  309.  
  310. bool operator < (L l0, L l1) {
  311.     if (sameDir(l0, l1)) {
  312.         return l1.include(l0[0]);
  313.     } else {
  314.         return cmp( l0.dir(), l1.dir() );
  315.     }
  316. }
  317.  
  318. bool check(L u, L v, L w) {
  319.     return w.include(isLL(u,v));
  320. }
  321.  
  322. vector<P> halfPlaneIS(vector<L> &l) {
  323.     sort(l.begin(), l.end());
  324.     deque<L> q;
  325.     for (int i = 0; i < (int)l.size(); ++i) {
  326.         if (i && sameDir(l[i], l[i - 1])) continue;
  327.         while (q.size() > 1 && !check(q[q.size() - 2], q[q.size() - 1], l[i])) q.pop_back();
  328.         while (q.size() > 1 && !check(q[1], q[0], l[i])) q.pop_front();
  329.         q.push_back(l[i]);
  330.     }
  331.     while (q.size() > 2 && !check(q[q.size() - 2], q[q.size() - 1], q[0])) q.pop_back();
  332.     while (q.size() > 2 && !check(q[1], q[0], q[q.size() - 1])) q.pop_front();
  333.     vector<P> ret;
  334.     for (int i = 0; i < (int)q.size(); ++i) ret.push_back(isLL(q[i], q[(i + 1) % q.size()]));
  335.     return ret;
  336. }
  337.  
  338. P inCenter(P A, P B, P C) {
  339.     double a = (B - C).abs(), b = (C - A).abs(), c = (A - B).abs();
  340.     return (A * a + B * b + C * c) / (a + b + c);
  341. }
  342.  
  343. P circumCenter(P a, P b, P c) {
  344.     P bb = b - a, cc = c - a;
  345.     double db = bb.abs2(), dc = cc.abs2(), d = 2 * bb.det(cc);
  346.     return a - P(bb.y * dc - cc.y * db, cc.x * db - bb.x * dc) / d;
  347. }
  348.  
  349. P othroCenter(P a, P b, P c) {
  350.     P ba = b - a, ca = c - a, bc = b - c;
  351.     double Y = ba.y * ca.y * bc.y,
  352.     A = ca.x * ba.y - ba.x * ca.y,
  353.     x0 = (Y + ca.x * ba.y * b.x - ba.x * ca.y * c.x) / A,
  354.     y0 = -ba.x * (x0 - c.x) / ba.y + ca.y;
  355.     return {x0, y0};
  356. }
  357.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!