CodeChef TWOROADS（计算几何+拉格朗日乘数法）

谁践踏了优雅 2022-01-12 09:23 172阅读 0赞

## 题面 ##

[传送门][Link 1]

简要题意：给出\\(n\\)个点，请求出两条直线，并最小化每个点到离它最近的那条直线的距离的平方和，\\(n\\leq 100\\)

orz [Shinbokuow][]

## 前置芝士 ##

给出\\(n\\)个点，请求出一条直线，使所有点到它距离的平方和最小，点带插入和删除

如果我们设\\(y=kx+b\\)，设点\\(i\\)为\\((x\_i,y\_i)\\)，那么就是要求我们最小化

\\\[ \\begin\{aligned\} Ans &=\\sum\_\{i=1\}^n\{(kx\_i-y\_i+b)^2\\over k^2+1\}\\\\ &=\\sum\_\{i=1\}^n\{k^2x\_i^2+y\_i^2+b^2-2kx\_iy\_i+2kbx\_i-2by\_i\\over k^2+1\}\\\\ &=\{1\\over k^2+1\}\\left(k^2\\sum\_\{i=1\}^nx\_i^2+\\sum\_\{i=1\}^ny\_i^2+nb^2-2k\\sum\_\{i=1\}^nx\_iy\_i+2kb\\sum\_\{i=1\}^nx\_i-2b\\sum\_\{i=1\}^ny\_i\\right) \\end\{aligned\} \\\]

我们令这个柿子为\\(f(k,b)\\)，也就是说这是一个以\\(k,b\\)为自变量的函数

接下来，我们用拉格朗日乘数法对\\(b\\)求偏导数

上面这话说的简直不是人话，翻译一下的话，可以理解为，我们假设\\(k\\)是一个定值\\(k\_0\\)，并认为\\(b\\)是变量，对它求导。导数为\\(0\\)的点就是当\\(k=k\_0\\)时的极大或极小值（因为这里函数可以无限大所以肯定是极小值），那么显然对于任意一个\\(k\\)都有一个最优的\\(b\\)，且\\(b\\)是一个关于\\(k\\)的一次函数，我们把\\(b\\)代入就行了

也就是说

\\\[\{\\partial f\\over \\partial b\}=\{1\\over k^2+1\}\\left(2nb+2k\\sum\_\{i=1\}^nx\_i-2\\sum\_\{i=1\}^ny\_i\\right)=0\\\]

如果我们令

\\\[\\overline\{x\}=\{1\\over n\}\\sum\_\{i=1\}^nx\_i,\\overline\{y\}=\{1\\over n\}\\sum\_\{i=1\}^ny\_i\\\]

那么\\(b\\)就可以表示为

\\\[b=\\overline\{y\}-k\\overline\{x\}\\\]

然后我们把代入原式，可以解得

\\\[f(k,b)=\{Ak^2+Bk+C\\over k^2+1\}\\\]

其中

\\\[A=\\sum x\_i^2-\{(\\sum x\_i)^2\\over n\}\\\]

\\\[B=\{2\\sum x\_i\\sum y\_i\\over n\}-2\\sum x\_iy\_i\\\]

\\\[C=\\sum y\_i^2-\{(\\sum y\_i)^2\\over n\}\\\]

然后我们要求\\(f(k,b)\\)的最小值，移项之后可以得到

\\\[(A-f(k,b))k^2+Bk+C-f(k,b)=0\\\]

以下我们设\\(\\alpha=f(k,b)\\)。因为这里需要保证\\(k\\)有解，也就是说

\\\[B^2-4(A-\\alpha)(C-\\alpha)\\geq 0\\\]

整理之后有

\\\[-4\\alpha^2+4(A+C)\\alpha+B^2-4AC\\geq 0\\\]

数形结合一下发现左边是个开口向下的二次函数，与\\(x\\)轴有一个或两个交点，那么我们取小一点的那个就好了

等会儿？万一这个二次函数最大值小于\\(0\\)呢？那不就无解了么？

实际上是不会的，可以用两种方法考虑

感性理解：想一想它代表的意义，再看看我手里的锤子，您说有没有解？

数学方法：因为我们考虑二次函数\\(ax^2+bx+c=0\\)，它的最值为\\(\{4ac-b^2\\over 4a\}\\)，那么代入之后可以发现哪个二次函数的最大值为

