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

题面

传送门

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

orz Shinbokuow

前置芝士

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

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

[egin{aligned} Ans &=sum_{i=1}^n{(kx_i-y_i+b)^2over k^2+1}\ &=sum_{i=1}^n{k^2x_i^2+y_i^2+b^2-2kx_iy_i+2kbx_i-2by_iover k^2+1}\ &={1over k^2+1}left(k^2sum_{i=1}^nx_i^2+sum_{i=1}^ny_i^2+nb^2-2ksum_{i=1}^nx_iy_i+2kbsum_{i=1}^nx_i-2bsum_{i=1}^ny_i ight) end{aligned} ]

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

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

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

也就是说

[{partial fover partial b}={1over k^2+1}left(2nb+2ksum_{i=1}^nx_i-2sum_{i=1}^ny_i ight)=0 ]

如果我们令

[overline{x}={1over n}sum_{i=1}^nx_i,overline{y}={1over n}sum_{i=1}^ny_i ]

那么(b)就可以表示为

[b=overline{y}-koverline{x} ]

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

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

其中

[A=sum x_i^2-{(sum x_i)^2over n} ]

[B={2sum x_isum y_iover n}-2sum x_iy_i ]

[C=sum y_i^2-{(sum y_i)^2over 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 ]

整理之后有

[-4alpha^2+4(A+C)alpha+B^2-4ACgeq 0 ]

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

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

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

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

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

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

所以肯定有解啦！

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

题解

然后是题解了……

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

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

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

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

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

复杂度(O(n^3log n))，其中(O(nlog 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
",res/n);
	return 0;
}