RUPC Day2 I - One - Lilliput Steps

問題文 : One

解法 :

放物線の長さは,

$\displaystyle \int_{\alpha}^{\beta} \sqrt{1+4a^2(x-p)^2} dx = [ \displaystyle \frac{1}{4a}(\log |u^2 + \sqrt{u^2 + 1} | + u \sqrt{u^2 + 1}) ]_{\alpha}^{\beta}$

で求まる. (ここで, $u = 2a(x-p)$ とおいた.)

交点を列挙してこの関数を適用していけば良い. 誤差に気をつけて計算すべし...

コード :

#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>

using namespace std;

const long double EPS = 1e-9;

bool eq(long double a, long double b)
{
	return (fabs(b - a) <= EPS);
}

struct two {
	long double a, p, q;
	two(long double a, long double p, long double q) : a(a), p(p), q(q) {}
};

struct Point {
	long double x, y;
	Point(long double x, long double y) : x(x), y(y) {}
};

bool operator < (const Point &a, const Point &b)
{
	return (a.x < b.x);
}

vector<Point> get(two u, two v)
{
	vector<Point> ret;
	long double a = u.a - v.a,
			b = -2 * u.a * u.p - (-2 * v.a * v.p),
			c = (u.a * u.p * u.p + u.q) - (v.a * v.p * v.p + v.q);
	
	if (eq(a, 0)){
		if (!eq(b, 0)){
			long double x = -c / b;
			Point p(x, v.a * (x - v.p) * (x - v.p) + v.q);
			ret.push_back(p);
		}
	}
	else {
		long double D = b * b - 4 * a * c;
		
		if (D + EPS >= 0){
			if (eq(D, 0)){
				long double x = -b / (2 * a);
				ret.push_back(Point(x, v.a * (x - v.p) * (x - v.p) + v.q));
			}
			else {
				long double x1 = (-b + sqrt(D)) / (2 * a),
							 x2 = (-b - sqrt(D)) / (2 * a);
				ret.push_back(Point(x1, v.a * (x1 - v.p) * (x1 - v.p) + v.q));
				ret.push_back(Point(x2, v.a * (x2 - v.p) * (x2 - v.p) + v.q));
			}
		}
	}
	
	return (ret);	
}

long double integral(two t, long double x)
{
	long double u = 2 * t.a * (x - t.p);
	return ((log(fabs(u + sqrt(1 + u * u))) + u * sqrt(1 + u * u)) / (4 * t.a));
}

int main()
{
	int W, H, N;
	vector<two> v;
	vector<Point> pt;
	
	scanf("%d %d %d", &W, &H, &N);
	
	for (int i = 0; i < N; i++){
		int a, p, q;
		scanf("%d %d %d", &a, &p, &q);
		v.push_back(two(a, p, q));
	}
	
	for (int i = 0; i < N; i++){
		for (int j = i + 1; j < N; j++){
			vector<Point> a = get(v[i], v[j]);
			for (int s = 0; s < a.size(); s++){
				pt.push_back(a[s]);
			}
		}
		pt.push_back(Point(0, v[i].a * v[i].p * v[i].p + v[i].q));
		pt.push_back(Point(W, v[i].a * (W - v[i].p) * (W - v[i].p) + v[i].q));
		vector<Point> a = get(v[i], two(0, 0, 0));
		for (int s = 0; s < a.size(); s++)
			pt.push_back(a[s]);
	}
	
	sort(pt.begin(), pt.end());
	long double len = 0.0;
	for (int i = 0; i < (int)pt.size() - 1; i++){
		if (0 <= pt[i].x + EPS && pt[i + 1].x - EPS <= W && !eq(pt[i].x, pt[i + 1].x)){
			int cand = -1;
			long double val = -1e20;
			long double piv = (pt[i].x + pt[i + 1].x) / 2.0;
			for (int j = 0; j < N; j++){
				if (v[j].a * (piv - v[j].p) * (piv - v[j].p) + v[j].q > val){
					cand = j;
					val = v[j].a * (piv - v[j].p) * (piv - v[j].p) + v[j].q;
				}
			}
			if (0 <= val + EPS){
				//printf("%Lf %Lf %d\n", pt[i].x, pt[i + 1].x, cand);
				//printf("%Lf %Lf\n", integral(v[cand], pt[i + 1].x), integral(v[cand], pt[i].x));
				len += fabs(integral(v[cand], pt[i + 1].x) - integral(v[cand], pt[i].x));
			}
		}
	}
	
	printf("%.15Lf\n", len);
	
	return (0);
}