\\\[ \\begin\{aligned\} \{4\\times -4\\times (B^2-4AC)-16(A+C)^2\\over -16\} &=B^2-4AC+(A+C)^2\\\\ &=B^2+(A-C)^2\\geq 0 \\end\{aligned\} \\\]

所以肯定有解啦！

于是我们就可以做到\\(O(1)\\)插入，\\(O(1)\\)删除了

## 题解 ##

然后是题解了……

对于两条直线，它们两个的两条垂直的角平分线会把平面分成四个部分，其中两个相对的区域离一个平面比较近，另外两个相对的区域离另一条平面比较近

但是枚举角平分线这件事情就会变得非常辣手……

我们可以枚举点对\\(a,b\\)，并令直线\\(ab\\)为第一条角平分线，钦定\\(a\\)在其中一端，\\(b\\)在另一端，那么不难发现我们这样枚举其实等价于枚举完了所有的角平分线

然后第一条角平分线就解决了，接下来是和它垂直的第二条

我们计算出每个点在第一条角平分线上的投影长度，然后按投影长度排个序，那么第二条角平分线对平面的分割只有\\(O(n)\\)种情况，直接找过去就可以了

复杂度\\(O(n^3\\log n)\\)，其中\\(O(n\\log n)\\)是排序的复杂度

//minamoto #include<bits/stdc++.h> #define R register #define inline __inline__ __attribute__((always_inline)) #define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i) #define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i) #define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v) template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;} using namespace std; const int N=105;const double eps=1e-10; inline int sgn(R double x){return x<-eps?-1:x>eps;} struct node{ double x,y; inline node(){} inline node(R double xx,R double yy):x(xx),y(yy){} inline node operator +(const node &b)const{return node(x+b.x,y+b.y);} inline node operator -(const node &b)const{return node(x-b.x,y-b.y);} inline double operator *(const node &b)const{return x*b.y-y*b.x;} inline node operator *(const double &b)const{return node(x*b,y*b);} inline double operator ^(const node &b)const{return x*b.x+y*b.y;} }p[N],v; struct qwq{ node p;double d;bool in; inline qwq(){} inline qwq(R node pp,R double dd,R bool ii):p(pp),d(dd),in(ii){} inline bool operator <(const qwq &b)const{return d<b.d;} }st[N]; struct Line{ double x,y,xx,yy,xy;int sz; inline void ins(R node p){x+=p.x,y+=p.y,xx+=p.x*p.x,yy+=p.y*p.y,xy+=p.x*p.y,++sz;}; inline void del(R node p){x-=p.x,y-=p.y,xx-=p.x*p.x,yy-=p.y*p.y,xy-=p.x*p.y,--sz;}; inline void clr(){x=y=xx=yy=xy=sz=0;} double calc(){ if(!sz)return 0; double xa=x/sz,ya=y/sz,A=xx-sz*xa*xa; double B=2*xa*ya*sz-2*xy,C=yy-sz*ya*ya; double a=4,b=-4*(A+C),c=4*A*C-B*B; return (-b-sqrt(b*b-4*a*c))/(a*2); } }l1,l2; double res=1e18;int n; int main(){ // freopen("testdata.in","r",stdin); scanf("%d",&n); fp(i,1,n)scanf("%lf%lf",&p[i].x,&p[i].y); fp(a,1,n)fp(b,1,n)if(a!=b){ v=p[b]-p[a]; fp(i,1,n)st[i]=qwq(p[i],v^(p[i]-p[a]),v*(p[i]-p[a])>eps); st[a].in=1; sort(st+1,st+1+n),l1.clr(),l2.clr(); fp(i,1,n)st[i].in?l1.ins(st[i].p):l2.ins(st[i].p); cmin(res,l1.calc()+l2.calc()); fp(i,1,n){ if(st[i].in)l1.del(st[i].p),l2.ins(st[i].p); else l2.del(st[i].p),l1.ins(st[i].p); cmin(res,l1.calc()+l2.calc()); } } printf("%.10lf\n",res/n); return 0; }

转载于:https://www.cnblogs.com/bztMinamoto/p/10700369.html

[Link 1]: https://www.codechef.com/problems/TWOROADS
[Shinbokuow]: http://wyfcyx.is-programmer.com/posts/190179.html