http://hihocoder.com/contest/hiho40/problem/1
描述
这一次我们就简单一点了,题目在此:
在之前的几周中我们了解到二分法作为分治中最常见的方法,适用于单调函数,逼近求解某点的值。
但当函数是凸形函数时,二分法就无法适用,这时就需要用到三分法。
从三分法的名字中我们可以猜到,三分法是对于需要逼近的区间做三等分:
我们发现lm这个点比rm要低,那么我们要找的最小点一定在[left,rm]之间。如果最低点在[rm,right]之间,就会出现在rm左右都有比他低的点,这显然是不可能的。 同理,当rm比lm低时,最低点一定在[lm,right]的区间内。
利用这个性质,我们就可以在缩小区间的同时向目标点逼近,从而得到极值。
接下来我们回到题目上,抛物线和点之间的距离可以简单的用直线公式计算:
即d = min{sqrt((X - x)^2+(aX^2+bX+c-y)^2)}
该公式展开后为4次,需要采用求导等方法来求极值。对于计算机编程来说是很麻烦的一件事。
进一步观察题目,我们可以发现根据带入的X值不同,d的长度恰好满足凸形函数。
而我们要求的最短距离d,正好就是这个凸形函数的极值。
那么三分法不就正好可以用来解决这道题目了么?
需要注意在解题过程中一定要想清楚如何划分区间,我们求的各个变量到底是什么含义。
另外,这道题还有一个小小的trick,在解决的时候请一定要小心。
输入
第1行:5个整数a,b,c,x,y。前三个数构成抛物线的参数,后两个数x,y表示P点坐标。-200≤a,b,c,x,y≤200
输出
第1行:1个实数d,保留3位小数(四舍五入)
样例输入
2 8 2 -2 6
样例输出
2.437
1 #include <iostream> 2 #include <math.h> 3 using namespace::std; 4 int a; 5 double Xp,Yp; 6 const double precision=1E-6; 7 8 inline double d(double x) 9 { 10 return sqrt((Xp-x)*(Xp-x) + (Yp-a*x*x)*(Yp-a*x*x)); 11 } 12 13 14 int main(int argc, const char* argv[]) 15 { 16 int b,c; 17 int X0,Y0; 18 cin>>a>>b>>c>>X0>>Y0; 19 double D; 20 if (a==0) { 21 D=fabs(1.*(b*X0-Y0+c)/sqrt(b*b+1)); 22 }else{ 23 double Xc=-1.*b/(2*a); 24 double Yc=c-1.*b*b/(4*a); 25 Xp=X0-Xc; 26 Yp=Y0-Yc; 27 double x1,x2; 28 if (a<0){ 29 a=-a; 30 Yp=-Yp; 31 } 32 if (Xp<0)Xp=-Xp; 33 if (Yp<0) { 34 x1=0.; 35 x2=Xp; 36 }else{ 37 x1=Xp; 38 x2=sqrt(Yp/a); 39 if (x1>x2)swap(x1, x2); 40 } 41 double lmx,rmx,lm,rm; 42 while (fabs(d(x2)-d(x1))>precision) { 43 lmx=x1+1./3*(x2-x1); 44 lm=d(lmx); 45 rmx=x1+2./3*(x2-x1); 46 rm=d(rmx); 47 if (lm<rm){ 48 x2=rmx; 49 }else{ 50 x1=lmx; 51 } 52 53 } 54 D=d((x1+x2)/2); 55 } 56 printf("%.3lf\n",D); 57 return 0; 58 }
以前只知道怎么求零点,求函数极值一般是求导后求零点。
现在大开眼界三分法可以求凸函数的极值,想想以后可能会学到巧妙方法求鞍点之类的问题就兴奋不已。
好像没什么可说的,后面计算某些参数要除以a,所以a判断是不是零很重要,如果是a==0的话就用点到直线的距离公式就可以了。
不是零的话就把整个体系变到第一象限去,这样方便点。关于区间的选择(x1,x2)暂时没想到更好的办法,就是代码中的笨办法。
要保留3位小数,当时precision=1E-4的时候d通过不了,后来索性就改到1e-6了,个人比较喜欢6这个数字。