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

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!