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 long
int 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;
}