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

 

BZOJ3564:[SHOI2014]信号增幅仪 坐标变换+计算几何+最小圆覆盖

思路:

首先我们将所有的点都顺时针旋转\(a\)度,这样就相当于坐标系逆时针旋转了\(a\)度,那么现在椭圆的长轴就与\(x\)轴平行了.

那么还要考虑椭圆现在是横向伸长了\(p\)倍,于是我们就将所有点目前的横坐标缩小\(p\)倍.

然后就是随机增量法的最小圆覆盖了.

 

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
 
typedef double f2;
 
static const f2 eps=1e-7;
inline int dcmp(const f2&x){return fabs(x)<eps?0:(x<0?-1:1);}
 
#define N 50010
 
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
     
    inline void read(){scanf("%lf%lf",&x,&y);}
     
    bool operator<(const Point&B)const{return dcmp(x-B.x)==0?y<B.y:x<B.x;}
     
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p)const{return Point(p*x,p*y);}
    Point operator/(const f2&p)const{return Point(x/p,y/p);}
}P[N];
 
f2 getdis(const Point&a,const Point&b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
f2 cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
 
Point getmid(const Point&a,const Point&b){return Point((a.x+b.x)/2,(a.y+b.y)/2);}
Point getvec(const Point&v){return Point(-v.y,v.x);}
Point rot(const Point&p,f2 alp){return Point(p.x*cos(alp)-p.y*sin(alp),p.x*sin(alp)+p.y*cos(alp));}
 
struct Line{
    Point p,v;
    Line(){}
    Line(const Point&p1,const Point&p2):p(p1),v(p1-p2){}
};
Line getvecline(const Point&p1,const Point&p2){
    Line re;
    re.p=getmid(p1,p2),re.v=getvec(p1-p2);
    return re;
}
Point getlineintersection(const Line&l1,const Line&l2){
    return l1.p+l1.v*(cross(l1.p-l2.p,l2.v)/cross(l1.v,l2.v));
}
Point getpoint(Point p1,Point p2,Point p3){
    static Point P[3];P[0]=p1,P[1]=p2,P[2]=p3;sort(P,P+3);p1=P[0],p2=P[1],p3=P[2];
    if(dcmp(cross(p1-p2,p2-p3))==0)return getmid(p1,p3);
    else return getlineintersection(getvecline(p1,p2),getvecline(p2,p3));
}
 
int main(){
    int n;scanf("%d",&n);register int i,j,k;
    for(i=1;i<=n;++i)P[i].read();
     
    f2 alp,p;scanf("%lf%lf",&alp,&p);
     
    for(i=1;i<=n;++i)P[i]=rot(P[i],-alp/180*acos(-1.0)),P[i].x/=p;
     
    for(i=2;i<=n;++i)swap(P[i],P[rand()%i+1]);
    Point o=P[1];f2 r=0;
    for(i=2;i<=n;++i)if(getdis(o,P[i])>r){
        o=P[i],r=0;
        for(j=1;j<i;++j)if(getdis(o,P[j])>r){
            o=getmid(P[i],P[j]),r=getdis(o,P[i]);
            for(k=1;k<j;++k)if(getdis(o,P[k])>r)o=getpoint(P[i],P[j],P[k]),r=getdis(o,P[i]);
        }
    }
     
    printf("%.3lf",r);
     
    return 0;
}

BZOJ2829:信用卡凸包 计算几何+凸包+坐标变换

思路:

我们发现最终的答案就是缩小一圈后的矩形的顶点们形成的凸包的周长再加上一个半径为\(r\)的圆的周长.

因为多边形外角和为\(360\)度.这个不用多说.

 

此外还有一点点问题就是旋转后的坐标变换.

如果是\((x,y)\)以原点为旋转中心逆时针旋转\(\alpha\)弧度,那么新的坐标为\((xcos\alpha-ysin\alpha,xsin\alpha+ycos\alpha)\).

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
 
typedef double f2;
 
#define N 200010
 
static const f2 eps=1e-8;
inline int dcmp(const f2&x){return (fabs(x)<eps)?0:(x<0?-1:1);}
 
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
     
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p)const{return Point(p*x,p*y);}
    Point operator/(const f2&p)const{return Point(x/p,y/p);}
     
    bool operator<(const Point&B)const{return dcmp(x-B.x)==0?(dcmp(y-B.y)<0):(dcmp(x-B.x)<0);}
}P[N<<2],PP[N<<2],stk[N<<2];
 
