一个完全不会计算几何的蒟蒻的自我拯救……
辛普森积分
有的时候会有一些毒瘤计算几何题，要求的图形面积边缘是一段函数，而这个函数解析式通常非常繁琐，没办法直接用公式积分，所以就需要用辛普森积分求近似值。
辛普森积分的用途就是在精度要求不高的时候（通常是求图形面积）求函数积分的近似值，大概步骤就是在积分区间$[a,b]$中不断二分，每次找$a,b,frac{a+b}{2}$三点用一条抛物线来拟合原函数，当误差小于eps的时候就返回答案。
辛普森积分公式
$$int_{a}^{b}f(x)dxapproxfrac{b-a}{6} imesleft[f(a)+f(b)+4fleft(frac{a+b}{2} ight) ight]$$
下面推导一波：
设$g(x)$是一个关于$x$的二次函数，且$g(x)=Ax^2+Bx+C$，对于定积分$int_{0}^{a}g(x)dx$求积得$frac{Ax^3}{3}+frac{Bx^2}{2}+Cx+D$，其中$D$为常数可以忽略；
令$H(x)=int_{0}^{a}g(x)dx$，则求其中的一段定积分就有$int_{a}^{b}g(x)dx=H(b)-H(a)$；
设$g(x)$是当前的拟合函数，那么有：
$$int_{a}^{b}f(x)dxapproxint_{a}^{b}g(x)dx$$
$$=H(b)-H(a)$$
$$=frac{A}{3}(b^3-a^3)+frac{B}{2}(b^2-a^2)+C(b-a)$$
$$=frac{b-a}{6} imesleft[2A(a^2+ab+b^2)+3B(a+b)+6C ight]$$
然后大力拆括号配方：
$$=frac{b-a}{6} imesleft[(Aa^2+Ba+C)+(Ab^2+Bb+C)+A(a^2+2ab+b^2)+2B(a+b)+4C ight]$$
$$=frac{b-a}{6} imesleft[g(a)+g(b)+4Aleft(frac{a+b}{2} ight)^2+4Bleft(frac{a+b}{2} ight)+4C ight]$$
$$=frac{b-a}{6} imesleft[g(a)+g(b)+4gleft(frac{a+b}{2} ight) ight]$$
在实际使用的时候，$g(x)$拟合的曲线可以用$f(x)$代替，因此：
$$int_{a}^{b}f(x)dxapproxfrac{b-a}{6} imesleft[f(a)+f(b)+4fleft(frac{a+b}{2} ight) ight]$$
实现：
具体实现就是不断二分区间，用辛普森积分公式求出近似值，当精度达到要求后就返回；
写起来很方便，核心代码就十几行左右……
注意用辛普森的时候要大力调eps，否则可能会在WA和TLE之间徘徊……
代码：
模板：洛咕P4525 自适应辛普森法1
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 using namespace std;
10 typedef long long ll;
11 double a,b,c,d,l,r;
12 double f(double x){
13 return (c*x+d)/(a*x+b);
14 }
15 double simpson(double l,double r){
16 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
17 }
18 double sint(double l,double r,double v){
19 double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
20 if(fabs(lv+rv-v)<=eps){
21 return v;
22 }
23 return sint(l,mid,lv)+sint(mid,r,rv);
24 }
25 double getint(double l,double r){
26 return sint(l,r,simpson(l,r));
27 }
28 int main(){
29 scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
30 printf("%.6lf",getint(l,r));
31 return 0;
32 }
一个精度更高的做法：
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 //#define eps 1e-7
9 using namespace std;
10 typedef long long ll;
11 const double eps=1e-7;
12 double a,b,c,d,l,r;
13 double f(double x){
14 return (c*x+d)/(a*x+b);
15 }
16 double simpson(double l,double r){
17 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
18 }
19 double sint(double l,double r,double ep,double v){
20 double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
21 if(fabs(lv+rv-v)<=ep*15){
22 return lv+rv+(lv+rv-v)/15;
23 }
24 return sint(l,mid,ep/2,lv)+sint(mid,r,ep/2,rv);
25 }
26 double getint(double l,double r,double ep){
27 return sint(l,r,ep,simpson(l,r));
28 }
29 int main(){
30 scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
31 printf("%.6lf",getint(l,r,eps));
32 return 0;
33 }
例题
【BZOJ1502】【NOI2005】月下柠檬树
若干个圆台投影之后形状大概是若干个圆加外公切线围成的图形，然后最顶上是个圆锥，整个图形是轴对称的，大力辛普森积分就好了；
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 using namespace std;
10 typedef long long ll;
11 typedef double db;
12 struct cir{
13 db x,y,xx,yy;
14 }a[1001];
15 int n,tot=0;
16 db al,mi,mx,t,h[1001],s[1001],r[1001];
17 db f(db x){
18 db ret=0;
19 for(int i=0;i<=n;i++){
20 if(fabs(s[i]-x)<r[i]){
21 ret=max(ret,sqrt(r[i]*r[i]-(s[i]-x)*(s[i]-x)));
22 }
23 }
24 for(int i=1;i<=tot;i++){
25 if(x>a[i].x&&x<a[i].xx){
26 ret=max(ret,a[i].y+(a[i].yy-a[i].y)*(x-a[i].x)/(a[i].xx-a[i].x));
27 }
28 }
29 return ret;
30 }
31 db sim(db l,db r,db fl,db fm,db fr){
32 return (fl+fr+4*fm)*(r-l)/6;
33 }
34 db sint(db l,db r,db fl,db fm,db fr){
35 db mid=(l+r)/2,L=f((l+mid)/2),R=f((mid+r)/2),s=sim(l,r,fl,fm,fr),sl=sim(l,mid,fl,L,fm),sr=sim(mid,r,fm,R,fr);
36 if(fabs(sl+sr-s)<eps)return sl+sr;
37 return sint(l,mid,fl,L,fm)+sint(mid,r,fm,R,fr);
38 }
39 int main(){
40 scanf("%d%lf",&n,&al);
41 for(int i=0;i<=n;i++){
42 scanf("%lf",&h[i]);
43 if(i)h[i]+=h[i-1];
44 s[i]=h[i]/tan(al);
45 }
46 for(int i=0;i<n;i++){
47 scanf("%lf",&r[i]);
48 mi=min(mi,s[i]-r[i]);
49 mx=max(mx,s[i]+r[i]);
50 }
51 mx=max(mx,s[n]);
52 for(int i=0;i<n;i++){
53 t=s[i+1]-s[i];
54 if(t>fabs(r[i+1]-r[i])){
55 a[++tot].x=s[i]-r[i]*(r[i+1]-r[i])/t;
56 a[tot].y=sqrt(r[i]*r[i]-(a[tot].x-s[i])*(a[tot].x-s[i]));
57 a[tot].xx=s[i+1]-r[i+1]*(r[i+1]-r[i])/t;
58 a[tot].yy=sqrt(r[i+1]*r[i+1]-(a[tot].xx-s[i+1])*(a[tot].xx-s[i+1]));
59 }
60 }
61 printf("%.2lf",sint(mi,mx,0,f((mi+mx)/2),0)*2);
62 return 0;
63 }
【BZOJ2178】圆的面积并
先把被包含的圆去掉，然后大力切圆就好了
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 using namespace std;
10 typedef long long ll;
11 typedef double db;
12 struct node{
13 db x,y;
14 node(){}
15 node(db _x,db _y){x=_x,y=_y;}
16 friend bool operator <(node a,node b){
17 return (fabs(a.x-b.x)<eps)?a.y<b.y:a.x<b.x;
18 }
19 }p[1001];
20 struct cir{
21 db x,y,r;
22 bool ch;
23 friend bool operator <(cir a,cir b){
24 return a.x<b.x;
25 }
26 node getf(db xx){
27 if(fabs(x-xx)>=r)return node(0,0);
28 return node(y-sqrt(r*r-(x-xx)*(x-xx)),y+sqrt(r*r-(x-xx)*(x-xx)));
29 }
30 }c[1001];
31 int n;
32 db mi=inf,mx=-inf;
33 db dis(cir a,cir b){
34 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
35 }
36 db f(db x){
37 db ret=0,ls=-inf;
38 int top=0;
39 for(int i=1;i<=n;i++){
40 p[++top]=c[i].getf(x);
41 if(p[top].x==0&&p[top].y==0)top--;
42 }
43 sort(p+1,p+top+1);
44 for(int i=1;i<=top;i++){
45 if(p[i].x>ls){
46 ret+=p[i].y-p[i].x;
47 ls=p[i].y;
48 }else if(p[i].y>ls){
49 ret+=p[i].y-ls;
50 ls=p[i].y;
51 }
52 }
53 return ret;
54 }
55 db sim(db l,db r,db fl,db fm,db fr){
56 return (fl+4*fm+fr)*(r-l)/6;
57 }
58 db sint(db l,db r,db fl,db fm,db fr){
59 db mid=(l+r)/2,L=f((l+mid)/2),R=f((mid+r)/2),s=sim(l,r,fl,fm,fr),sl=sim(l,mid,fl,L,fm),sr=sim(mid,r,fm,R,fr);
60 if(fabs(sl+sr-s)<eps)return sl+sr;
61 return sint(l,mid,fl,L,fm)+sint(mid,r,fm,R,fr);
62 }
63 int main(){
64 scanf("%d",&n);
65 for(int i=1;i<=n;i++){
66 scanf("%lf%lf%lf",&c[i].x,&c[i].y,&c[i].r);
67 mi=min(mi,c[i].x-c[i].r);
68 mx=max(mx,c[i].x+c[i].r);
69 }
70 sort(c+1,c+n+1);
71 for(int i=1;i<=n;i++){
72 if(!c[i].ch){
73 for(int j=i+1;j<=n;j++){
74 if(!c[j].ch&&dis(c[i],c[j])<=c[i].r-c[j].r){
75 c[j].ch=true;
76 }
77 }
78 }
79 }
80 for(int i=1;i<=n;i++){
81 if(c[i].ch)swap(c[i--],c[n--]);
82 }
83 printf("%.3lf",sint(mi,mx,0,f((mi+mx)/2),0));
84 return 0;
85 }