计算几何模板


#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
#include<stdlib.h>
using namespace std;
const int MAXN=550;
const double eps=1e-8;

struct Point
{
    double x,y,z;
    Point(){}

    Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}

    //两向量之差
    Point operator -(const Point p1)
    {
        return Point(x-p1.x,y-p1.y,z-p1.z);
    }

    //两向量之和
    Point operator +(const Point p1)
    {
        return Point(x+p1.x,y+p1.y,z+p1.z);
    }

    //叉乘
    Point operator *(const Point p)
    {
        return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
    }

    Point operator *(double d)
    {
        return Point(x*d,y*d,z*d);
    }

    Point operator / (double d)
    {
        return Point(x/d,y/d,z/d);
    }

    //点乘
    double  operator ^(Point p)
    {
        return (x*p.x+y*p.y+z*p.z);
    }
};

struct CH3D
{
    struct face
    {
        //表示凸包一个面上的三个点的编号
        int a,b,c;
        //表示该面是否属于最终凸包上的面
        bool ok;
    };
    //初始顶点数
    int n;
    //初始顶点
    Point P[MAXN];
    //凸包表面的三角形数
    int num;
    //凸包表面的三角形
    face F[8*MAXN];
    //凸包表面的三角形
    int g[MAXN][MAXN];
    //向量长度
    double vlen(Point a)
    {
        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
    }
    //叉乘
    Point cross(const Point &a,const Point &b,const Point &c)
    {
        return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),
                     (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),
                     (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)
                     );
    }
    //三角形面积*2
    double area(Point a,Point b,Point c)
    {
        return vlen((b-a)*(c-a));
    }
    //四面体有向体积*6
    double volume(Point a,Point b,Point c,Point d)
    {
        return (b-a)*(c-a)^(d-a);
    }
    //正:点在面同向
    double dblcmp(Point &p,face &f)
    {
        Point m=P[f.b]-P[f.a];
        Point n=P[f.c]-P[f.a];
        Point t=p-P[f.a];
        return (m*n)^t;
    }
    void deal(int p,int a,int b)
    {
        int f=g[a][b];//搜索与该边相邻的另一个平面
        face add;
        if(F[f].ok)
        {
            if(dblcmp(P[p],F[f])>eps)
              dfs(p,f);
            else
            {
                add.a=b;
                add.b=a;
                add.c=p;//这里注意顺序,要成右手系
                add.ok=true;
                g[p][b]=g[a][p]=g[b][a]=num;
                F[num++]=add;
            }
        }
    }
    void dfs(int p,int now)//递归搜索所有应该从凸包内删除的面
    {
         F[now].ok=0;
         deal(p,F[now].b,F[now].a);
         deal(p,F[now].c,F[now].b);
         deal(p,F[now].a,F[now].c);
    }
    bool same(int s,int t)
    {
        Point &a=P[F[s].a];
        Point &b=P[F[s].b];
        Point &c=P[F[s].c];
        return fabs(volume(a,b,c,P[F[t].a]))<eps &&
               fabs(volume(a,b,c,P[F[t].b]))<eps &&
               fabs(volume(a,b,c,P[F[t].c]))<eps;
    }
    //构建三维凸包
    void create()
    {
        int i,j,tmp;
        face add;

        num=0;
        if(n<4)return;
    //**********************************************
        //此段是为了保证前四个点不共面
        bool flag=true;
        for(i=1;i<n;i++)
        {
            if(vlen(P[0]-P[i])>eps)
            {
                swap(P[1],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前三个点不共线
        for(i=2;i<n;i++)
        {
            if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps)
            {
                swap(P[2],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前四个点不共面
        for(int i=3;i<n;i++)
        {
            if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps)
            {
                swap(P[3],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
    //*****************************************
        for(i=0;i<4;i++)
        {
            add.a=(i+1)%4;
            add.b=(i+2)%4;
            add.c=(i+3)%4;
            add.ok=true;
            if(dblcmp(P[i],add)>0)swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;
            F[num++]=add;
        }
        for(i=4;i<n;i++)
        {
            for(j=0;j<num;j++)
            {
                if(F[j].ok&&dblcmp(P[i],F[j])>eps)
                {
                    dfs(i,j);
                    break;
                }
            }
        }
        tmp=num;
        for(i=num=0;i<tmp;i++)
          if(F[i].ok)
            F[num++]=F[i];

    }
    //表面积
    double area()
    {
        double res=0;
        if(n==3)
        {
            Point p=cross(P[0],P[1],P[2]);
            res=vlen(p)/2.0;
            return res;
        }
        for(int i=0;i<num;i++)
          res+=area(P[F[i].a],P[F[i].b],P[F[i].c]);
        return res/2.0;
    }
    double volume()
    {
        double res=0;
        Point tmp(0,0,0);
        for(int i=0;i<num;i++)
           res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(res/6.0);
    }
    //表面三角形个数
    int triangle()
    {
        return num;
    }
    //表面多边形个数
    int polygon()
    {
        int i,j,res,flag;
        for(i=res=0;i<num;i++)
        {
            flag=1;
            for(j=0;j<i;j++)
              if(same(i,j))
              {
                  flag=0;
                  break;
              }
            res+=flag;
        }
        return res;
    }
    //三维凸包重心
    Point barycenter()
    {
        Point ans(0,0,0),o(0,0,0);
        double all=0;
        for(int i=0;i<num;i++)
        {
            double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
            ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
            all+=vol;
        }
        ans=ans/all;
        return ans;
    }
    //点到面的距离
    double ptoface(Point p,int i)
    {
        return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
    }
};
CH3D hull;
int main()
{
   // freopen("in.txt","r",stdin);
   // freopen("out.txt","w",stdout);
    while(scanf("%d",&hull.n)==1)  //输入点
    {
        for(int i=0;i<hull.n;i++)
        {
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z)  //输入每个点;
        }
        hull.create();  //创建凸包
        printf("%d\n",hull.polygon());//求凸包表面多边形个数
}
//给一个三维凸包,求重心到表面的最短距离
//模板题:三维凸包+多边形重心+点面距离

 	 while(scanf("%d",&hull.n)==1)
    {
        for(int i=0;i<hull.n;i++)
        {
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        }
        hull.create(); //建立凸包
        Point p=hull.barycenter() //凸包重心; 
        double minn=1e20;
        for(int i=0;i<hull.num;i++)  //凸包面的个数
        {
            minn=min(minn,hull.ptoface(p,i));  //p点到第i个面的距离
        }
        printf("%.3lf\n",minn);
    }
while(scanf("%d",&hull.n)==1)
    {
        for(int i=0;i<hull.n;i++)
        {
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        }
        hull.create();
        printf("%.3f\n",hull.area());//凸包表面积
    }
    return 0;
}



给出n个点,用这些点构造三角形,然后找出相似三角形个数的最大值。
/*
数相似三角形个数
用map记录角度
注意:
1、去除重点
2、排除共线的情况
3、精度问题。

WA了8次了
*/

#include<stdio.h>
#include<math.h>
#include<map>
#include<iostream>
#include<string.h>
using namespace std;
const double eps=1e-6;
map<long long,int>mp;
struct Node
{
    int x,y;
}node[50];

double dis(Node a,Node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

void solve(Node a,Node b,Node c)
{
    double la=dis(b,c);
    double lb=dis(a,c);
    double lc=dis(a,b);
    if(la+lb<=lc+eps)return;
    if(lb+lc<=la+eps)return;
    if(la+lc<=lb+eps)return;
    double A=acos((lb*lb+lc*lc-la*la)/(2*lb*lc));
    double B=acos((la*la+lc*lc-lb*lb)/(2*la*lc));
    double C=acos((la*la+lb*lb-lc*lc)/(2*la*lb));
    if(A<eps||B<eps||C<eps)return;
    int t1=(int)(A*10000);
    int t2=(int)(B*10000);
    int t3=(int)(C*10000);
    if(t1>t2)swap(t1,t2);
    if(t1>t3)swap(t1,t3);
    if(t2>t3)swap(t2,t3);
    if(t1==0)return;
    long long t=t1*1000000*1000000+t2*1000000+t3;
    mp[t]++;
}
int hole[220][220];
int main ()
{
    int n;
    while(scanf("%d",&n),n)
    {
        int t=0;
        memset(hole,0,sizeof(hole));
        for(int i=0;i<n;i++)
        {
            scanf("%d%d",&node[t].x,&node[t].y);
            if(hole[node[t].x+100][node[t].y+100]==0)
            {
                hole[node[t].x+100][node[t].y+100]=1;
                t++;
            }
        }
        n=t;
        mp.clear();
        for(int i=0;i<n;i++)
          for(int j=i+1;j<n;j++)
            for(int k=j+1;k<n;k++)
              solve(node[i],node[j],node[k]);
        int ans=0;
        map<long long,int>::iterator it;
        for(it=mp.begin();it!=mp.end();it++)
        {
            int t=it->second;
            if(t>ans)ans=t;
        }
        printf("%d\n",ans);
    }
    return 0;
}

题意:
给定n个点,现在要求找一个点作为一个半径为1的圆的圆心,使得这个圆能覆盖的点个数最多,输出最多能覆盖多少个点。
题解:
我们枚举两个点,认为他们是在这个圆上,那么就能确定圆心。之后可以从所有点到圆心的距离判断覆盖情况。
证:假设现在确定了一个能覆盖最多点的圆,那么我们可以移动这个圆使得覆盖的点数不变,且至少有两个点在这个圆上。
每两个点可以确定两个圆心,但只要用一个就够了(都用两点构成向量的上方圆心或者下方圆心)。
证明:若圆可以覆盖三个以上的点(两个点不需要找圆心。。),那么任意取覆盖的最外围三个点连线得到三角形ABC。我们枚举的时候枚举AB,AC,无论是先向量上方还是下方,都必定存在一个枚举的圆心,在三角形ABC内部。也就是说只要枚举一个圆心就能确定答案了(可以减少一半的时间消耗)。

注意:两点距离大于2的可以不枚举;初始化ans=1,否则一个点的时候会出错;距离判断时可以用距离的平方判断,因为sqrt很耗时间。
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<math.h>
using namespace std;
const int MAXN=330;
const double eps=1e-4;

double p[MAXN][2];
double xx1,yy1,xx2,yy2;

double dis(int i,int j)
{
    return sqrt((p[i][0]-p[j][0])*(p[i][0]-p[j][0])+(p[i][1]-p[j][1])*(p[i][1]-p[j][1]));
}

void get_center_point(int a,int b)
{
    double x0=(p[a][0]+p[b][0])/2;
    double y0=(p[a][1]+p[b][1])/2;
    double d=sqrt(1-((p[a][0]-p[b][0])*(p[a][0]-p[b][0])+(p[a][1]-p[b][1])*(p[a][1]-p[b][1]))/4);
    if(fabs(p[a][1]-p[b][1])<1e-6)
    {
        xx1=x0;
        yy1=y0+d;
        xx2=x0;
        yy2=y0-d;
    }
    else
    {
        double tmp=atan(-(p[a][0]-p[b][0])/(p[a][1]-p[b][1]));
        double dx=d*cos(tmp);
        double dy=d*sin(tmp);
        xx1=x0+dx;
        yy1=y0+dy;
        xx2=x0-dx;
        yy2=y0-dy;
    }
}

int main()
{
    int T;
    int n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        for(int i=0;i<n;i++)
           scanf("%lf%lf",&p[i][0],&p[i][1]);
        int ans=1;
        for(int i=0;i<n;i++)
          for(int j=i+1;j<n;j++)
          {
              if(dis(i,j)>2.0)continue;
              get_center_point(i,j);
              int cnt=0;
              for(int i=0;i<n;i++)
                if(sqrt((p[i][0]-xx1)*(p[i][0]-xx1)+(p[i][1]-yy1)*(p[i][1]-yy1))<1+eps)
                  cnt++;
              if(cnt>ans)ans=cnt;
              cnt=0;
              for(int i=0;i<n;i++)
                if(sqrt((p[i][0]-xx2)*(p[i][0]-xx2)+(p[i][1]-yy2)*(p[i][1]-yy2))<1+eps)
                  cnt++;
              if(cnt>ans)ans=cnt;
          }
        printf("%d\n",ans);
    }
    return 0;
}
题目:在三维空间中,给出太阳的位置(看成点光源),空间中存在一个凸物体,问物体在一个给定的平面内的影子大小。
分析:计算几何、凸包、三维坐标转换。求解分为四个步骤:1.物体顶点的位置判断;2.求物体顶点在平面上的投影;3.将投影点转化到x-y平面上;4.构造凸包、求面积。
            1.物体顶点的位置判断。过光源做平行于已知平面的平面。有三种情况:
               (1)所有点都在光源平面的上方或在光源平面内,此时没有影子;
               (2)所有点都在光源下方,此时有一个可以计算的影子;
               (3)其他情况,影子的面积为无限大。
               至于判断方法,把点带入平面方程,根据结果的符号可以确定点在平面的哪一侧。注意给定的法向量方向,需要先转化到指向点的方向(如果点带入已知平面结果为负,则把abcd均乘以-1即可,完成平面翻转)。
            2.求解顶点在平面上的投影。根据直线的上的已知点(x,y,z)和法向量(dx,dy,dz)可以得到直线点集:(x,y,z)+t(dx,dy,dz)。然后将点带入平面ax+by+cz=d即可求出t,确定直线与平面交点(x+tdx,y+tdy,z+tdz)。在求解之前先进点的位置判断,可以避免通过t的正负来确定交点是否存在(因为实际上是射线,t>=0)。
            3.三维坐标转换。先绕z轴旋转(让y轴落在(a,b,0)方向上),然后绕x轴旋转即可(让z周落在(a,b,c)方向上)。下面给出坐标转换公式:
               绕z轴旋转:x`=xcosC-ysinC;y`=xsinC+ycosC;z`=z;
               绕x轴旋转:x`=x;y`=ycosA-zsinA;z`=ysinA+zcosA;
               绕z轴旋转:x`=zsinB+xcosB;y`=y;z`=zcosB-xsinB;
               推导过程并不复杂,因为有一维是不变的,可以直接看成平面的推导。
            4.构造凸包求面积。这里直接看成是二维的即可(z全相等,可以忽略)。然后将凸包分割成三角形求利用叉乘面积即可。
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
#include <vector>
using namespace std;
const int MAXN=110;
const double eps=1e-8;
const double PI=acos(-1.0);
inline int sign(double d)
{
    if(d > eps)return 1;
    else if(d < -eps)return -1;
    else return 0;
}
struct Point  //三维点坐标
{
    double x,y,z;
    Point(double xx=0,double yy=0,double zz=0):x(xx),y(yy),z(zz){}
    Point operator +(const Point p1)
    {
        return Point(x+p1.x,y+p1.y,z+p1.z);
    }
    Point operator -(const Point p1)
    {
        return Point(x-p1.x,y-p1.y,z-p1.z);
    }
    Point operator *(const Point p1)
    {
        return Point(y*p1.z-z*p1.y,z*p1.x-x*p1.z,x*p1.y-y*p1.x);
    }
    Point operator *(double d)
    {
        return Point(x*d,y*d,z*d);
    }
    Point operator /(double d)
    {
        return Point(x/d,y/d,z/d);
    }
    double operator ^(const Point p1)
    {
        return x*p1.x+y*p1.y+z*p1.z;
    }
    double getLen()
    {
        return sqrt(x*x+y*y+z*z);
    }
    void read()
    {
        scanf("%lf%lf%lf",&x,&y,&z);
    }
};
Point ps[MAXN];

inline Point get_point(Point st,Point ed,Point tp)//求点tp在直线st-ed上的垂足
{
    double t1=(tp-st)^(ed-st);
    double t2=(ed-st)^(ed-st);
    double t=t1/t2;
    Point ans=st + (ed-st)*t;
    return ans;
}
inline double dist(Point st,Point ed)
{
    return sqrt((ed-st)^(ed-st));
}
Point rotate(Point st,Point ed,Point tp,double A)//沿着直线st-ed旋转角度A,从ed往st看是逆时针
{
    Point root=get_point(st,ed,tp);
    Point e=(ed-st)/dist(ed,st);//单位向量
    Point r=tp-root;
    Point tmp=e*r;
    Point ans=r*cos(A)+tmp*sin(A)+root;
    return ans;
}
double a,b,c,d;
int n;
bool input()
{
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
    if(a==0&&b==0&&c==0&&d==0)return false;
    scanf("%d",&n);
    for(int i=0;i<=n;i++)
        ps[i].read();
    return true;
}


//****************************************************
//二维凸包部分
struct point
{
    double x,y;
};
point list[MAXN];
int stack[MAXN],top;
double cross(point p0,point p1,point p2)//p0p1 X p0p2
{
    return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
double dis(point p1,point p2)
{
    return sqrt((p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y));
}
bool cmp(point p1,point p2)
{
    double tmp=cross(list[0],p1,p2);
    if(sign(tmp)>0)return true;
    else if(sign(tmp)==0&&dis(list[0],p1)<dis(list[0],p2))return true;
    else return false;
}
void init(int n)
{
    point p0;
    int k;
    p0.x=list[0].x;
    p0.y=list[0].y;
    k=0;
    for(int i=1;i<n;i++)
    {
        if(p0.y>list[i].y || ((p0.y==list[i].y)&&(p0.x>list[i].x)))
        {
            p0.x=list[i].x;
            p0.y=list[i].y;
            k=i;
        }
    }
    list[k]=list[0];
    list[0]=p0;
    sort(list+1,list+n,cmp);
}
void graham(int n)
{
    int i;
    if(n==1)
    {
        top=0;
        stack[0]=0;
        return;
    }
    if(n==2)
    {
        top=1;
        stack[0]=0;
        stack[1]=1;
        return;
    }
    for(int i=0;i<n;i++)stack[i]=i;
    top=1;
    for(i=2;i<n;i++)
    {
        while(top>0 && cross(list[stack[top-1]],list[stack[top]],list[i])<=0)top--;
        top++;
        stack[top]=i;
    }
}
//****************************************************
void solve()
{
    if(n<=2)
    {
        printf("0.00\n");
        return;
    }
    Point p1(a,b,c);//法向量
    Point tp(0,0,0);
    Point e(0,0,1);
    if(sign(a)!=0)tp.x=d/a;
    else if(sign(b)!=0)tp.y=d/b;
    else if(sign(c)!=0)tp.x=d/c;
    ps[n+1]=tp;
    Point vec=p1*e;
    double A=(e^p1)/p1.getLen();
    A=acos(A);
    if(sign(A)!=0 && sign(A-PI)!=0)
    {
        Point p0(0,0,0);
        for(int i=0;i<=n+1;i++)
             ps[i]=rotate(p0,vec,ps[i],A);
    }
    for(int i=0;i<=n+1;i++)
        ps[i].z-=ps[n+1].z;
    if(sign(ps[n].z)<0)
    {
        for(int i=0;i<=n;i++)
            ps[i].z=-ps[i].z;
    }
    int cnt1=0;//记录ps[n]上方的点的个数
    int cnt2=0;//记录和ps[n]高度一样的点的个数
    for(int i=0;i<n;i++)
    {
        int k=sign(ps[i].z-ps[n].z);
        if(k>0)cnt1++;
        else if(k==0)cnt2++;
    }
    if(cnt1+cnt2==n)
    {
        printf("0.00\n");
        return;
    }
    if(cnt1>0 || cnt2>0)
    {
        printf("Infi\n");
        return;
    }
    for(int i=0;i<n;i++)
    {
        double u=ps[n].z/(ps[i].z-ps[n].z);
        list[i].x=ps[n].x+u*(ps[i].x-ps[n].x);
        list[i].y=ps[n].y+u*(ps[i].y-ps[n].y);
    }
    init(n);
    graham(n);
    double ans=0;
    for(int i=0;i<top;i++)
        ans+=(list[stack[i]].x * list[stack[i+1]].y - list[stack[i+1]].x*list[stack[i]].y);
    ans+=(list[stack[top]].x * list[stack[0]].y - list[stack[0]].x*list[stack[top]].y);
    printf("%.2lf\n",ans/2);
}
int main()
{
    while(input())solve();
    return 0;
}
两园面积交


#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;

const double eps = 1e-8;
const double PI = acos(-1.0);
struct Point
{
    double x,y;
    void input()
    {
        scanf("%lf%lf",&x,&y);
    }
};
double dist(Point a,Point b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
//两个圆的公共部分面积
double Area_of_overlap(Point c1,double r1,Point c2,double r2)
{
    double d = dist(c1,c2);
    if(r1 + r2 < d + eps)return 0;
    if(d < fabs(r1 - r2) + eps)
    {
        double r = min(r1,r2);
        return PI*r*r;
    }
    double x = (d*d + r1*r1 - r2*r2)/(2*d);
    double t1 = acos(x / r1);
    double t2 = acos((d - x)/r2);
    return r1*r1*t1 + r2*r2*t2 - d*r1*sin(t1);
}

int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    Point c1,c2;
    double r1,r2;
    int T;
    scanf("%d",&T);
    int iCase = 0;
    while(T--)
    {
        iCase ++;
        c1.input();
        scanf("%lf",&r1);
        c2.input();
        scanf("%lf",&r2);
        printf("Case %d: %.6lf\n",iCase,Area_of_overlap(c1,r1,c2,r2));
    }
    return 0;
}
给你两条直线,问你最短距离是多少,并且求出最短距离的点,需要求得点是唯一的。

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
struct Point
{
    double x,y,z;
    Point(double _x = 0,double _y = 0,double _z = 0)
    {
        x = _x;
        y = _y;
        z = _z;
    }
    Point operator +(const Point &b)const
    {
        return Point(x + b.x, y + b.y,z + b.z);
    }
    Point operator -(const Point &b)const
    {
        return Point(x - b.x, y - b.y,z - b.z);
    }
    Point operator /(double k)
    {
        return Point(x/k,y/k,z/k);
    }
    Point operator *(double k)
    {
        return Point(x*k,y*k,z*k);
    }
    double operator *(const Point &b)const
    {
        return x*b.x + y*b.y + z*b.z;
    }
    Point operator ^(const Point &b)const
    {
        return Point(y*b.z - z *b.y,z*b.x - x*b.z, x*b.y - y * b.x);
    }
    void input()
    {
        scanf("%lf%lf%lf",&x,&y,&z);
    }
    void output()
    {
        printf("%.6lf %.6lf %.6lf",x,y,z);
    }
};
double dis(Point a,Point b)
{
    return sqrt((a.x - b.x) * (a.x- b.x)  + (a.y - b.y) *(a.y - b.y) + (a.z - b.z)*(a.z - b.z));
}
double norm(Point a)
{
    return sqrt(a.x *a.x + a.y *a.y + a.z * a.z);
}
Point A,B,C,D;
Point mid1,mid2;
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    int T;
    scanf("%d",&T);
    while(T--)
    {
        A.input();
        B.input();
        C.input();
        D.input();
        Point p1 = (B-A);
        Point p2 = (D-C);
        Point p = p1^p2;

        double dd = (p*(A-C))/norm(p);
        dd = fabs(dd);
        printf("%.6lf\n",dd);

        double t1 = ( (C-A)^p2 )*(p1^p2);
        t1 /= norm(p1^p2)*norm(p1^p2);
        double t2 =  ( (C-A)^p1 )*(p1^p2);
        t2 /= norm(p1^p2)*norm(p1^p2);
        mid1 = A + (p1 * t1);
        mid2 = C + (p2 * t2);
        mid1.output();
        printf(" ");
        mid2.output();
        printf("\n");

    }
    return 0;
} /* 
题目意思:
给了平面上的n条线段:
让你在x轴上[0,L]的范围内找一个点作为圆心,画一个圆,这个圆不能把线段包含在里面去。
求最大的半径。
求解:
二分最终的答案,求解。
对于二分的半径值,对每条线段找出对于x位置上满足半径要求的段。
看n条线段有没有交集就可以了
 */

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;

const double eps = 1e-8;
const double inf = 1e5;
const double pi = acos(-1.0);
int sgn(double x)
{
    if(fabs(x) < eps)return 0;
    else if(x < 0)return -1;
    return 1;
}
struct Point
{
    double x,y;
    Point(){}
    Point(double _x,double _y)
    {
        x = _x;
        y = _y;
    }
    void input()
    {
        scanf("%lf%lf",&x,&y);
    }
    Point operator -(const Point &b)const
    {
        return Point(x-b.x,y-b.y);
    }
    double operator *(const Point &b)const
    {
        return x*b.x + y*b.y;
    }
    double operator ^(const Point &b)const
    {
        return x*b.y - y*b.x;
    }
    double len()
    {
        return hypot(x,y);
    }
    double distance(Point p)
    {
        return hypot(x-p.x,y-p.y);
    }
};
struct Line
{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e)
    {
        s = _s;
        e = _e;
    }
    void input()
    {
        s.input();
        e.input();
    }
    double length()
    {
        return s.distance(e);
    }
    double dispointtoline(Point p)
    {
        return fabs((p-s)^(e-s))/length();
    }
    double dispointtoseg(Point p)
    {
        if(sgn((p-s)*(e-s)) < 0 || sgn((p-e)*(s-e)) < 0)
            return min(p.distance(s),p.distance(e));
        return dispointtoline(p);
    }
};
const int MAXN = 2020;
Line line[MAXN];

double get(int i,double L)
{
    double l = 0, r = L;
    while(r - l >= eps)
    {
        double mid = (l + r)/2;
        double midmid = (r + mid)/2;
        if(line[i].dispointtoseg(Point(mid,0)) < line[i].dispointtoseg(Point(midmid,0)))
            r = midmid - eps;
        else l = mid + eps;
    }
    return l;
}
int n;
double L;

struct Node
{
    double a;
    int c;
    Node(double _a = 0,int _c = 0)
    {
        a = _a;
        c = _c;
    }
};
bool cmp(Node a,Node b)
{
    if(sgn(a.a - b.a) != 0)return a.a < b.a;
    else return a.c > b.c;
}

double bin1(int i,double l,double r,double R)
{
    while(r-l >= eps)
    {
        double mid = (l+r)/2;
        if(line[i].dispointtoseg(Point(mid,0)) <= R)
            r = mid - eps;
        else l = mid + eps;
    }
    return l;
}
double bin2(int i,double l,double r,double R)
{
    while(r-l >= eps)
    {
        double mid = (l+r)/2;
        if(line[i].dispointtoseg(Point(mid,0)) <= R)
            l = mid + eps;
        else r = mid - eps;
    }
    return l;
}

bool gao(double r)
{
    vector<Node>vec;
    for(int i = 0;i < n;i++)
    {
        double tmp = get(i,L);
        if(sgn(line[i].dispointtoseg(Point(tmp,0)) - r) >= 0)
        {
            vec.push_back(Node(0,1));
            vec.push_back(Node(L,-1));
            continue;
        }
        if(sgn(line[i].dispointtoseg(Point(0,0)) - r) >= 0)
        {
            double tt = bin1(i,0,tmp,r);
            vec.push_back(Node(0,1));
            vec.push_back(Node(tt,-1));
        }
        if(sgn(line[i].dispointtoseg(Point(L,0)) - r) >= 0)
        {
            double tt = bin2(i,tmp,L,r);
            vec.push_back(Node(tt,1));
            vec.push_back(Node(L,-1));
        }
    }
    sort(vec.begin(),vec.end(),cmp);
    int cnt = 0;
    for(int i = 0;i < vec.size();i++)
    {
        cnt += vec[i].c;
        if(cnt >= n)return true;
    }
    return false;
}

double solve()
{
    double l = 0, r = inf;
    while(r-l >= eps)
    {
        double mid = (l+r)/2;
        if(gao(mid))l = mid+eps;
        else r = mid - eps;
    }
    return l;
}

int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%lf",&n,&L);
        for(int i = 0;i < n;i++)
            line[i].input();
        printf("%.3lf\n",solve());
    }
    return 0;
}


  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值