题解-洛谷P2571 [SCOI2010]传送带
本文主要讲述了某蒟蒻历经一天一夜尝试 \(O(1)\) 做法的一点点思考，以及对梯度下降法及各种优化算法较为详细的介绍
如有不足之处，或有更好的想法，望不吝赐教～
题目传送门
\(O(1)\) 做法の尝试
推柿子
容易得到，最优路径一定是从 \(A\) 移动到 \(AB\) 上的某点 \(M\)，再由 \(M\) 经平面移动到 \(CD\) 上的某点 \(N\)，最后由 \(N\) 移动到 \(D\)。
所花时间：
\[t=\frac{AM}P+\frac{MN}R+\frac{ND}Q
\]
找到最优的 \(M\)、\(N\) 即可。
直接设坐标比较麻烦，可以设线段的比值「取 \(0\) ~ \(1\)」
设\(AM=xAB,ND=yCD(x,y\in[0,1])\)
则：
\[x_M-x_A=x(x_B-x_A),y_M-y_A=x(y_B-y_A)
\\
x_M=x_A+x(x_B-x_A),y_M=y_A+x(y_B-y_A)
\]
同理：
\[x_N-x_D=y(x_C-x_D),y_N-y_D=y(y_C-y_D)
\\
x_N=x_D+y(x_C-x_D),y_N=y_D+y(y_C-y_D)
\]
\[t=\frac{x\cdot AB}P+\frac{MN}R+\frac{y\cdot CD}Q
\]
\[\begin{aligned}
MN&=\sqrt{(x_M-x_N)^2+(y_M-y_N)^2}
\\&=\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D))^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D))^2)}
\end{aligned}
\]
这一大串式子看起来很难受，不妨简化一下。
设：
\[c_1=\frac{AB}P\qquad c_2=\frac{CD}Q\qquad c_3=\frac1R
\]
则
\[t=c_1\cdot x+c_2\cdot y+c_3\cdot MN
\]
设：
\[\begin{aligned}
&c_4=x_B-x_A\quad&c_5=x_D-x_C\qquad\quad&c_6=x_A-x_D
\\
&c_7=y_B-y_A\quad&c_8=y_D-y_C\qquad\quad&c_9=y_A-y_D
\end{aligned}
\]
则
\[MN=\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}
\]
\[t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}
\]
c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9))
答案即为求解 \(f(x,y)「x,y\in[0,1]」\) 的最小值。
对 \(x\)、\(y\) 求偏导：
\[\begin{aligned}
\frac{\partial f}{\partial x}&=c_1+c_3\cdot\frac
{2\cdot(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+2\cdot(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}
{2\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}
\\
&=c_1+c_3\cdot\frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}
{\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}
\\
\frac{\partial f}{\partial y}&=c_2+c_3\cdot\frac
{2\cdot(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+2\cdot(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}
{2\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}
\\
&=c_2+c_3\cdot\frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}
{\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}
\end{aligned}
\]
令 \(\frac{\partial f}{\partial x}=\frac{\partial f}{\partial y}=0\)，解出所有驻点：
\[\begin{cases}
\begin{aligned}
c_1+c_3\cdot\frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7}
{\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}=0
\\
c_2+c_3\cdot\frac{(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8}
{\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}}=0
\end{aligned}
\end{cases}
\]
将 \(c_1\)、\(c_2\) 移到等号右边，将 \(\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}\) 乘过去，两边开平方：
\[\begin{cases}
\begin{align}
c_3^2\cdot[(c_4\cdot x+c_5\cdot y+c_6)\cdot c_4+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_7]^2
=c_1^2\cdot[(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2]
\\
c_3^2\cdot[(c_4\cdot x+c_5\cdot y+c_6)\cdot c_5+(c_7\cdot x+c_8\cdot y+c_9)\cdot c_8]^2
=c_2^2\cdot[(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2]
\end{align}
\end{cases}
\]
这样一来，方程组的形式就优美了许多，我们可以轻松求解。
用正则表达式稍作整理：
\[c3^2*((c4*x+c5*y+c6)*c4+(c7*x+c8*y+c9)*c7)^2=c1^2*((c4*x+c5*y+c6)^2+(c7*x+c8*y+c9)^2),c3^2*((c4*x+c5*y+c6)*c5+(c7*x+c8*y+c9)*c8)^2=c2^2*((c4*x+c5*y+c6)^2+(c7*x+c8*y+c9)^2)
\]
塞到求解方程组的计算器中，即可求出 \(x\)、\(y\) 的值：
$$
\begin{cases}
\begin{aligned}
x&=\frac{c5\cdot c9-c6\cdot c8}{c4\cdot c8-c5\cdot c7}
\\y&=-\frac{c4\cdot c9-c6\cdot c7}{c4\cdot c8-c5\cdot c7}=\frac{c6\cdot c7-c4\cdot c9}{c4\cdot c8-c5\cdot c7}
\end{aligned}
\end{cases}
$$
x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7)
将 \(x\)、\(y\) 代入原方程，可以利用表达式化简工具验证一下：
\[(c3)^2*(((c4)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c5)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c6)*(c4)+((c7)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c8)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c9))*(c7))^2-((c1)^2*(((c4)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c5)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c6))^2+((c7)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c8)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c9))^2))
\]
\[(c3)^2*(((c4)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c5)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c6))*(c5)+((c7)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c8)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c9))*(c8))^2-((c2)^2*(((c4)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c5)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c6))^2+((c7)*(c5*c9-c6*c8)/(c4*c8-c5*c7)+(c8)*(c6*c7-c4*c9)/(c4*c8-c5*c7)+(c9))^2))
\]
上面两个式子的值都为 \(0\)，答案正确。
此时，\(x\)、\(y\) 的取值即为 \(f(x,y)\) 的极值点。
将 \(x\)、\(y\) 代入 \(f(x,y)\)，得：
\[c1*(c5*c9-c6*c8)/(c4*c8-c5*c7)+c2*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c3*sqrt((c4*(c5*c9-c6*c8)/(c4*c8-c5*c7)+c5*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c6)*(c4*(c5*c9-c6*c8)/(c4*c8-c5*c7)+c5*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c6)+(c7*(c5*c9-c6*c8)/(c4*c8-c5*c7)+c8*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c9)*(c7*(c5*c9-c6*c8)/(c4*c8-c5*c7)+c8*(c6*c7-c4*c9)/(c4*c8-c5*c7)+c9))
\]
化简，得：
\[((c1*c5-c2*c4)*c9-c1*c6*c8+c2*c6*c7)/(c4*c8-c5*c7)
\]
即
\[f(x,y)=\frac{(c_1\cdot c_5-c_2\cdot c_4)\cdot c_9-c_1\cdot c_6\cdot c_8+c_2\cdot c_6\cdot c_7}{c_4\cdot c_8-c_5\cdot c_7}
\]
((c1*c5-c2*c4)*c9-c1*c6*c8+c2*c6*c7)/(c4*c8-c5*c7)
这就完了？
当然没有。。。
首先，你不知道解出来的 \(f(x,y)\) 是极大值还是极小值。
其次，就算是极小值，你也别忘了 \(x,y\in[0,1]\) 的条件。
再者，边界值也有可能是合法区间内的最值。
因此，可以先判断 \(x\)，\(y\) 的范围，若均在 \([0,1]\) 内则用 \(f(x,y)\) 更新答案；然后用四个边界值 \(f(0,0),f(0,1),f(1,0),f(1,1)\) 更新答案。
code:
#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))
using namespace std;
long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;
struct node{
long double x,y;
}A,B,C,D,M,N;
inline long double dis(node a,node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline long double f(long double x,long double y)
{
return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}
signed main()
{
scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
printf("%.2Lf\n",ans);
return 0;
}
特判——「一元」
若只想到这里，你只能获得 \(60pts\)。
还有一点需要处理：分母为 \(0\) 的情况。
考虑 \(x、y、f(x,y)\) 在 \(c_4\cdot c_8-c_5\cdot c_7=0\) 时没有定义。
展开发现，此时 \((x_B-x_A)\cdot(y_D-y_C)=(x_D-x_C)\cdot(y_B-y_A)\)。
注意 \(AB\)、\(CD\) 可能为 \(0\)。
特判这两种情况：
设 \(CD=0\)，那么容易想到，最优路径一定是从 \(A\) 移动到 \(AB\) 上的某点 \(M\)，再由 \(M\) 经平面移动到 \(D\)。
显然这是原先的一种特殊情况：\(N\) 与 \(C、D\) 重合
此时只需对原式略作修改即可。
所花时间：
\[t=\frac{AM}P+\frac{MD}R=\frac{x\cdot AB}P+\frac{MD}R
\]
\[\begin{aligned}
MD=MN&=\sqrt{(x_M-x_N)^2+(y_M-y_N)^2}
\\&=\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D))^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D))^2)}
\\&=\sqrt{(x_A+x(x_B-x_A)-x_D)^2+(y_A+x(y_B-y_A)-y_D)^2)}
\end{aligned}
\]
容易发现，此时 \(y\) 的取值不再影响结果，\(t\) 变成了关于 \(x\) 的一元函数。
重新用一下前面定义的几个参数：
\[c_1=\frac{AB}P\qquad c_2=\frac{CD}Q=0\qquad c_3=\frac1R
\]
\[\begin{aligned}
&c_4=x_B-x_A\quad&c_5=x_D-x_C=0\qquad\quad&c_6=x_A-x_D
\\
&c_7=y_B-y_A\quad&c_8=y_D-y_C=0\qquad\quad&c_9=y_A-y_D
\end{aligned}
\]
\[MD=\sqrt{(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2}
\]
\[t=f(x)=c_1\cdot x+c_3\cdot\sqrt{(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2}
\]
答案即为求解 \(f(x)「x\in[0,1]」\)的最小值
直接求导：
\[\begin{aligned}
\frac{\mathrm df}{\mathrm dx}&=c_1+c_3\cdot\frac
{2\cdot(c_4\cdot x+c_6)\cdot c_4+2\cdot(c_7\cdot x+c_9)\cdot c_7}
{2\sqrt{(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2}}
\\
&=c_1+c_3\cdot\frac
{(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7}
{\sqrt{(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2}}
\end{aligned}
\]
令 \(\frac{\mathrm df}{\mathrm dx}=0\)，解出所有驻点：
\[c_1+c_3\cdot\frac
{(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7}
{\sqrt{(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2}}=0
\]
稍作化简：
\[c_3^2\cdot[(c_4\cdot x+c_6)\cdot c_4+(c_7\cdot x+c_9)\cdot c_7]^2
=c_1^2\cdot[(c_4\cdot x+c_6)^2+(c_7\cdot x+c_9)^2]
\]
不就是个一元二次方程嘛，手解也能解出来「就是懒得解」
再次用正则表达式整理：
\[c3^2*((c4*x+c6)*c4+(c7*x+c9)*c7)^2=c1^2*((c4*x+c6)^2+(c7*x+c9)^2)
\]
塞进那个工具里，然后就得到了两个解：
\[\begin{aligned}
x_1&={{\sqrt{c_{3}^2\,c_{7}^2+c_{3}^2\,c_{4}^2-c_{1}^2}\,\left(
c_{1}\,c_{4}\,c_{9}-c_{1}\,c_{6}\,c_{7}\right)+\left(\left(c_{1}^2-
c_{3}^2\,c_{4}^2\right)\,c_{7}-c_{3}^2\,c_{7}^3\right)\,c_{9}-c_{3}^
2\,c_{4}\,c_{6}\,c_{7}^2+\left(c_{1}^2\,c_{4}-c_{3}^2\,c_{4}^3
\right)\,c_{6}}\over{c_{3}^2\,c_{7}^4+\left(2\,c_{3}^2\,c_{4}^2-
c_{1}^2\right)\,c_{7}^2+c_{3}^2\,c_{4}^4-c_{1}^2\,c_{4}^2}}
\\
x_2&=-{{
\sqrt{c_{3}^2\,c_{7}^2+c_{3}^2\,c_{4}^2-c_{1}^2}\,\left(c_{1}\,c_{4}
\,c_{9}-c_{1}\,c_{6}\,c_{7}\right)+\left(c_{3}^2\,c_{7}^3+\left(
c_{3}^2\,c_{4}^2-c_{1}^2\right)\,c_{7}\right)\,c_{9}+c_{3}^2\,c_{4}
\,c_{6}\,c_{7}^2+\left(c_{3}^2\,c_{4}^3-c_{1}^2\,c_{4}\right)\,c_{6}
}\over{c_{3}^2\,c_{7}^4+\left(2\,c_{3}^2\,c_{4}^2-c_{1}^2\right)\,
c_{7}^2+c_{3}^2\,c_{4}^4-c_{1}^2\,c_{4}^2}}
\end{aligned}
\]
x1=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4)
x2=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4)
同理，当 \(AB=0\) 时，\(M\) 与 \(A、B\) 重合：
\[c_1=\frac{AB}P=0\qquad c_2=\frac{CD}Q\qquad c_3=\frac1R
\]
\[\begin{aligned}
&c_4=x_B-x_A=0\quad&c_5=x_D-x_C\qquad\quad&c_6=x_A-x_D
\\
&c_7=y_B-y_A=0\quad&c_8=y_D-y_C\qquad\quad&c_9=y_A-y_D
\end{aligned}
\]
\[t=\frac{AN}R+\frac{ND}Q=\frac{AN}R+\frac{y\cdot CD}Q=c_2\cdot y+c_3\cdot AN
\]
\[\begin{aligned}
AN=MN&=\sqrt{(x_M-x_N)^2+(y_M-y_N)^2}
\\&=\sqrt{(x_A+x(x_B-x_A)-x_D-y(x_C-x_D))^2+(y_A+x(y_B-y_A)-y_D-y(y_C-y_D))^2)}
\\&=\sqrt{(x_A-x_D-y(x_C-x_D))^2+(y_A-y_D-y(y_C-y_D))^2)}
\\&=\sqrt{(c_5\cdot y+c_6)^2+(c_8\cdot y+c_9)^2)}
\end{aligned}
\]
\[t=f(y)=c_2\cdot y+c_3\cdot\sqrt{(c_5\cdot y+c_6)^2+(c_8\cdot y+c_9)^2)}
\]
后面步骤相似，最后化简为：
\[c_3^2\cdot[(c_5\cdot y+c_6)\cdot c_5+(c_8\cdot y+c_9)\cdot c_8]^2
=c_2^2\cdot[(c_5\cdot y+c_6)^2+(c_8\cdot y+c_9)^2]
\]
\[c3^2*((c5*y+c6)*c5+(c8*y+c9)*c8)^2=c2^2*((c5*y+c6)^2+(c8*y+c9)^2)
\]
\[\begin{aligned}
y_1&={{\sqrt{c_{3}^2\,c_{8}^2+c_{3}^2\,c_{5}^2-c_{2}^2}\,\left(
c_{2}\,c_{5}\,c_{9}-c_{2}\,c_{6}\,c_{8}\right)+\left(\left(c_{2}^2-
c_{3}^2\,c_{5}^2\right)\,c_{8}-c_{3}^2\,c_{8}^3\right)\,c_{9}-c_{3}^
2\,c_{5}\,c_{6}\,c_{8}^2+\left(c_{2}^2\,c_{5}-c_{3}^2\,c_{5}^3
\right)\,c_{6}}\over{c_{3}^2\,c_{8}^4+\left(2\,c_{3}^2\,c_{5}^2-
c_{2}^2\right)\,c_{8}^2+c_{3}^2\,c_{5}^4-c_{2}^2\,c_{5}^2}}
\\
y_2&=-{{
\sqrt{c_{3}^2\,c_{8}^2+c_{3}^2\,c_{5}^2-c_{2}^2}\,\left(c_{2}\,c_{5}
\,c_{9}-c_{2}\,c_{6}\,c_{8}\right)+\left(c_{3}^2\,c_{8}^3+\left(
c_{3}^2\,c_{5}^2-c_{2}^2\right)\,c_{8}\right)\,c_{9}+c_{3}^2\,c_{5}
\,c_{6}\,c_{8}^2+\left(c_{3}^2\,c_{5}^3-c_{2}^2\,c_{5}\right)\,c_{6}
}\over{c_{3}^2\,c_{8}^4+\left(2\,c_{3}^2\,c_{5}^2-c_{2}^2\right)\,
c_{8}^2+c_{3}^2\,c_{5}^4-c_{2}^2\,c_{5}^2}}
\end{aligned}
\]
y1=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5)
y2=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5)
至于 \(AB\)、\(CD\) 同时为 \(0\) 的情况，只需计算 \(A、D\) 间距离即可。
code:
#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))
using namespace std;
long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;
struct node{
long double x,y;
}A,B,C,D,M,N,t;
inline long double dis(node a,node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline long double f(long double x,long double y)
{
return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}
signed main()
{
scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
if(c4*c8-c5*c7)
{
x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
}
else
{
if(C.x==D.x&&C.y==D.y)
{
x=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
ans=min(ans,f(x,y));
x=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
ans=min(ans,f(x,y));
}
else if(A.x==B.x&&A.y==B.y)
{
y=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
ans=min(ans,f(x,y));
y=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
ans=min(ans,f(x,y));
}
else if(C.x==D.x&&C.y==D.y&&A.x==B.x&&A.y==B.y) ans=min(ans,dis(A,D));
}
ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
printf("%.2Lf\n",ans);
return 0;
}
特判——「平行」
若只想到这里，你还是只能获得\(60pts\)。。。
还有最后一种情况：即是 \(AB、CD\) 不为 \(0\)，\(AB//CD\) 也会使 \(c_4\cdot c_8-c_5\cdot c_7=0\)「因为\((x_B-x_A)\cdot(y_D-y_C)=(x_D-x_C)\cdot(y_B-y_A)\)」
回到一开始的式子：
\[t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}
\]
似乎很难下手了……
或许我们可以换一种思路「乱搞」：充分利用几何知识
两条直线平行，我们可以直接考虑以下几种方案：
- \(A\rightarrow D\)
- 过 \(C\) 向 \(AB\) 作垂线，若与 \(AB\) 有交点 \(M\)，则：\(A\rightarrow M\rightarrow C\rightarrow D\)
- 过 \(D\) 向 \(AB\) 作垂线，若与 \(AB\) 有交点 \(M\)，则：\(A\rightarrow M\rightarrow D\)
- 过 \(B\) 向 \(CD\) 作垂线，若与 \(CD\) 有交点 \(N\)，则：\(A\rightarrow B\rightarrow N\rightarrow D\)
- 过 \(A\) 向 \(CD\) 作垂线，若与 \(CD\) 有交点 \(N\)，则：\(A\rightarrow N\rightarrow D\)
还有可能两条线段在同一条直线上，需要特判一下。
对这几种方案分别取最小值，得到正确答案的概率就大了许多。
code:
#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))
using namespace std;
long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;
struct node{
long double x,y;
}A,B,C,D,M,N,t;
inline long double dis(node a,node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline long double f(long double x,long double y)
{
return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}
inline bool eq_direction(node A1,node B1,node A2,node B2)//判断共线向量A1B1、A2B2是否同向
{
if(B1.x-A1.x&&B2.x-A2.x) return (B1.x-A1.x)*(B2.x-A2.x)>0;
else return (B1.y-A1.y)*(B2.y-A2.y)>0;
}
inline node foot_point(node P,node A,node B)//过P点作直线AB的垂线,返回垂足
{
node res;
//计算AB的解析式
long double a=B.y-A.y,b=A.x-B.x,c=-a*A.x-b*A.y;
res.x=(b*b*P.x-a*b*P.y-a*c)/(a*a+b*b);
res.y=(a*a*P.y-a*b*P.x-b*c)/(a*a+b*b);
return res;
}
inline bool on_line(node P,node A,node B)//判断在直线AB上的点P是否在线段AB上,只需判断P是否在A、B围成的矩形中即可
{
return min(A.x,B.x)<=P.x&&P.x<=max(A.x,B.x)&&min(A.y,B.y)<=P.y&&P.y<=max(A.y,B.y);
}
signed main()
{
scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
if(c4*c8-c5*c7)
{
x=(c5*c9-c6*c8)/(c4*c8-c5*c7),y=(c6*c7-c4*c9)/(c4*c8-c5*c7);
if(x>=0&&x<=1&&y>=0&&y<=1) ans=min(ans,f(x,y));
}
else
{
if(C.x==D.x&&C.y==D.y)
{
x=(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+((c1*c1-c3*c3*c4*c4)*c7-c3*c3*c7*c7*c7)*c9-c3*c3*c4*c6*c7*c7+(c1*c1*c4-c3*c3*c4*c4*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
ans=min(ans,f(x,y));
x=-(sqrt(c3*c3*c7*c7+c3*c3*c4*c4-c1*c1)*(c1*c4*c9-c1*c6*c7)+(c3*c3*c7*c7*c7+(c3*c3*c4*c4-c1*c1)*c7)*c9+c3*c3*c4*c6*c7*c7+(c3*c3*c4*c4*c4-c1*c1*c4)*c6)/(c3*c3*c7*c7*c7*c7+(2*c3*c3*c4*c4-c1*c1)*c7*c7+c3*c3*c4*c4*c4*c4-c1*c1*c4*c4);
ans=min(ans,f(x,y));
}
else if(A.x==B.x&&A.y==B.y)
{
y=(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+((c2*c2-c3*c3*c5*c5)*c8-c3*c3*c8*c8*c8)*c9-c3*c3*c5*c6*c8*c8+(c2*c2*c5-c3*c3*c5*c5*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
ans=min(ans,f(x,y));
y=-(sqrt(c3*c3*c8*c8+c3*c3*c5*c5-c2*c2)*(c2*c5*c9-c2*c6*c8)+(c3*c3*c8*c8*c8+(c3*c3*c5*c5-c2*c2)*c8)*c9+c3*c3*c5*c6*c8*c8+(c3*c3*c5*c5*c5-c2*c2*c5)*c6)/(c3*c3*c8*c8*c8*c8+(2*c3*c3*c5*c5-c2*c2)*c8*c8+c3*c3*c5*c5*c5*c5-c2*c2*c5*c5);
ans=min(ans,f(x,y));
}
//else if(C.x==D.x&&C.y==D.y&&A.x==B.x&&A.y==B.y) ans=min(ans,dis(A,D));
else
{
ans=min(ans,dis(A,D)/R);
if(c4*(C.y-B.y)==c7*(C.x-B.x))//已知AB//CD,若AB//BC,则A、B、C、D四点共线
{
ans=min(ans,dis(A,B)/P+dis(B,C)/R+dis(C,D)/Q);//ABCD
ans=min(ans,dis(A,B)/P+dis(B,D)/R);//ABD
ans=min(ans,dis(A,C)/R+dis(C,D)/Q);//ACD
if(!eq_direction(A,C,B,C)&&!eq_direction(A,D,B,D))
{
if(eq_direction(A,C,C,D)) ans=min(ans,dis(A,C)/P+dis(C,D)/max(P,Q));//ACDB
else ans=min(ans,dis(A,D)/P);//ADCB
}
if(!eq_direction(A,C,A,D)&&!eq_direction(B,C,B,D))
{
if(eq_direction(C,A,A,B)) ans=min(ans,dis(A,B)/max(P,Q)+dis(B,D)/Q);//CABD
else ans=min(ans,dis(A,D)/Q);//CBAD
}
}
else//不共线
{
t=foot_point(C,A,B);
if(on_line(t,A,B)) ans=min(ans,dis(A,t)/P+dis(t,C)/R+dis(C,D)/Q);//AMCD
t=foot_point(D,A,B);
if(on_line(t,A,B)) ans=min(ans,dis(A,t)/P+dis(t,D)/R);//AMD
t=foot_point(B,C,D);
if(on_line(t,C,D)) ans=min(ans,dis(A,B)/P+dis(B,t)/R+dis(t,D)/Q);//ABND
t=foot_point(A,C,D);
if(on_line(t,C,D)) ans=min(ans,dis(A,t)/R+dis(t,D)/Q);//AND
}
}
}
ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
printf("%.2Lf\n",ans);
return 0;
}
瓶颈
emmm...
仍是\(60pts...\)
乱搞终究不是正道。。。
不难发现，\(O(1)\)做法的瓶颈卡在了\(AB//CD\)使得求导后分母为\(0\)的情况。
若有大佬可以通过调整参数、用更好的方法求最值的方式来绕过这一特殊情况，请不吝赐教～
梯度下降法
\(O(1)\)：我死了QAQ……
再次回到起点：
\[t=f(x,y)=c_1\cdot x+c_2\cdot y+c_3\cdot\sqrt{(c_4\cdot x+c_5\cdot y+c_6)^2+(c_7\cdot x+c_8\cdot y+c_9)^2}
\]
非 \(O(1)\) 复杂度求解多元函数最值的方法有很多：粒子群算法、模拟退火、三分套三分、牛顿迭代法……
在此介绍梯度下降法。
梯度
了解多元微积分的各位大佬们都知道，梯度是一个向量，指向多元函数的值变化最快的方向，大小即为变化率。
「不了解的可以看看 3Blue1Brown 的作者 Grant Sanderson 讲的多元微积分课程，在b站上已有熟肉【链接】」
百度百科对梯度的定义：
设二元函数\(z=f(x,y)\) 在平面区域\(D\)上具有一阶连续偏导数，则对于每一个点\(P(x,y)\)都可定出一个向量
\[\left\{ \frac{\partial f}{\partial x},\frac{\partial f}{\partial y} \right\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j} \]该
函数「也不是不对，但称为向量更好理解」就称为函数\(z=f(x,y)\)在点\(P(x,y)\)的梯度，记作\(\text{grad}f(x,y)\)或\(\nabla f(x,y)\)。于是有：
\[\text{grad}f(x,y)=\nabla f(x,y)=\left\{ \frac{\partial f}{\partial x},\frac{\partial f}{\partial y} \right\}=f_x(x,y)\hat{i}+f_y(x,y)\hat{j} \]其中
\[\nabla=\frac{\partial}{\partial x}\hat{i}+\frac{\partial}{\partial y}\hat{j} \]称为（二维的）向量微分算子或Nabla算子，且有
\[\nabla f=\frac{\partial f}{\partial x}\hat{i}+\frac{\partial f}{\partial y}\hat{j} \]
梯度下降法的主要思想
再来介绍梯度下降法。
正如我们先前对 \(f(x,y)\) 求偏导，令偏导值为 \(0\)，解出所有驻点的操作，实际上就是试图直接令梯度值为 \(0\)。
但是也正如我们遇到的问题，有时直接求解比较困难，我们不妨通过迭代，让梯度值逐渐下降，即逐渐接近极值点。
例如，对于函数 \(f(x)=x^2\)，求导得 \(f'(x)=2x\)。
不妨设 \(x_0=1\)，则从点 \((1,1)\) 开始，计算其梯度值：\(f'(1)=2>0\)。
梯度值 \(>0\)，说明在该点沿使 \(x\) 增大的方向「即 \(x\) 轴正半轴方向」，\(f(x)\) 函数值会增大，增大的速率为 \(2\)。
容易发现，在 \((3,9)\) 处，\(f'(3)=6\)，增大的速率为 \(6\)，观察图像也可得知，沿 \(x\) 轴正半轴方向，\(f(x)\) 的图像变陡，上升速率变快。
而在 \((0,0)\) 处，\(f'(0)=0\)，函数图像在这一点的切线斜率是平的，增大的速率无限接近 \(0\)，此时梯度值为 \(0\)，\((0,0)\) 是函数 \(f(x)\) 的一个驻点，同时 \(f(x)\) 也在此取到最小值。
因此，我们有如下策略：
梯度值 \(>0\)，向左移动一点；梯度值 \(<0\)，向右移动一点。
即每次让 \(x\) 朝与梯度值符号相反的方向移动，使梯度值逐渐下降，最终趋于 \(0\)。
用数学语言来描述，即：
\[x_{k+1}=x_k-\nabla f(x_k)
\]
其中 \(x_k\) 为第 \(k\) 次迭代时点的横坐标，\(x_{k+1}\) 为第 \(k+1\) 次迭代移动到的点的横坐标，\(x_0\) 表示初始横坐标。
\(\nabla f(x)\) 表示在 \(x\) 处的梯度值，\(\nabla f(x)=\frac{\text df}{\text dx}\)。
\(x_k-\nabla f(x_k)\) 表示向与梯度值符号相反的方向移动，符合之前的策略
优化——「学习率」
但这样就有一个问题：
例如 \(f(x)=x^2,x_0=10\)，则 \(\nabla f(x)=2x\)：
\[\begin{aligned}
x_1&=x_0-2\cdot x_0=-10
\\
x_2&=x_1-2\cdot x_1=10
\\
x_3&=x_2-2\cdot x_2=-10
\\
...
\end{aligned}
\]
可以看到，虽然横坐标一直在变化，但一直在 \(10、-10\) 之间振荡，函数值始终没有降到最低点。
为此，我们需要一个参数控制移动的距离，这个参数被称作学习率，用 \(\eta\) 表示：
\[x_{k+1}=x_k-\eta\cdot\nabla f(x_k)
\]
显然，在刚刚的计算中，\(\eta=1\)。
当 \(\eta>1\) 时，函数值不仅不会降到最低点，甚至会越来越大：
\(e.g.\qquad\eta=1.05\)
\[\begin{aligned}
x_1&=x_0-2\cdot x_0=-11
\\
x_2&=x_1-2\cdot x_1=12.1
\\
x_3&=x_2-2\cdot x_2=-13.31
\\
...
\end{aligned}
\]
而若将 \(\eta\) 调低到一个较小值，如 \(0.02\)：
\[x_1=x_0-2\cdot x_0=9.6
\\
x_2=x_1-2\cdot x_1=9.2
\\
x_3=x_2-2\cdot x_2=8.8
\\
...
\]
降低的速度太慢，容易超时。
当 \(\eta=0.2\) 时，迭代 \(10\) 次左右即可降入谷底：
\[\begin{aligned}
x_1&=x_0-2\cdot x_0=6&||\nabla f(x_0)||=12
\\
x_2&=x_1-2\cdot x_1=3.6&||\nabla f(x_1)||=7.2
\\
x_3&=x_2-2\cdot x_2=2.16&||\nabla f(x_2)||=4.32
\\
...
\end{aligned}
\]
可以看到，右边的梯度值在每次迭代后都会下降，故称为梯度下降法「Gradient Descent」
只要选择合适的学习率，梯度就可以下降到任意小：
\[\lim_{k\to\infty}f(x_k)=\min f(x)
\]
可以用泰勒公式进行严格证明，此处不再赘述（逃
因此可以通过梯度值的大小作为终止条件「也可以直接用迭代次数控制精度」
对于二元函数，同样可以用梯度下降法求解极值：
\[f(x_{k+1},y_{k+1})=f(x_k,y_k)-\eta\cdot\nabla f(x_k,y_k)
\]
\(e.g.\qquad f(x,y)=x^2+2y^2,(x_0,y_0)=(-3.5,-3.5),\eta=0.1\)，则\(\nabla f(x,y)=(2x,4y)\)
\[(x_1,y_1)=(x_0,y_0)-\eta\cdot\nabla f(x_0,y_0)=(-2.8,-2.1)
\\
(x_2,y_2)=(x_1,y_1)-\eta\cdot\nabla f(x_1,y_1)=(-2.24,-1.26)
\\
...
\]
代码——「梯度下降法」
code：
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
inline long double f(long double x)
{
return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}
inline long double numerical_diff(long double x)//数值微分法估计一阶导数
{
long double dx=1e-6;
return (f(x+dx)-f(x-dx))/(dx*2);
}
inline long double gradient_descent(long double x,long double eta)//梯度下降法
{
long double eps=1e-6;
//int cnt=0;
while(abs(numerical_diff(x))>eps)
{
x-=eta*numerical_diff(x);
//cnt++;
}
//cout<<cnt<<endl;
return f(x);
}
signed main()
{
cout<<gradient_descent(0,0.02)<<endl;
return 0;
}
\(e.g.\qquad f(x)=3x^4-x^3+2x^2-9x+5\sqrt{(x+3)^2+(5x+6)^2}-25\)
经过\(9\)次迭代达到目标精度。
优化——动量梯度下降法「MGD」
需要指出的是，梯度值为 \(0\) 的驻点不一定是函数的极值点，如 \(f(x)=x^3\) 在 \(x=0\) 处梯度值为 \(0\)，但并不是函数的极值点：
「实际上，\((0,0)\) 是 \(f(x)\) 的拐点，此处函数的凹凸性发生改变」
同时，对于非单峰函数来说，梯度下降法的结果易受到初始值的影响，也就是陷入局部最优解。
现实生活中，一个小球从高处落下，大概率会越过比较低的坎继续下降。
我们同样可以引入惯性来优化：
用梯度模拟受力，使之不直接控制移动距离，而是给小球一个速度 \(v\)；同时可以引入阻力，也就是速度衰减率 \(\beta\) 使小球减速。
于是，整个过程就像一个有动量的小球在空间中来回滚动，故称动量梯度下降法「Gradient Descent with Momentum，简称MGD」
根据动量定理
\[F\Delta t=m\Delta v
\]
移项，得：
\[\Delta v=\frac{\Delta t}{m}\cdot F
\]
直接令 \(\eta=\frac{\Delta t}{m}\)，则：
\[\Delta v=\eta\cdot F
\]
但在实际编程中，不必那么精细地刻画阻力，直接用一个速度衰减率 \(\beta\) 代替即可：
\[v_{k+1}=\beta\cdot v_k-\eta\cdot\nabla f(x_k)
\]
code：
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
inline long double f(long double x)
{
return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}
inline long double numerical_diff(long double x)
{
long double dx=1e-6;
return (f(x+dx)-f(x-dx))/(dx*2);
}
inline long double mgd(long double x,long double eta,long double beta)//动量梯度下降法
{
long double v=0,eps=1e-6;
int cnt=0;
while(abs(numerical_diff(x))>eps)
{
v=beta*v-eta*numerical_diff(x);
x+=v;
cnt++;
}
cout<<cnt<<endl;
return f(x);
}
signed main()
{
cout<<mgd(0,0.03,0.05)<<endl;
return 0;
}
迭代次数：\(13\)。
「多峰函数的优化效果较明显」
优化——自适应梯度算法「AdaGrad」
对比一下两种算法的真实移动距离 \(\Delta x_i\)：
| 梯度下降法 | 动量梯度下降法 | |
|---|---|---|
| \(\Delta x_i\) | \(\eta\cdot\frac{\partial f}{\partial x_i}\) | \(\eta\cdot v_i\) |
动量梯度下降法只优化了 \(\frac{\partial f}{\partial x_i}\) 的部分。
很自然的想到，是否可以优化学习率 \(\eta\)？
答案是肯定的，这种算法称为自适应梯度算法「AdaGrad【Adapt+Gradient】」，基本思想是开始时快速移动接近目标，然后减速提高精度「梯度对其产生的动力效果逐渐减弱」。实现时记录每次梯度下降的梯度平方和 \(h\)，将 \(\eta\) 除以 \(\sqrt{h}\) 进行调整「为防止分母出现 \(0\)，往往加上一个较小的常量」
code：
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
inline long double f(long double x)
{
return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}
inline long double numerical_diff(long double x)
{
long double dx=1e-6;
return (f(x+dx)-f(x-dx))/(dx*2);
}
inline long double AdaGrad(long double x,long double eta,long double beta)//自适应梯度算法
{
long double v=0,h=0,grad=numerical_diff(x),eps=1e-6;
//int cnt=0;
while(abs(grad)>eps)
{
h+=grad*grad;
v=beta*v-eta*grad/(sqrt(h)+eps);
x+=v;
grad=numerical_diff(x);
//cnt++;
}
//cout<<cnt<<endl;
return f(x);
}
signed main()
{
cout<<AdaGrad(0,0.5,0.05)<<endl;
return 0;
}
迭代次数：\(12\)。
「多峰函数的优化效果较明显」
优化——自适应动量算法「Adam」
自适应梯度算法的缺陷：若初始点离极值点较远，可能还没到极值点，小球就跑不动了「相应的有 RMSProp 优化算法，通过把 \(h\) 乘以一个衰减率 \(\beta_2\) 来逐步遗忘之前的梯度，专业一点说是“指数移动平均”，与动量梯度下降法有异曲同工之妙，此处不作详细介绍」
我们也可以将动量梯度下降法「MGD」与自适应梯度算法「AdaGrad」融合在一起，得到自适应动量算法「Adaptive Momentum Estimation，简称 Adam」
在梯度下降过程中，Adam 算法用梯度模拟小球受力改变速度 \(v\)，同时会增大小球的质量而改变小球移动的难易程度，二者的作用最终影响移动距离。
code：
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
inline long double f(long double x)
{
return 3*x*x*x*x-x*x*x+2*x*x-9*x+5*sqrt((x+3)*(x+3)+(5*x+6)*(5*x+6))-25;
}
inline long double numerical_diff(long double x)
{
long double dx=1e-6;
return (f(x+dx)-f(x-dx))/(dx*2);
}
inline long double Adam(long double x,long double eta,long double beta1,long double beta2)
{
long double v=0,m=0,grad=numerical_diff(x),t,t1=beta1,t2=beta2,eps=1e-6;
//int cnt=0;
while(abs(grad)>eps)
{
t=eta*sqrt(1-t2)/(1-t1);
v+=(1-beta1)*(grad-v);
m+=(1-beta2)*(grad*grad-m);
x-=t*v/(sqrt(m)+eps);
t1*=beta1,t2*=beta2;
grad=numerical_diff(x);
//cnt++;
}
//cout<<cnt<<endl;
return f(x);
}
signed main()
{
cout<<Adam(0,0.5,0.6,0.9999)<<endl;
return 0;
}
迭代次数：\(49\)。
综合来看，AdaGrad 的算法似乎更具优势，但毕竟各个算法特点不同，适用的函数也不同，用哪种算法还是应根据实际应用而定。
二维写法
code：
struct node{
long double x,y;
};
inline long double f(long double x,long double y)
{
return ...;
}
inline long double part_x(long double x,long double y)
{
long double dx=1e-6;
return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}
inline long double part_y(long double x,long double y)
{
long double dy=1e-6;
return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}
inline node gd(long double x,long double y,long double eta,long double beta)
{
long double eps=1e-6;
node res;
int cnt=0,T=100;
while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
{
x-=eta*part_x(x,y);
y-=eta*part_y(x,y);
//cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
}
res.x=x,res.y=y;
return res;
}
inline node mgd(long double x,long double y,long double eta,long double beta)
{
long double vx=0,vy=0,eps=1e-6;
node res;
int cnt=0,T=100;
while(T--&&abs(part_x(x,y))+abs(part_y(x,y))>eps)
{
vx=beta*vx-eta*part_x(x,y);
vy=beta*vy-eta*part_y(x,y);
x+=vx,y+=vy;
//cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
}
res.x=x,res.y=y;
return res;
}
inline node AdaGrad(long double x,long double y,long double eta,long double beta)
{
long double vx=0,vy=0,hx=0,hy=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),eps=1e-6;
node res;
int cnt=0,T=100;
while(T--&&grad>1e-5)
{
hx+=gradx*gradx;
hy+=grady*grady;
vx=beta*vx-eta*gradx/(sqrt(hx)+eps);
vy=beta*vy-eta*grady/(sqrt(hy)+eps);
x+=vx,y+=vy;
gradx=part_x(x,y);
grady=part_y(x,y);
grad=abs(gradx)+abs(grady);
//cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
}
res.x=x,res.y=y;
return res;
}
inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-7;
node res;
int cnt=0,T=100;
while(T--&&grad>1e-5)
{
t=eta*sqrt(1-t2)/(1-t1);
vx+=(1-beta1)*(gradx-vx);
vy+=(1-beta1)*(grady-vy);
mx+=(1-beta2)*(gradx*gradx-mx);
my+=(1-beta2)*(grady*grady-my);
x-=t*vx/(sqrt(mx)+eps);
y-=t*vy/(sqrt(my)+eps);
t1*=beta1,t2*=beta2;
gradx=part_x(x,y);
grady=part_y(x,y);
grad=abs(gradx)+abs(grady);
//cout<<x<<' '<<y<<' '<<f(x,y)<<' '<<++cnt<<endl;
}
res.x=x,res.y=y;
return res;
}
本题代码
回到此题，直接用 Adam 交了一发，AC～「可能运气较好……」
code：
#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#include<cstdlib>
#include<ctime>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))
using namespace std;
long double P,Q,R,ans=INT_MAX;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;
struct node{
long double x,y;
}A,B,C,D,M,N,t;
inline long double dis(node a,node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline long double f(long double x,long double y)
{
return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}
inline long double part_x(long double x,long double y)
{
long double dx=1e-6;
return (f(x+dx,y)-f(x-dx,y))/(dx*2);
}
inline long double part_y(long double x,long double y)
{
long double dy=1e-6;
return (f(x,y+dy)-f(x,y-dy))/(dy*2);
}
inline node Adam(long double x,long double y,long double eta,long double beta1,long double beta2)
{
long double vx=0,vy=0,mx=0,my=0,gradx=part_x(x,y),grady=part_y(x,y),grad=abs(gradx)+abs(grady),t,t1=beta1,t2=beta2,eps=1e-6;
node res;
int T=20;
while(T--&&grad>1e-4)
{
t=eta*sqrt(1-t2)/(1-t1);
vx+=(1-beta1)*(gradx-vx);
vy+=(1-beta1)*(grady-vy);
mx+=(1-beta2)*(gradx*gradx-mx);
my+=(1-beta2)*(grady*grady-my);
x-=t*vx/(sqrt(mx)+eps);
y-=t*vy/(sqrt(my)+eps);
t1*=beta1,t2*=beta2;
gradx=part_x(x,y);
grady=part_y(x,y);
grad=abs(gradx)+abs(grady);
}
res.x=x,res.y=y;
return res;
}
inline long double Rand()
{
return abs(rand()*rand()*rand()%(10000000)/(long double)(10000000));
}
inline void gradient()
{
for(int T=1;T<=200000;T++)
{
t=Adam(Rand(),Rand(),Rand(),Rand(),Rand());
if(t.x>=0&&t.x<=1&&t.y>=0&&t.y<=1) ans=min(ans,f(t.x,t.y));
}
}
signed main()
{
srand(time(0));
scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
gradient();
ans=Min(f(0,0),f(0,1),f(1,0),f(1,1),ans);
printf("%.2Lf\n",ans);
return 0;
}
三分套三分法
突然发现三分套三分跑的贼快，于是也打算简单实现一下。
模板就不讲了，直接贴一下之前的博客：「二分法&三分法模板」
上代码：
#include<iostream>
#include<cstdio>
#include<cmath>
#include<climits>
#define Min(a,b,c,d,e) (min(a,min(b,min(c,min(d,e)))))
using namespace std;
long double P,Q,R,ans=INT_MAX,eps=1e-6;
long double c1,c2,c3,c4,c5,c6,c7,c8,c9;
long double x,y;
struct node{
long double x,y;
}A,B,C,D,M,N,t;
inline long double dis(node a,node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline long double f(long double x,long double y)
{
return c1*x+c2*y+c3*sqrt((c4*x+c5*y+c6)*(c4*x+c5*y+c6)+(c7*x+c8*y+c9)*(c7*x+c8*y+c9));
}
inline long double g(long double x)
{
long double l=0,r=1,mid;
while(r-l>eps)
{
mid=(l+r)/2;
if(f(x,mid-eps)<f(x,mid+eps)) r=mid;
else l=mid;
}
return f(x,mid);
}
signed main()
{
scanf("%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf%Lf",&A.x,&A.y,&B.x,&B.y,&C.x,&C.y,&D.x,&D.y,&P,&Q,&R);
c1=dis(A,B)/P,c2=dis(C,D)/Q,c3=1/R;
c4=B.x-A.x,c5=D.x-C.x,c6=A.x-D.x,c7=B.y-A.y,c8=D.y-C.y,c9=A.y-D.y;
long double l=0,r=1,mid;
while(r-l>eps)
{
mid=(l+r)/2;
if(g(mid-eps)<g(mid+eps)) r=mid;
else l=mid;
}
printf("%.2Lf\n",g(mid));
return 0;
}
$$ Thanks\quad for\quad reading. $$
————THE——END————