f2 getdis(const Point&a,const Point&b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point move(const Point&p,const f2&Ang,const f2&x,const f2&y){
    return Point(p.x*cos(Ang)-p.y*sin(Ang)+x,p.x*sin(Ang)+p.y*cos(Ang)+y);
}
 
int main(){
    int n;scanf("%d",&n);
     
    f2 a,b,r,x,y,Ang;scanf("%lf%lf%lf",&a,&b,&r);a-=2*r,b-=2*r;
     
    register int i,j;
    int cnt=0;
    for(i=1;i<=n;++i){
        scanf("%lf%lf%lf",&x,&y,&Ang);
         
        P[++cnt]=move(Point(b/2,a/2),Ang,x,y);
        P[++cnt]=move(Point(-b/2,a/2),Ang,x,y);
        P[++cnt]=move(Point(b/2,-a/2),Ang,x,y);
        P[++cnt]=move(Point(-b/2,-a/2),Ang,x,y);
    }
     
    int top=0;
    sort(P+1,P+cnt+1);
    int revins;for(revins=cnt;revins!=1&&P[revins].x==P[revins-1].x;--revins);
    for(i=revins;i<=(n+revins)/2;++i)swap(P[i],P[n+revins-i]);
     
    for(i=1;i<=cnt;++i){
        if(!top)stk[++top]=P[i];
        else{while(top>1&&(P[i].y-stk[top].y)*(stk[top].x-stk[top-1].x)<=(stk[top].y-stk[top-1].y)*(P[i].x-stk[top].x))--top;stk[++top]=P[i];}
    }
    int nowtop=top;
    for(i=revins-1;i>=1;--i){
        if(top==nowtop)stk[++top]=P[i];
        else{while(top>=nowtop+1&&(P[i].y-stk[top].y)*(stk[top].x-stk[top-1].x)<=(stk[top].y-stk[top-1].y)*(P[i].x-stk[top].x))--top;stk[++top]=P[i];}
    }--top;
     
    f2 re=0;
    for(i=1;i<top;++i)re+=getdis(stk[i],stk[i+1]);re+=getdis(stk[top],stk[1]);
     
    printf("%.2lf",re+2*acos(-1.0)*r);
     
    return 0;
}

凸多边形面积并 计算几何+扫描线

思路:

(本题现在可以在JDFZOJ1005找到)

我们首先来围观一篇神题解.

http://www.cnblogs.com/ch3656468/archive/2011/10/17/2215551.html

(果然\(O(n^3)\)不是最优的么...)

首先,将所有的凸多边形上的点按照逆时针排列,随后将所有的向量分为向左和向右两种,竖直方向的向量直接无视.

对于某些向量,它们可能在一条直线上,因此我们就将这些向量一起处理.

对于所有向左的向量(他们在上边界上),我们处理出所有在凸多边形边界上,但却不在凸多边形内部的部分与\(x\)轴形成的梯形面积的和.

然后对于向右的向量们也用同样的方法计算.

(然后看图不知道怎么样就发现)答案就是两者相减.时间复杂度仅为\(O(n^2logn)\).

我的代码:

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<climits>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef double f2;
 
static const f2 eps=1e-8;
static const f2 INF=1e50;
inline int dcmp(const f2&x){return fabs(x)<eps?0:(x<0?-1:1);}
 
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
    bool operator<(const Point&B)const{return x<B.x||(dcmp(x-B.x)==0&&y<B.y);}
    bool operator!=(const Point&B)const{return dcmp(x-B.x)!=0||dcmp(y-B.y)!=0;}
    bool operator==(const Point&B)const{return dcmp(x-B.x)==0&&dcmp(y-B.y)==0;}
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p)const{return Point(p*x,p*y);}
    Point operator/(const f2&p)const{return Point(x/p,y/p);}
};
f2 Cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
f2 Dot(const Point&a,const Point&b){return a.x*b.x+a.y*b.y;}
 
