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

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

       如图，类似二分的定义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));
	}
}