2014-2015 ACM-ICPC Northeastern European Regional Contest (NEERC 14) 題解
Tag
A（構造）
B（貪心）
D（定積分，座標系）
E（構造，圖） 但是不會
F（bitset模擬）
I（思維，動態規劃）
J（搜尋）
K（算術）
D. Damage Assessment
題意：兩個球缺和一個圓柱組成了一個傾斜擺放的幾何體，求其所盛水的體積。
思路
將原問題分為圓柱體和球缺兩個獨立部分，分別對水的截面面積做定積分得到體積。
求定積分的數值解用到 自適應Simpson演算法
先考慮一個圓上面積問題：
（1）子問題
如圖，已知 \(R\) 和 \(D(0\leq D \leq R)\)，求左邊弓形面積。
顯然面積是 \(F_1(R,D)=r^2\arccos(1-\frac{D}{R})-\sqrt{2RD-D^2}(R-D)\)
（2）圓柱內水的體積
設傾斜角 \(\omega =\arcsin(\frac{t}{l})\)，先定義圓柱側面長方形的左上端點的高度 \(h_1=\frac{d}{\cos \omega}\)，
則答案為 \(\int_0^{\min(h,t+h_1)}F_2(x)dx\)，其中 \(F_2(x)\) 是圖中藍色水平面截圓柱體的面積，\(x\) 是水平面的高度。
具體來說，考慮將一個 \(F_2(x)\) 截出來的橢圓的長軸乘以一個 \(\sin \omega\) ，得到以 \(d\) 為直徑的圓（其面積也乘了個 \(\sin \omega\)
（3）球缺內水的體積
以左邊的球缺為例，不妨設水介面 \(h<h_1\)，換一個與（2）不同的座標系：
紫色是水介面 \(h\) ，紅色箭頭的指向就是積分方向（\(x\)正半軸方向），紅箭頭的根部就是原點，藍色是其中一個截面，對於固定的一個截面（知道這條藍線的\(x\)座標），可以根據座標系的知識（比如算一下直線和1/4圓交點）算出這條藍線的高度，進而轉化為 （1） 的情形。
（4）一些細節
本題的資料是帶量綱的，注意讀入和輸出的單位換算。（可以直接把所有輸入數除以100）
注意這兩個邊界：
-
當 \(t=0\)，\(\sin \omega=0\)，計算（2） \(F_2(x)\) 的時候直接返回一個長方形面積就行
-
當 \(t=l\)，\(\cos \omega=0\)，此時左邊的球缺是完整的，右邊的球缺（\(h>t\)時）的截面是個圓，需特判
更多實現細節可以看程式碼
#include<bits/stdc++.h>
using namespace std;
const double PI = 3.1415926535897932384626433;
struct Integration{ // 自適應Simpson演算法求定積分
typedef double Func(double);
const Func*f;Integration(const Func&g):f(&g){}
typedef pair<double,double> pdd;
pdd add(pdd a,pdd b)const{
double mid=(a.first+b.first)/2;
return (pdd){mid,f(mid)};
}
#define simpson(p1,p2,p3) (((p3).first-(p1).first)*((p1).second+4*(p2).second+(p3).second)/6)
double find(pdd p1,pdd p3,pdd p5,double EPS,double res,int dep)const{
pdd p2=add(p1,p3),p4=add(p3,p5);
double fl=simpson(p1,p2,p3),fr=simpson(p3,p4,p5),d=(fl+fr-res)/15;
if(abs(d)<=EPS&&dep<0)return fl+fr+d;
return find(p1,p2,p3,EPS/2,fl,dep-1)+find(p3,p4,p5,EPS/2,fr,dep-1);
}
double operator()(double l,double r,double EPS=1e-7/*精度*/,int dep=16/*最小遞迴深度*/)const{
pdd p1(l,f(l)),p3(r,f(r)),p2=add(p1,p3);
return find(p1,p2,p3,EPS,simpson(p1,p2,p3),dep);
}
#undef simpson
};
double d, l, r, t, h;
//(1) 弓形面積
double F1(double r, double d) {
auto F0 = [](double r, double d) -> double { // D<=R 的情形
return r * r * acos(1 - d / r) - sqrt(d * (2 * r - d)) * (r - d);
};
if (d == r)
return PI * r * r / 2;
else if (d < r)
return F0(r, d);
else
return PI * r * r - F0(r, 2 * r - d);
}
//(2) 柱體截面積
double h1, cs, b, sc;
double F2(double x) {
if (t < 1e-6) return 2 * sqrt(x * (d - x)) * l; // 特判 直接返回一個長方形
double res = sc;
if (x < h1) res -= F1(b, (h1 - x) / cs);
if (x > t) res -= F1(b, (x - t) / cs);
return res / t * l;
}
//(3) 的左球缺
double k, y3, c;
double F3(double x) {
double ym = sqrt(x * (2 * r - x));
if (t + 1e-6 > l) return PI * ym * ym; // 特判 直接返回一個圓
return F1(ym, ym + max(min(ym, y3 - k * x), -ym));
}
//(3) 的右球缺
double y4;
double F4(double x) {
double ym = sqrt(x * (2 * r - x));
return F1(ym, ym + max(min(ym, y4 + k * x), -ym));
}
signed main() {
freopen("damage.in", "r", stdin);freopen("damage.out", "w", stdout); // 牛客上提交似乎要去掉這句？
scanf("%lf%lf%lf%lf%lf", &d, &l, &r, &t, &h);
d /= 100; l /= 100; r /= 100; t /= 100; h /= 100;// 所有輸入數除以100以和輸出的單位同步
// 求解全部輔助引數
cs = sqrt(1 - (t / l) * (t / l)); // cos w
h1 = d * cs; // 圓柱側面長方形的左上端點的高度
b = d / 2; // 圓柱側面半徑
sc = b * b * PI; // 圓柱側面(直徑為d的圓)積
c = r - sqrt(r * r - b * b); // 球缺的厚度
k = t / sqrt(l * l - t * t);
y3 = b - (h1 - h) / cs + k * c;
y4 = (h - t) / cs - b - k * c;
// 根據幾何情況統計答案
double ans = (Integration(F2))(0, min(h, t + h1)) + (Integration(F3))(0, c);
if (t + 1e-6 > l) {
if (h > t) ans += (Integration(F3))(0, h - t);
} else {
ans += (Integration(F4))(0, c);
}
printf("%.9lf\n", ans);
return 0;
}