struct Line{
    Point p,v;
    Line(){}
    Line(const Point&p1,const Point&p2):p(p1),v(p1-p2){}
    f2 asky(f2 x){
        return dcmp(x-p.x)==0?p.y:p.y+(x-p.x)/v.x*v.y;
    }
    bool operator==(const Line&B)const{
        f2 b=dcmp(0-p.x)==0?p.y:p.y+(0-p.x)/v.x*v.y;
        f2 _b=dcmp(0-B.p.x)==0?B.p.y:B.p.y+(0-B.p.x)/B.v.x*B.v.y;
        return dcmp(Cross(v,B.v))==0&&dcmp(b-_b)==0;
    }
    bool operator<(const Line&B)const{
        f2 k1=dcmp(v.x)==0?INF:v.y/v.x,k2=dcmp(B.v.x)==0?INF:B.v.y/B.v.x;
        return k1<k2;
    }
 
};
Point GetLineIntersection(const Line&l1,const Line&l2){
    if(dcmp(Cross(l1.v,l2.v))==0)return Point(-INF,-INF);
    return l1.p+l1.v*(Cross(l1.p-l2.p,l2.v)/Cross(l1.v,l2.v));
}
 
#define N 110
int n;
struct Polygen{
    Point P[N];int siz;
    void Init(){scanf("%d",&siz);for(int i=1;i<=siz;++i)scanf("%lf%lf",&P[i].x,&P[i].y);P[siz+1]=P[1];}
}Poly[11];
 
struct Quest{
    f2 x;int type;
    Quest(){}
    Quest(f2 _x,int _type):x(_x),type(_type){}
    bool operator<(const Quest&B)const{return x<B.x;}
};
Line Lines[N*N],reallines[N*N],q1[N*N],q2[N*N],seq[N*N];int cnt,id,n1,n2;
 
inline f2 Solve(Line*L,int m){
    sort(L+1,L+m+1);
    register int i,j,k;
    for(id=0,i=1;i<=m;i=j+1){for(j=i;j!=m&&L[j]==L[j+1];++j);seq[++id]=L[i];}
     
    f2 res=0;
     
    vector<Quest>q;
    for(i=1;i<=id;++i){
        q.clear();
        for(j=1;j<=n;++j){
            bool touch=0;
            for(k=1;k<=Poly[j].siz&&!touch;++k)if(seq[i]==Line(Poly[j].P[k],Poly[j].P[k+1]))touch=1;
            if(touch)continue;
             
            vector<Point>get;
            for(k=1;k<=Poly[j].siz;++k){
                Point x=GetLineIntersection(seq[i],Line(Poly[j].P[k],Poly[j].P[k+1]));
                if(x!=Point(-INF,-INF)&&Dot(Poly[j].P[k]-x,Poly[j].P[k+1]-x)<=0)get.push_back(x);
            }
            sort(get.begin(),get.end());
            vector<Point>realget;
            for(k=0;k<get.size();++k)if(!k||get[k]!=get[k-1])realget.push_back(get[k]);
             
            if(realget.size()==2){
                q.push_back(Quest(realget[0].x,0)),q.push_back(Quest(realget[1].x,1));
            }
        }
         
        f2 Mx,Mn;
        for(j=1;j<=m;++j)if(dcmp(Cross(seq[i].v,L[j].v))==0&&dcmp(seq[i].asky(0)-L[j].asky(0))==0){
            Mx=max(L[j].p.x,(L[j].p+L[j].v).x),Mn=min(L[j].p.x,(L[j].p+L[j].v).x);
            q.push_back(Quest(Mn,2)),q.push_back(Quest(Mx,3));
        }
         
        sort(q.begin(),q.end());
         
        int in=0,cnt=0;
        for(j=0;j<q.size();++j){
            f2 y0=seq[i].asky(q[j==0?0:j-1].x),y1=seq[i].asky(q[j].x);
            if(!in&&cnt)res+=(y0+y1)*(j==0?0:q[j].x-q[j-1].x)*.5;
            if(q[j].type==0)++in;if(q[j].type==1)--in;
            if(q[j].type==2)++cnt;if(q[j].type==3)--cnt;
        }
    }
     
    return res;
}
int main(){
    //freopen("tt.in","r",stdin);
    scanf("%d",&n);register int i,j;
    for(i=1;i<=n;++i){
        Poly[i].Init();
        for(j=1;j<=Poly[i].siz;++j){
            if(dcmp(Poly[i].P[j].x-Poly[i].P[j+1].x)==0)continue;
            if(dcmp(Poly[i].P[j].x-Poly[i].P[j+1].x)<0)q2[++n2]=Line(Poly[i].P[j],Poly[i].P[j+1]);
            else q1[++n1]=Line(Poly[i].P[j],Poly[i].P[j+1]);
        }
    }
     
    printf("%.3lf",Solve(q1,n1)-Solve(q2,n2));
     
    return 0;
}

