计算几何模板

#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;}