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
#include
#include
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


						
							

						
							
							
			题解-洛谷计算几何-计算几何

P4207 [NOI2005] 月下柠檬树

相关

标签