[POI2007]对称轴Axes of Symmetry-osi 计算几何+KMP

思路:

网上的题解也有很多了,主要介绍一下我得到的姿势.

我们抽象出边和点的信息,将多边形成一个长度为\(2n\)的字符串.

然后我们考虑有多少个这个字符串的循环同构串与这个串的反串相同.每出现一个就有一个对称轴.

想想好像是挺显然的.

#include<cstdio>
#include<cctype>
#include<cstring>
#include<iostream>
#include<algorithm>
 
typedef long long LL;
 
#define N 100010
 
struct Point{
    int x,y;
    inline void read(){scanf("%d%d",&x,&y);}
    Point(){}
    Point(int _x,int _y):x(_x),y(_y){}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
}P[N];
inline LL Cross(const Point&a,const Point&b){return(LL)a.x*b.y-(LL)a.y*b.x;}
inline LL sqr(const int&x){return(LL)x*x;}
inline LL dis2(const Point&a,const Point&b){return sqr(a.x-b.x)+sqr(a.y-b.y);}
inline LL abs(const LL&x){return x<0?-x:x;}
LL seq1[N<<2],seq2[N<<1];
 
int n;
inline int f(int x){if(x<=0)return n;else if(x>n)return 1;else return x;}
 
int pre[N<<1];
int main(){
    int cas;scanf("%d",&cas);register int i,j;
    while(cas--){
        scanf("%d",&n);for(i=1;i<=n;++i)P[i].read();
        int cnt=0;for(i=1;i<=n;++i)seq1[++cnt]=Cross(P[f(i-1)]-P[i],P[i]-P[f(i+1)])|(1LL<<62),seq1[++cnt]=dis2(P[i],P[f(i+1)]);
         
        for(i=1;i<=cnt;++i)seq2[i]=seq1[cnt-i+1];
        for(i=cnt+1;i<=2*cnt;++i)seq1[i]=seq1[i-cnt];
         
        for(j=pre[1]=0,i=2;i<=cnt;++i){
            while(j&&seq2[j+1]!=seq2[i])j=pre[j];
            if(seq2[j+1]==seq2[i])++j;pre[i]=j;
        }
        int res=0;
        for(i=1,j=0;i<2*cnt;++i){
            while(j&&seq2[j+1]!=seq1[i])j=pre[j];
            if(seq2[j+1]==seq1[i])++j;if(j==cnt)++res;
        }
        printf("%d\n",res);
    }return 0;
}

APIO2010信号覆盖 单调性+计算几何+容斥原理

思路:

一开始只是YY出了一个\(O(n^3logn)\)的算法.

我们枚举每两个点,同时再枚举另外一个点,就能得到这三个点的圆心前在两个点垂直平分线的位置.

那么现在我们再枚举每个点出现在这种圆中有多少种情况:这个点已经给定时,我们可以考虑圆心在垂直平分线上的可行位置,通过解不等式发现仅仅满足一个限制就是可行的圆心.于是就可以用数据结构维护快速得到有多少个可行的圆心.

然而正解的思路非常巧妙:

考虑每四个点,要么形成一个凸多边形,要么形成一个凹多边形.

