圆的相关计算(刘汝佳版)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;

/*Point的模板*/
struct Point{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
//向量+向量=向量，点+向量=点
Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
//点-点=向量
Vector operator - (Vector A, Vector B) { return Vector(A.x-B.x, A.y-B.y); }
//向量*数=向量
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
//向量/数=向量
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }

bool operator < (const Point& a, const Point& b) {
return a.x < b.x || (a.x == b.x && a.y < b.y);
}

//比较
const double eps = 1e-10;
int dcmp(double x){
if(fabs(x) < eps) return 0;
else return x < 0 ? -1: 1;
}

bool operator == (const Point& a,const Point& b) {
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}

//基本计算
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A,A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B))/Length(A)/Length(B); } //两向量的夹角
double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; } //叉积
double Area2(Point A, Point B, Point C) { return Cross(B-A,C-A); } //有向面积

//A向量逆时针旋转α rad
//x'=xcosα-ysinα;
//y'=xsinα+ycosα;
Vector Rotate(Vector A, double rad){
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

//A的单位法线,也就是逆时针90°，长度变为1，注意A要非零向量
Vector Normal(Vector A){
double L=Length(A);
return Vector(-A.y/L,A.x/L);
}
/*Point的模板*/

/*Circle的模板*/
const double PI=acos(-1);
struct Circle
{
Point c;
double r;
Circle(Point c, double r):c(c), r(r) {}
Point point(double a) {      //根据极角给出圆上的点
return Point(c.x + cos(a)*r, c.y + sin(a)*r);
}
};
struct Line
{
Vector v; //参数方程L=p+v*t
Point p;
Line(Point p,Vector v):v(v), p(p) {}
Point point(double t){     //根据参数方程L=v+(p-v)*t中t求出线上的点
return p+v*t;
}
};
//求直线L与圆C的交点，两个交点参数方程值t1,t2，sol保存交点本身，函数返回交点个数
int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){
double a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
double e = a*a+c*c, f = 2*(a*b+c*d), g = b*b+d*d-C.r*C.r;
double delta = f*f-4*e*g;
if(dcmp(delta) < 0) return 0;
if(dcmp(delta) == 0) {
t1 = t2 = -f / (2*e) ; sol.push_back(L.point(t1));
return 1;
}

t1 = (-f-sqrt(delta)) / (2*e); sol.push_back(L.point(t1));
t2 = (-f+sqrt(delta)) / (2*e); sol.push_back(L.point(t2));
return 2;
}
double Angle(Vector v) { return atan2(v.y,v.x); }   //求v的极角
//求圆与圆的交点，函数返回交点数，sol保存交点
int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol){
double d = Length(C1.c-C2.c);
if(dcmp(d)==0){
if(dcmp(C1.r-C2.r)==0) return -1;  //两圆重合
return 0;
}
if(dcmp(C1.r+C2.r-d)<0) return 0;  //相离
if(dcmp(fabs(C1.r-C2.r)-d)>0) return 0;  //内含

double a=Angle(C2.c-C1.c);     //求极角
double da=acos((C1.r*C1.r+d*d-C2.r*C2.r) / (2*C1.r*d));
Point p1=C1.point(a-da), p2=C1.point(a+da);

sol.push_back(p1);
if(p1==p2) return 1;
sol.push_back(p2);
return 2;
}
//过点p到圆C做切线，函数返回切线条数，v[i]保存第i条切线的向量
int getTangents(Point p, Circle C, Vector* v){
Vector u=C.c-p;
double dist=Length(u);
if(dist<C.r) return 0;
else if(dcmp(dist-C.r)==0) {
v[0]=Rotate(u,PI/2);
return 1;
}
else {
double ang=asin(C.r/dist);
v[0]=Rotate(u, -ang);
v[1]=Rotate(u, +ang);
return 2;
}
}
//求两圆的公切线，返回切线条数，-1表示无穷条，a[i]和b[i]表示第i切线在圆A和圆B上的切点
int getTangents(Circle A, Circle B, Point* a, Point* b){
int cnt=0;
if(A.r<B.r) { swap(A, B); swap(a, b); }  //保证圆A比圆B大
int d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);
int rdiff=A.r-B.r;
int rsum=A.r+B.r;
if(d2<rdiff*rdiff) return 0;  //内含

double base=Angle(B.c-A.c);	// atan2(B.c.y-A.c.y,B.c.x-A.c.x);

if(d2==0&&A.r==B.r) return -1;   //两圆重合，有无数条
if(d2==rdiff*rdiff){            //内切，一条切线
a[cnt]=A.point(base); b[cnt]=B.point(base); cnt++;
return 1;
}
//有外共切线
double ang=acos((A.r-B.r)/sqrt(d2));
a[cnt]=A.point(base+ang); b[cnt]=B.point(base+ang); cnt++;
a[cnt]=A.point(base-ang); b[cnt]=B.point(base-ang); cnt++;
if(d2==rsum*rsum){      //存在一条
a[cnt]=A.point(base); b[cnt]=B.point(PI+base); cnt++;
}
else if(d2>rsum*rsum){  //存在两条
double ang=acos((A.r+B.r)/sqrt(d2));
a[cnt]=A.point(base+ang); b[cnt]=B.point(PI+base+ang); cnt++;
a[cnt]=A.point(base-ang); b[cnt]=B.point(PI+base-ang); cnt++;
}
return cnt;
}
/*Circle的模板*/