# 自适应Simpson积分

Simpson公式的话是一个对于三次及以下函数成立的连续区间求定积分的公式.(好像还有许多其他的公式……)

http://blog.csdn.net/greatwall1995/article/details/8639135

http://blog.csdn.net/xl2015190026/article/details/53518077

hdu1724:Ellipse *真尼玛裸

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
typedef double db;
const db eps=1e-4;
db K,B;
inline db f(db x){
return std::sqrt(K*x*x+B);
}
inline db calc(db l,db r){
return (f(l)+4*f((l+r)*0.5)+f(r))*(r-l)/6;
}
inline db simpson(db l,db r){
db mid=(l+r)*0.5;
if(std::fabs(calc(l,mid)+calc(mid,r)-calc(l,r))<eps)return calc(l,r);
return simpson(l,mid)+simpson(mid,r);
}
int main(){
int T,a,b,l,r;scanf("%d",&T);
while(T--){
scanf("%d%d%d%d",&a,&b,&l,&r);
K=-(db)(b*b)/(a*a),B=b*b;
printf("%.3f\n",simpson(l,r)*2);
}
}

bzoj2178:圆的面积并 *也真尼玛裸(如果你不怕被hack的话)

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ft first
#define sd second
#define mmp(a,b) (std::make_pair(a,b))
typedef double db;
typedef std::pair<db,db> pdd;
const int N=1010;
const db eps=1e-13;
pdd seg[N];
struct Cir{
int x,y,r;
inline db dis(Cir a){
return std::sqrt((a.x-x)*(a.x-x)+(a.y-y)*(a.y-y));
}
inline bool in(db pos){return x-r<pos&&x+r>pos;}
inline bool in(Cir a){
return r>=dis(a)+a.r;
}
inline pdd get(db pos){
db len=std::sqrt(r*r-(pos-x)*(pos-x));
return mmp(y-len,y+len);
}
}cir[N];
bool die[N];
int n;
inline db f(db pos){
int i,len=0;db ret=0,l,r;
for(i=1;i<=n;++i)
if(cir[i].in(pos))
seg[++len]=cir[i].get(pos);
std::sort(seg+1,seg+(len+1));
for(i=1;i<=len;++i){
l=seg[i].ft,r=seg[i].sd;
if(i!=1&&seg[i-1].sd>l)l=seg[i-1].sd;
if(l<r)ret+=r-l;
else seg[i].sd=l;
}
return ret;
}
inline db calc(db l,db r,db fl,db fr,db fm){
return (r-l)/6*(fl+4*fm+fr);
}
inline db Simpson(db l,db r,db mid,db fl,db fr,db fm,db A){
db m1=(l+mid)/2,m2=(mid+r)/2;
db fm1=f(m1),fm2=f(m2);
db L=calc(l,mid,fl,fm,fm1),R=calc(mid,r,fm,fr,fm2);
if(std::fabs(L+R-A)<=eps)return A;
return Simpson(l,mid,m1,fl,fm,fm1,L)+Simpson(mid,r,m2,fm,fr,fm2,R);
}
int main(){
scanf("%d",&n);
int i,j,l=2000,r=-2000;
for(i=1;i<=n;++i){
l=std::min(l,cir[i].x-cir[i].r);
r=std::max(r,cir[i].x+cir[i].r);
}
for(i=1;i<=n;++i)
if(!die[i])
for(j=1;j<=n;++j)
if(j!=i&&!die[j])
if(cir[i].in(cir[j]))
die[j]=true;
j=0;
for(i=1;i<=n;++i)
if(!die[i])
cir[++j]=cir[i];
n=j;
db fl=f(l),fr=f(r),fm=f((l+r)/2.);
printf("%.3f\n",Simpson(l,r,(l+r)/2.,fl,fr,fm,calc(l,r,fl,fr,fm)));
return 0;
}

