求凸函数的极值的一般方法是三分
三分的思想大概是这样的：
例如我们要求下凸函数的极值
在区间[L,R]上，
我们定义m1为区间的第一个三等分点
定义m2为区间的第二个三等分点
设函数值为F(x)
则若F(m1)<F(m2)，证明解在[L,m2]中
否则解在[m1,R]中
一般三分的写法是迭代，注意控制精度和时间的平衡
UVa 1476
很容易发现一堆二次函数求max之后还是一个凸函数
之后三分即可
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std; const int maxn=10010;
int T,n;
int a[maxn],b[maxn],c[maxn]; double F(double x){
double ans=a[1]*x*x+b[1]*x+c[1];
for(int i=2;i<=n;++i)ans=max(ans,a[i]*x*x+b[i]*x+c[i]);
return ans;
} int main(){
scanf("%d",&T);
while(T--){
scanf("%d",&n);
for(int i=1;i<=n;++i)scanf("%d%d%d",&a[i],&b[i],&c[i]);
double L=0.0,R=1000.0;
for(int i=1;i<=100;++i){
double m1=L+(R-L)/3,m2=m1+(R-L)/3;
if(F(m1)<F(m2))R=m2;
else L=m1;
}
printf("%.4lf\n",F(L));
}return 0;
}
ZOJ 3203
很坑的题目，首先我们先暴力导出公式发现是个凸函数
之后分类讨论影子会不会射到墙上(其实我的代码trick了一下，实在懒得分类讨论了）
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<cmath>
#define eps 1e-6
using namespace std; int T;
double H,h,D; double F(double M){
return h*(H*(D-M)-h*D)/(h*(D-M)-h*D)+(D-M);
//return (h*(H*M-h*D)+M*(h*M-h*D))/(h*M-h*D);
} int main(){
scanf("%d",&T);
while(T--){
scanf("%lf%lf%lf",&H,&h,&D);
double L=0.0,R=min(H,D);
for(int i=1;i<=10000;++i){
double len=(R-L)/3;
double m1=L+len,m2=m1+len;
if(F(m1)>F(m2))R=m2;
else L=m1;
}
printf("%.3lf\n",F(L));
}return 0;
}
BZOJ 1857 SCOI 传送带
首先我们设在第一个传送带上走x距离，第二个传送带上走y距离
因为中间走的是欧几里得距离
所以不难发现对于x整体是一堆类二次函数，由于UVa的题目的解决经验
之后我们尝试三分x
然后对于y，现在的问题就转换成了给定一个类二次函数，最小化G(y)
这也是可以三分的
所以我们就可以三分求出F(x)的极值了
第一次提交TLE是因为三分的迭代次数设的太多了
第二次提交WA是因为计算的时候出现了除0错误(并不知道为什么不会RE）
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#define eps 1e-8
using namespace std; double L1,L2;
double P,Q,R;
struct Point{
double x,y;
void read(){scanf("%lf%lf",&x,&y);}
Point(double x=0,double y=0):x(x),y(y){}
}A,B,C,D,u,v,nowA,nowD;
typedef Point Vector;
Vector operator -(const Point &A,const Point &B){return Vector(A.x-B.x,A.y-B.y);}
Vector operator +(const Point &A,const Point &B){return Vector(A.x+B.x,A.y+B.y);}
Vector operator *(const Point &A,const double &p){return Vector(A.x*p,A.y*p);}
double Dot(const Point &A,const Point &B){return A.x*B.x+A.y*B.y;}
double Len(const Vector &A){return sqrt(Dot(A,A));} double G(double x){
double t=0.0;
if(fabs(Dot(v,v))>eps)t=sqrt(x*x/Dot(v,v));
nowD=D+v*t;
return Len(nowD-nowA)/R+x/Q;
}
double F(double x){
double t=0.0;
if(fabs(Dot(u,u))>eps)t=sqrt(x*x/(Dot(u,u)));
nowA=A+u*t;
double L=0.0,R=L2;
for(int i=1;i<=200;++i){
double len=(R-L)/3;
double m1=L+len,m2=m1+len;
if(G(m1)<G(m2))R=m2;
else L=m1;
}return G(L)+x/P;
} int main(){
A.read();B.read();u=B-A;L1=Len(u);
C.read();D.read();v=C-D;L2=Len(v);
scanf("%lf%lf%lf",&P,&Q,&R);
double L=0.0,R=L1;
for(int i=1;i<=200;++i){
double len=(R-L)/3;
double m1=L+len,m2=m1+len;
if(F(m1)<F(m2))R=m2;
else L=m1;
}
printf("%.2lf\n",F(L));
return 0;
}
POJ 3301
完全想不到的神思路
首先我们注意到平行于坐标轴的最大正方形非常好算
最后的正方形一定是某个平行于坐标轴的正方形旋转了多少的角度之后得到的
所以我们三分这个旋转的角度，就可以得到答案啦
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<cmath>
using namespace std; const int maxn=52;
const double pi=acos(-1.0);
int T,n;
struct Point{
double x,y;
void read(){scanf("%lf%lf",&x,&y);}
}A[maxn],B[maxn]; double F(double x){
for(int i=1;i<=n;++i){
B[i].x=A[i].x*cos(x)-A[i].y*sin(x);
B[i].y=A[i].x*sin(x)+A[i].y*cos(x);
}
double mxy=-1e12,mny=1e12,mxx=-1e12,mnx=1e12;
for(int i=1;i<=n;++i){
mxy=max(mxy,B[i].y);
mny=min(mny,B[i].y);
mxx=max(mxx,B[i].x);
mnx=min(mnx,B[i].x);
}
double len=max(mxy-mny,mxx-mnx);
return len*len;
} int main(){
scanf("%d",&T);
while(T--){
scanf("%d",&n);
for(int i=1;i<=n;++i)A[i].read();
double L=0,R=2*pi;
for(int i=1;i<=300;++i){
double len=(R-L)/3;
double m1=L+len,m2=m1+len;
if(F(m1)<F(m2))R=m2;
else L=m1;
}printf("%.2lf\n",F(L));
}return 0;
}
总结：感觉三分这个算法需要注意的东西很多的
首先，我们需要对题目建立模型，并发现它是个凸函数(几何画板胡乱证明之）
其次，在写三分的过程中，要注意迭代次数
迭代次数过少会导致精度不够
过多会导致TLE
最后，在计算F(x)也就是函数值的时候要尽量避免精度误差，并且防止除0错误