读书人

[计算几何][SHOI2008]保险的航线flight

发布时间: 2012-12-20 09:53:21 作者: rapoo

[计算几何][SHOI2008]安全的航线flight

#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <algorithm>using namespace std;#define eps 0.000000001#define oo 999999999999.struct point {double x, y;};struct seg {point x, y;};struct line {point f; double c;};int sig(double a) {return a < - eps ? - 1 : a > eps;}double len(point a) {return sqrt(a.x * a.x + a.y * a.y);}double operator *(point a, point b) {return a.x * b.y - a.y * b.x;}double operator %(point a, point b) {return a.x * b.x + a.y * b.y;}point operator +(point a, point b) {return (point) {a.x + b.x, a.y + b.y};}point operator -(point a, point b) {return (point) {a.x - b.x, a.y - b.y};}point operator *(line a, line b){  double c = a.f * b.f;return (point) {(b.f.y * a.c - a.f.y * b.c) / c,     (a.f.x * b.c - b.f.x * a.c) / c};}point make_point(point a, double dist){  double x = a.x * a.x;  double y = a.y * a.y, z = x + y;  return (point) {sqrt(x * dist * dist / z) * sig(a.x),    sqrt(y * dist * dist / z) * sig(a.y)};}line make_line(seg a){  point f = (point) {a.y.y - a.x.y, a.x.x - a.y.x};  return (line) {f, f.x * a.x.x + f.y * a.x.y};}point select(seg a, point b){if (sig(a.x.x - a.y.x)){if ((a.x.x - b.x) * (a.y.x - b.x) < eps) return b;if (a.x.x > a.y.x) swap(a.x, a.y);return b.x < a.x.x ? a.x : a.y;}else{ if ((a.x.y - b.y) * (a.y.y - b.y) < eps) return b;if (a.x.y > a.y.y) swap(a.x, a.y);return b.y < a.x.y ? a.x : a.y;}}point operator *(seg a, line b){  point p = make_line(a) * b;  if (sig(a.x.x - a.y.x) && (a.x.x - p.x) * (a.y.x - p.x) > - eps) return (point) {oo};  if (sig(a.x.y - a.y.y) && (a.x.y - p.y) * (a.y.y - p.y) > - eps) return (point) {oo};  return p;}point cross_circle(line a, seg b){line l = make_line(b);line r = (line) {(point) {l.f.y, - l.f.x}, a.f.x * l.f.y - a.f.y * l.f.x};point p = l * r, q;double dist = len(p - a.f), A, B;if (dist - a.c > - eps) return (point) {oo};dist = sqrt(a.c * a.c - dist * dist);q = make_point(b.x - b.y, dist);A = len(select(b, p + q) - b.x), B = len(select(b, p - q) - b.x);return A < B ? (point) {A, B} : (point) {B, A};}void cross_poly(int k, point poly[], seg b, int & tot, point u[]){int i, j; double a[31];seg s; point p; line l = make_line(b); j = 0;for (i = 1; i <= k; ++ i){s = (seg) {poly[i - 1], poly[i]}, p = s * l;if (p.x != oo) a[++ j] = len(select(b, p) - b.x);if (! sig((b.x - poly[i]) * (b.y - poly[i])))a[++ j] = len(select(b, poly[i]) - b.x);}sort(a + 1, a + j + 1);for (i = 1; i <= j; i += 2)u[++ tot] = (point) {a[i], a[i + 1]};}int C, n;int num[35], total[35];point fli[35];point lan[35][35], cover[35][1000];bool cmp(point A, point B) {return A.x < B.x;}bool okay(double dist){int i, j, k, tot; double le;point p; point cov[1000], poly[6]; seg s;for (i = 2; i <= n; ++ i){tot = total[i], memcpy(cov, cover[i], sizeof cover[i]);for (j = 1; j <= C; ++ j){s = (seg) {fli[i - 1], fli[i]};for (k = 1; k <= num[j]; ++ k){p = cross_circle((line) {lan[j][k], dist}, s);if (p.x != oo) cov[++ tot] = p;p = lan[j][k] - lan[j][k - 1];p = make_point((point) {p.y, - p.x}, dist);poly[1] = lan[j][k - 1], poly[2] = lan[j][k - 1] + p;poly[3] = lan[j][k] + p, poly[0] = poly[4] = lan[j][k];cross_poly(4, poly, s, tot, cov);}}sort(cov + 1, cov + tot + 1, cmp);p = cov[1], le = 0;for (j = 2; j <= tot; ++ j)if (cov[j].x > p.y) le += p.y - p.x, p = cov[j];else if (cov[j].y > p.y) p.y = cov[j].y;le += p.y - p.x;if (len(s.y - s.x) - le > eps) return 0;}return 1;}int main(){freopen("1020.in", "r", stdin);freopen("1020.out", "w", stdout);int i, j; double l, r, s;scanf("%d%d", & C, & n);for (i = 1; i <= n; ++ i)scanf("%lf%lf", & fli[i].x, & fli[i].y);for (i = 1; i <= C; ++ i){scanf("%d", num + i);for (j = 1; j <= num[i]; ++ j)scanf("%lf%lf", & lan[i][j].x, & lan[i][j].y);s = lan[i][num[i]] * lan[i][1];for (j = 2; j <= num[i]; ++ j)s += lan[i][j - 1] * lan[i][j];if (s < - eps)for (j = 1; j << 1 <= num[i]; ++ j)swap(lan[i][j], lan[i][num[i] - j + 1]);lan[i][0] = lan[i][num[i]];for (j = 2; j <= n; ++ j)cross_poly(num[i], lan[i], (seg) {fli[j - 1], fli[j]}, total[j], cover[j]);}for (l = 0, r = 10000; r - l > 0.001; )okay(s = (l + r) / 2) ? r = s : l = s;printf("%.2lf", s);return 0;}


1楼emoizhang前天 22:50
桂兰Jerry盾
Re: JerryDung昨天 15:17
回复emoizhangnOrz

读书人网 >其他相关

热点推荐