POJ 1263 简单几何？！

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <vector>
using namespace std;
const double DINF = 1E100;
struct Point
{
	double x, y;
	Point(double x = 0, double y = 0): x(x), y(y) {}
};
typedef Point Vector;
typedef vector<Point> Polygon;
Vector operator +(Vector A, Vector B)//
{
	return Vector(A.x + B.x, A.y + B.y);
}
Vector operator -(Point A, Point 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 Cross(Vector A, Vector B)//
{
	return A.x * B.y - A.y * B.x;
}
double Angle(Vector A, Vector B)//
{
	return acos(Dot(A, B) / Length(A) / Length(B));
}
Vector Normal(Vector A)//
{
	double L = Length(A);
	return Vector(-A.y / L, A.x / L);
}
Vector Rotate(Vector A, double rad) //
{
	return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}
struct Line
{
	Point p;
	Vector v;
	double ang;
	Line() {};
	Line(Point P, Vector V)
	{
		p = P; v = V;
		ang = atan2(v.y, v.x);
	}
	bool operator < (const Line& L) const
	{
		return ang < L.ang;
	}
	Point point(double t)
	{
		return p + v * t;
	}
	Line move(double d)
	{
		return Line(p + Normal(v) * d, v);
	}
};
bool OnLeft(Line L, Point K)
{
	return Cross(L.v, K - L.p) > 0;
}
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);
	}
};
int GetLineCircleIntersection(Line L, Circle C, int num, vector<pair<Point, int> >& sol)
{
	double t1, t2;
	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(make_pair(L.point(t1), num));
		return 1;
	}
	// 相交
	t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(make_pair(L.point(t1), num));
	t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(make_pair(L.point(t2), num));
	return 2;
}
Point Reflect(Line K, Point X)
{
	double A = K.v.y;
	double B = -K.v.x;
	double C = K.p.y * K.v.x - K.p.x * K.v.y;
	Point Q;
	Q.x = ((B * B - A * A) * X.x - 2 * A * B * X.y - 2 * A * C) / (A * A + B * B);
	Q.y = ((A * A - B * B) * X.y - 2 * A * B * X.x - 2 * B * C) / (A * A + B * B);
	return Q;
}
std::vector<Circle> res;
Line L;
double X, Y, R, PX, PY;
int n, kase;
int main(int argc, char const *argv[])
{
	while (~scanf("%d", &n) && n)
	{
		res.clear();
		for (int i = 0; i < n; i++)
		{
			scanf("%lf%lf%lf", &X, &Y, &R);
			res.push_back(Circle(Point(X, Y), R));
		}
		scanf("%lf%lf%lf%lf", &X, &Y, &PX, &PY);
		L = Line(Point(X, Y), Point(PX, PY));
		printf("Scene %d\n", ++kase);
		int cnt = 0;
		while (true)
		{
			vector<pair<Point, int> >temP;
			Point initP;
			double initPx = DINF, temPx;
			int CirNum = 0;
			for (int i = 0; i < res.size(); i++)
				GetLineCircleIntersection(L, res[i], i, temP);
			for (int i = 0; i < temP.size(); i++)
			{
				temPx = (temP[i].first.x - L.p.x) / L.v.x;
				if (dcmp(temPx) > 0 && temPx < initPx) initPx = temPx, CirNum = temP[i].second;
			}
			if (dcmp(fabs(initPx) - DINF) >= 0) break;
			cnt++;
			if (cnt > 10) break;
			initP = L.point(initPx);
			printf("%d ", CirNum + 1);
			Point NextP = Reflect(Line(initP, res[CirNum].c - initP), L.p);
			L = Line(initP, NextP - initP);
		}
		if (cnt <= 10) printf("inf\n\n");
		else printf("...\n\n");
	}
	return 0;
}

这题卡了一天了。。。确实是一道思路简洁的几何，一道直线打在球上，发生反射。

要准备好一个反射的模板，

if (dcmp(fabs(initPx) - DINF) >= 0) break;
			cnt++;
			if (cnt > 10) break;

就这玩意儿，第十个打空就应该inf，否则。。。 还以为自己模板错了。