题解：P10624 [ICPC 2013 WF] Pollution Solution

ChrisWangZi · 2026-01-26 23:13:18 · 题解

P10624 [ICPC 2013 WF] Pollution Solution 题解

题目传送门

前言：本文将尽量从头梳理本题。有一定数学基础的读者可直接跳至Part 1。

Part 0 前置知识：积分入门及自适应辛普森法

导数和积分

对于任意连续且可导函数 f(x)，记其导函数为 f'(x)。

我们用 \mathrm{d} 来表示其后面变量的“微小变化量”（微元），因此 \mathrm{d}x 表示 x 的微小变化量。

（原函数与导函数的关系常写作 \displaystyle\frac{\mathrm{d}}{\mathrm{d}x}f(x)=f'(x)，但本文将尽量不使用该符号。）

则有 \displaystyle\int_l^r f'(x)\mathrm{d}x=f(r)-f(l)（牛顿-莱布尼茨公式）。

我们将上述在 \displaystyle\int 上下有数（上式 l 与 r）的积分称为定积分，将没有上下界的积分称为不定积分。不定积分的求解过程可以理解为求导的逆运算（再此不过多讨论）。

定积分的几何意义是该函数在 [l,r] 上与 x 轴围成的面积（可依此理解牛顿-莱布尼茨公式）。

基本求导与不定积分运算（本题解涉及到的）

任意函数 f(x)、g(x) 与 h(x)，满足 h(x)=f(x)+g(x)，则有 h'(x)=f'(x)+g'(x)。

任意函数 f(x)、g(x) 与 h(x)，满足 h(x)=f(x)+g(x)，则有 \displaystyle\int h(x)\mathrm{d}x=\int f(x)\mathrm{d}x+\int g(x)\mathrm{d}x+C。

（由于常函数 f(x)=c 的导函数 f'(x)=0，因此求不定积分时都要加上任意常数 C。）

任意函数 f(x) 与 g(x)，满足 g(x)=cf(x)（ c 为常数），则有 g'(x)=cf'(x)。

任意函数 f(x) 与 g(x)，满足 g(x)=cf(x)（ c 为常数），则有 \displaystyle\int g(x)\mathrm{d}x=c\int f(x)\mathrm{d}x+C。

任意幂函数 f(x)=x^c（ c 为常数），有 f'(x)=cx^{c-1}。

任意幂函数 f(x)=x^c（ c 为常数），有 \displaystyle\int f(x)\mathrm{d}x=\frac{1}{c+1}x^{c+1}+C。

对于任意二次函数 f(x)=ax^2+bx+c，由幂函数求导公式，有 \displaystyle\int f(x)\mathrm{d}x=\frac{a}{3}x^3+\frac{b}{2}x^2+cx+C。

辛普森公式

注：该推导强烈建议初学者自己动笔推一遍。

令 f(x)=ax^2+bx+c。

由二次函数求导公式，有

:::align{center}

:::

移项，将系数一样的进行合并，得

:::align{center}

:::

因式分解，得

:::align{center}

:::

提取公因式 \displaystyle\frac{r-l}{6}，得

:::align{center}

:::

拆项，移项，得

:::align{center}

:::

由 f(x)=ax^2+bx+c，得

:::align{center}

:::

这便是辛普森公式——对于二次函数的定积分求法。

自适应辛普森法

由于上述推导是基于 f(x) 为二次函数（抛物线）而进行的，所以对于大多数情况（非二次函数），该公式不成立。

接下来就轮到自适应辛普森法登场了——一种保证精度又有效控制时间复杂度的方法。

考虑任意函数 f(x)上，若 f(x) 在 [l,r] 上“长得很像”抛物线，那么我们就用辛普森公式进行近似求值。

具体地，采用分治思想。我们要求 \displaystyle\int_l^r f(x)\mathrm{d}x。先求出 [l,r] 段的值 s，再分别求出 [l,mid) 与 [mid,r] 的值 s1 与 s2，若 |s-(s1+s2)|<\epsilon（\epsilon 根据题目精度要求确定），那么我们认为 f(x) 在 [l,r] 上近似为抛物线了，否则继续对 [l,mid) 和 [mid,r] 段进行分治处理。

习题：

P4525 【模板】自适应辛普森法 1

P4526 【模板】自适应辛普森法 2

参考代码：

double fr(double x);
double sim(double l,double r){
    return (r-l)*(fr(l)+fr(r)+4*fr((l+r)/2))/6;
}
const double eps;
double sol(double l,double r){
    double mid=(l+r)/2;
    double s1=sim(l,r);
    double s2=sim(l,mid)+sim(mid,r);
    if(abs(s1-s2)<=eps){
        return s2;
    }
    else{
        return sol(l,mid)+sol(mid,r);
    }
}

