並不會計算幾何的並不對勁的人想必是非常之菜的。

很對勁的太刀流在這裏->

辛普森積分本身是把一段函數f(x)當成二次函數算，用(f(l)+4*f((l+r)/2)+f(r))*(r-l)/6算出這段曲線在l到r之間與x軸所夾面積。

但是這聽上去非常不精確，所以有人發明了自適應的辛普森積分。

每一次分別算出這段左一半、右一半和整個的辛普森積分，如果左+右與整個之差大於eps就繼續遞歸處理左、右半部分。

至於為什麽不把一段函數當成一次函數算，是因為會被f(m)恰好在直線上這種數據卡：

而二次函數只能用這種手形數據卡：

那為什麽不用三次函數呢？想必是因為三次函數復雜的讓人手不健康而且手形數據不是那麽好構造吧。

————————例題分界線————————
bzoj2178: 圓的面積並
Time Limit: 20 Sec Memory Limit: 259 MB
Submit: 2415 Solved: 617
[Submit][Status][Discuss]
Description
給出N個圓，求其面積並
Input
先給一個數字N ,N< = 1000 接下來是N行是圓的圓心，半徑，其絕對值均為小於1000的整數
Output
面積並，保留三位小數
————————題解分界線—————————
可以把直線x=x₁ 在圓中的部分的長度看成是關於x₁的函數。對這個函數求自適應的辛普森積分，但是復雜度非常之玄學。

經某些人的實驗，發現eps設為1e-14會TLE，設為1e-12會WA，所以好像沒什麽可選的…

註意要刪去所有被包含的圓、加一堆常數優化才能過。

#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#define fir first
#define sec second
#define eps (1e-13)
#define maxn 1010
#define four 4.0
#define six 6.0
#define
 
 two 2.0
#define re register 
using namespace std;
typedef pair<double,double>Pdd; 
typedef struct point{double x,y;}P;
typedef struct circle
{
    P o;
    double r;
    bool operator <(const circle &z)const{ return o.x<z.o.x;}
    Pdd getd(double x)  
    {  
        if(r<=fabs(o.x-x))   return Pdd(0,0);//相離 
        double dis=sqrt(r*r-(o.x-x)*(o.x-x));  
        return Pdd(o.y-dis,o.y+dis); //下、上交點 
    }  
}C;
int n,no[maxn],cntt;
double minl=10101010.0,maxr=-10101010.0;
C cir[maxn],circ[maxn];
Pdd pd[maxn];
double D(P _1,P _2){return sqrt((_1.x-_2.x)*(_1.x-_2.x)+(_1.y-_2.y)*(_1.y-_2.y));}

double getf(double x)
{//此線在多個圓中的部分的並 
    double ans=0,last=-10101010;  
    int cnt=0;  
    for(re int i=1;i<=cntt;i++)  
    {  
        pd[++cnt]=circ[i].getd(x);  
        cnt=(pd[cnt]==Pdd(0,0))?cnt-1:cnt;  
    }  
    sort(pd+1,pd+cnt+1);  
    for(re int i=1;i<=cnt;i++)  
    {  
        if(pd[i].fir>last)  
            ans+=pd[i].sec-pd[i].fir,last=pd[i].sec;   
        else if(pd[i].sec>last)  
            ans+=pd[i].sec-last,last=pd[i].sec;  
    }  
    //printf("***%.8lf %.8lf***",x,ans);
    return ans;  
} 
inline double getsim(double l,double r,double fl,double fm,double fr){return ((r-l)/six)*(fl+four*fm+fr);}
double simpson(double l,double m,double r,double fl,double fm,double fr)
{//是對f(x)表示直線x在圓內的部分的總長度這個函數求積分 
    double lm=(l+m)/two,rm=(m+r)/two,flm=getf(lm),frm=getf(rm);
    double lans=getsim(l,m,fl,flm,fm),rans=getsim(m,r,fm,frm,fr),totans=getsim(l,r,fl,fm,fr);
    //printf("%.8lf %.8lf %.8lf %.8lf\n",l,r,lans+rans,totans);system("pause");
    if(fabs(lans+rans-totans)<=eps)return totans;
    else return simpson(l,lm,m,fl,flm,fm)+simpson(m,rm,r,fm,frm,fr);
}
int main()
{
    memset(no,0,sizeof(no));
    scanf("%d",&n);
    for(re int i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&cir[i].o.x,&cir[i].o.y,&cir[i].r);
        if(minl>cir[i].o.x-cir[i].r)minl=cir[i].o.x-cir[i].r;
        if(maxr<cir[i].o.x+cir[i].r)maxr=cir[i].o.x+cir[i].r;
    }
    sort(cir+1,cir+n+1);
    for(re int i=1;i<=n;i++)
    {
        if(no[i])continue;    
        for(int j=i+1;j<=n;j++)
        {
            if(no[j])continue;
            if(D(cir[i].o,cir[j].o)+cir[j].r<=cir[i].r)  
                no[j]=1;  
        }
    }//去除被包含的圓 
    for(re int i=1;i<=n;i++)
        if(!no[i])circ[++cntt]=cir[i];
    double ans=simpson(minl,(minl+maxr)/two,maxr,0,getf((minl+maxr)/two),0);
    printf("%.3lf",ans);
    return 0;
}

宣傳一波廣義嘔吐之光，歡迎加入。

並不對勁的辛普森積分