sgu233：The Greatest Angle(计算几何)

题目大意：
$~~~~~~$给出一个圆$O$和严格在圆内的两点$A,B$，在圆上求一点$C$使得$\angle ACB$最大。

分析：
$~~~~~~$首先我们得知道$A,B,C$所构成的外接圆与圆$O$相切。如果$O$到$AB$的垂线垂足位于$AB$中点使得有两个切点的话我们考虑离$AB$近的那一个。(不知道请走隔壁数学组)
$~~~~~~$知道这个之后就是乱搞了…(本蒟蒻参考的刘汝佳的做法)
http://wenku.baidu.com/link?url=DYfmAczQnMSS2aXUkQdLQ7bkNl7yLq05nfLElMDu249JWIqHhytV3kQg_T36tw3dCOBJeYYIS3Yj-9Y-QfSYQuTGQT9ic9g0FPeYcAtxQ8O

AC code：

#include <cstdio>
#include <cmath>
#include <algorithm>
typedef double DB;
#define ONLINE_JUDGE
using namespace std;

const DB eps = 1e-8;
const DB pi = acos(-1.0);

struct pot
{
    DB x, y;
    pot(DB x = 0, DB y = 0):x(x),y(y){}
    void read() {scanf("%lf%lf", &x, &y);}
    void print() {printf("%.6lf %.6lf", x, y);}
    DB length() {return sqrt(x*x+y*y);}
    friend pot operator + (const pot &a, const pot &b) {return pot(a.x+b.x, a.y+b.y);} 
    friend pot operator - (const pot &a, const pot &b) {return pot(a.x-b.x, a.y-b.y);}
    friend pot operator * (const pot &a, DB k) {return pot(a.x*k, a.y*k);}
    friend pot operator / (const pot &a, DB k) {return pot(a.x/k, a.y/k);}
    pot unit() {return pot(x, y)/length();}
}O, A, B, C, vCP, vAC, oP, D;
DB R, r;
DB AC, CO, cosPCO, OCA;
DB P, Q;
DB a, b, c, deita;
DB PC;
DB p, q;

DB sqr(DB x) {return x*x;}
DB dis2(const pot &a, const pot &b) {return sqr(a.x-b.x)+sqr(a.y-b.y);}
DB dis(const pot &a, const pot &b) {return sqrt(dis2(a, b));}
bool is_zero(DB x)
{
    if(x >= -eps && x <= eps) return true;
    return false;
}

DB calc_angle(const pot &a, const pot &b, const pot &c)
{
    DB ac2 = dis2(a, c), bc2 = dis2(b, c), ab2 = dis2(a, b);
    if(is_zero(bc2) || is_zero(ab2)) return 0;
    return acos((ac2-bc2-ab2)/2/sqrt(bc2)/sqrt(ab2));
}

bool is_correct(const pot &P)
{
    if(is_zero(dis(P, O)+r-R)) return true;
    return false;   
}

int main()
{
    #ifndef ONLINE_JUDGE
    freopen("sgu233.in", "r", stdin);
    freopen("sgu233.out", "w", stdout);
    #endif

    int test;
    scanf("%d", &test);
    while(test--)
    {
        O.read(), scanf("%lf", &R);
        A.read(), B.read(), C = (A+B)/2;
        AC = dis(A, C), CO = dis(C, O);
        OCA = calc_angle(A, C, O);
        cosPCO = cos(OCA+pi/2);
        P = sqr(R)-dis2(O, C)+dis2(A, C);
        Q = 2*CO*cosPCO;
        a = sqr(Q)-4*sqr(R), b = 2*P*Q, c = sqr(P)-4*dis2(A, C)*sqr(R);
        if(is_zero(a)) PC = -c/b;
        else
        {
            deita = sqr(b)-4*a*c;
            if(is_zero(deita)) deita = 0;
            p = -b/2/a, q = sqrt(deita)/2/a;
            PC = p+q;
            if(sqr(p+q) > sqr(p-q)) PC = p-q;
        }
        r = sqrt(sqr(PC)+dis2(A, C));
        vAC = C-A, vAC = vAC/vAC.length();
        vCP = pot(-vAC.y, vAC.x), PC = fabs(PC);
        if(is_correct(C+vCP*PC)) oP = C+vCP*PC;
        else oP = C-vCP*PC;
        D = (oP-O).unit()*R+O;
        D.print();puts("");
    }

    #ifndef ONLINE_JUDGE
    fclose(stdin);
    fclose(stdout);
    #endif
    return 0;
}