两个圆的切线与切点

两个圆的切线与切点
https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_G

int sgn(double x) {
	if(fabs(x)<eps)return 0;
	return x<0?-1:1;
}
struct Point { //定义点和基本运算
	double x,y;
	double ang;
	Point() {}
	Point(double x,double y):x(x),y(y) {}
	Point operator + (Point B) {
		return Point(x+B.x,y+B.y);
	}
	Point operator - (Point B) {
		return Point(x-B.x,y-B.y);
	}
	Point operator * (double k) {
		return Point(x*k,y*k);   //长度增大k倍
	}
	Point operator / (double k) {
		return Point(x/k,y/k);   //长度缩小k倍
	}
	bool operator == (Point B) {
		return sgn(x-B.x)==0 && sgn(y-B.y)==0;
	}
	double operator ^(Point B) {
		return x*B.y-y*B.x;
	}
	double distance(Point p) {
		return hypot(x-p.x,y-p.y);
	}
};

typedef Point Vector;
double Cross(Vector A,Vector B) {
	return A.x*B.y - A.y*B.x;   //叉积
}

struct Line {
	Point p1,p2;//线上的两个点
	Line() {}
	Line(Point p1,Point p2):p1(p1),p2(p2) {}
};

struct Circle {
	Point c;//圆心
	double r;//半径
	Circle() {}
	Circle(Point c,double r):c(c),r(r) {}
	Circle(double x,double y,double _r) {
		c=Point(x,y);
		r = _r;
	}
	Point point(double ang) { //圆上与圆心极坐标为ang的点-----------------
		return Point(c.x+cos(ang)*r,c.y+sin(ang)*r);
	}
};
//0,-1:没有切线 a[]是c1上的切点,b[]是c2
int getTangents(Circle A, Circle B, Point *a, Point *b) {
	int cnt = 0;        //存切点用
	if(sgn(A.r - B.r) < 0) {
		swap(A, B);
		swap(a, b);
	}
	double d = sqrt((A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y));     //圆心距
	double rdiff = A.r - B.r;      //两圆半径差
	double rsum = A.r + B.r;       //两圆半径和
	if(sgn(d - rdiff) < 0) return 0;        //1.内含
	double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);      //向量AB的极角
	if(sgn(d) == 0) return -1;        //2.重合
	if(sgn(d - rdiff) == 0) {      //3.内切
		a[cnt] = b[cnt] = A.point(base);
		cnt++;
		return 1;
	}
	double ang = acos((A.r - B.r) / d);
	a[cnt] = A.point(base + ang);
	b[cnt] = B.point(base + ang);
	cnt++;      //4.相交（外切、外离的外公切线也在此求出）
	a[cnt] = A.point(base - ang);
	b[cnt] = B.point(base - ang);
	cnt++;      //两条外公切线的切点
	if(sgn(d - rsum) == 0) {       //5.外切
		a[cnt] = b[cnt] = A.point(base);
		cnt++;
	} else if(sgn(d - rsum) > 0) {   //6.外离
		double ang = acos((A.r + B.r) / d);
		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;

}

bool cmp(Point a,Point b) {
	if(sgn(a.x-b.x)!=0)return a.x < b.x;
	else if(sgn(a.x-b.x)==0)return a.y < b.y;
}

void work() {
	Circle a,b;
	scanf("%lf%lf%lf",&a.c.x,&a.c.y,&a.r);
	scanf("%lf%lf%lf",&b.c.x,&b.c.y,&b.r);
	Point p1[5],p2[5];
	int cnt = getTangents(a,b,p1,p2);
//	cout << cnt << endl;
	if(cnt == -1 ||cnt==0)return ;
	vector<Point>ans;
	for(int i=0; i<cnt; i++) {
		ans.push_back(p1[i]);
	}
	sort(ans.begin(),ans.end(),cmp);
	for(int i=0; i<cnt; i++) {
		printf("%.10f %.10f
",ans[i].x,ans[i].y);
	}

}