【SHTSC2014】信号增幅仪

题目大意

给定平面直角坐标系中的n个点，有一个在某一特定的方向上伸长p倍的椭圆。问如何定下将椭圆的圆心，使所有点被覆盖？求这个椭圆的最小半短轴。

题解

引入

题目有两个难点：①这个椭圆在某一个特定的方向伸长。②这是个椭圆，且半长轴的长度为半短轴的p倍。
如果是一个正圆，那么很好求最小圆覆盖。那么能不能够利用题目给的条件来转换成最小圆覆盖问题呢？
这需要我们掌握图形变换的一些常识。
考虑一个任意的图形，如果图形上的所有点的横坐标扩大p倍，那么整个图形横向拉长p倍。
如果一个点x与原点O的连线l与x轴的方位角为α，那么将图逆时针旋转α后，x在新的x轴上。
根据以上两点，将n个点顺时针旋转α，横坐标除以p，就转化为最小圆覆盖问题了。

那么最小圆覆盖问题怎么做？

使用随机增量法。
也就是说，通过让圆的半径不断变大而使所有点被圆覆盖。
①随机地加入某一个点i，如果这个点在1~i-1的最小覆盖圆内/上，那么不管。
②如果在圆外，那么使i在新的圆上。
再寻找一个点j，如果j在圆外，那么圆必经点j，那么新圆就是以ij为直径的圆，现在要把前j-1个点加入回这个圆中。否则不管。
如果点k在当前的最小圆外，说明k一定在包括i,j和前k个点的最小覆盖圆上，我们可以直接得出这个圆了。否则不管。
③当点完全随机时，第i个点在圆外的概率为3/i，所以可以证明此算法的期望复杂度为$O(n)$。
这说明跑这个算法的时候，点必须随机。

常识

旋转公式：将$(x,y)$逆时针旋转$α$得到$(x',y')$
$x'=x*cos(α)-y*sin(α)$
$y'=x*sin(α)+y*cos(α)$
三点定圆：已知点$(x_1,y_1),(x_2,y_2),(x_3,y_3)$，求圆心$(x,y)$
根据圆心到顶点的距离相等，可列出：
$(x_1-x)*(x_1-x)+(y_1-y)*(y_1-y)=(x2_-x)*(x_2-x)+(y_2-y)*(y_2-y)$
$(x_2-x)*(x_2-x)+(y_2-y)*(y_2-y)=(x_3-x)*(x_3-x)+(y_3-y)*(y_3-y)$
化简，得：
$2x(x_2-x_1)+2y(y_2-y_1)=x_2^2+y_2^2-x_1^2-y_1^2$
$2x(x_3-x_2)+2y(y_3-y_2)=x_3^2+y_3^2-x_2^2-y_2^2$
我们需要用到已知的来推未知的，所以
令$a1=2(x_2-x_1),b1=2(y_2-y_1),c1=x_2^2+y_2^2-x_1^2-y_1^2$
$\ \ \ a2=2(x_3-x_2),b2=2(y_3-y_2),c2=x_3^2+y_3^2-x_2^2-y_2^2$
则$a1*x+b1*y=c1,a2*x+b2*y=c2$。
解这个二元一次方程组，得
$x=\frac{c1*b2-c2*b1}{a1*b2-a2*b1},y=\frac{a1*c2-a2*c1}{a1*b2-a2*b1}$

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 50010
#define DB double
#define eps 1e-9
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
struct note{
    DB x,y;
};note p[N],c;
int i,n;
DB r,ang,pw,X,Y,pi;
DB sgn(DB x){
    if(fabs(x)<eps)return 0;
    return x>0?1:-1; 
}
DB dis(note A,note B){
    return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)); 
}
note new_circle(note A,note B,note C){
    DB a1,a2,b1,b2,c1,c2;
    a1=2*(B.x-A.x);
    b1=2*(B.y-A.y);
    c1=B.x*B.x+B.y*B.y-A.x*A.x-A.y*A.y;
    a2=2*(C.x-B.x);
    b2=2*(C.y-B.y);
    c2=C.x*C.x+C.y*C.y-B.x*B.x-B.y*B.y;
    note circle;
    circle.x=(c1*b2-c2*b1)/(a1*b2-a2*b1);
    circle.y=(a1*c2-a2*c1)/(a1*b2-a2*b1);
    return circle;
}
void min_circle(){
    int i,j,k;
    c.x=c.y=r=0;
    fo(i,2,n)
        if(sgn(dis(p[i],c)-r)>0){
            c=p[i];
            r=0;
            fo(j,1,i-1)
                if(sgn(dis(p[j],c)-r)>0){
                    c.x=(p[i].x+p[j].x)/2.0;
                    c.y=(p[i].y+p[j].y)/2.0;
                    r=dis(p[i],p[j])/2.0;
                    fo(k,1,j-1)
                        if(sgn(dis(p[k],c)-r)>0){
                            c=new_circle(p[i],p[j],p[k]);
                            r=dis(p[i],c);
                        }
                }
        }
}
int main(){
    scanf("%d",&n);
    fo(i,1,n)scanf("%lf%lf",&p[i].x,&p[i].y);
    scanf("%lf%lf",&ang,&pw);
    pi=acos(-1); 
    ang=(360-ang)*pi/180.0;
    fo(i,1,n){
        X=p[i].x*cos(ang)-p[i].y*sin(ang);
        Y=p[i].x*sin(ang)+p[i].y*cos(ang);
        p[i].x=X/(DB)pw;
        p[i].y=Y;
    }
    random_shuffle(p+1,p+n+1);
    min_circle();
    printf("%.3lf",r);
    return 0;
}