Codechef 13.9 TWOROADS 数学+计算几何+拉格朗日乘数法

 

BZOJ2508:简单题 拉格朗日乘数法

思路:考虑若给定一个定点\((x_0,y_0)\),以及某条直线\(ax+by+c=0\),那么距离的平方为:

\[\frac{(ax_0+by_0+c)^2}{a^2+b^2}\]

那么对于若干条直线,我们总能将距离的平方和表示为如下的式子:

\[Ax_0^2+By_0^2+Cx_0y_0+Dx_0+Ey_0+F\]

这是一个无限制求函数最值的问题,我们利用拉格朗日乘数法,搞出两个求导后的式子,令他们均为0,然后带入这个方程就行了.计算的复杂度为\(O(1)\).

容易发现增删直线对于系数的更改都是\(O(1)\)的.

于是我们就在\(O(m)\)的时间内解决了.

(对于系数的判断很"简单")

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
    
typedef double f2;
    
#define N 120010
f2 a[N],b[N],c[N],A,B,C,D,E,F;
   
void print(f2 x){
    if(fabs(x)<1e-6)puts("0.00");else printf("%.2lf\n",fabs(x));
}
 
f2 G[3][4],ansx,ansy;
inline void Solve(){
    register int i,j,k;f2 t;int n=2;
    if(fabs(G[1][1])<1e-6){
        if(fabs(G[1][2])<1e-6)ansx=0,ansy=G[2][3]/G[2][2];
        else ansy=G[1][3]/G[1][2],ansx=(G[2][3]-G[2][2]*ansy)/G[2][1];
        return;
    }
    if(fabs(G[2][2])<1e-6){
        if(fabs(G[2][1])<1e-6)ansx=G[1][3]/G[1][1],ansy=0;
        else ansx=G[2][3]/G[2][1],ansy=(G[1][3]-G[1][1]*ansx)/G[1][2];
        return;
    }
    t=-G[2][1]/G[1][1];G[2][1]=0,G[2][2]+=t*G[1][2],G[2][3]+=t*G[1][3];
    G[2][3]=(fabs(G[2][2])<1e-6)?0:G[2][3]/G[2][2],G[1][3]-=G[1][2]*G[2][3],G[1][3]/=G[1][1];
    ansx=G[1][3],ansy=G[2][3];
}
 
int main(){
    f2 x0,y0,x1,y1,t;int num=0,cnt,now=0;
    int q,qte,del;scanf("%d",&q);
    while(q--){
        scanf("%d",&qte);
        if(qte==0){
            scanf("%lf%lf%lf%lf",&x0,&y0,&x1,&y1);cnt=++num;++now;
            a[cnt]=y0-y1,b[cnt]=x1-x0,c[cnt]=y0*(x0-x1)-x0*(y0-y1),t=a[cnt]*a[cnt]+b[cnt]*b[cnt];
            A+=a[cnt]*a[cnt]/t,B+=b[cnt]*b[cnt]/t;C+=2*a[cnt]*b[cnt]/t;
            D+=2*a[cnt]*c[cnt]/t,E+=2*b[cnt]*c[cnt]/t,F+=c[cnt]*c[cnt]/t;
        }
        else if(qte==1){
            scanf("%d",&del);cnt=del,t=a[cnt]*a[cnt]+b[cnt]*b[cnt],--now;
            A-=a[cnt]*a[cnt]/t,B-=b[cnt]*b[cnt]/t,C-=2*a[cnt]*b[cnt]/t;
            D-=2*a[cnt]*c[cnt]/t,E-=2*b[cnt]*c[cnt]/t,F-=c[cnt]*c[cnt]/t;
        }
        else{
            if(now==0){puts("0.00");continue;}
            G[1][1]=2*A,G[1][2]=C,G[1][3]=-D,G[2][1]=C,G[2][2]=2*B,G[2][3]=-E;
            Solve();
            print(A*ansx*ansx+B*ansy*ansy+C*ansx*ansy+D*ansx+E*ansy+F);
        }
    }
    return 0;
}

BZOJ2876:[Noi2012]骑行川藏 拉格朗日乘数法

思路:套用拉格朗日乘数法在限制下函数最值的方法,添加未知数\(\lambda\),不过限制是什么呢?

显然能量应该满足限制.在最优意义下不难发现能量应该全部耗尽.于是需要我们最优化的函数如下:

\[min{\sum_{i=1}^{n}\frac{s_i}{v_i}+\lambda(\sum_{i=1}^{n}k_i(v_i-v_i')^2{s_i}-E)}\]

对于每个\(v_i\)求导,得到如下等式:

\[2\lambda{k_i}v_i^2(v_i-v_i')=1\]

对于\(\lambda\)求导,有:

\[\sum_{i=1}^{n}k_i(v_i-v_i')^2{s_i}=E\]

首先显然合理的\(v_i\)都应该满足\(v_i>v_i'\),我们考虑给定\(\lambda\)时,我们能够通过二分又第一种限制解出每个\(v_i\).

那么\(\lambda\)是什么?考虑第二种限制,显然左面的值是关于\(\lambda\)单调变化的,我们解出所有\(v_i\)check一下就可以通过二分解出\(\lambda\)的值了.

总之,二分套二分,时间复杂度\(O(nlog^2{n})\).

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef double f2;
 
#define N 10010
f2 s[N],k[N],_v[N];
 
f2 solve(f2 lam,int i){
    f2 L=_v[i],R=1e10,mid;
    while(R-L>1e-12){
        mid=(L+R)/2;
        if(2*lam*k[i]*mid*mid*(mid-_v[i])<1)L=mid;else R=mid;
    }return L;
}
int main(){
    int n;f2 E;scanf("%d%lf",&n,&E);
    register int i;for(i=1;i<=n;++i)scanf("%lf%lf%lf",&s[i],&k[i],&_v[i]);
     
    f2 llam=0,rlam=1e5,mid,add,tmp;
    while(rlam-llam>1e-12){
        mid=(llam+rlam)/2;
        for(add=0,i=1;i<=n;++i)tmp=solve(mid,i),add+=k[i]*s[i]*(tmp-_v[i])*(tmp-_v[i]);
        if(add>E)llam=mid;else rlam=mid;
    }
     
    f2 ans=0;for(i=1;i<=n;++i)ans+=s[i]/solve(llam,i);
     
    printf("%.10lf",ans);
     
    return 0;
}