三分法——求解凸性函数的极值问题

二分法作为分治中最常见的方法，适用于单调函数，逼近求解某点的值。但当函数是凸性函数时，二分法就无法适用，这时三分法就可以“大显身手”～～



如图，类似二分的定义Left和Right，mid = (Left + Right) / 2，midmid = (mid + Right) / 2; 如果mid靠近极值点，则Right = midmid；否则(即midmid靠近极值点)，则Left = mid;

程序模版如下：

double Calc(Type a)
{
/* 根据题目的意思计算 */
}

void Solve(void)
{
double Left, Right;
double mid, midmid;
double mid_value, midmid_value;
Left = MIN; Right = MAX;
while (Left + EPS < Right)
{
mid = (Left + Right) / 2;
midmid = (mid + Right) / 2;
mid_area = Calc(mid);
midmid_area = Calc(midmid);
// 假设求解最大极值.
if (mid_area >= midmid_area) Right = midmid;
else Left = mid;
}
}

现根据几道的OJ题目来分析三分法的具体实现。

描述

In this problem, you're to calculate the distance between a point P(xp, yp, zp) and a segment (x1, y1, z1) ? (x2, y2, z2), in a 3D space, i.e. the minimal distance from P to any point Q(xq, yq, zq) on the segment (a segment is part of a line).

输入

The first line contains a single integer T (1 ≤ T ≤ 1000), the number of test cases. Each test case is a single line containing 9 integers xp, yp, zp, x1, y1, z1, x2, y2, z2. These integers are all in [-1000,1000].

输出

For each test case, print the case number and the minimal distance, to two decimal places.

样例输入

3

0 0 0 0 1 0 1 1 0

1 0 0 1 0 1 1 1 0

-1 -1 -1 0 1 0 -1 0 -1

样例输出

Case 1: 1.00

Case 2: 0.71

Case 3: 1.00题意为在一条线段上找到一点，与给定的Ｐ点距离最小。很明显的凸性函数，用三分法来解决。

Calc函数即为求某点到P点的距离。

#include<iostream>
#include<cstdlib>
#include<stdio.h>
#include<memory.h>
#include<math.h>
#define eps 1e-8
using namespace std;

