1 const double eps=1e-10; 2 const double PI=acos(-1.0); 3 using namespace std; 4 struct Point{ 5 double x; 6 double y; 7 Point(double x=0,double y=0):x(x),y(y){} 8 void operator<<(Point &A) {cout<<A.x<<' '<<A.y<<endl;} 9 }; 10 11 int dcmp(double x) {return (x>eps)-(x<-eps); } 12 int sgn(double x) {return (x>eps)-(x<-eps); } 13 typedef Point Vector; 14 15 Vector operator +(Vector A,Vector B) { return Vector(A.x+B.x,A.y+B.y);} 16 Vector operator -(Vector A,Vector B) { return Vector(A.x-B.x,A.y-B.y); } 17 Vector operator *(Vector A,double p) { return Vector(A.x*p,A.y*p); } 18 Vector operator /(Vector A,double p) {return Vector(A.x/p,A.y/p);} 19 ostream &operator<<(ostream & out,Point & P) { out<<P.x<<' '<<P.y<<endl; return out;} 20 bool operator< (const Point &A,const Point &B) { return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0); } 21 bool operator== ( const Point &A,const Point &B) { return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;} 22 23 double Dot(Vector A,Vector B) {return A.x*B.x+A.y*B.y;} 24 double Cross(Vector A,Vector B) {return A.x*B.y-B.x*A.y; } 25 double Length(Vector A) { return sqrt(Dot(A, A));} 26 double Angle(Vector A,Vector B) {return acos(Dot(A,B)/Length(A)/Length(B));} 27 double Area2(Point A,Point B,Point C ) {return fabs(Cross(B-A, C-A));} 28 29 Vector Rotate(Vector A,double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));} 30 Vector Normal(Vector A) {double L=Length(A);return Vector(-A.y/L,A.x/L);} 31 32 Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){ 33 Vector u=P-Q; 34 double t=Cross(w, u)/Cross(v,w); 35 return P+v*t; 36 }//输入两个点斜式方程输出交点 37 38 double DistanceToLine(Point P,Point A,Point B){ 39 Vector v1=P-A; Vector v2=B-A; 40 return fabs(Cross(v1,v2))/Length(v2); 41 }//输入三个点，输出P到直线AB的距离 42 43 double DistanceToSegment(Point P,Point A,Point B){ 44 if(A==B) return Length(P-A); 45 Vector v1=B-A; 46 Vector v2=P-A; 47 Vector v3=P-B; 48 if(dcmp(Dot(v1,v2))==-1) return Length(v2); 49 else if(Dot(v1,v3)>0) return Length(v3); 50 else return DistanceToLine(P, A, B); 51 }//输入三个点，输出P到线段AB的距离 52 53 Point GetLineProjection(Point P,Point A,Point B){ 54 Vector v=B-A; 55 Vector v1=P-A; 56 double t=Dot(v,v1)/Dot(v,v); 57 return A+v*t; 58 }//输入三个点，输出P在直线AB上的投影点 59 60 bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){ 61 double c1=Cross(b1-a1, a2-a1); 62 double c2=Cross(b2-a1, a2-a1); 63 double c3=Cross(a1-b1, b2-b1); 64 double c4=Cross(a2-b1, b2-b1); 65 return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0; 66 }//输入四个点，判断两条线段是否有交点 67 68 bool OnSegment(Point P,Point A,Point B){ 69 return dcmp(Cross(P-A, P-B))==0&&dcmp(Dot(P-A,P-B))<=0; 70 }//输入三个点，判断P是否在线段AB上 71 72 double PolygonArea(Point *p,int n){ 73 double area=0; 74 for(int i=1;i<n-1;i++){ 75 area+=Cross(p[i]-p[0], p[i+1]-p[0]); 76 } 77 return area/2; 78 }//输入p[i]首地址和p[]的数量，（p[0]到p[n-1]），输出多边形面 79 80 Point read_point(){ 81 Point P; 82 scanf("%lf%lf",&P.x,&P.y); 83 return P; 84 }//输入一个点 85 86 // ---------------与圆有关的-------- 87 88 struct Circle{ 89 Point c; 90 double r; 91 Circle(Point c=Point(0,0),double r=0):c(c),r(r) {} 92 Point point(double a){ 93 return Point(c.x+r*cos(a),c.y+r*sin(a)); 94 } 95 }; 96 97 struct Line{ 98 Point p; 99 Vector v; 100 Line(Point p=Point(0,0),Vector v=Vector(0,1)):p(p),v(v) {} 101 Point point(double t){ 102 return Point(p+v*t); 103 } 104 }; 105 106 int getLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point> &sol){ 107 double a=L.v.x; 108 double b=L.p.x-C.c.x; 109 double c=L.v.y; 110 double d=L.p.y-C.c.y; 111 double e=a*a+c*c; 112 double f=2*(a*b+c*d); 113 double g=b*b+d*d-C.r*C.r; 114 double delta=f*f-4*e*g; 115 116 if(dcmp(delta)<0) return 0; 117 if(dcmp(delta)==0){ 118 t1=t2=-f/(2*e); 119 sol.push_back(L.point(t1)); 120 return 1; 121 } 122 123 else{ 124 t1=(-f-sqrt(delta))/(2*e); 125 t2=(-f+sqrt(delta))/(2*e); 126 127 sol.push_back(L.point(t1)); 128 sol.push_back(L.point(t2)); 129 130 return 2; 131 } 132 } 133 134 // 向量极角公式 135 136 double angle(Vector v) {return atan2(v.y,v.x);} 137 138 int getCircleCircleIntersection(Circle C1,Circle C2,vector<Point> &sol) 139 { 140 double d=Length(C1.c-C2.c); 141 142 if(dcmp(d)==0) 143 { 144 if(dcmp(C1.r-C2.r)==0) return -1; // 重合 145 else return 0; // 内含 0 个公共点 146 } 147 148 if(dcmp(C1.r+C2.r-d)<0) return 0; // 外离 149 if(dcmp(fabs(C1.r-C2.r)-d)>0) return 0; // 内含 150 151 double a=angle(C2.c-C1.c); 152 double da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d)); 153 154 Point p1=C1.point(a-da); 155 Point p2=C1.point(a+da); 156 157 sol.push_back(p1); 158 159 if(p1==p2) return 1; // 相切 160 else{ 161 sol.push_back(p2); 162 return 2; 163 } 164 } 165 166 167 // 求点到圆的切线 168 169 int getTangents(Point p,Circle C,Vector *v){ 170 Vector u=C.c-p; 171 172 double dist=Length(u); 173 174 if(dcmp(dist-C.r)<0) return 0; 175 176 else if(dcmp(dist-C.r)==0){ 177 v[0]=Rotate(u,PI/2); 178 return 1; 179 } 180 181 else{ 182 double ang=asin(C.r/dist); 183 v[0]=Rotate(u,-ang); 184 v[1]=Rotate(u,+ang); 185 return 2; 186 } 187 } 188 189 // 求两圆公切线 190 int getTangents(Circle A,Circle B,Point *a,Point *b){ 191 int cnt=0; 192 193 if(A.r<B.r){ 194 swap(A,B); swap(a, b); // 有时需标记 195 } 196 197 double d=Length(A.c-B.c); 198 double rdiff=A.r-B.r; 199 double rsum=A.r+B.r; 200 201 if(dcmp(d-rdiff)<0) return 0; // 内含 202 203 double base=angle(B.c-A.c); 204 if(dcmp(d)==0&&dcmp(rdiff)==0) return -1 ; // 重合 无穷多条切线 205 206 if(dcmp(d-rdiff)==0) // 内切 外公切线 207 { 208 a[cnt]=A.point(base); 209 b[cnt]=B.point(base); 210 cnt++; 211 return 1; 212 } 213 214 // 有外公切线的情形 215 216 double ang=acos(rdiff/d); 217 a[cnt]=A.point(base+ang); 218 b[cnt]=B.point(base+ang); 219 cnt++; 220 a[cnt]=A.point(base-ang); 221 b[cnt]=B.point(base-ang); 222 cnt++; 223 224 if(dcmp(d-rsum)==0) // 外切 有内公切线 225 { 226 a[cnt]=A.point(base); 227 b[cnt]=B.point(base+PI); 228 cnt++; 229 } 230 231 else if(dcmp(d-rsum)>0) // 外离 又有两条外公切线 232 { 233 double ang_in=acos(rsum/d); 234 a[cnt]=A.point(base+ang_in); 235 b[cnt]=B.point(base+ang_in+PI); 236 cnt++; 237 a[cnt]=A.point(base-ang_in); 238 b[cnt]=B.point(base-ang_in+PI); 239 cnt++; 240 } 241 242 return cnt; 243 } 244 245 Point Zero=Point(0,0); 246 247 double TriAngleCircleInsection(Circle C, Point A, Point B){ 248 Vector OA = A-C.c, OB = B-C.c; 249 Vector BA = A-B, BC = C.c-B; 250 Vector AB = B-A, AC = C.c-A; 251 double DOA = Length(OA), DOB = Length(OB),DAB = Length(AB), r = C.r; 252 if(dcmp(Cross(OA,OB)) == 0) return 0; 253 if(dcmp(DOA-C.r) < 0 && dcmp(DOB-C.r) < 0) return Cross(OA,OB)*0.5; 254 else if(DOB < r && DOA >= r) { 255 double x = (Dot(BA,BC) + sqrt(r*r*DAB*DAB-Cross(BA,BC)*Cross(BA,BC)))/DAB; 256 double TS = Cross(OA,OB)*0.5; 257 return asin(TS*(1-x/DAB)*2/r/DOA)*r*r*0.5+TS*x/DAB; 258 } 259 else if(DOB >= r && DOA < r) { 260 double y = (Dot(AB,AC)+sqrt(r*r*DAB*DAB-Cross(AB,AC)*Cross(AB,AC)))/DAB; 261 double TS = Cross(OA,OB)*0.5; 262 return asin(TS*(1-y/DAB)*2/r/DOB)*r*r*0.5+TS*y/DAB; 263 } 264 else if(fabs(Cross(OA,OB)) >= r*DAB || Dot(AB,AC) <= 0 || Dot(BA,BC) <= 0) { 265 if(Dot(OA,OB) < 0){ 266 if(Cross(OA,OB) < 0) return (-acos(-1.0)-asin(Cross(OA,OB)/DOA/DOB))*r*r*0.5; 267 else return ( acos(-1.0)-asin(Cross(OA,OB)/DOA/DOB))*r*r*0.5; 268 } 269 else return asin(Cross(OA,OB)/DOA/DOB)*r*r*0.5; 270 } 271 else { 272 double x = (Dot(BA,BC)+sqrt(r*r*DAB*DAB-Cross(BA,BC)*Cross(BA,BC)))/DAB; 273 double y = (Dot(AB,AC)+sqrt(r*r*DAB*DAB-Cross(AB,AC)*Cross(AB,AC)))/DAB; 274 double TS = Cross(OA,OB)*0.5; 275 return 276 (asin(TS*(1-x/DAB)*2/r/DOA)+asin(TS*(1-y/DAB)*2/r/DOB))*r*r*0.5 + TS*((x+y)/DAB-1); 277 } 278 }