BZOJ3775: 点和直线(计算几何+拉格朗日乘数法)

题面

传送门

题解

劲啊……

没有和(Claris)一样推,用了类似于(Shinbokuow)推已知点求最短直线的方法,结果(WA)了好几个小时,拿(Claris)代码拍了几个小时都没找到(bug)在哪儿,最后发现是我一个除法的地方忘记除数为(0)的情况了……甘霖娘……

公式恐惧症患者可以直接转去结论了

设直线为(ax+by+c=0),点为((x,y)),记(d_i=a_i^2+b_i^2),那么就是要我们最小化

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

以下为了方便,记(A^2=sum{a_i^2over d_i})(B^2,C^2)同理,以及(AB=sum{a_ib_iover 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 fover partial y}=2yB^2+2xAB+2BC=0 ]

解得

[y={-xAB-BCover B^2} ]

代入原式可以化为

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

其中

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

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

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

易知(alpha geq 0)(证明下面有)

不过这里其实还有一个尴尬的情况就是有可能(B^2=0),也就是说所有直线的(b_i=0),不过我们转过头去看会发现这种情况下(y)(f(x,y))完全没有影响,而且(alpha,eta,gamma)的值分别就是(A^2,2AC,C^2)。所以这种情况其实并不会有影响

如果(alpha eq 0),我们要最小化(f(x,y)),同时还需要满足方程

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

有解

代入根的判别式,可知需要满足

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

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

最小值显然了

如果(alpha=0),则

[f(x,y)=eta x+gamma ]

(Claris)说这种情况下答案就等于(gamma)……然而我实在看不出为啥……我怎么感觉可以无限小呢……然而它要是变成负数显然不符合常理啊……有哪位鸽鸽知道为什么的么可以在下面留言哦qwq

然后就做完了

ps:关于(alphageq 0)的证明

因为有

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

首先显然(B^2geq 0),如果(B^2=0),那么根据上面所说(alpha=A^2geq 0),所以假设(B^2>0),我们需要证明

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

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

代入原来的值

[left(sum {a_i^2over d_i} ight)left(sum {b_i^2over d_i} ight)geq left(sum {a_ib_iover d_i} ight)^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]='',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
",res);
                break;
            }
        }
    }
    return 0;
}
原文地址:https://www.cnblogs.com/bztMinamoto/p/10705106.html