由于没有任意四个点在一个圆上,那么对于一个凸多边形对于答案的贡献是2,一个凹多边形对于答案的贡献是1.

考虑如何计算凹多边形的数目.(然后我们也可以知道凸多边形的数目)

我们枚举中间的凹点,那么发现形成凹多边形的另外三个点必然不在凹点的同侧.因此我们就转化为计算总数减去在凹点同侧的三元组数目,这个利用极角序排序以及单调性可以\(O(nlogn)\)得到.

因此就在\(O(n^2logn)\)的时间内做出来了.

#include<cmath>
#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef long long LL;
typedef double f2;
 
#define N 1510
int x[N],y[N];
 
struct Point{
    int x,y;
    f2 Ang;
    inline void read(){scanf("%d%d",&x,&y);}
    Point(){}
    Point(int _x,int _y):x(_x),y(_y){}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    bool operator<(const Point&B)const{return Ang<B.Ang;}
    inline void print(){printf("%d %d\n",x,y);}
}P[N],seq[N<<1];int cnt;
inline LL Cross(const Point&a,const Point&b){return(LL)a.x*b.y-(LL)a.y*b.x;}
 
int main(){
    //freopen("tt.in","r",stdin);
    int n;scanf("%d",&n);int i,j,k;
    for(i=1;i<=n;++i)P[i].read();
     
    LL res=0;
    for(i=1;i<=n;++i){
        for(cnt=0,j=1;j<=n;++j)if(i^j)seq[++cnt]=P[j],seq[cnt].Ang=atan2(P[j].y-P[i].y,P[j].x-P[i].x);
        sort(seq+1,seq+cnt+1);
        for(j=cnt+1;j<=2*cnt;++j)seq[j]=seq[j-cnt];
         
        res+=(LL)(n-1)*(n-2)*(n-3)/6;
        int get;k=1;
        for(j=1;j<=cnt;++j){
            while(k+1!=j+cnt&&Cross(seq[j]-P[i],P[i]-seq[k+1])<=0)++k;
            //printf("j=%d ",j);seq[j].print();printf("k=%d ",k),seq[k].print();puts("");
            get=k-j;
            res-=(LL)get*(get-1)/2;
        }
    }
     
    LL all=(LL)n*(n-1)*(n-2)*(n-3)/24;
    LL ret=res+(all-res)*2;
    printf("%.3lf",3.0+(f2)ret*6/n/(n-1)/(n-2));
     
    return 0;
}

POJ2187 Beauty Contest 计算几何+旋转卡壳

思路:

本人自己写了一个傻逼的不能再傻逼的御用凸包模板,谁要是扒的话简直就是作死...

然后这道题目主要就是利用旋转卡壳寻找对踵点,然后更新答案.

然后旋转卡壳非常牛逼是\(O(n)\)的.这是为什么呢?

对于凸包上的每一条边,对踵点与这条边形成的三角形的面积必定最大.因此,这条边之外的点是成一个单峰函数分布的.

此外我们还很容易从凸包上得到单调性.

利用这两个性质,就很容易做到凸包上的线性扫描了.

#include<set>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;

#define N 50010

struct Point{
	int x,y;
	Point(){}
	Point(int _x,int _y):x(_x),y(_y){}
	bool operator<(const Point&B)const{return(x<B.x)||(x==B.x&&y>B.y);}
	Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
	Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
}P[N];int cnt;
inline int Cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
inline int sqr(int x){return x*x;}
inline int dis2(const Point&a,const Point&b){return sqr(a.x-b.x)+sqr(a.y-b.y);}
set<Point>S;

Point stk[N];int top;

