P4288 [SHOI2014]信号增幅仪最小圆覆盖

传送门

题意：给出一些点，让后让你用一个 $ba=1p\frac{b}{a}=\frac{1}{p}$ 的且绕x正方向逆时针转动了a°的椭圆覆盖，求最小的b。

由于题目中椭圆是绕x正方向逆时针旋转了a°，为了方便写出方程，我们可以先把它顺时针转回来，这样就得到椭圆方程 $x2a2+y2b2=1\frac {x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1$ ，因为 $ba=1p\frac{b}{a}=\frac{1}{p}$ ，我们可以让 $a = p$ , $b = 1$ ，这样方程就变成了 $x2p2+y21=1\frac {x^{2}}{p^{2}}+\frac{y^{2}}{1}=1$ ，用高中的知识我们可以做如下换元 $x′=xpx'=\frac{x}{p}$ ， $y = y$ 。这样就变成了一个圆 $x^2+y^2=1$ ，问题转换成用半径最小的圆覆盖掉这些点，这样直接套板子就好啦。

我这里random_shuffle了两下，因为我一下过不去，应该是被卡精度了。

//#pragma GCC optimize(2)
#include<cstdio>
#include<iostream>
#include<string>
#include<cstring>
#include<map>
#include<cmath>
#include<cctype>
#include<vector>
#include<set>
#include<queue>
#include<algorithm>
#include<sstream>
#include<ctime>
#include<cstdlib>
#define X first
#define Y second
#define L (u<<1)
#define R (u<<1|1)
#define pb push_back
#define mk make_pair
#define Mid (tr[u].l+tr[u].r>>1)
#define Len(u) (tr[u].r-tr[u].l+1)
#define random(a,b) ((a)+rand()%((b)-(a)+1))
#define db puts("---")
using namespace std;//void rd_cre() { freopen("d://dp//data.txt","w",stdout); srand(time(NULL)); }
//void rd_ac() { freopen("d://dp//data.txt","r",stdin); freopen("d://dp//AC.txt","w",stdout); }
//void rd_wa() { freopen("d://dp//data.txt","r",stdin); freopen("d://dp//WA.txt","w",stdout); }typedef long long LL;
typedef unsigned long long ULL;
typedef pair<int,int> PII;const int N=100010,mod=1e9+7,INF=0x3f3f3f3f;// `计算几何模板`
const double eps = 1e-15;
const double inf = 1e15;
const double pi = acos(-1.0);
const int maxp = 100010;
//`Compares a double to zero`
int sgn(double x){if(fabs(x) < eps)return 0;if(x < 0)return -1;else return 1;
}
//square of a double
inline double sqr(double x){return x*x;}
/** Point* Point()               - Empty constructor* Point(double _x,double _y)  - constructor* input()             - double input* output()            - %.2f output* operator ==         - compares x and y* operator <          - compares first by x, then by y* operator -          - return new Point after subtracting curresponging x and y* operator ^          - cross product of 2d points* operator *          - dot product* len()               - gives length from origin* len2()              - gives square of length from origin* distance(Point p)   - gives distance from p* operator + Point b  - returns new Point after adding curresponging x and y* operator * double k - returns new Point after multiplieing x and y by k* operator / double k - returns new Point after divideing x and y by k* rad(Point a,Point b)- returns the angle of Point a and Point b from this Point* trunc(double r)     - return Point that if truncated the distance from center to r* rotleft()           - returns 90 degree ccw rotated point* rotright()          - returns 90 degree cw rotated point* rotate(Point p,double angle) - returns Point after rotateing the Point centering at p by angle radian ccw*/
struct Point{double x,y;Point(){}Point(double _x,double _y){x = _x;y = _y;}void input(){scanf("%lf%lf",&x,&y);}void output(){printf("%.2f %.2f\n",x,y);}bool operator == (Point b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;}bool operator < (Point b)const{return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;}Point operator -(const Point &b)const{return Point(x-b.x,y-b.y);}//叉积double operator ^(const Point &b)const{return x*b.y - y*b.x;}//点积double operator *(const Point &b)const{return x*b.x + y*b.y;}//返回长度double len(){return hypot(x,y);//库函数}//返回长度的平方double len2(){return x*x + y*y;}//返回两点的距离double distance(Point p){return hypot(x-p.x,y-p.y);}Point operator +(const Point &b)const{return Point(x+b.x,y+b.y);}Point operator *(const double &k)const{return Point(x*k,y*k);}Point operator /(const double &k)const{return Point(x/k,y/k);}//`计算pa  和  pb 的夹角`//`就是求这个点看a,b 所成的夹角`//`测试 LightOJ`double rad(Point a,Point b){Point p = *this;return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));}//`化为长度为r的向量`Point trunc(double r){double l = len();if(!sgn(l))return *this;r /= l;return Point(x*r,y*r);}//`逆时针旋转90度`Point rotleft(){return Point(-y,x);}//`顺时针旋转90度`Point rotright(){return Point(y,-x);}//`绕着p点逆时针旋转angle`Point rotate(Point p,double angle){Point v = (*this) - p;double c = cos(angle), s = sin(angle);return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);}
};
/** Stores two points* Line()                         - Empty constructor* Line(Point _s,Point _e)        - Line through _s and _e* operator ==                    - checks if two points are same* Line(Point p,double angle)     - one end p , another end at angle degree* Line(double a,double b,double c) - Line of equation ax + by + c = 0* input()                        - inputs s and e* adjust()                       - orders in such a way that s < e* length()                       - distance of se* angle()                        - return 0 <= angle < pi* relation(Point p)              - 3 if point is on line*                                  1 if point on the left of line*                                  2 if point on the right of line* pointonseg(double p)           - return true if point on segment* parallel(Line v)               - return true if they are parallel* segcrossseg(Line v)            - returns 0 if does not intersect*                                  returns 1 if non-standard intersection*                                  returns 2 if intersects* linecrossseg(Line v)           - line and seg* linecrossline(Line v)          - 0 if parallel*                                  1 if coincides*                                  2 if intersects* crosspoint(Line v)             - returns intersection point* dispointtoline(Point p)        - distance from point p to the line* dispointtoseg(Point p)         - distance from p to the segment* dissegtoseg(Line v)            - distance of two segment* lineprog(Point p)              - returns projected point p on se line* symmetrypoint(Point p)         - returns reflection point of p over se**/
struct Line{Point s,e;Line(){}Line(Point _s,Point _e){s = _s;e = _e;}bool operator ==(Line v){return (s == v.s)&&(e == v.e);}//`根据一个点和倾斜角angle确定直线,0<=angle<pi`Line(Point p,double angle){s = p;if(sgn(angle-pi/2) == 0){e = (s + Point(0,1));}else{e = (s + Point(1,tan(angle)));}}//ax+by+c=0Line(double a,double b,double c){if(sgn(a) == 0){s = Point(0,-c/b);e = Point(1,-c/b);}else if(sgn(b) == 0){s = Point(-c/a,0);e = Point(-c/a,1);}else{s = Point(0,-c/b);e = Point(1,(-c-a)/b);}}void input(){s.input();e.input();}void adjust(){if(e < s)swap(s,e);}//求线段长度double length(){return s.distance(e);}//`返回直线倾斜角 0<=angle<pi`double angle(){double k = atan2(e.y-s.y,e.x-s.x);if(sgn(k) < 0)k += pi;if(sgn(k-pi) == 0)k -= pi;return k;}//`点和直线关系`//`1  在左侧`//`2  在右侧`//`3  在直线上`int relation(Point p){int c = sgn((p-s)^(e-s));if(c < 0)return 1;else if(c > 0)return 2;else return 3;}// 点在线段上的判断bool pointonseg(Point p){return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;}//`两向量平行(对应直线平行或重合)`bool parallel(Line v){return sgn((e-s)^(v.e-v.s)) == 0;}//`两线段相交判断`//`2 规范相交`//`1 非规范相交`//`0 不相交`int segcrossseg(Line v){int d1 = sgn((e-s)^(v.s-s));int d2 = sgn((e-s)^(v.e-s));int d3 = sgn((v.e-v.s)^(s-v.s));int d4 = sgn((v.e-v.s)^(e-v.s));if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||(d4==0 && sgn((e-v.s)*(e-v.e))<=0);}//`直线和线段相交判断`//`-*this line   -v seg`//`2 规范相交`//`1 非规范相交`//`0 不相交`int linecrossseg(Line v){int d1 = sgn((e-s)^(v.s-s));int d2 = sgn((e-s)^(v.e-s));if((d1^d2)==-2) return 2;return (d1==0||d2==0);}//`两直线关系`//`0 平行`//`1 重合`//`2 相交`int linecrossline(Line v){if((*this).parallel(v))return v.relation(s)==3;return 2;}//`求两直线的交点`//`要保证两直线不平行或重合`Point crosspoint(Line v){double a1 = (v.e-v.s)^(s-v.s);double a2 = (v.e-v.s)^(e-v.s);return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));}//点到直线的距离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);}//`返回线段到线段的距离`//`前提是两线段不相交，相交距离就是0了`double dissegtoseg(Line v){return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));}//`返回点p在直线上的投影`Point lineprog(Point p){return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );}//`返回点p关于直线的对称点`Point symmetrypoint(Point p){Point q = lineprog(p);return Point(2*q.x-p.x,2*q.y-p.y);}
};
//圆
struct circle{Point p;//圆心double r;//半径circle(){}circle(Point _p,double _r){p = _p;r = _r;}circle(double x,double y,double _r){p = Point(x,y);r = _r;}//`三角形的外接圆`//`需要Point的+ /  rotate()  以及Line的crosspoint()`//`利用两条边的中垂线得到圆心`//`测试：UVA12304`circle(Point a,Point b,Point c){Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));p = u.crosspoint(v);r = p.distance(a);}//`三角形的内切圆`//`参数bool t没有作用，只是为了和上面外接圆函数区别`//`测试：UVA12304`circle(Point a,Point b,Point c,bool t){Line u,v;double m = atan2(b.y-a.y,b.x-a.x), n = atan2(c.y-a.y,c.x-a.x);u.s = a;u.e = u.s + Point(cos((n+m)/2),sin((n+m)/2));v.s = b;m = atan2(a.y-b.y,a.x-b.x) , n = atan2(c.y-b.y,c.x-b.x);v.e = v.s + Point(cos((n+m)/2),sin((n+m)/2));p = u.crosspoint(v);r = Line(a,b).dispointtoseg(p);}//输入void input(){p.input();scanf("%lf",&r);}//输出void output(){printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r);}bool operator == (circle v){return (p==v.p) && sgn(r-v.r)==0;}bool operator < (circle v)const{return ((p<v.p)||((p==v.p)&&sgn(r-v.r)<0));}//面积double area(){return pi*r*r;}//周长double circumference(){return 2*pi*r;}//`点和圆的关系`//`0 圆外`//`1 圆上`//`2 圆内`int relation(Point b){double dst = b.distance(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r)==0)return 1;return 0;}//`线段和圆的关系`//`比较的是圆心到线段的距离和半径的关系`int relationseg(Line v){double dst = v.dispointtoseg(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r) == 0)return 1;return 0;}//`直线和圆的关系`//`比较的是圆心到直线的距离和半径的关系`int relationline(Line v){double dst = v.dispointtoline(p);if(sgn(dst-r) < 0)return 2;else if(sgn(dst-r) == 0)return 1;return 0;}//`两圆的关系`//`5 相离`//`4 外切`//`3 相交`//`2 内切`//`1 内含`//`需要Point的distance`//`测试：UVA12304`int relationcircle(circle v){double d = p.distance(v.p);if(sgn(d-r-v.r) > 0)return 5;if(sgn(d-r-v.r) == 0)return 4;double l = fabs(r-v.r);if(sgn(d-r-v.r)<0 && sgn(d-l)>0)return 3;if(sgn(d-l)==0)return 2;if(sgn(d-l)<0)return 1;}//`求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点`//`需要relationcircle`//`测试：UVA12304`int pointcrosscircle(circle v,Point &p1,Point &p2){int rel = relationcircle(v);if(rel == 1 || rel == 5)return 0;double d = p.distance(v.p);double l = (d*d+r*r-v.r*v.r)/(2*d);double h = sqrt(r*r-l*l);Point tmp = p + (v.p-p).trunc(l);p1 = tmp + ((v.p-p).rotleft().trunc(h));p2 = tmp + ((v.p-p).rotright().trunc(h));if(rel == 2 || rel == 4)return 1;return 2;}//`求直线和圆的交点，返回交点个数`int pointcrossline(Line v,Point &p1,Point &p2){if(!(*this).relationline(v))return 0;Point a = v.lineprog(p);double d = v.dispointtoline(p);d = sqrt(r*r-d*d);if(sgn(d) == 0){p1 = a;p2 = a;return 1;}p1 = a + (v.e-v.s).trunc(d);p2 = a - (v.e-v.s).trunc(d);return 2;}//`得到过a,b两点，半径为r1的两个圆`int gercircle(Point a,Point b,double r1,circle &c1,circle &c2){circle x(a,r1),y(b,r1);int t = x.pointcrosscircle(y,c1.p,c2.p);if(!t)return 0;c1.r = c2.r = r;return t;}//`得到与直线u相切，过点q,半径为r1的圆`//`测试：UVA12304`int getcircle(Line u,Point q,double r1,circle &c1,circle &c2){double dis = u.dispointtoline(q);if(sgn(dis-r1*2)>0)return 0;if(sgn(dis) == 0){c1.p = q + ((u.e-u.s).rotleft().trunc(r1));c2.p = q + ((u.e-u.s).rotright().trunc(r1));c1.r = c2.r = r1;return 2;}Line u1 = Line((u.s + (u.e-u.s).rotleft().trunc(r1)),(u.e + (u.e-u.s).rotleft().trunc(r1)));Line u2 = Line((u.s + (u.e-u.s).rotright().trunc(r1)),(u.e + (u.e-u.s).rotright().trunc(r1)));circle cc = circle(q,r1);Point p1,p2;if(!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2);c1 = circle(p1,r1);if(p1 == p2){c2 = c1;return 1;}c2 = circle(p2,r1);return 2;}//`同时与直线u,v相切，半径为r1的圆`//`测试：UVA12304`int getcircle(Line u,Line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4){if(u.parallel(v))return 0;//两直线平行Line u1 = Line(u.s + (u.e-u.s).rotleft().trunc(r1),u.e + (u.e-u.s).rotleft().trunc(r1));Line u2 = Line(u.s + (u.e-u.s).rotright().trunc(r1),u.e + (u.e-u.s).rotright().trunc(r1));Line v1 = Line(v.s + (v.e-v.s).rotleft().trunc(r1),v.e + (v.e-v.s).rotleft().trunc(r1));Line v2 = Line(v.s + (v.e-v.s).rotright().trunc(r1),v.e + (v.e-v.s).rotright().trunc(r1));c1.r = c2.r = c3.r = c4.r = r1;c1.p = u1.crosspoint(v1);c2.p = u1.crosspoint(v2);c3.p = u2.crosspoint(v1);c4.p = u2.crosspoint(v2);return 4;}//`同时与不相交圆cx,cy相切，半径为r1的圆`//`测试：UVA12304`int getcircle(circle cx,circle cy,double r1,circle &c1,circle &c2){circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r);int t = x.pointcrosscircle(y,c1.p,c2.p);if(!t)return 0;c1.r = c2.r = r1;return t;}//`过一点作圆的切线(先判断点和圆的关系)`//`测试：UVA12304`int tangentline(Point q,Line &u,Line &v){int x = relation(q);if(x == 2)return 0;if(x == 1){u = Line(q,q + (q-p).rotleft());v = u;return 1;}double d = p.distance(q);double l = r*r/d;double h = sqrt(r*r-l*l);u = Line(q,p + ((q-p).trunc(l) + (q-p).rotleft().trunc(h)));v = Line(q,p + ((q-p).trunc(l) + (q-p).rotright().trunc(h)));return 2;}//`求两圆相交的面积`double areacircle(circle v){int rel = relationcircle(v);if(rel >= 4)return 0.0;if(rel <= 2)return min(area(),v.area());double d = p.distance(v.p);double hf = (r+v.r+d)/2.0;double ss = 2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));double a1 = acos((r*r+d*d-v.r*v.r)/(2.0*r*d));a1 = a1*r*r;double a2 = acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d));a2 = a2*v.r*v.r;return a1+a2-ss;}//`求圆和三角形pab的相交面积`//`测试：POJ3675 HDU3982 HDU2892`double areatriangle(Point a,Point b){if(sgn((p-a)^(p-b)) == 0)return 0.0;Point q[5];int len = 0;q[len++] = a;Line l(a,b);Point p1,p2;if(pointcrossline(l,q[1],q[2])==2){if(sgn((a-q[1])*(b-q[1]))<0)q[len++] = q[1];if(sgn((a-q[2])*(b-q[2]))<0)q[len++] = q[2];}q[len++] = b;if(len == 4 && sgn((q[0]-q[1])*(q[2]-q[1]))>0)swap(q[1],q[2]);double res = 0;for(int i = 0;i < len-1;i++){if(relation(q[i])==0||relation(q[i+1])==0){double arg = p.rad(q[i],q[i+1]);res += r*r*arg/2.0;}else{res += fabs((q[i]-p)^(q[i+1]-p))/2.0;}}return res;}
};double get_dis(Point a,Point b)
{double dx=a.x-b.x;double dy=a.y-b.y;return sqrt(dx*dx+dy*dy);
}int n;
Point q[N];Point rotate(Point a, double b)
{return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}int main()
{
//	ios::sync_with_stdio(false);
//	cin.tie(0);scanf("%d",&n);for(int i=0;i<n;i++) scanf("%lf%lf",&q[i].x,&q[i].y);double p,a; scanf("%lf%lf",&a,&p);a=a/180*pi;for(int i=0;i<n;i++){Point cmp(0,0);//q[i].rotate(cmp,a);q[i]=q[i].rotate(cmp,-a);q[i].x/=p;}random_shuffle(q,q+n);random_shuffle(q,q+n);circle c(q[0],0);for(int i=1;i<n;i++)if(sgn(get_dis(c.p,q[i])-c.r)>0){c={q[i],0};for(int j=0;j<i;j++)if(sgn(get_dis(c.p,q[j])-c.r)>0){c={(q[i]+q[j])/2,get_dis(q[i],q[j])/2};for(int k=0;k<j;k++)if(sgn(get_dis(c.p,q[k])-c.r)>0){circle cc(q[i],q[j],q[k]);c=cc;}}}printf("%.3f\n",c.r);return 0;
}
/**/