首页 > 其他分享> > 2014-2015 ACM-ICPC Northeastern European Regional Contest (NEERC 14) 题解

2014-2015 ACM-ICPC Northeastern European Regional Contest (NEERC 14) 题解

2022-05-16 17:33:42 作者：互联网

gym100553

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)\)，\(R<D\leq 2R\) 的情况同理。

（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\)），对于特定的 \(x\) 计算出损失的部分的“长度”（例如图中红线），即根据（1）算出左右两边损失了弓形后剩下的面积，最后再除以 \(\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;
}