自適應Simpson積分

Simpson公式的話是一個對於三次及如下函數成立的連續區間求定積分的公式.(好像還有許多其餘的公式……)
雖然他是一個普通的定積分公式,可是他不只能夠用來求不規則函數的定積分還能夠用來瞎搞求面積,你能夠把用這個式子求面積當作是對於定積分本質的運用,也能夠直接把圖像當作函數,我習慣於按照後者來理解.
http://blog.csdn.net/greatwall1995/article/details/8639135
上面那個博客已經把Simpson積分說的已經比較明白了.
http://blog.csdn.net/xl2015190026/article/details/53518077
上面這個博客裏提到了對於Simpson積分的修正——I.乘15 II.把精度劈半
可是據個人經驗,不修正才穩啊.
對於自適應Simpson積分的代碼實現,在遞歸的過程當中,除了左右端點之外,其餘的參量能夠根據我的喜愛與具體狀況選擇性地傳入.
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 void read(){scanf("%d%d%d",&x,&y,&r);}
  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){
    cir[i].read();
    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;
}

個人代碼

這我的比我打得要漂亮得多http://blog.csdn.net/liutian429073576/article/details/51154077