#define f(x) (x>top?1:x)
int main(){
	//freopen("tt.in","r",stdin);
	int n;scanf("%d",&n);register int i,j;
	
	if(n<=1){puts("0");return 0;}
	
	int x,y;
	while(n--){scanf("%d%d",&x,&y);if(S.find(Point(x,y))==S.end())S.insert(Point(x,y)),P[++cnt]=Point(x,y);}
	
	sort(P+1,P+cnt+1);
	
	bool allline=1;
	for(i=1;i+2<=cnt;++i)if(Cross(P[i]-P[i+1],P[i+1]-P[i+2])!=0)allline=0;
	if(allline){printf("%d",dis2(P[1],P[cnt]));return 0;}
	
	int revins;for(revins=cnt;revins!=1&&P[revins].x==P[revins-1].x;--revins);
	for(i=revins;i<=(n+revins)/2;++i)swap(P[i],P[n+revins-i]);
	
	for(i=1;i<=cnt;++i){
		if(!top)stk[++top]=P[i];
		else{while(top>1&&(P[i].y-stk[top].y)*(stk[top].x-stk[top-1].x)<=(stk[top].y-stk[top-1].y)*(P[i].x-stk[top].x))--top;stk[++top]=P[i];}
	}
	int nowtop=top;
	for(i=revins-1;i>=1;--i){
		if(top==nowtop)stk[++top]=P[i];
		else{while(top>nowtop+1&&(P[i].y-stk[top].y)*(stk[top].x-stk[top-1].x)<=(stk[top].y-stk[top-1].y)*(P[i].x-stk[top].x))--top;stk[++top]=P[i];}
	}--top;
	
	int r=2,res=0;
	for(i=1;i<=top;++i){
		while(Cross(stk[i]-stk[f(i+1)],stk[i]-stk[r])<=Cross(stk[i]-stk[f(i+1)],stk[i]-stk[f(r+1)]))r=f(r+1);
		res=max(res,dis2(stk[i],stk[r]));
		res=max(res,dis2(stk[f(i+1)],stk[r]));
	}
	
	printf("%d",res);
	return 0;
}

hdu3932 Groundhog Build Home 计算几何+最小圆覆盖

思路:

本来想写的是求出凸包然后在凸包上\(O(n^3)\)进行暴力的算法.虽然平面上的随机点形成的凸包点数很少,不过显然如果直接给出一个凸包就直接卡掉了.

于是我学习的是随机增量法求最小圆覆盖.

网上题解很多,就不解释了(还不太懂),如果随机打乱一下的话,好像是能够做到期望\(O(n)\),相当厉害的算法.

一个麻烦就是求出三角形的外心,用了一大堆叉积来算好像很傻逼...

#include<cstdio>
#include<cstring>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;

#define N 1010

#define eps 1e-8
typedef double f2;
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
    bool operator<(const Point&B)const{return x<B.x;}
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
    Point operator*(const f2&p){return Point(p*x,p*y);}
    Point operator/(const f2&p){return Point(x/p,y/p);}
}P[N];
inline f2 sqr(const f2&x){return x*x;}
f2 Cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;} 
f2 getdist(const Point&a,const Point&b){return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}

struct Line{
    Point p,v;
    Line(){}
    Line(Point _p,Point _v):p(_p),v(_v){}
};
Point GetLineIntersection(Line L1,Line L2){
    return L1.p+L1.v*(Cross(L1.p-L2.p,L2.v)/Cross(L1.v,L2.v));
}
Point Getvert(const Point&a){return Point(-a.y,a.x);}
Point Getmid(const Point&a,const Point&b){return Point((a.x+b.x)/2,(a.y+b.y)/2);}

Line GetLine(const Point&a,const Point&b){
    return Line(Getmid(a,b),Getvert(a-b));
}
Point Getcenter(Point a,Point b,Point c){
    Point P[3]={a,b,c};sort(P,P+3);a=P[0],b=P[1],c=P[2];
    if(fabs(Cross(a-b,b-c))<eps)return Getmid(a,c);
    else return GetLineIntersection(GetLine(a,b),GetLine(b,c));
}
    
int main(){
    int X,Y,n;register int i,j,k;
    while(scanf("%d%d%d",&X,&Y,&n)!=EOF){
        for(i=1;i<=n;++i)scanf("%lf%lf",&P[i].x,&P[i].y);
        random_shuffle(P+1,P+n+1);
        Point center=P[1];f2 r=0;
        for(i=2;i<=n;++i){
            if(getdist(center,P[i])>r){
                center=P[i],r=0;
                for(j=1;j<i;++j){
                    if(getdist(center,P[j])>r){
                        center=Getmid(P[i],P[j]),r=getdist(P[i],P[j])*.5;
                        for(k=1;k<j;++k)
                            if(getdist(center,P[k])>r)center=Getcenter(P[i],P[j],P[k]),r=getdist(center,P[i]);
                    }
                }
            }
        }
        printf("(%.1lf,%.1lf).\n%.1lf\n",center.x,center.y,r);
    }
    return 0;
}