Part 1 关于本题

要求一个多边形与半圆相交的面积。

对于所有非竖直（与 x 轴平行）的边，其比在多边形所围成面积的下方或上方。通过 \Omicron(n^2) 的时间暴力判断（判断两线段位置关系，统计位于其上或下方边的数量）。

对于所有边，求其上方与半圆围成的面积（将两个定积分相减），若不存在或积出来负值（如圆弧完全位于线段下方）则跳过这条边。对于处在下方的边，将该面积与答案相加，因为其上方必然存在相交面积与一个上边（有可能是半圆的上圆弧）；对于处在上方的边，将答案减去该面积，因为其上方已经没有相交面积了。

（特别地，计算面积须计算线段（所在直线）与半圆的交点，不能直接使用线段端点积分。）

最终复杂度 \Omicron(n^2)。

附上代码：

#include<iostream>
#include<cmath>
using namespace std;
double k,b;
int n;
double r; 
double fr(double x){ 
    return sqrt(r*r-x*x)-(k*x+b);
}
double sim(double l,double r){
    return (r-l)*(fr(l)+fr(r)+4*fr((l+r)/2))/6;
}
const double eps=1e-4;
double sol(double l,double r){
    double mid=(l+r)/2;
    double s1=sim(l,r);
    double s2=sim(l,mid)+sim(mid,r);
    if(abs(s1-s2)<=eps){
        return s2;
    }
    else{
        return sol(l,mid)+sol(mid,r);
    }
}
struct pt{
    double x,y;
    double mid;
    int t1,t2;
    int dr;
}a[105];
double py(double X,int i,int j){//判断i,j点所连直线在x=X时y的值 
    return a[i].y+(a[i].y-a[i+1].y)/(a[i].x-a[i+1].x)*(X-a[i].x);
}
bool down(int i,int j,int k,int l){//判断两直线位置关系 
    return py(a[k].mid,i,j)<py(a[k].mid,k,l);
}
int main(){
    ios::sync_with_stdio(0);
    cin>>n>>r;
    for(int i=1;i<=n;i++){
        cin>>a[i].t1>>a[i].t2; 
        a[i].x=a[i].t1,a[i].y=a[i].t2;
    }
    a[n+1]=a[1];
    for(int i=1;i<=n;i++){
        //防止竖直线段造成干扰，故选择一非整数横坐标 
        if((a[i].t1+a[i+1].t1)%2==0){
            a[i].mid=(a[i].x+a[i+1].x)/2+0.5; 
        }
        else{
            a[i].mid=(a[i].x+a[i+1].x)/2;
        }
    }
    for(int i=1;i<=n;i++){
        if(a[i].t1==a[i+1].t1){
            continue;
        }
        int cnt=0;
        for(int j=1;j<=n;j++){
            if(i==j){
                continue;
            }
            if(a[j].t1==a[j+1].t1){
                continue;
            }
            if((a[j].x<=a[i].mid&&a[i].mid<=a[j+1].x)||(a[j+1].x<=a[i].mid&&a[i].mid<=a[j].x)){
                if(down(j,j+1,i,i+1)){
                    cnt++;
                }
            }
        }
        //统计其下方线段数量 
        if(cnt%2==0){
            a[i].dr=1;//下边 
        }
        else{
            a[i].dr=-1;//上边 
        }
    }
    double ans=0;
    for(int i=1;i<=n;i++){
        if(a[i].x==a[i+1].x){
            continue;
        }
        double mL=min(a[i].x,a[i+1].x);
        double mR=max(a[i].x,a[i+1].x);
        double L,R;
        k=(a[i].y-a[i+1].y)/(a[i].x-a[i+1].x);
        b=a[i].y-a[i].x*k;
        //计算判别式 
        double del=pow(2*b*k,2)-4*(k*k+1)*(b*b-r*r);
        if(del<0){
            //无交点，跳过 
            L=1,R=0;
        }
        else{
            //计算交点，与线段端点比较 
            double x1=(-2*b*k-sqrt(del))/(2*(k*k+1));
            double x2=(-2*b*k+sqrt(del))/(2*(k*k+1));
            L=max(x1,mL);
            R=min(x2,mR);
        }
        if(L>R){
            continue;
        }
        if(a[i].dr==1){
            ans+=sol(L,R);
        }
        else{
            ans-=sol(L,R);
        }
    }
    printf("%.3lf",ans);
    return 0;
}