CodeChef TWOROADS（計算幾何+拉格朗日乘數法）

題面

傳送門

簡要題意：給出\(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;
}