struct point
{
double x,y,z;
};
point l,r,p;
point MID(point a,point b)
{
point k;
k.x=(a.x+b.x)*0.5;
k.y=(a.y+b.y)*0.5;
k.z=(a.z+b.z)*0.5;
return k;
}
double dist(point a,point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
int sgn(double a)
{
return (a>eps)-(a<-eps);
}
point work()
{
point mid,midmid;
while(sgn(dist(l,r))>0)
{
mid=MID(l,r);
midmid=MID(mid,r);
if(dist(mid,p)<dist(midmid,p))
r=midmid;
else
l=mid;
}
return r;
}
int main()
{
int t,i,j,k,cas=1;
scanf("%d",&t);
while(t--)
{
scanf("%lf%lf%lf",&p.x,&p.y,&p.z);
scanf("%lf%lf%lf",&l.x,&l.y,&l.z);
scanf("%lf%lf%lf",&r.x,&r.y,&r.z);
l=work();
printf("Case %d: %.2lf\n",cas++,dist(l,p));
}
return 0;
}

ZOJ 3203 Light Bulb



如图，人左右走动，求影子L的最长长度。

根据图，很容易发现当灯，人的头部和墙角成一条直线时(假设此时人站在A点)，此时的长度是影子全在地上的最长长度。当人再向右走时，影子开始投影到墙上，当人贴着墙，影子长度即为人的高度。所以当人从A点走到墙，函数是先递增再递减，为凸性函数，所以我们可以用三分法来求解。

下面只给出Calc函数，其他直接套模版即可。

double Calc(double x)

{

return (h * D - H * x) / (D - x) + x;

}

//-------common header---------------
#include <stdio.h>
#include <vector>
#include <stack>
#include <math.h>
#include <algorithm>
typedef unsigned long u32;
using namespace std;
int main()
{
int caseNum = 0;
scanf("%d",&caseNum);
while(caseNum--)
{
double H,h,D;
scanf("%lf %lf %lf",&H,&h,&D);
double result = 0.0;
if(H-h>=D)
{
result = h;
}
else if(H*H<=D*(H-h))
{
result = h*D/H;
}
else
{
double a = -sqrt((H-h)*D)+H;
result = D*(h-a)/(H-a)+a;
}

printf("%.3f/n",(float)result);
}
return 1;
}

heru 5081 Turn the corner 08年哈尔滨regional网络赛



汽车拐弯问题，给定X, Y, l, d判断是否能够拐弯。首先当X或者Y小于d，那么一定不能。

其次我们发现随着角度θ的增大，最大高度ｈ先增长后减小，即为凸性函数，可以用三分法来求解。

这里的Calc函数需要比较繁琐的推倒公式：

s = l * cos(θ) + w * sin(θ) - x;

h = s * tan(θ) + w * cos(θ);

其中s为汽车最右边的点离拐角的水平距离, h为里拐点最高的距离, θ范围从0到90。

#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
double pi = acos(-1.0);
double x,y,l,w,s,h;
double cal(double a)
{
s = l*cos(a)+w*sin(a)-x;
h = s*tan(a)+w*cos(a);
return h;
}
int main()
{
double left,right,mid,midmid;
while(scanf("%lf%lf%lf%lf",&x,&y,&l,&w)!=EOF)
{
left = 0.0;
right = pi/2;
while(fabs(right-left)>1e-8)
{
mid = (left+right)/2;
midmid = (mid+right)/2;
if(cal(mid)>=cal(midmid))right = midmid;
else left = mid;
}
if(cal(mid)<=y)printf("yes\n");
else printf("no\n");
}
return 0;
}

POJ 3301 Texas Trip

题意为给定n(n <= 30)个点，求出饱含这些点的面积最小的正方形。

有两种解法，一种为逼近法，就是每次m分角度，求出最符合的角度，再继续m分，如此进行times次，即可求出较为精确的解。(m 大概取10, times取30即可)

第二种解法即为三分法，首先旋转的角度只要在0到180度即可，超过180度跟前面的相同的。坐标轴旋转后，坐标变换为：

X’ = x * cosa - y * sina;

y’ = y * cosa + x * sina;

至于这题的函数是否是凸性的，为什么是凸性的，我也无法给出准确的证明，希望哪位路过的大牛指点一下～～

#include <cmath>
#include <iostream>
using namespace std;

#define MAX 33
#define eps 0.00000005
#define max(a,b) (a>b?a:b)
struct { double x, y; } p[MAX];

double cal1 ( int n, double d )
{
double dis1 = 0.0, temp;
for ( int i = 1; i < n; i++ )
{
for ( int j = i + 1; j <= n; j++ )
{
temp = fabs( (p[i].y-p[j].y) * sin(d) + (p[i].x-p[j].x) * cos(d));
dis1 = max ( dis1, temp );
}
}
return dis1;
}

double cal2 ( int n, double d )
{
double dis2 = 0.0, temp;
for ( int i = 1; i < n; i++ )
{
for ( int j = i + 1; j <= n; j++ )
{
temp = fabs( (p[i].y-p[j].y) * cos(d) - (p[i].x-p[j].x) * sin(d));
dis2 = max ( dis2, temp );
}
}
return dis2;
}

int main()
{
int cs, n;
scanf("%d",&cs);
while ( cs-- )
{
scanf("%d",&n);
for ( int i = 1; i <= n; i++ )
scanf("%lf%lf", &p[i].x, &p[i].y);

double ans = 1000, low = 0, high = asin(1.0);
while ( high - low > eps )
{
double mid1 = low + ( high - low ) / 3;
double mid2 = high - ( high - low ) / 3;
double len1 = max ( cal1 ( n, mid1 ), cal2 ( n, mid1 ));
double len2 = max ( cal1 ( n, mid2 ), cal2 ( n, mid2 ));

if ( len1 < len2 )
{
ans = len1;
high = mid2;
}
else
{
ans = len2;
low = mid1;
}
}
printf("%.2lf\n",ans*ans);
}
return 0;
}

hdu 3400 Line belt

典型的三分法，先三分第一条线段，找到一个点，然后根据这个点再三分第二条线段即可，想出三分的思路基本就可以过了。

对于求解一些实际问题，当公式难以推导出来时，二分、三分法可以较为精确地求解出一些临界值，且效率也是令人满意的。

#include<iostream>
#include<cmath>
using namespace std;
#define eps 1e-6

struct POINT
{
double x,y;
}a,b,c,d;
double p,q,r;

double dis(POINT A,POINT B)
{
return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+eps);
}

double ftime(POINT t,double distd)//time between t and D
{
POINT temp;
temp.x=d.x+(c.x-d.x)*((q*distd)/dis(c,d));
temp.y=d.y+(c.y-d.y)*((q*distd)/dis(c,d));
return dis(t,temp)/r+distd;
}

double time(double posa)//posa->time between t and A
{
POINT t;
t.x=a.x+(b.x-a.x)*((posa*p)/dis(a,b));
t.y=a.y+(b.y-a.y)*((posa*p)/dis(a,b));
//test distance on CD
double low=0,high=dis(c,d)/q,mid,m;
while(high-low>eps)
{
m=(high+low)/2;
mid=(high+m)/2;
if(ftime(t,m)<=ftime(t,mid))
high=mid;
else
low=m;
}
return ftime(t,low)+posa;
}

int main()
{
freopen("1.txt","w",stdout);
int T;
cin>>T;
while(T--)
{
cin>>a.x>>a.y>>b.x>>b.y>>c.x>>c.y>>d.x>>d.y;
cin>>p>>q>>r;
double low=0,high=dis(a,b)/p,mid,m;
while(high-low>eps)
{
m=(low+high)/2;
mid=(low+m)/2;
if(time(mid)>=time(m))
low=mid;
else
high=m;
}
printf("%.2lf/n",time(low));
}
}