P4207 [NOI2005] 月下柠檬树

https://www.luogu.com.cn/problem/P4207

想象一下，投出的影子中圆的半径不会变，而圆的距离要除以 \(\tan(\alpha)\)

于是就变成了一堆圆放在平面上，然后每相邻的两个画一条公切线，求围出来的图形的面积
显然是对称的，于是只求一边

考虑怎么求公切线（这里显然需要的是外公切线），设一大一小两个圆直径分别为 \(R,r\)，圆心分别是 \(O_1,O_2\)，现在求的是 \(\vec{O_1O_2}\) 左边的公切线，从 \(O_1\) 向公切线做垂线设为 \(l\)，再从 \(O_2\) 向 \(l\) 做垂线
则此时得到了一个直角三角形，由于上面是个矩形，所以直角顶点（就是刚才的垂足）和 \(O_1\) 距离为 \(R-r\)
设顶点在 \(O_1\) 的那个角为 \(\alpha\)，则 \(\cos(\alpha)=\frac{R-r}{d}\)，其中 \(d\) 是两圆的圆心距
得到 \(\alpha\) 以后，用 \(\vec{O_1O_2}\) 逆时针旋转 \(\alpha-\frac{\pi}{2}\) 后得到的向量，与两圆分别求切点即可得到共切点

向量和圆求切点就随便按切点的角度求就行了
然后你发现只要将刚才的 \(\alpha-\frac{\pi}{2}\) 改成 \(\alpha+\frac{\pi}{2}\) 就能得到另一条外公切线，不过这里不需要

由于这些公切线可能会相交，直接分割求面积可能需要割割补补很麻烦，所以直接辛普森即可
求某点函数值的时候枚举每个圆、每条公切线分别求一编纵坐标，返回最大值
树的最上面那个点当作半径为 \(0\) 的圆处理；相邻的圆可能出现包含的情况，这样求出来交点坐标是 nan，那么没有点会在两个 nan 构成的线段上，所以没有影响不用特判
精度不能卡着 \(10^{-2}\) 开

写完这题以后想到的：其实求公切线直接设 \(l:ax+by+c=0\)，然后用点到直线距离公式 \(\dfrac{|ax_1+by_1+c|}{\sqrt(a^2+b^2)}=r_1\) 解方程组，比这简单的多。。。
关于那个绝对值直接枚举正负号就行

