1020: [SHOI2008]安全的航线flight

Description

在设计航线的时候，安全是一个很重要的问题。首先，最重要的是应采取一切措施确保飞行不会发生任何事故，但同时也需要做好最坏的打算，一旦事故发生，就要确保乘客有尽量高的生还几率。当飞机迫降到海上的时候，最近的陆地就是一个关键的因素。航线中最危险的地方就是距离最近的陆地最远的地方，我们称这种点为这条航线“孤地点”。孤地点到最近陆地的距离被称为“孤地距离”。作为航空公司的高级顾问，你接受的第一个任务就是尽量找出一条航线的孤地点，并计算这条航线的孤地距离。为了简化问题，我们认为地图是一个二维平面，陆地可以用多边形近似，飞行线路为一条折线。航线的起点和终点都在陆地上，但中间的转折点是可能在海上（如下图所示，方格标示出了孤地点）。
Input

输入的第一行包括两个整数C和N（1≤C≤20，2≤N≤20），分别代表陆地的数目的航线的转折点的数目。接下来有N行，每行有两个整数x,y。(x,y)表示一个航线转折点的坐标，第一个转折点为航线的起点，最后一个转折点为航线的终点。接下来的输入将用来描述C块大陆。每块输入由一个正整数M开始（M≤30），M表示多边形的顶点个数，接下来的M行，每行会包含两个整数x,y，(x,y)表示多边形的一个顶点坐标，我们保证这些顶点以顺时针或逆时针给出了该多边形的闭包，不会出现某些边相交的情况。此外我们也保证输入数据中任何两块大陆不会相交。输入的所有坐标将保证在-10000到10000的范围之间。
Output

输出一个浮点数，表示航线的孤地距离，数据保留2位小数。
Sample Input
1 2
-9 -6
5 1
3
0 16
-16 -12
17 -6
Sample Output
0.00

原来的解法是二分答案，然后把陆地扩展，再判断是否覆盖了航线，但是太繁琐，我根本写不出来

所以我用的是莫涛的那种解法

没想到我竟然错在求垂足上，囧..........

但是答案竟然正确了，时间变长了好多，查了好久才查出来

1.初始化孤地点可能位于的线段集合为整条航线。
2.对于长L的某条线段，左端点与陆地的最近点为P1，右端点与陆地的最近点为P2，那么该线段上的孤地距离将受P1与P2影响。具体来说，利用二分求出该线段上的点P使得dis(p1,p)=dis(p2,p)
令r=dis(p1,p),若r小于已有的最优答案，那么可以删除该线段。（当然，这个只要r-0.001<ans就行了，只要精度达到了就行，不然就是死循环了）
4.取所有线段的中点更新答案。
5.将所有线段从中点分成左右两条线段。
6.不断进行2，3，4直到线段的集合为空。

这个做法效率很高，（哇，好开心，在BZOJ上是pascal第一）

我们现在要讨论的就是这个算法的正确性

现在我们有一条线段和对应的p1和p2，分别是左端点最近的点和右端点最近的点

有三种情况

然后我们发现线段上的点到自己最近点的距离不会超过max(dis(p1,p),dis(a,p1),(b,p2))（a，b分别为线段的左右端点）

所以我们的删除操作是对的，我们删除的都是不会更新答案的线段，就像ydc说的一样，这个就像是搜索剪枝

时间效率也很好，只是细节要注意，不要像我一样傻×，连求垂足都求错

/**************************************************************
    Problem: 1020
    User: 1997cb
    Language: Pascal
    Result: Accepted
    Time:96 ms
    Memory:6488 kb
****************************************************************/
 
const
    maxn=22;
    maxm=33;
    maxq=100000;
    eps=1e-16;
type
    point=record
      x,y:double;
    end;
    polygon=record
      tot:longint;
      a:array[0..maxm]of point;
    end;
    seg=record
      a,b,neara,nearb:point;
    end;
var
    c,head,tail:longint;
    ans:double;
    s:array[0..maxq]of seg;
    p:array[0..maxn]of polygon;
  
function max(x,y:double):double;
begin
    if x>y then exit(x);
    exit(y);
end;
  
function cj(x1,y1,x2,y2:double):double;
begin
    exit(x1*y2-y1*x2);
end;
  
function on(var a,b,c:point):boolean;
begin
    exit((abs(cj(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y))<=eps) and ((a.x-c.x)*(b.x-c.x)<=eps) and ((a.y-c.y)*(b.y-c.y)<=eps));
end;
  
