[HDU 4855] Goddess （极角排序+三分）

HDU - 4855
借这道题学了下极角排序和三分求凸函数最大值
按所有左右切线，圆心的弧度值排序，从而分割出角度上的若干个区间
射线的角度在这些区间里变动时，每个割线的长度都是凸函数，叠加起来也是凸函数
所以可以对每个区间利用三分法求得最大值，再对这些最大值取 max 即为答案

1) 利用三分法可以求得单峰函数的最大值
2) 利用 atan2(y,x)可以比较方便地求出向量与 x轴正方向的夹角，取值范围 [-Pi,Pi]，从逆时针旋转 Pi 到顺时针旋转 Pi
但是排序后要插入两个额外的角， -Pi和 Pi，这是为了防止有圆穿过 x负半轴，而无法进行比较
3) 计算圆是否与射线相交时，可以利用旋转矩阵旋转射线和圆心，将射线与 x正半轴对齐
此时距离即为旋转后圆心的 y坐标，x坐标若为负数，则表明圆在射线另一侧

#if _WIN32||_WIN64#define lld I64d#define llu I64u#endif#include <cstdio>#include <iostream>#include <cstdlib>#include <cstring>#include <algorithm>#include <cmath>#include <map>#include <set>using namespace std;typedef pair<int,int> Pii;#define MST(a,b) memset(a,b,sizeof(a))#define CLR(a) MST(a,0)#define LL long long#define ULL unsigned long longint maxx(int a,int b){return a>b?a:b;}int minn(int a,int b){return a<b?a:b;}int abss(int a){return a<0?(-a):a;}const double eps=1e-11,PI=acos(-1.0);struct Circle{    double x,y,r;};struct Vector{    double x,y;    Vector(double tx=0, double ty=0):x(tx),y(ty){};};double VLen(Vector u){return sqrt(u.x*u.x+u.y*u.y);}double operator ^ (Vector u, Vector v){return u.x*v.y-u.y*v.x;}int N,ap;Circle inpt[250];double angl[800];double calc(double);int main(){    while(~scanf("%d", &N))    {        ap=0;        angl[++ap]=PI;        angl[++ap]=-PI; // to calculate the circle go across the x-asis        for(int i=1; i<=N; i++)        {            scanf("%lf%lf%lf", &inpt[i].x, &inpt[i].y, &inpt[i].r);            double acent=atan2(inpt[i].y,inpt[i].x);            double distc=VLen(Vector(inpt[i].x,inpt[i].y));            double delta=asin(inpt[i].r/distc);            double ang1=acent-delta,ang2=acent+delta;            if(ang1<-PI&&fabs(ang1+PI)>eps) ang1+=2*PI;            if(ang2>PI&&fabs(ang2-PI)>eps) ang2-=2*PI;            angl[++ap]=acent;            angl[++ap]=ang1;            angl[++ap]=ang2;        }        sort(angl+1,angl+1+ap);        ap=unique(angl+1,angl+1+ap)-angl-1;        double ans=0.0;//      for(int i=1; i<=ap; i++) printf("%lf\n", angl[i]);        for(int i=1; i<ap; i++)        {            double l=angl[i],r=angl[i+1];            double res=0;            while(l<r&&fabs(l-r)>eps)            {//              printf("l:%f r:%f res:%f\n", l, r, res);system("pause");                double ll=l+(r-l)/3;                double rr=l+2*(r-l)/3;                double resl=calc(ll);                double resr=calc(rr);                if(resl>resr)                {                    res=resl;                    r=rr-eps;                }                else                {                    res=resr;                    l=ll+eps;                }                if(res>ans&&fabs(res-ans)>eps) ans=res;            }        }        printf("%.8f\n", ans);    }    return 0;}double calc(double ang){    double m11=cos(ang),m12=sin(ang),m21=-sin(ang),m22=cos(ang);// use matrix to rotate the center of circles    double temp=0.0;    for(int i=1; i<=N; i++)    {        double x=m11*inpt[i].x+m12*inpt[i].y; // has the same direction with the laser        double y=fabs(m21*inpt[i].x+m22*inpt[i].y); // the distance to the laser        if((y<inpt[i].r&&fabs(y-inpt[i].r)>eps)&&(x>0))        {            temp+=2*sqrt(inpt[i].r*inpt[i].r-y*y);        }    }    return temp;}