三分法小结

昨天晚上快要回宿舍的时候偶然发现一道看起来简单的计算几何题，搜了一下题解，遇见一个新的方法，三分法，以前只是听过二分法，于是今天刷了几道三分法的题目。

 （上图来自http://hi.baidu.com/czyuan_acm/item/81b21d1910ea729c99ce33db）

     二分法主要用来搜索单调函数的某个值，而三分法主要用于凸（凹）函数的极值。以凸函数为例，如果mid=(left+right)/2,midmid=(mid+right)/2,fun(mid)>=fun(midmid),则right=midmid;

例1：hdu(4454)  stealing a cake

给你一个点，一个圆，一个矩形，问从这个点到圆，然后到矩形的最短距离。

       就是从这一题注意三分法的。注意到在（0，PI），（PI,2*PI)内距离先减小后增加（sad,不会证明T^T……），分别三分角度即可

#include <cmath>#include <cstdio>#include <cstring>#include <iostream>using namespace std;const double eps=1e-8;const double pi=acos(-1);struct Point{    double x,y;    Point(double x=0,double y=0)        :x(x),y(y) { }};Point operator-(const Point &lhs,const Point &rhs){    return Point(lhs.x-rhs.x,lhs.y-rhs.y);}int dcmp(double a){    if(fabs(a)<=eps)    return 0;    return a<0?-1:1;}bool operator==(const Point &a,const Point &b){    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}double cross(const Point &a,const Point &b){    return a.x*b.y-a.y*b.x;}double dot(const Point &a,const Point &b){    return a.x*b.x+a.y*b.y;}double length(const Point &a){    return sqrt(dot(a,a));}double dist2line(const Point &p,const Point &a,const Point &b){    if(a==b)    return length(p-a);    Point v1=b-a,v2=p-a,v3=p-b;    if(dcmp(dot(v1,v2))<0)  return length(v2);    if(dcmp(dot(v1,v3))>0)  return length(v3);    return fabs(cross(v1,v2))/length(v1);}double dist2rec(const Point &p,Point *rec){    double ans=dist2line(p,rec[0],rec[1]);    ans=min(dist2line(p,rec[1],rec[2]),ans);    ans=min(dist2line(p,rec[2],rec[3]),ans);    ans=min(dist2line(p,rec[3],rec[0]),ans);    return ans;}Point getPoint(const Point &c,int r,double a){    return Point(c.x+cos(a)*r,c.y+sin(a)*r);}double solve(Point *rec,const Point &c,const Point &p,int R){    double l=0,r=pi;    while(l+eps<=r)    {        double mid=(l+r)/2;        double midmid=(mid+r)/2;        Point p1=getPoint(c,R,mid),p2=getPoint(c,R,midmid);        double dis1=dist2rec(p1,rec)+length(p-p1);        double dis2=dist2rec(p2,rec)+length(p-p2);        if(dis1>dis2)   l=mid;        else r=midmid;    }    Point p0=getPoint(c,R,l);    double ans1=dist2rec(p0,rec)+length(p-p0);    l=pi,r=2*pi;    while(l+eps<=r)    {        double mid=(l+r)/2;        double midmid=(mid+r)/2;        Point p1=getPoint(c,R,mid),p2=getPoint(c,R,midmid);        double dis1=dist2rec(p1,rec)+length(p-p1);        double dis2=dist2rec(p2,rec)+length(p-p2);        if(dis1>dis2)   l=mid;        else r=midmid;    }    p0=getPoint(c,R,l);    double ans2=dist2rec(p0,rec)+length(p-p0);    return min(ans1,ans2);}int main(){    double x,y,r;    //freopen("4544.in","r",stdin);    while(scanf("%lf%lf",&x,&y)!=EOF)    {        if(x==0.0&&y==0.0)  break;        Point p(x,y);        scanf("%lf%lf%lf",&x,&y,&r);        Point c(x,y);        Point rec[4];        double x1,y1;        scanf("%lf%lf%lf%lf",&x,&y,&x1,&y1);        rec[0]=Point(x,y),rec[1]=Point(x1,y);        rec[2]=Point(x1,y1),rec[3]=Point(x,y1);        printf("%.2lf\n",solve(rec,c,p,r));    }    return 0;}