bzoj1502:[NOI2005]月下柠檬树 ***

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
typedef double db;
const db eps=1e-5;
const db Pi=std::acos(-1.);
const int N=510;
int n,m;
db alpha;
struct Cir{
db x,r,l;
bool in(Cir a){return r>a.r+std::fabs(a.x-x);}
bool in(db pos){return x-r<pos&&pos<x+r;}
db h(db pos){return std::sqrt(r*r-(x-pos)*(x-pos));}
}cir[N];
inline bool comp1(Cir a,Cir b){return a.l<b.l;}
struct Line{
db l,r,k,b;
bool in(db pos){return l<pos&&pos<r;}
db h(db pos){return k*pos+b;}
}line[N];
inline bool comp2(Line a,Line b){return a.l<b.l;}
inline Line get(db aks1,db w1,db aks2,db w2){
Line ret;
ret.l=aks1,ret.r=aks2;
ret.k=(w1-w2)/(aks1-aks2),ret.b=w1-aks1*ret.k;
return ret;
}
inline db f(db pos){
db ret=0.;int i;
for(i=1;i<=n;++i){
if(cir[i].in(pos))
ret=std::max(ret,cir[i].h(pos));
if(pos<cir[i].l)break;;
}
for(i=1;i<=m;++i){
if(line[i].in(pos))
ret=std::max(ret,line[i].h(pos));
if(pos<line[i].l)break;
}
return ret*2.;
}
inline db calc(db l,db r,db fl,db fr,db fm){
return (r-l)*(fl+4*fm+fr)/6;
}
inline db Simpson(db l,db r,db mid,db fl,db fr,db fm,db A){
db m1=(l+mid)*0.5,m2=(mid+r)*0.5;
db fm1=f(m1),fm2=f(m2);
db L=calc(l,mid,fl,fm,fm1),R=calc(mid,r,fm,fr,fm2);
if(std::fabs(L+R-A)<=eps)return L+R;
return Simpson(l,mid,m1,fl,fm,fm1,L)+Simpson(mid,r,m2,fm,fr,fm2,R);
}
int main(){
scanf("%d%lf",&n,&alpha);
++n,alpha=std::tan(alpha);
int i;db h,s=0.,l=1e18,r=-1e18;
for(i=1;i<=n;++i)
scanf("%lf",&h),cir[i].x=(s+=h)/alpha;
for(i=1;i<=n;++i){
if(i==n)h=0.;
else scanf("%lf",&h);
l=std::min(l,cir[i].x-h);
r=std::max(r,cir[i].x+h);
cir[i].r=h;
cir[i].l=(cir[i].x-h);
}
for(i=1;i<n;++i){
if(cir[i].in(cir[i+1])||cir[i+1].in(cir[i]))continue;
if(cir[i].r==cir[i+1].r)line[++m]=get(cir[i].x,cir[i+1].x,cir[i].r,cir[i].r);
else{
h=std::fabs(cir[i].x-cir[i+1].x);
h=h+h*(std::min(cir[i].r,cir[i+1].r))/(std::fabs(cir[i].r-cir[i+1].r));
if(cir[i].r>cir[i+1].r)s=std::acos(cir[i].r/h);
else s=std::acos(-cir[i+1].r/h);
line[++m]=get(cir[i].x+cir[i].r*std::cos(s),cir[i].r*std::sin(s),cir[i+1].x+cir[i+1].r*std::cos(s),cir[i+1].r*std::sin(s));
}
}
std::sort(cir+1,cir+(n+1),comp1);
std::sort(line+1,line+(m+1),comp2);
db mid=(l+r)*0.5,fl=f(l),fr=f(r),fm=f(m);
printf("%.2f\n",Simpson(l,r,mid,fl,fr,fm,calc(l,r,fl,fr,fm)));
return 0;
}

posted @ 2018-03-11 11:53  TS_Hugh  阅读(369)  评论(1编辑  收藏