HDU4752 Polygon

点击打开链接

题目大意是给定一个多边形，一条抛物线，求抛物线在多边形里面部分的长度的总和。

基本思想：对抛物线和多边形求交点，求出来交点之后判断那几段在多边形内，然后对那几段进行曲线积分，然后加和即可。

我的最初思想是求出来所有的交点之后，然后判断两个交点的中点是不是在多边形内，然后确定区间之后直接曲线积分，但是如果这么做的话，如果多边形的边在抛物线上面交错相交，时间复杂度为O(nn)，稳稳地不可以，于是乎这道题被卡住了……

然后几何大神适牛大大放出了这道题的解题报告，我便去拜读了一下，发现不止姿势优美，而且非常好写，尤其是判断抛物线恰好落在多边形顶点上面的时候，那个交点是否选择，他的判断方法真的是神思想……详情参见适牛大大的解题报告

至于code……几乎是仿写了适牛的写法……TAT

（同时学会了自适应辛普森积分法……）

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <vector>
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("%.2lf %.2lf", 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 bool operator == (const Point2D &a, const Point2D &b){
		return dblcmp(a.x - b.x) == 0 && dblcmp(a.y - b.y) == 0;
	}
	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);
	}
	double norm(){
		return hypot(x, y);
	}
	Point2D trunc(double r){
		double l = norm();
		if (!dblcmp(l)) return *this;
		r /= l;
		return Point2D(x * r, y * r);
	}
};
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).norm();
}
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));
}
double a, b, c, l, r;
double f0(double x){
	return sqrt(1 + sqr(2 * a * x + b));
}
double adaptive_simpsons_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 = (c + b) / 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 s2 + (s2 - s) / 15;
	return adaptive_simpsons_aux(f, a, c, eps / 2, sl, fa, fc, fd, depth - 1) + 
		   adaptive_simpsons_aux(f, c, b, eps / 2, sr, fc, fb, fe, depth - 1);
}
double adaptive_simpsons(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_simpsons_aux(f, a, b, eps, s, fa, fb, fc, depth);
}
double f(double x){
	return a * x * x + b * x + c;
}
int n;
const int maxn = 2e4 + 10;
Point2D p[maxn];
vector<double> h;
bool inSeg(Point2D l, Point2D r, Point2D mid){
	return dblcmp(mid.x - fmin(l.x, r.x)) >= 0
		&& dblcmp(fmax(l.x, r.x) - mid.x) >= 0
		&& dblcmp(mid.y - fmin(l.y, r.y)) >= 0
		&& dblcmp(fmax(l.y, r.y) - mid.y) >= 0;
}
bool check(Point2D pre, Point2D now, Point2D nxt){
	pre = pre - now;
	nxt = nxt - now;
	pre = pre.trunc(0.1);
	nxt = nxt.trunc(0.1);
	pre = pre + now;
	nxt = nxt + now;
	if (dblcmp(pre.y - f(pre.x)) * dblcmp(nxt.y - f(nxt.x)) < 0) return true;
	return false;
}
void solve(){
	scanf("%lf%lf%lf%lf%lf", &a, &b, &c, &l, &r);
	for (int i = 0; i < n; i++)
		p[i].input();
	h.resize(0);
	if (n <= 2){
		puts("0.00");
		return;
	}
	for (int i = 0; i < n; ++i){
		Point2D now, pre, nxt, Inter, Inter2;
		pre = p[(i + n - 1) % n];
		now = p[i];
		nxt = p[(i + 1) % n];
		if (dblcmp(now.x - nxt.x) == 0){
			Inter = Point2D(now.x, f(now.x));
			if (inSeg(now, nxt, Inter)){
				h.push_back(Inter.x);
				if (Inter == now)
					if (!check(pre, now, nxt))
						h.push_back(Inter.x);
			}
		}else{
			double kk = (nxt.y - now.y) / (nxt.x - now.x);
			double bb = nxt.y - nxt.x * kk;
			double A = sqr(b - kk) - 4.0 * (c - bb) * a;
			if (dblcmp(A) > 0){
				double x1 = (-(b - kk) + sqrt(A)) / 2.0 / a;
				double x2 = (-(b - kk) - sqrt(A)) / 2.0 / a;
				Inter = Point2D(x1, f(x1));
				Inter2 = Point2D(x2, f(x2));
				if (inSeg(now, nxt, Inter)){
					if (!(Inter == pre) && !(Inter == nxt)) h.push_back(Inter.x);
					if (Inter == now && !check(pre, now, nxt)) h.push_back(Inter.x);
				}
				if (inSeg(now, nxt, Inter2)){
					if (!(Inter2 == pre) && !(Inter2 == nxt)) h.push_back(Inter2.x);
					if (Inter2 == now && !check(pre, now, nxt)) h.push_back(Inter2.x);
				}
			}
		}
	}
	sort(h.begin(), h.end());
	double ans = 0;
	for (int i = 0; i + 1 < h.size(); i += 2){
		double L = fmax(l, h[i]), R = fmin(r, h[i + 1]);
		if (dblcmp(R - L) > 0)
			ans += adaptive_simpsons(f0, L, R, eps, 10);
		if (dblcmp(R - L) < 0) break;
	}
	double out = ans;
	printf("%.2lf\n", out + eps);
}
int main(){
	#ifndef ONLINE_JUDGE
	freopen("in.txt", "rt", stdin);
	#endif
	while(scanf("%d", &n) == 1)
		solve();
	return 0;
}