#include<cstdio>
#include<cmath>
#include<algorithm>
const double PI=acos(-1);
#define EPS 1e-10
struct Vector{
	double x,y;
	inline Vector(){x=y=0;}
	inline Vector(double _x,double _y){x=_x;y=_y;}
	inline double len()const{return __builtin_sqrt(x*x+y*y);}
	inline double len2()const{return x*x+y*y;}
	inline double atan()const{return __builtin_atan2(y,x);}
	inline void operator += (const Vector &a){x+=a.x;y+=a.y;}
	inline void operator -= (const Vector &a){x-=a.x;y-=a.y;}
	inline void operator *= (const double &a){x*=a;y*=a;}
	inline void operator /= (const double &a){x/=a;y/=a;}
	inline Vector operator + (const Vector &a)const{return Vector(x+a.x,y+a.y);}
	inline Vector operator - (const Vector &a)const{return Vector(x-a.x,y-a.y);}
	inline Vector operator * (const double &a)const{return Vector(x*a,y*a);}
	inline Vector operator / (const double &a)const{return Vector(x/a,y/a);}
	inline double operator * (const Vector &a)const{return x*a.x+y*a.y;}
	inline double operator ^ (const Vector &a)const{return x*a.y-y*a.x;}
};
struct Circle{
	Vector o;
	double r;
};
struct Line{
	Vector way,o;
	inline Line(){}
	inline Line(double a,double b,double c){//ax+by+c=0
		if(std::abs(a)>=EPS) o=Vector(-c/a,0);
		else o=Vector(0,-c/b);
		way=Vector(b,-a);
	}
	inline Line(const Vector &a,const Vector &b){o=a;way=b-a;}
};
inline Vector rotate(const Vector &v,double a){//counterclockwise, will not change the len
	double x=v.x,y=v.y,cos=__builtin_cos(a),sin=__builtin_sin(a);
	return Vector(x*cos-y*sin,x*sin+y*cos);
}
inline double getAngle(const Vector &a,const Vector &b){
	return b.atan()-a.atan();
}
inline Vector getIntersection(const Vector &v,const Circle &c,int type){
	double a=v.atan();
	return Vector(c.o.x+__builtin_cos(a+PI/2)*c.r*type,c.o.y+__builtin_sin(a+PI/2)*c.r*type);
}
inline Line getCommonTangent(const Circle &a,const Circle &b,Vector *pa=NULL,Vector *pb=NULL,int type=1){
//type: 1= on left of the vector (b-a) -1= on right
	Vector d=b.o-a.o;
	double ang=__builtin_acos((a.r-b.r)/d.len());
	d=rotate(d,type*(ang-PI/2));
	Vector _pa=getIntersection(d,a,type),_pb=getIntersection(d,b,type);
	if(pa) *pa=_pa;
	if(pb) *pb=_pb;
	return Line(_pa,_pb);
}
inline double getYByX(const Circle &c,double x){
	if(x>c.o.x+c.r) return 0;
	return __builtin_sqrt(c.r*c.r-(x-c.o.x)*(x-c.o.x));
}
inline double getYByX(const Line &l,double x){
	return l.o.y+l.way.y/l.way.x*(x-l.o.x);
}
#define N 100006
int n;
double alpha;
Circle c[N];
Line tangent[N];
Vector va[N],vb[N];
double h[N];
inline void init(){
	n++;
	double cot=1/__builtin_tan(alpha);
	c[1].o=Vector(0,0);
	for(int i=2;i<=n;i++) c[i].o=c[i-1].o+Vector(h[i-1]*cot,0);
	for(int i=1;i<n;i++) tangent[i]=getCommonTangent(c[i],c[i+1],vb+i,va+i+1);
}
inline double calc(double x){
	double ans=0;
	for(int i=1;i<=n;i++) ans=std::max(ans,getYByX(c[i],x));
	for(int i=1;i<n;i++)if(vb[i].x<=x&&x<=va[i+1].x) ans=std::max(ans,getYByX(tangent[i],x));
//		printf("x= %lf   ans= %lf\n",x,ans);
	return ans;
}
inline double simpson(double l,double r){
	return (r-l)*(calc(l)+4*calc((l+r)/2)+calc(r))/6;
}
double asr(double l,double r,double eps,double ans){
	double mid=(l+r)/2;
	double fl=simpson(l,mid),fr=simpson(mid,r);
	if(std::abs(fl+fr-ans)<=15*eps) return fl+fr+(fl+fr-ans)/15;
	return asr(l,mid,eps/2,fl)+asr(mid,r,eps/2,fr);
}
int main(){
	scanf("%d%lf",&n,&alpha);
	for(int i=0;i<=n;i++) scanf("%lf",h+i);
	for(int i=1;i<=n;i++) scanf("%lf",&c[i].r);
	init();
		#ifdef LOCAL
		puts("\n\n");
		for(int i=1;i<=n;i++) printf("%lf %lf r=%lf    %lf %lf  %lf %lf\n",c[i].o.x,c[i].o.y,c[i].r,va[i].x,va[i].y,vb[i].x,vb[i].y);
		puts("\n\n");
		#endif
	double l=1e10,r=-1e10;
	for(int i=1;i<=n;i++) l=std::min(l,c[i].o.x-c[i].r),r=std::max(r,c[i].o.x+c[i].r);
		#ifdef LOCAL
		printf("l= %lf   r= %lf\n\n\n",l,r);
		#endif
	printf("%.2lf\n",asr(l,r,1e-3,simpson(l,r))*2);
	return 0;
}