Changsha F Curve in Circle

怎么说呢，这道题也是受了适牛的指导，看解题报告看了三四天才能A掉……

转向适牛解题报告

最初的思想是定小圆转大圆，然后通过球长度求出来我们想要的轨迹方程，对于多边形的每一条边，我们求出来每一个端点的边界值，套用辛普森积分法……

然后和适牛说了我的想法，对此适牛直接把这个问题转化成了三角形内部的问题，化为计较坐标，用神一般的三角函数直接予以解决……

code还是仿写适牛……代码功力太渣……TAT，有些东西真的实现不了

（适牛的code真的好清爽……仰慕……）

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>
using namespace std;
const double eps = 1e-8;
int dblcmp(double x){
	if (fabs(x) < eps) return 0;
	return x > eps ? 1 : -1;
}
const double pi = acos(-1.0);
inline double sqr(double x){
	return x * x;
}
struct Point2D{
	double x, y;
	Point2D(){}
	Point2D(double a, double b) : x(a), y(b){}
	void input(){
		scanf("%lf%lf", &x, &y);
	}
	void output(){
		printf("%lf %lf\n", x, y);
	}
	friend Point2D operator + (const Point2D &a, const Point2D &b){
		return Point2D(a.x + b.x, a.y + b.y);
	}
	friend Point2D operator - (const Point2D &a, const Point2D &b){
		return Point2D(a.x - b.x, a.y - b.y);
	}
	friend Point2D operator * (const double a, const Point2D &b){
		return Point2D(a * b.x, a * b.y);
	}
	friend Point2D operator / (const Point2D &a, const double &b){
		return Point2D(a.x / b, a.y / b);
	}
	friend bool operator == (const Point2D &a, const Point2D &b){
		return dblcmp(a.x - b.x) == 0 && dblcmp(a.y - b.y) == 0;
	}
	friend bool operator < (const Point2D &a, const Point2D &b){
		return dblcmp(a.x - b.x) < 0 || dblcmp(a.x - b.x) == 0 && dblcmp(a.y - b.y) < 0;
	}
	double len(){
		return hypot(x, y);
	}
};
double det(const Point2D &a, const Point2D &b){
	return a.x * b.y - a.y * b.x;
}
double dot(const Point2D &a, const Point2D &b){
	return a.x * b.x + a.y * b.y;
}
double dist(const Point2D &a, const Point2D &b){
	return (a - b).len();
}
Point2D rotate_point(const Point2D &p, double A){
	double tx = p.x, ty = p.y;
	return Point2D(tx * cos(A) - ty * sin(A), tx * sin(A) + ty * cos(A));
}
const int maxp = 100 + 10;
struct Polygon2D{
	int n;
	Point2D p[maxp];
	void copy(Polygon2D &A){
		A.n = n;
		for (int i = 0; i < n; i++){
			A.p[i] = p[i];
		}
	}
	void input(){
		for (int i = 0; i < n; i++)
			p[i].input();
	}
	struct cmp{
		Point2D p;
		cmp(const Point2D &p0){p = p0;}
		bool operator ()(const Point2D &aa, const Point2D &bb){
			Point2D a = aa, b = bb;
			int d = dblcmp(det(a - p, b - p));
			if (d == 0)
				return dblcmp(dist(a, p) - dist(b, p)) < 0;
			return d > 0;
		}
	};
	void norm(){
		Point2D mi = p[0];
		for (int i = 0; i < n; i++)
			mi = min(p[i], mi);
		sort(p, p + n, cmp(mi));
	}
	void getconvex(Polygon2D &convex){
		int i, j, k;
		sort(p, p + n);
		convex.n = n;
		for (i = 0; i < min(2, n); i++)
			convex.p[i] = p[i];
		if (n <= 2) return;
		int &top = convex.n;
		top = 1;
		for (int i = 2; i < n; i++){
			while(top && dblcmp(det(convex.p[top] - p[i], convex.p[top - 1] - p[i])) <= 0)
				top -= 1;
			convex.p[++top] = p[i];
		}
		int tmp = top;
		convex.p[++top] = p[n - 2];
		for (int i = n - 3; i >= 0; i--){
			while(top != tmp && dblcmp(det(convex.p[top] - p[i], convex.p[top - 1] - p[i])) <= 0)
				top -= 1;
			convex.p[++top] = p[i];
		}
	}
};
double R, r, A, B, C, k, la;
double T(double x){
	double fm = sin(k * x) / sin(C) + sin(A - k * x) / sin(B);
	double ret = la / fm + R - r;
	return ret;
}
double dT(double x){
	 double fz = la * sin(B) * sin(C) * (-k * sin(C) * cos(A - k * x) + k * sin(B) * cos(k * x));
	 double fm = sqr(sin(C) * sin(A - k * x) + sin(B) * sin(k * x));
	 return -fz / fm;
}
double f0(double x){
	double dx = dT(x) * cos(x) - T(x) * sin(x);
	double dy = dT(x) * sin(x) + T(x) * cos(x);
	return sqrt(sqr(dx) + sqr(dy));
}
double adaptive_simpspons_aux(double (*f)(double), double a, double b, double eps, double s, 
	double fa, double fb, double fc, int depth){
	double c = (a + b) / 2, h = b - a;
	double d = (a + c) / 2, e = (b + c) / 2;
	double fd = f(d), fe = f(e);
	double sl = (fa + 4 * fd + fc) * h / 12;
	double sr = (fc + 4 * fe + fb) * h / 12;
	double s2 = sl + sr;
	if (depth <= 0 || fabs(s2 - s) <= 15 * eps)
		return s + (s2 - s) / 15;
	return adaptive_simpspons_aux(f, a, c, eps / 2, sl, fa, fc, fd, depth - 1) + 
		   adaptive_simpspons_aux(f, c, b, eps / 2, sr, fc, fb, fe, depth - 1);
}
double adaptive_simpspons(double (*f)(double), double a, double b, double eps, int depth){
	double c = (a + b) / 2, h = b - a;
	double fa = f(a), fb = f(b), fc = f(c);
	double s = (fa + 4 * fc + fb) * h / 6;
	return adaptive_simpspons_aux(f, a, b, eps, s, fa, fb, fc, depth);
}
Polygon2D p, q;
void solve(){
	k = R / r;
	p.input();
	p.getconvex(q);
	double ans = 0;
	q.norm();
	for (int i = 0; i < q.n; i++){
		int j = (i + 1) % q.n;
		double t = atan2(q.p[j].y, q.p[j].x)- atan2(q.p[i].y , q.p[i].x); 
		if (dblcmp(t) < 0) t += pi * 2.0;
		la = (q.p[i] - q.p[j]).len();
		double lb = q.p[j].len(), lc = q.p[i].len();
		A = (sqr(lb) + sqr(lc) - sqr(la)) / 2.0 / lb / lc;
		A = acos(A);
		B = (sqr(la) + sqr(lc) - sqr(lb)) / 2.0 / la / lc;
		B = acos(B);
		C = (sqr(la) + sqr(lb) - sqr(lc)) / 2.0 / la / lb;
		C = acos(C);
		ans += adaptive_simpspons(f0, 0, t / k, eps, 20);
	}
	printf("%.3lf\n", ans + eps);
}
int main(){
	#ifndef ONLINE_JUDGE
	freopen("in", "rt", stdin);
	#endif
	while(~scanf("%lf%lf%d", &R, &r, &p.n))
		solve();
	return 0;
}