</pre><p></p><p>例2 poj(3301) <a target=_blank href="http://http://poj.org/problem?id=3301">Texas Trip</a></p><p>     给你不超过30个点，问这些点的最小包围正方形。角度的枚举范围在0-PI,然后根据坐标旋转公式求出旋转后的点，然后统计一下。。。</p><p><pre name="code" class="cpp">#include <cmath>#include <cstdio>#include <cstring>#include <iostream>using namespace std;const int maxn=34;const double pi=acos(-1.0);const double inf=1<<30;const double eps=1e-12;struct Point{    double x,y;    Point(double x=0,double y=0)        :x(x),y(y)  { }}p[maxn];double check(double angle,int n){    double minx=inf,miny=inf,maxy=-inf,maxx=-inf;    double tcos=cos(angle),tsin=sin(angle);    for(int i=0;i<n;i++)    {        double tx=tcos*p[i].x-tsin*p[i].y;        double ty=tsin*p[i].x+tcos*p[i].y;        minx=min(minx,tx),miny=min(miny,ty);        maxx=max(maxx,tx),maxy=max(maxy,ty);    }    return max(maxx-minx,maxy-miny);}int main(){    int t,n;    //freopen("3301.in","r",stdin);    scanf("%d",&t);    while(t--)    {        scanf("%d",&n);        for(int i=0;i<n;i++)            scanf("%lf%lf",&p[i].x,&p[i].y);        double l=0,r=pi;        while(l+eps<r)        {            double mid=(l+r)/2;            double midmid=(mid+r)/2;            double d1=check(mid,n),d2=check(midmid,n);            if(d1<=d2)                r=midmid;            else                l=mid;        }        double ans=check(l,n);        printf("%.2f\n",ans*ans);    }    return 0;}

例3：hdu(3400) Line belt

     给出两条线段，ab,cd,问从a点到d点的最短时间。不难想到，最优策略是在ab上走一段，然后走到cd上某一点，然后走向d点。ab上的时间函数是单调递增，从ab上某点到d点的时间函数是凹函数，因此可以三分ab上的点，确定了这个点之后，三分cd上的那一点，三分套三分。

    

#include <cmath>#include <cstdio>#include <cstring>#include <iostream>using namespace std;const double eps=1e-8;double inf=1<<26;struct Point{    double x,y;    Point(double x=0,double y=0)        :x(x),y(y) { }};Point getPoint(){    double x,y;    scanf("%lf%lf",&x,&y);    return Point(x,y);}Point a,b,c,d;double vab,vcd,vo;double distan(const Point &a,const Point &b){    return sqrt( (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) );}double getLength(double fir,double sec){    Point p1,p2;    p1.x=a.x+(b.x-a.x)*fir;    p1.y=a.y+(b.y-a.y)*fir;    p2.x=c.x+(d.x-c.x)*sec;    p2.y=c.y+(d.y-c.y)*sec;    return distan(a,p1)/vab+distan(p1,p2)/vo+distan(p2,d)/vcd;}double solve_cd(double fir){    double l=0,r=1;    while(l+eps<=r)    {        double mid=(l+r)/2;        double midmid=(mid+r)/2;        double ans1=getLength(fir,mid);        double ans2=getLength(fir,midmid);        if(ans1<=ans2)            r=midmid;        else            l=mid;    }    return getLength(fir,l);}double solve_abcd(){    double l=0,r=1;    while(l+eps<=r)    {        double mid=(l+r)/2;        double midmid=(mid+r)/2;        double ans1=solve_cd(mid);        double ans2=solve_cd(midmid);        if(ans1<=ans2)            r=midmid;        else            l=mid;    }    return solve_cd(l);}int main(){    int t;    //freopen("3400.in","r",stdin);    scanf("%d",&t);    while(t--)    {        a=getPoint(),b=getPoint();        c=getPoint(),d=getPoint();        scanf("%lf%lf%lf",&vab,&vcd,&vo);        printf("%.2lf\n",solve_abcd());    }    return 0;}