自适应Simpson积分
首先介绍一下Simpson积分的公式$$\int_{a}^{b}f(x)dx \approx \frac {b-a} {6} (f(a)+4f(\frac{a+b} {2})+f(b))$$
可以用这个求比较平滑的曲线的积分的近似...
但是直接用这个公式来做会误差很大
具体用法是：
用这个方法积$[l,mid]$，积$[mid,r]$，再积$[l,r]$判断是不是在误差范围内
如果不在就递归地细化$[l,mid]$和$[mid,r]$
比如说SDOI2017龙与地下城这道题
就用这种方法积了一个正态分布曲线
#include<bits/stdc++.h> using namespace std; inline int read() { int x=0;char ch=getchar(); while(ch<'0' || '9'<ch)ch=getchar(); while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar(); return x; } typedef double db; db x,y; namespace Fast_Fouriet_Transform { typedef complex<db> cp; const int M=800009; const db pi=acos(-1); cp a[M],b[M]; db per,ans[M]; int m,l,rev[M]; inline void FFT(cp *a,int n,int f) { for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int h=2;h<=n;h<<=1) { cp w(cos(2*pi/h),f*sin(2*pi/h)); for(int j=0;j<n;j+=h) { cp wn(1,0),x,y; for(int k=j;k<j+(h>>1);k++) { x=a[k],y=wn*a[k+(h>>1)]; a[k]=x+y; a[k+(h>>1)]=x-y; wn*=w; } } } if(f==-1) for(int i=0;i<n;i++) a[i].real()/=(db)n; } int mian() { per=1.0/(db)x; for(m=1,l=0;m<(y*x);m<<=1)l++; for(int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1)); memset(a,0,sizeof(a)); memset(b,0,sizeof(b)); a[0].real()=1.0; for(int i=0;i<x;i++) b[i].real()=per; FFT(a,m,1);FFT(b,m,1); for(int i=0;i<m;i++) { int cnt=y; while(cnt) { if(cnt&1) a[i]=a[i]*b[i]; b[i]=b[i]*b[i]; cnt>>=1; } } FFT(a,m,-1); ans[0]=a[0].real(); for(int i=1;i<m;i++) ans[i]=a[i].real()+ans[i-1]; for(int i=1;i<=10;i++) { int a=read(),b=read(); if(a==0) printf("%.8lf\n",ans[b]); else printf("%.8lf\n",ans[b]-ans[a-1]); } return 0; } } namespace simpsons { const db pi2=sqrt(acos(-1)*2.0); const db eps=1e-12; db sigma,mu; inline db sqr(db x){return x*x;} inline db f(db x){return exp(-sqr(x)/2)/pi2;} inline bool dcmp(db x){return fabs(x)<15*eps;} inline void init() { sigma=sqrt((sqr(x)-1.0)/12.0); mu=(x-1.0)/2.0; } inline db simpson(db a,db b,db fl,db fmid,db fr) { return (b-a)*(fl+4.0*fmid+fr)/6.0; } inline db solve(db l,db r) { db mid=(l+r)/2.0,lm=(l+mid)/2.0,mr=(mid+r)/2.0; db fl=f(l),fmid=f(mid),fr=f(r),flm=f(lm),fmr=f(mr); db sipl=simpson(l,mid,fl,flm,fmid); db sipr=simpson(mid,r,fmid,fmr,fr); db sipm=simpson(l,r,fl,fmid,fr); if(dcmp(sipl+sipr-sipm)) return sipl+sipr+(sipl+sipr-sipm)/15; else return solve(l,mid)+solve(mid,r); } int mian() { init(); db base=solve(0,(x-1)*y); for(int i=1;i<=10;i++) { db a=(((db)read()-0.5)/y-mu)*sqrt(y)/sigma; db b=(((db)read()+0.5)/y-mu)*sqrt(y)/sigma; printf("%.8lf\n",solve(0,b)-solve(0,a)); } } } int main() { int T=read(); while(T--) { scanf("%lf%lf",&x,&y); if(x*y<=400000) Fast_Fouriet_Transform::mian(); else simpsons::mian(); } return 0; }
还有经典的NOI2005月下柠檬树
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<cmath> using namespace std; #define MAXN 1010 double alpha; int N,num; #define INF 1e12 #define eps 1e-5 struct Point { double x,y; Point (double X=0,double Y=0) {x=X; y=Y;} }; struct Circle { double r; Point c; Circle(Point C=(Point){0,0},double R=0) {c=C; r=R;} }C[MAXN]; struct Line { Point s,t; double k,b; Line(Point S=(Point){0,0},Point T=(Point){0,0}) { s=S,t=T; if (s.x>t.x) swap(s,t); k=(s.y-t.y)/(s.x-t.x); b=s.y-k*s.x; } double f(double x) {return k*x+b;} }l[MAXN]; int dcmp(double x) {if (fabs(x)<eps) return 0; return x<0? -1:1;} double F(double x) { double re=0; for (int i=1; i<=N; i++) //枚举圆是否与扫描线有交 { double d=fabs(x-C[i].c.x); if (dcmp(d-C[i].r)>0) continue; double len=2*sqrt(C[i].r*C[i].r-d*d); re=max(re,len); } for (int i=1; i<=num; i++) //枚举公切线 if (x>=l[i].s.x && x<=l[i].t.x) re=max(re,2*l[i].f(x)); return re; } //利用扫描线去判断 double Calc(double l,double r) {double mid=(l+r)/2; return (F(l)+F(r)+F(mid)*4)*(r-l)/6;} double Simpson(double l,double r,double now) { double mid=(l+r)/2; double x=Calc(l,mid),y=Calc(mid,r); if (!dcmp(now-x-y)) return now; else return Simpson(l,mid,x)+Simpson(mid,r,y); } void Solve() { double L=INF,R=-INF; for (int i=1; i<=N+1; i++) L=min(L,C[i].c.x-C[i].r),R=max(R,C[i].c.x+C[i].r); // printf("%lf\n%lf\n",L,R); for (int i=1; i<=N; i++) { double d=C[i+1].c.x-C[i].c.x; if (dcmp(d-fabs(C[i].r-C[i+1].r))<0) continue; //特判小圆被大圆覆盖的情况 double sina=(C[i].r-C[i+1].r)/d,cosa=sqrt(1-sina*sina); l[++num]=(Line){(Point){C[i].c.x+C[i].r*sina,C[i].r*cosa},(Point){C[i+1].c.x+C[i+1].r*sina,C[i+1].r*cosa}}; } printf("%.2lf\n",Simpson(L,R,Calc(L,R))); } int main() { scanf("%d%lf",&N,&alpha); double h,r; for (int i=1; i<=N+1; i++) scanf("%lf",&h), C[i]=(Circle){((Point){(h/tan(alpha))+C[i-1].c.x,0}),0}; for (int i=1; i<=N; i++) scanf("%lf",&r),C[i].r=r; // for (int i=1; i<=N+1; i++) // printf("%d %.2lf %.2lf\n",i,C[i].c.x,C[i].r); Solve(); return 0; }
感觉...算个黑科技吧
攒着以后可能会有点小用