BZOJ1822:[JSOI2010]Frozen Nova 冷冻波 二分+计算几何+网络流

思路:

一眼网络流,套一个二分就行了.

关键是如何判断巫妖和精灵之间能否有圆阻挡.

一开始我是进行了复杂的判断:如果线段平行于x轴或者平行于y轴直接特判;

否则求出圆心到直线的距离,如果距离不超过半径,则用一个简单的公式求出圆心到直线最近点的横坐标,利用这个横坐标与线段两个端点的横坐标比较,得出结果.

可是数据跟我想的根本不是一回事啊啊啊!

根本只需要判断距离啊啊啊!!而且相切的情况不算相交啊啊啊!!

什么坑爹数据,我是领教到了.

#include<cstdio>
#include<cstring>
#include<climits>
#include<iostream>
#include<algorithm>
using namespace std;

typedef long long LL;

#define INF 0x3f3f3f3f

#define N 210
#define M 210
#define P 210

struct MaxflowSolver{
	static const int V=410;
	static const int E=100010;
	int head[V],next[E],end[E],flow[E],ind,d[V];
	inline void reset(){ind=0,memset(head,-1,sizeof head);}
	inline void addedge(int a,int b,int f){int q=ind++;end[q]=b,next[q]=head[a],head[a]=q,flow[q]=f;}
	inline void make(int a,int b,int f){addedge(a,b,f),addedge(b,a,0);}
	inline bool bfs(int S,int T){
		static int q[V];int fr=0,ta=0;
		memset(d,-1,sizeof d),d[S]=0,q[ta++]=S;
		while(fr^ta){
			int i=q[fr++];for(int j=head[i];j!=-1;j=next[j])if(flow[j]&&d[end[j]]==-1)d[end[j]]=d[i]+1,q[ta++]=end[j];
		}return d[T]!=-1;
	}
	inline int dinic(int p,int T,int Maxflow){
		if(p==T)return Maxflow;
		int last=Maxflow;
		for(int j=head[p];j!=-1;j=next[j])if(flow[j]&&d[end[j]]==d[p]+1){
			int to=dinic(end[j],T,last>flow[j]?flow[j]:last);
			if(to){flow[j]-=to,flow[j^1]+=to;if(!(last-=to))return Maxflow;}
		}d[p]=-1;return Maxflow-last;
	}
	inline int Maxflow(int S,int T){
		int res=0;while(bfs(S,T))res+=dinic(S,T,INF);return res;
	}
}Stg;
int x1[N],y1[N],x2[M],y2[M],lim[N],cd[N],ox[P],oy[P],r[P];

int G[N][M];

inline int sqr(int x){return x*x;}



int main(){
	int n,m,p;scanf("%d%d%d",&n,&m,&p);register int i,j,k;
	for(i=1;i<=n;++i)scanf("%d%d%d%d",&x1[i],&y1[i],&lim[i],&cd[i]);
	for(i=1;i<=m;++i)scanf("%d%d",&x2[i],&y2[i]);
	for(i=1;i<=p;++i)scanf("%d%d%d",&ox[i],&oy[i],&r[i]);
	
	int a,b,c;
	for(i=1;i<=n;++i)for(j=1;j<=m;++j)if(sqr(x1[i]-x2[j])+sqr(y1[i]-y2[j])<=sqr(lim[i])){
		bool ok=1;
		for(k=1;k<=p&&ok;++k){
			a=y1[i]-y2[j],b=x2[j]-x1[i],c=(x1[i]-x2[j])*y1[i]-(y1[i]-y2[j])*x1[i];
			if((LL)(a*ox[k]+b*oy[k]+c)*(a*ox[k]+b*oy[k]+c)<(LL)r[k]*r[k]*(a*a+b*b))ok=0;
		}
		if(ok)G[i][j]=1;
	}
	
	int L=0,R=1<<30,mid;
	while(L<R){
		mid=(L+R)>>1;
		Stg.reset();
		for(i=1;i<=n;++i)Stg.make(0,i,mid/cd[i]+1);
		for(i=1;i<=m;++i)Stg.make(n+i,n+m+1,1);
		for(i=1;i<=n;++i)for(j=1;j<=m;++j)if(G[i][j])Stg.make(i,n+j,1);
		if(Stg.Maxflow(0,n+m+1)==m)R=mid;else L=mid+1;
	}
	
	if(L==1<<30)puts("-1");else printf("%d",L);
	return 0;
}

