BZOJ3775: 點和直線（計算幾何+拉格朗日乘數法）

題面

傳送門

題解

勁啊……

沒有和\(Claris\)同樣推，用了相似於\(Shinbokuow\)推已知點求最短直線的方法，結果\(WA\)了好幾個小時，拿\(Claris\)代碼拍了幾個小時都沒找到\(bug\)在哪兒，最後發現是我一個除法的地方忘記除數爲\(0\)的狀況了……甘霖娘……

公式恐懼症患者能夠直接轉去結論了

設直線爲\(ax+by+c=0\)，點爲\((x,y)\)，記\(d_i=a_i^2+b_i^2\)，那麼就是要咱們最小化

\[ \begin{aligned} f(x,y) &=\sum {(a_ix+b_iy+c_i)^2\over d_i}\\ &=\sum {a_i^2x^2+b_i^2y^2+c_i^2+2a_ib_ixy+2a_ic_ix+2b_ic_iy\over d_i} \end{aligned} \]

如下爲了方便，記\(A^2=\sum{a_i^2\over d_i}\)，\(B^2,C^2\)同理，以及\(AB=\sum{a_ib_i\over d_i}\)，\(BC,AC\)同理，那麼原式能夠表示成

\[f(x,y)=x^2A^2+y^2B^2+C^2+2xyAB+2xAC+2yBC\]

用拉格朗日乘數法對\(y\)求偏導數（這句話的意思大概就是，咱們認爲\(x\)是一個常數，那麼對於每個\(x=x_0\)，\(y\)都會有一個極值點，而這個極值點就是它導數爲\(0\)的點，因此咱們把\(y\)看作變量求導）

\[{\partial f\over \partial y}=2yB^2+2xAB+2BC=0\]

解得

\[y={-xAB-BC\over B^2}\]

代入原式能夠化爲

\[f(x,y)=\alpha x^2+\beta x+\gamma\]

其中

\[\alpha=A^2-{(AB)^2\over B^2}\]

\[\beta=2AC-{2(AB)(BC)\over B^2}\]

\[\gamma=C^2-{(BC)^2\over B^2}\]

易知\(\alpha \geq 0\)（證實下面有）

不過這裏其實還有一個尷尬的狀況就是有可能\(B^2=0\)，也就是說全部直線的\(b_i=0\)，不過咱們轉過頭去看會發現這種狀況下\(y\)對\(f(x,y)\)徹底沒有影響，並且\(\alpha,\beta,\gamma\)的值分別就是\(A^2,2AC,C^2\)。因此這種狀況其實並不會有影響

若是\(\alpha\neq 0\)，咱們要最小化\(f(x,y)\)，同時還須要知足方程

\[\alpha x^2+\beta x+\gamma-f(x,y)=0\]

有解

代入根的判別式，可知須要知足

\[\beta^2-4\alpha(\gamma-f(x,y))\geq 0\]

\[f(x,y)\geq \gamma-{\beta^2\over 4\alpha}\]

最小值顯然了

若是\(\alpha=0\)，則

\[f(x,y)=\beta x+\gamma\]

\(Claris\)說這種狀況下答案就等於\(\gamma\)……然而我實在看不出爲啥……我怎麼感受能夠無限小呢……然而它要是變成負數顯然不符合常理啊……有哪位鴿鴿知道爲何的麼能夠在下面留言哦qwq

而後就作完了

ps：關於\(\alpha\geq 0\)的證實

由於有

\[\alpha=A^2-{(AB)^2\over B^2}\]

首先顯然\(B^2\geq 0\)，若是\(B^2=0\)，那麼根據上面所說\(\alpha=A^2\geq 0\)，因此假設\(B^2>0\)，咱們須要證實

\[A^2-{(AB)^2\over B^2}\geq 0\]

即

\[A^2B^2\geq (AB)^2\]

代入原來的值

\[\left(\sum {a_i^2\over d_i}\right)\left(\sum {b_i^2\over d_i}\right)\geq \left(\sum {a_ib_i\over d_i}\right)^2\]

這就是柯西不等式啊……顯然成立

而後沒有而後了

//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)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
int read(char *s){
    R int len=0;R char ch;while(((ch=getc())>'9'||ch<'0'));
    for(s[++len]=ch;(ch=getc())>='0'&&ch<='9';s[++len]=ch);
    return s[len+1]='\0',len;
}
double readdb()
{
    R double x=0,y=0.1,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(x=ch-'0';(ch=getc())>='0'&&ch<='9';x=x*10+ch-'0');
    for(ch=='.'&&(ch=getc());ch>='0'&&ch<='9';x+=(ch-'0')*y,y*=0.1,ch=getc());
    return x*f;
}
inline int getop(){R char ch;while((ch=getc())>'9'||ch<'0');return ch-'0';}
const int N=2e5+5;const double eps=1e-7;
inline int sgn(R double x){return x<-eps?-1:x>eps;}
struct node{
    double aa,bb,cc,ab,bc,ac;int sz;
    inline void ins(R double a,R double b,R double c,R double d){
        ++sz,aa+=a*a*d,bb+=b*b*d,cc+=c*c*d,ab+=a*b*d,ac+=a*c*d,bc+=b*c*d;
    }
    inline void del(R double a,R double b,R double c,R double d){
        --sz,aa-=a*a*d,bb-=b*b*d,cc-=c*c*d,ab-=a*b*d,ac-=a*c*d,bc-=b*c*d;
    }
    double calc(){
        if(!sz)return 0;
        double invb=sgn(bb)?1.0/bb:0;
        double a=aa-ab*ab*invb,b=2*ac-2*ab*bc*invb,c=cc-bc*bc*invb;
        return !sgn(a)?c:c-b*b*0.25/a;
    }
}q;
struct Line{
    double a,b,c,d;
    inline Line(){}
    inline Line(R double x,R double y,R double xx,R double yy){
        !sgn(x-xx)?(a=1,b=0,c=-x):(a=(yy-y)/(xx-x),b=-1,c=y-a*x);
        d=1.0/(a*a+b*b);
    }
}L[N];
int top,op,i;double x,y,xx,yy,res;
int main(){
//  freopen("testdata.in","r",stdin);
//  freopen("testdata.out","w",stdout);
    for(int T=read();T;--T){
        op=getop();
        switch(op){
            case 0:{
                x=readdb(),y=readdb(),xx=readdb(),yy=readdb();
                L[++top]=Line(x,y,xx,yy),q.ins(L[top].a,L[top].b,L[top].c,L[top].d);
                break;
            }
            case 1:{
                i=read(),q.del(L[i].a,L[i].b,L[i].c,L[i].d);
                break;
            }
            case 2:{
                res=q.calc();
                if(res<1e-3&&res>-1e-3)res=0;
                printf("%.2lf\n",res);
                break;
            }
        }
    }
    return 0;
}