SGU 233 The Greatest Angle（求圆上一点使得ACB角度最大）

题目链接：http://acm.sgu.ru/problem.php?contest=0&problem=233

题意：给出一个圆O以及圆内两点AB。在圆上求一点C使得交ACB最大？

思路：求一个经过AB且与C相切的圆Q，切点就是C。设Q的圆心为P，则AP=CP,O、P、C共线，所以OP+PC=R。所以OP+AP=R。所以P在一个椭圆上。

int DB(double x)
{
    if(x>1e-8) return 1;
    if(x<-1e-8) return -1;
    return 0;
}


class point
{
public:
    double x,y;

    point(){}
    point(double _x,double _y)
    {
        x=_x;
        y=_y;
    }

    void get()
    {
        RD(x,y);
    }

    point operator+(point a)
    {
        return point(x+a.x,y+a.y);
    }

    point operator-(point a)
    {
        return point(x-a.x,y-a.y);
    }

    double operator*(point a)
    {
        return x*a.y-y*a.x;
    }

    point operator*(double t)
    {
        return point(x*t,y*t);
    }

    double operator^(point a)
    {
        return x*a.x+y*a.y;
    }

    double getLen()
    {
        return sqrt(x*x+y*y);
    }

    point operator/(double t)
    {
        return point(x/t,y/t);
    }

    double getAng(point a)
    {
        return atan2(*this*a,*this^a);
    }

    int operator==(const point &a) const
    {
        return DB(x-a.x)==0&&DB(y-a.y)==0;
    }

    int operator<(const point &a) const
    {
        if(DB(x-a.x)) return x<a.x;
        return y<a.y;
    }

    point adjust(double p)
    {
        point t=*this;
        double x=p/t.getLen();
        return point(x*t.x,x*t.y);
    }

    point zhuanNi(double ang)
    {
        return point(x*cos(ang)-y*sin(ang),x*sin(ang)+y*cos(ang));
    }

    point zhuanShun(double ang)
    {
        return point(x*cos(ang)+y*sin(ang),-x*sin(ang)+y*cos(ang));
    }

    point vertical()
    {
        return point(y,-x);
    }

    void print()
    {
        if(DB(x)==0) x=0;
        if(DB(y)==0) y=0;
        printf("%.6lf %.6lf\n",x,y);
    }
};


int C;
point O,A,B,ans,o,tempA,tempB,tempO;
double R,ang;

double test(point ans)
{
    ans=ans.zhuanNi(ang)+o;
    ans=(ans-tempO).adjust(R);
    ans=ans+tempO;
    return fabs((tempA-ans).getAng(tempB-ans));
}

void calCross(double a,double b,point A,point B)
{
    point M=(A+B)/2,N=M+(A-B).vertical();
    double aa=N.y-M.y,bb=M.x-N.x,cc=N.y*(N.x-M.x)-N.x*(N.y-M.y);
    point temp1,temp2;
    double x,y,det;
    if(DB(bb)==0)
    {
        x=-cc/aa;
        y=b*sqrt(1-x*x/(a*a));
        temp1=point(x,y);
        temp2=point(x,-y);
        det=test(temp1)-test(temp2);
        if(DB(det)==1) ans=temp1;
        else ans=temp2;
    }
    else
    {
        double a1=a*a*aa*aa+b*b*bb*bb;
        double b1=2*aa*cc*a*a;
        double c1=a*a*cc*cc-a*a*b*b*bb*bb;
        det=b1*b1-4*a1*c1;
        if(DB(det)<=0)
        {
            x=-b1/(2*a1);
            y=(-cc-aa*x)/bb;
            ans=point(x,y);
        }
        else
        {
            x=(-b1-sqrt(det))/(2*a1);
            y=(-cc-aa*x)/bb;
            temp1=point(x,y);
            x=(-b1+sqrt(det))/(2*a1);
            y=(-cc-aa*x)/bb;
            temp2=point(x,y);
            det=test(temp1)-test(temp2);
            if(DB(det)==1) ans=temp1;
            else ans=temp2;
        }
    }
}

int main()
{
    RD(C);
    while(C--)
    {
        O.get();
        RD(R);
        A.get();
        B.get();
        o=(A+O)/2;
        tempA=A;
        tempB=B;
        tempO=O;
        A=A-o;
        B=B-o;
        O=O-o;
        ang=point(100,0).getAng(O);
        A=A.zhuanShun(ang);
        B=B.zhuanShun(ang);
        O=O.zhuanShun(ang);
        double c=(A-O).getLen()/2;
        double a=R/2;
        double b=sqrt(a*a-c*c);
        calCross(a,b,A,B);
        ans=ans.zhuanNi(ang)+o;
        ans=(ans-tempO).adjust(R);
        ans=ans+tempO;
        ans.print();
    }
    return 0;
}