function jiao(var a,b,c,d:point):boolean;
begin
    exit((cj(b.x-a.x,b.y-a.y,d.x-a.x,d.y-a.y)*cj(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y)<=eps) and (cj(d.x-c.x,d.y-c.y,a.x-c.x,a.y-c.y)*cj(d.x-c.x,d.y-c.y,b.x-c.x,b.y-c.y)<=eps));
end;
  
function include(var a:polygon;var b:point):boolean;
var
    i,tot:longint;
    k:point;
begin
    for i:=1 to a.tot do
      if on(a.a[i],a.a[i mod a.tot+1],b) then exit(true);
    tot:=0;
    k.x:=-100000;
    k.y:=100000;
    for i:=1 to a.tot do
      if jiao(a.a[i],a.a[i mod a.tot+1],k,b) then inc(tot);
    if tot and 1=1 then exit(true);
    exit(false);
end;
  
procedure get(var near,a:point;var dis:double;var b,c:point);
var
    d:double;
begin
    if (c.x-b.x)*(a.x-b.x)+(c.y-b.y)*(a.y-b.y)<=eps then
    begin
      d:=sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
      if dis>d+eps then
      begin
        near:=b;
        dis:=d;
      end;
      exit;
    end;
    if (c.x-b.x)*(a.x-c.x)+(c.y-b.y)*(a.y-c.y)>=-eps then
    begin
      d:=sqrt(sqr(a.x-c.x)+sqr(a.y-c.y));
      if dis>d+eps then
      begin
        near:=c;
        dis:=d;
      end;
      exit;
    end;
    d:=cj(c.x-b.x,c.y-b.y,a.x-b.x,a.y-b.y)/sqrt(sqr(b.x-c.x)+sqr(b.y-c.y));
    if dis>abs(d)+eps then
    begin
      dis:=abs(d);
      near.x:=a.x+(c.y-b.y)*d/sqrt(sqr(b.x-c.x)+sqr(b.y-c.y));
      near.y:=a.y-(c.x-b.x)*d/sqrt(sqr(b.x-c.x)+sqr(b.y-c.y));
    end;
end;
  
procedure find(var a,b:point);
var
    i,j:longint;
    dis:double;
begin
    for i:=1 to c do
      if include(p[i],b) then
      begin
        a:=b;
        exit;
      end;
    dis:=1<<30;
    for i:=1 to c do
      for j:=1 to p[i].tot do
        get(a,b,dis,p[i].a[j],p[i].a[j mod p[i].tot+1]);
    ans:=max(ans,dis);
end;
  
procedure init;
var
    i,j:longint;
begin
    read(c,tail);
    head:=1;
    for i:=1 to tail do
      with s[i].a do
        read(x,y);
    for i:=1 to c do
      with p[i] do
        begin
          read(tot);
          for j:=1 to tot do
            with a[j] do
              read(x,y);
        end;
    for i:=1 to tail do
      with s[i] do
        find(neara,a);
    for i:=1 to tail-1 do
      begin
        s[i].b:=s[i+1].a;
        s[i].nearb:=s[i+1].neara;
      end;
end;
  
procedure work;
var
    l,r,mid:point;
    d:double;
begin
    while head<>tail do
      begin
        l:=s[head].a;
        r:=s[head].b;
        while (sqrt(sqr(l.x-r.x)+sqr(l.y-r.y))>1e-4) do
          begin
            mid.x:=(l.x+r.x)/2;
            mid.y:=(l.y+r.y)/2;
            with s[head] do
            if sqrt(sqr(mid.x-neara.x)+sqr(mid.y-neara.y))<sqrt(sqr(mid.x-nearb.x)+sqr(mid.y-nearb.y)) then l:=mid
            else r:=mid;
          end;
        with s[head] do
        d:=max(sqrt(sqr(l.x-neara.x)+sqr(l.y-neara.y)),sqrt(sqr(l.x-nearb.x)+sqr(l.y-nearb.y)));
        mid.x:=(s[head].a.x+s[head].b.x)/2;
        mid.y:=(s[head].a.y+s[head].b.y)/2;
        find(l,mid);
        if d>ans+0.001 then
        begin
          s[tail].a:=s[head].a;
          s[tail].b:=mid;
          s[tail].nearb:=l;
          s[tail].neara:=s[head].neara;
          tail:=tail mod maxq+1;
          s[tail].a:=mid;
          s[tail].b:=s[head].b;
          s[tail].neara:=l;
          s[tail].nearb:=s[head].nearb;
          tail:=tail mod maxq+1;
        end;
        head:=head mod maxq+1;
      end;
    write(ans:0:2);
end;
  
begin
    init;
    work;
end.