白书上也有讲,就是一个在精度要求不高的时候计算积分的算法。
推荐博客:https://blog.csdn.net/VictoryCzt/article/details/80660113
我感觉其实功利一点的话,不用完全理解它的原理,大概明白然后会套板子做题即可(蒟蒻口胡)。
洛谷P4525
模板题,直接上模板。
#include<bits/stdc++.h> using namespace std; const double eps=1e-7; double a,b,c,d,L,R; double f(double x) { //具体计算公式 return (c*x+d)/(a*x+b); } double simpson(double l,double r) { //辛普森积分公式 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6; } double asr(double l,double r,double exps,double val) { //自适应辛普森积分 double mid=(l+r)/2; double lval=simpson(l,mid),rval=simpson(mid,r); if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15; else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval); } int main() { scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R); printf("%.6lf ",asr(L,R,eps,simpson(L,R))); return 0; }
洛谷P4526
另一道模板题,不过辛普森积分好像不适用于发散函数和反常积分(积分区间无穷大)。所以我们要想办法处理这两种情况。①可以证明a<0函数发散,a>0函数收敛。②x到底大于某个数之后变得非常非常小,于是之后的就可以不积分了。于是这题就可以做了。
#include<bits/stdc++.h> using namespace std; const double eps=1e-7; double a; double f(double x) { //具体计算公式 return pow(x,a/x-x); } double simpson(double l,double r) { //辛普森积分公式 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6; } double asr(double l,double r,double exps,double val) { //自适应辛普森积分 double mid=(l+r)/2; double lval=simpson(l,mid),rval=simpson(mid,r); if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15; else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval); } int main() { scanf("%lf",&a); if (a<0) puts("orz"); else printf("%.5lf ",asr(1e-10,20,eps,simpson(1e-10,20))); return 0; }
HDU-1724
求椭圆面积一部分。精度要求很低,用辛普森积分。
#include<bits/stdc++.h> using namespace std; const double eps=1e-7; double a,b,l,r; double f(double x) { //具体计算公式 return 2*b*sqrt(1-x*x/a/a); } double simpson(double l,double r) { //辛普森积分公式 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6; } double asr(double l,double r,double exps,double val) { //自适应辛普森积分 double mid=(l+r)/2; double lval=simpson(l,mid),rval=simpson(mid,r); if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15; else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval); } int main() { int T; cin>>T; while (T--) { scanf("%lf%lf%lf%lf",&a,&b,&l,&r); printf("%.3lf ",asr(l,r,eps,simpson(l,r))); } return 0; }
BZOJ-2178
求圆的面积并。直接在x轴上积分,然后f函数就是直线经过的线段长度。其他还是套辛普森积分即可。
但是这题有非常多细节:第一,要先把被包含的圆去掉。第二,可以看出每次计算f函数都要耗费大量时间,所以要想办法通过尽量复用f函数值来减少计算次数。第三,这也是蒟蒻没想懂的一点,好像用我上面板子的控制精度方法是会超时的(即每次二分将精度也除以2的方法),这题可以用一个确切的精度1e-13然后一直用这个精度控制即可。
#include<bits/stdc++.h> using namespace std; const double eps=1e-13; const int N=1e3+10; int n; double x[N],y[N],r[N]; bool ban[N]; struct dat{ double x,y; bool operator < (const dat &rhs) const { return x<rhs.x; } }A[N<<1]; double f(double P) { //计算直线P经过线段长度 int num=0; for (int i=1;i<=n;i++) { if (fabs(x[i]-P)>=r[i]) continue; double d=fabs(x[i]-P); double l=sqrt(r[i]*r[i]-d*d); A[++num].x=y[i]-l; A[num].y=y[i]+l; } sort(A+1,A+num+1); double ret=0,lst=-0x3f3f3f3f; for (int i=1;i<=num;i++) if (A[i].x>lst) ret+=A[i].y-A[i].x,lst=A[i].y; else if (A[i].y>lst) ret+=A[i].y-lst,lst=A[i].y; return ret; } double simpson(double l,double r) { //辛普森积分公式 return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6; } double cal(double l,double r,double fl,double fr,double fm){ return (fl+fr+4*fm)*(r-l)/6; } double asr(double l,double r,double val,double L,double R,double M) { //自适应辛普森积分 double mid=(l+r)/2,LM=f((l+mid)/2),RM=f((mid+r)/2); double lval=cal(l,mid,L,M,LM),rval=cal(mid,r,M,R,RM); if (fabs(lval+rval-val)<=eps) return lval+rval+(lval+rval-val)/15; return asr(l,mid,lval,L,M,LM)+asr(mid,r,rval,M,R,RM); } int main() { cin>>n; double L=0x3f3f3f3f,R=-0x3f3f3f3f; for (int i=1;i<=n;i++) { scanf("%lf%lf%lf",&x[i],&y[i],&r[i]); L=min(L,x[i]-r[i]); R=max(R,x[i]+r[i]); } for (int i=1;i<=n;i++) for (int j=i+1;j<=n;j++) if (!ban[i] && !ban[j]) { double dist=sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j])); if (dist>fabs(r[i]-r[j])) continue; if (r[j]<=r[i]) ban[j]=1; else ban[i]=1; } int t=0; for (int i=1;i<=n;i++) if (!ban[i]) { t++; x[t]=x[i]; y[t]=y[i]; r[t]=r[i]; } n=t; double M=(L+R)/2; printf("%.3lf ",asr(L,R,simpson(L,R),f(L),f(R),f(M))); return 0; }
洛谷P4207/BZOJ-1502