POJ 1375 Intervals 解析几何 求圆的切线

题目大意：给出一个点，再给出都处于这个点之下的一些圆，求这个点光源照到这些圆上之后所得到的阴影的并集。

思路：求出每一个圆关于那个点的切线，每一个圆可以处理出来两个切线，这两个切线在x轴上交点的中间部分就是要求的阴影。最后将所有的阴影部分取并输出。

关于求切线，我是利用方向向量解方程做的。应该有更简洁的方法吧。。

CODE：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define MAX 510
#define EPS 1e-8
#define DCMP(a) (fabs(a) < EPS)
using namespace std;

struct Point{
	double x,y;

	Point(double _ = .0,double __ = .0):x(_),y(__) {}
	Point operator +(const Point &a)const {
		return Point(x + a.x,y + a.y);
	}
	Point operator -(const Point &a)const {
		return Point(x - a.x,y - a.y);
	}
	Point operator *(double a)const {
		return Point(x * a,y * a);
	}
	Point operator /(double a)const {
		return Point(x / a,y / a);
	}
	void Read() {
		scanf("%lf%lf",&x,&y);
	}
}start;
struct Circle{
	Point o;
	double r;

	void Read() {
		o.Read();
		scanf("%lf",&r);
	}
}circle;

struct Complex{
	double x,y;

	Complex(double _ = .0,double __ = .0):x(_),y(__) {}
	bool operator <(const Complex &a)const {
		if(x == a.x)	return y < a.y;
		return x < a.x;
	}
}ans[MAX];

int cnt;

inline double Calc(const Point &a,const Point &b)
{
	return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

inline Point Solve(double a1,double b1,double c1,double a2,double b2,double c2)
{
	double y = ((a1 * c2) - (a2 * c1)) / ((a1 * b2) - (a2 * b1));
	double x = (c1 - b1 * y) / a1;
	return Point(x,y);
}

inline Complex Work(const Point &p,const Circle &circle)
{
	double d = Calc(circle.o,p),_d = sqrt(d * d - circle.r * circle.r);

	Point v = Solve(circle.r,-_d,p.x - circle.o.x,_d,circle.r,p.y - circle.o.y);
	Point tangency = circle.o + v * circle.r;
	Point u = tangency - p;
	Point center = p + ((u / (p.y - tangency.y)) * p.y);
	double first = center.x;

	v = Solve(circle.r,_d,p.x - circle.o.x,-_d,circle.r,p.y - circle.o.y);
	tangency = circle.o + v * circle.r;
	u = tangency - p;
	center = p + ((u / (p.y - tangency.y)) * p.y);
	double second = center.x;

	return Complex(second,first);
}

int main()
{
	while(scanf("%d",&cnt),cnt) {
		start.Read();
		for(int i = 1; i <= cnt; ++i) {
			circle.Read();
			ans[i] = Work(start,circle);
		}
		sort(ans + 1,ans + cnt + 1);
		double l = ans[1].x,r = ans[1].y;
		for(int i = 2; i <= cnt; ++i) {
			if(ans[i].x > r) {
				printf("%.2f %.2f\n",l,r);
				l = ans[i].x,r = ans[i].y;
			}
			else	r = max(r,ans[i].y);
		}
		printf("%.2f %.2f\n\n",l,r);
	}
	return 0;
}