[BZOJ1502] [NOI2005] 月下柠檬树 [自适应Simpson]

Link
Luogu - https://www.luogu.org/problemnew/show/P4207
BZOJ - https://www.lydsy.com/JudgeOnline/problem.php?id=1502

---

很容易脑补出实际上要求的东西
具体一点就是一堆半圆和一堆直角梯形的面积并的两倍
有没有感觉很像CIRU或者2178？

---

Q：两圆的公切线？
A：
比方现在求 $(x-x_1{)}^2+y^2=r_1^2$ 和 $(x-x_2{)}^2+y^2=r_2^2$ 的公切线 $Ax+By+C=0$
则 $\frac{|A{x}_1+C|}{\sqrt{A^2+B^2}}=r_1$ 且 $\frac{|A{x}_2+C|}{\sqrt{A^2+B^2}}=r_2$
那么 $(A{x}_1+C{)}^2=r_1{A}^2+r_1{B}^2$ 且 $(A{x}_2+C{)}^2=r_2{A}^2+r_2{B}^2$
随便化一化
得到 $A^2(x_1^2-r_1-x_2^2+r_2)+B^2(r_1-r_2)=0$
并且
如果把 $Ax+By+C=0$ 写成 $y=kx+b$ 则 $A=k,B=-1,C=b$
由此
令 $B=-1$ 则 $A=\pm \sqrt{\frac{r_2-r_1}{x_{{}_1}^2-x_{{}_2}^2+r_2-r_1}}$
所以 $C=\mp \sqrt{r_1-\frac{(r_2-r_1)(x_{{}_1}^2-r_1)}{x_{{}_1}^2-x_{{}_2}^2+r_2-r_1}}$ （算的时候把 $r_1$ 也放进分子里会好一点大概吧）
然后直线方程就是 $y=Ax+C$

---

Q：求交点……
A：…………
对 $y=kx+b$ 和 $(x-x_0{)}^2+y^2=r^2$
$x^2-2x{x}_0+x_0^2+k^2{x}^2+2kbx+b^2=r^2$
$(k^2+1)x^2+2(kb-x_0)x+x_0^2+b^2-r^2=0$
直接套求根公式，反正已经知道 $\Delta =0$
得到 $x=\frac{-\text{−}-b^{\prime }\pm \sqrt{\Delta }}{2a^{\prime }}=\frac{x_0-\text{−}-kb}{k^2+1}$
$y=kx+b$

（求交点和求公切线都是我口胡的，如果感觉有不对的地方那肯定是我数学太差了
因为我最后并没有写这种方法
所以不保证正确性啊。）

---

然后暴力Simpson积分。这个我就不讲啦。

---

Q：细节很多。
A：计算几何是这样子的。是吗
实际上上面的计算过程完全可以被简化：注意到那个直角梯形可以被凑成矩形
而凑成矩形要凑上的另一半显然是个直角三角形 那就可以方便地……向量
然后就会变得巨肝胆
真好啊

---

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
#define db double
#define cdb const double&
const int MAXN = 504;
const double eps = 1e-8;
int n, tot;
double alpha, Mn=1145141919, Mx=-1145141919;
struct Circle
{
	double x, y, r;
	inline Circle(cdb a=0, cdb b=0, cdb c=0)
	{
		x=a,y=b,r=c;
	}
} C[MAXN];
struct Line
{
	double xa, ya, xb, yb, k, b;
	inline Line(cdb xaa=0, cdb yaa=0, cdb xbb=0, cdb ybb=0, cdb kk=0, cdb bb=0)
	{
		xa=xaa,ya=yaa,xb=xbb,yb=ybb,k=kk,b=bb;
	}
	inline bool GetEq()
	{
		if(xa==xb)return 1;
		k = (yb-ya)/(xb-xa);
		b = yb-k*xb;
		return 0;
	}
	inline double CalcY(cdb x)
	{
		return k*x+b;
	}
} L[MAXN];
inline double f(cdb x)
{
	db Ret = 0, Dist;
	for (register int i = 1; i <= n; ++i)
	{
		Dist = fabs(x - C[i].x);
		if (Dist - C[i].r > -eps) continue;
		Ret = max(Ret, 2.0 * sqrt(C[i].r*C[i].r-Dist*Dist));
	}
	for (register int i = 1; i <= tot; ++i)
	{
		if (x >= L[i].xa && x <= L[i].xb)
		{
			Ret = max(Ret, 2 * L[i].CalcY(x));
		}
	}
	return Ret;
}
inline double simpson(cdb l, cdb r)
{
	return (r-l)*(f(r)+f(l)+f((l+r)/2)*4)/6;
}
double divide(cdb L, cdb R, cdb Ans)
{
	db Mid = (L + R) / 2, l = simpson(L, Mid), r = simpson(Mid, R);
	if(fabs(l+r-Ans)<=eps*15) return l+r+(l+r-Ans)/15;
	return divide(L, Mid, l) + divide(Mid, R, r);
}
int main()
{
	scanf("%d%lf", &n, &alpha);
	{
		++n;
		register double h;
		for (register int i = 1; i <= n; ++i)
		{
			scanf("%lf", &h);
			C[i].x=C[i-1].x+h/tan(alpha);
			C[i].y=0;
		}
	}
	for (register int i = 1; i < n; ++i)
	{
		scanf("%lf", &C[i].r);
	}
	for (register int i = 1; i <= n; ++i)
	{
		Mn = min(Mn, C[i].x-C[i].r);
		Mx = max(Mx, C[i].x+C[i].r);
	}
	{
		register double Dist, Sn, Cs;
		for (register int i = 1; i < n; ++i)
		{
			Dist = C[i+1].x - C[i].x;
			if (Dist - fabs(C[i+1].r - C[i].r) < -eps) continue;
			Sn = (C[i].r - C[i+1].r) / Dist;
			Cs = sqrt(1.0 - Sn * Sn);
			L[++tot] = Line(C[i].x+Sn*C[i].r,Cs*C[i].r,C[i+1].x+Sn*C[i+1].r,Cs*C[i+1].r);
			L[tot].GetEq();
		}
	}
	printf("%.2lf", divide(Mn, Mx, simpson(Mn, Mx)));
	return 0;
}