BZOJ2338:[HNOI2011]数矩形 暴力+计算几何

思路:

首先搞出所有的线段,然后中点相同的放在一起,并且长度相同的放在一起.

那么,中点相同并且长度相同,这两个线段就能够构成一个矩形.

这样的线段不会超过\(O(n)\)条,所以我们可以直接暴力枚举.

(其实我觉得可以把所有的线段以中点为原点进行极角序排序,然后我们就是要找出两条线段使得它们之间的夹角尽可能接近90,这大概就可以线性扫描了.)

本代码完全避免了精度问题!速来点赞!

#include<cstdio>
#include<cstring>
#include<cctype>
#include<climits>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef long long LL;
 
struct Point{
    int x,y;
    Point(){}
    Point(int _x,int _y):x(_x),y(_y){}
 
    bool operator<(const Point&B)const{return(x!=B.x)?x<B.x:y<B.y;}
    bool operator!=(const Point&B)const{return(x!=B.x)||(y!=B.y);}
    Point operator+(const Point&B)const{return Point(x+B.x,y+B.y);}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
};
LL Cross(const Point&a,const Point&b){return(LL)a.x*b.y-(LL)a.y*b.x;}
LL Getans(Point p1,Point p2,Point q1,Point q2){
    LL res=Cross(p1-q1,p1-q2);
    return res<0?-res:res;
}
 
int gcd(int a,int b){return(!b)?a:gcd(b,a%b);}
 
struct Double{
    int x,y;
    Double(){}
    Double(int _x,int _y):x(_x),y(_y){}
 
    bool operator<(const Double&B)const{return(LL)x*B.y<(LL)y*B.x;}
};
 
 
#define N 1510
struct Segment{
    Point P1,P2;LL dis;
    bool insequal(const Segment&B)const{
        if(P1.x+P2.x!=B.P1.x+B.P2.x)return 0;
        if(P1.y+P2.y!=B.P1.y+B.P2.y)return 0;
        return 1;
    }
    bool operator<(const Segment&B)const{
        if(P1.x+P2.x!=B.P1.x+B.P2.x)return P1.x+P2.x<B.P1.x+B.P2.x;
        if(P1.y+P2.y!=B.P1.y+B.P2.y)return P1.y+P2.y<B.P1.y+B.P2.y;
        return dis<B.dis;
    }
    bool operator==(const Segment&B)const{return insequal(B)&&dis==B.dis;}
}S[N*N>>1];int id;
 
int x[N],y[N],n;
 
inline LL sqr(int x){
    return(LL)x*x;
}
 
int main(){
    scanf("%d",&n);register int i,j,k,p;for(i=1;i<=n;++i)scanf("%d%d",&x[i],&y[i]);
 
    for(i=1;i<=n;++i)for(j=i+1;j<=n;++j)++id,S[id].P1=Point(x[i],y[i]),S[id].P2=Point(x[j],y[j]),S[id].dis=sqr(x[i]-x[j])+sqr(y[i]-y[j]);
     
    sort(S+1,S+id+1);
 
    LL res=0;
    for(i=1;i<=id;i=j+1){
        for(j=i;j!=id&&S[j]==S[j+1];++j);
        for(k=i;k<=j;++k)for(p=k+1;p<=j;++p)res=max(res,Getans(S[k].P1,S[k].P2,S[p].P1,S[p].P2));
    }
 
    cout<<res<<endl;
 
    return 0;
}