POJ 3675 三角剖分

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
using namespace std;
struct Point
{
	double x, y;
	Point(double x = 0, double y = 0): x(x), y(y) {}
};
typedef Point Vector;
typedef vector<Point> Polygon;
Point read_point()
{
	double X, Y;
	scanf("%lf%lf", &X, &Y);
	return Point(X, Y);
}
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));
}
double Area2(Point A, Point B, Point C) //
{
	return Cross(B - A, C - A);
}
double TriArea(Point A, Point B, Point C) //
{
	return fabs(Cross(B - A, C - A)) / 2;
}
double Distance(Point A, Point B)
{
	return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
double Distance2(Point A, Point B)
{
	return (A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y);
}
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));
}
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) //
{
	Vector u = P - Q;
	double t = Cross(w, u) / Cross(v, w);
	return P + v * t;
}
double DistanceToLine(Point p, Point A, Point B) //
{
	Vector v1 = B - A, v2 = p - A;
	return fabs(Cross(v1, v2)) / Length(v1);
}
double DistanceToSegment(Point P, Point A, Point B) //
{
	if (A == B) return Length(P - A);
	Vector v1 = B - A, v2 = P - A, v3 = P - B;
	if (dcmp(Dot(v1, v2)) < 0) return Length(v2);
	else if (dcmp(Dot(v1, v3) > 0)) return Length(v3);
	else return fabs(Cross(v1, v2)) / Length(v1);
}
Point GetLineProjection(Point P, Point A, Point B) //
{
	Vector v = B - A;
	return A + v * (Dot(v, P - A) / Dot(v, v));
}
bool OnSegment(Point p, Point a1, Point a2) //
{
	return dcmp(Cross(a1 - p, a2 - p)) == 0 && dcmp(Dot(a1 - p, a2 - p)) < 0;
}
double PolygonArea(Point *p, int n) //
{
	double area = 0;
	for (int i = 1; i < n - 1; i++)
		area += Cross(p[i] - p[0], p[i + 1] - p[0]);
	return area / 2;
}
vector<Point> ConvexHull(vector<Point>& p)
{
	sort(p.begin(), p.end());
	p.erase(unique(p.begin(), p.end()), p.end());

	int n = p.size();
	int m = 0;
	vector<Point> ch(n + 1);
	for (int i = 0; i < n; i++)
	{
		while (m > 1 && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0) m--;
		ch[m++] = p[i];
	}
	int k = m;
	for (int i = n - 2; i >= 0; i--)
	{
		while (m > k && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0) m--;
		ch[m++] = p[i];
	}
	if (n > 1) m--;
	ch.resize(m);
	return ch;
}
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;
	}
};
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, 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 common_area(Circle C, Point A, Point B)
{
	if (A == B) return 0;
	if (A == C.c || B == C.c) return 0;
	double OA = Length(A - C.c), OB = Length(B - C.c);
	double d = DistanceToLine(Point(0, 0), A, B);
	int sg = dcmp(Cross(A, B));
	if (sg == 0) return 0;
	double angle = Angle(A, B);
	if (dcmp(OA - C.r) <= 0 && dcmp(OB - C.r) <= 0)
		return Cross(A, B) / 2;
	else if (dcmp(OA - C.r) >= 0 && dcmp(OB - C.r) >= 0 && dcmp(d - C.r) >= 0)
		return sg * C.r * C.r * angle / 2;
	else if (dcmp(OA - C.r) >= 0 && dcmp(OB - C.r) >= 0 && dcmp(d - C.r) < 0)
	{
		Point prj = GetLineProjection(Point(0, 0), A, B);
		if (OnSegment(prj, A, B))
		{
			vector<Point> p;
			Line L = Line(A, B - A);
			double t1, t2;
			GetLineCircleIntersection(L, C, t1, t2, p);
			double s1 = C.r * C.r * angle / 2;
			double s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
			s2 -= fabs(Cross(p[0], p[1]) / 2);
			s1 -= s2;
			return sg * s1;
		}
		else
			return sg * C.r * C.r * angle / 2;
	}
	else
	{
		if (dcmp(OB - C.r) < 0)
		{
			Point temp = A;
			A = B;
			B = temp;
		}
		Point inter_point;
		double t1, t2;
		Line L = Line(A, B - A);
		vector<Point> inter;
		GetLineCircleIntersection(L, C, t1, t2, inter);
		if (OnSegment(inter[0], A, B))
			inter_point = inter[0];
		else
			inter_point = inter[1];
		double s = fabs(Cross(inter_point, A) / 2);
		s += C.r * C.r * Angle(inter_point, B) / 2;
		return s * sg;
	}
}
const int maxn = 55;
double R;
int n;
Point P[maxn];

int main(int argc, char const *argv[])
{
	while (~scanf("%lf", &R) && R)
	{
		Circle C(Point(0, 0), R);
		scanf("%d", &n);
		for (int i = 0; i < n; i++)
			P[i] = read_point();
		P[n] = P[0];
		double ans = 0;
		for (int i = 0; i < n; i++)
			ans += common_area(C, P[i], P[i + 1]);
		printf("%.2lf\n", fabs(ans) + eps);
	}
	return 0;
}

这其实不能算三角剖分，只是把多边形分解称为三角形，求有向面积就可以了。最后不要忘记取正值。另外模板要好好保存。