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_JUDGEfreopen("in", "rt", stdin);#endifwhile(~scanf("%lf%lf%d", &R, &r, &p.n))solve();return 0;}