辛普森积分

这两天看了看辛普森积分，觉得这是很好的骗分工具，正好也比较简单，就随便写一写。这个东西对于三次以下的函数是正确的，但是对于三次以上的函数我们可以将其近似为低次函数套用Simpson公式，这个东西学名好像叫自适应Simpson积分。

先附上公式。

算法的大概流程就是不断的递归求对应区间的积分，当当前层和下一层的相对误差比较小时就退出，从而得到比较准确的值，但是这个其实是可以被卡掉的，比如一段波动的区间，我们取到的值有可能都是上顶点，就会直接退出，然而这并不是最后的答案。所以一般这个算法都是用来骗分的，2333。

附上代码

 1 double calc(double l,double r,double fl,double fr,double fm){
 2     return (r-l)/6.0*(fl+fr+4*fm);
 3 }
 4 double simpson(double l,double r,double mid,double fl,double fr,double fm,double s){
 5     double m1=(l+mid)/2.0,m2=(mid+r)/2.0;
 6     double f1=F(m1),f2=F(m2);
 7     double s1=calc(l,mid,fl,fm,f1),s2=calc(mid,r,fm,fr,f2);
 8     if(fabs(s-s1-s2)<eps)return s1+s2;
 9     return simpson(l,mid,m1,fl,fm,f1,s1)+simpson(mid,r,m2,fm,fr,f2,s2);
10 }

我一般在simpson函数里会将三个函数值一并传进去，否则在单次计算f的复杂度比较大时，会重复计算，浪费时间。

hdu1724，纯板子题。

这道题我写的还是不传函数值的版本，可以对照着看看

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <cmath>
 6 #define eps 1e-10
 7 using namespace std;
 8 double a,b;
 9 int T;
10 double f(double x){
11     return b*sqrt(1-(x*x)/(a*a));
12 }
13 double calc(double l,double r){
14     double mid=(l+r)/2.0;
15     return (r-l)/6*(f(l)+f(r)+4*f(mid));
16 }
17 double simpson(double l,double r,double now){
18     double mid=(l+r)/2.0;
19     double n1=calc(l,mid),n2=calc(mid,r);
20     if(fabs(now-n1-n2)<eps)return now;
21     return simpson(l,mid,n1)+simpson(mid,r,n2);
22 }
23 int main(){
24     scanf("%d",&T);
25     double l,r,ans;
26     while(T--){
27         scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
28         ans=simpson(l,r,calc(l,r))*2;
29         printf("%0.3f\n",ans);
30     }
31 }
View Code

bzoj2178，求圆的面积并

网上别人好像都是直接积，我是用了一个并查集，把连在一起的圆一起算，貌似这样被hack的概率比较小

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <cmath>
 6 #define eps 1e-13
 7 #define N 1005
 8 using namespace std;
 9 int v[N][N],cnt[N];
10 int n,fa[N],del[N],tot;
11 int le[N],ri[N];
12 int find(int x){return x==fa[x]?x:fa[x]=find(fa[x]);}
13 struct data{
14     double s,t;
15     data(){}
16     data(double a,double b){s=a,t=b;}
17 }L[N];
18 bool cmpdown(const data & a,const data & b){return a.s<b.s;}
19 struct cir{
20     double x,y,r;
21     void read(){scanf("%lf%lf%lf",&x,&y,&r);}
22 }p[N],d[N];
23 double dis(cir a,cir b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
24 bool cmpr(const cir & a,const cir & b){return a.r<b.r;}
25 bool cmpl(const int & a,const int & b){return p[a].x-p[a].r<p[b].x-p[b].r;}
26 bool check1(cir a,cir b){//包含
27     if(dis(a,b)+a.r<=b.r)return 1;
28     return 0;
29 }
30 bool check2(cir a,cir b){//相交
31     if(dis(a,b)<=a.r+b.r)return 1;
32     return 0;
33 }
34 double F(double x,int id){
35     tot=0;
36     for(int i=1;i<=cnt[id];i++){
37         int now=v[id][i];
38         if(x>p[now].x+p[now].r)continue;
39         if(x<p[now].x-p[now].r)break;
40         double l=sqrt(p[now].r*p[now].r-(p[now].x-x)*(p[now].x-x));
41         L[++tot]=data(p[now].y-l,p[now].y+l);
42     }
43     sort(L+1,L+tot+1,cmpdown);
44     double ans=0;
45     for(int i=1,j=1;i<=tot;i++){
46         j=i;double up=L[i].t;
47         while(j<tot&&L[j+1].s<=up){j++;up=max(up,L[j].t);}
48         ans+=up-L[i].s;i=j;
49     }
50     return ans;
51 }
52 double calc(double l,double r,double fl,double fr,double fm){
53     return (r-l)/6.0*(fl+fr+4*fm);
54 }
55 double simpson(double l,double r,double fl,double fr,double fm,double s,int id){
56     double mid=(l+r)/2.0,m1=(l+mid)/2.0,m2=(mid+r)/2.0;
57     double f1=F(m1,id),f2=F(m2,id),s1=calc(l,mid,fl,fm,f1),s2=calc(mid,r,fm,fr,f2);
58     if(fabs(s-s1-s2)<eps)return s1+s2;
59     return simpson(l,mid,fl,fm,f1,s1,id)+simpson(mid,r,fm,fr,f2,s2,id);
60 }
61 int main(){
62     memset(le,0x3f,sizeof le);
63     memset(ri,-0x3f,sizeof ri);
64     scanf("%d",&n);
65     for(int i=1;i<=n;i++)p[i].read(),fa[i]=i;
66     sort(p+1,p+n+1,cmpr);
67     for(int i=1;i<=n;i++)
68         for(int j=i+1;j<=n;j++)if(check1(p[i],p[j])){
69             del[i]=1;
70             break;
71         }
72     for(int i=1;i<=n;i++)if(!del[i])d[++tot]=p[i];
73     n=tot;
74     for(int i=1;i<=n;i++)p[i]=d[i];
75     for(int i=1;i<=n;i++)
76         for(int j=i+1;j<=n;j++)
77             if(check2(p[i],p[j]))
78                 fa[find(j)]=find(i);
79     for(int i=1;i<=n;i++){
80         int f=find(i);
81         v[f][++cnt[f]]=i;
82         le[f]=min(le[f],(int)(p[i].x-p[i].r));
83         ri[f]=max(ri[f],(int)(p[i].x+p[i].r));
84     }
85     double ans=0;
86     for(int i=1;i<=n;i++){
87         if(cnt[i]){
88             sort(v[i]+1,v[i]+cnt[i]+1,cmpl);
89             double fl=F(le[i],i),fr=F(ri[i],i),fm=F((le[i]+ri[i])/2.0,i);
90             ans+=simpson(le[i],ri[i],fl,fr,fm,calc(le[i],ri[i],fl,fr,fm),i);
91         }
92     }
93     printf("%0.3lf\n",ans);
94     return 0;
95 }
View Code

bzoj1502，转化完题面就是要求一些圆以及相邻圆公切线所包含的面积。

我被一个地方坑了一上午，最后发现是因为我先把包含的圆去掉再算的切线。身败名裂

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <iostream>
  4 #include <algorithm>
  5 #include <cmath>
  6 #define N 550
  7 #define eps 1e-5
  8 using namespace std;
  9 int tot,n,cnt;
 10 double alp;
 11 struct line{
 12     double x1,x2,y1,y2,K,B;
 13     bool operator < (const line & a)const{return x1<a.x1;}
 14 }l[N];
 15 struct cir{
 16     double x,r;
 17 }p[N],d[N];
 18 void solve(cir a,cir b){
 19     if(a.r==b.r){
 20         l[++cnt].x1=a.x;l[cnt].x2=b.x;
 21         l[cnt].y1=l[cnt].y2=a.r;
 22         l[cnt].K=0;
 23         l[cnt].B=a.r;
 24         return ;
 25     }
 26     if(a.r>b.r){
 27         double r1=a.r,r2=b.r,c=b.x-a.x;
 28         if(r2){
 29             double d=c*r2/(r1-r2);
 30             double e=sqrt(d*d-r2*r2);
 31             l[++cnt].x1=a.x+r1*r2/d;l[cnt].y1=r1*e/d;
 32             l[cnt].x2=b.x+r2*r2/d;l[cnt].y2=r2*e/d;
 33         }
 34         else{
 35             double d=sqrt(c*c-r1*r1);
 36             l[++cnt].x1=a.x+r1*r1/c;l[cnt].y1=r1*d/c;
 37             l[cnt].x2=b.x;l[cnt].y2=0;
 38         }
 39         l[cnt].K=(l[cnt].y2-l[cnt].y1)/(l[cnt].x2-l[cnt].x1);
 40         l[cnt].B=(l[cnt].y1)-l[cnt].K*l[cnt].x1;
 41         return ;
 42     }
 43     if(a.r<b.r){
 44         double r1=a.r,r2=b.r,c=b.x-a.x;
 45         double d=c*r1/(r2-r1);
 46         double e=sqrt(d*d-r1*r1);
 47         l[++cnt].x1=a.x-r1*r1/d;l[cnt].y1=r1*e/d;
 48         l[cnt].x2=b.x-r2*r1/d;l[cnt].y2=r2*e/d;
 49         l[cnt].K=(l[cnt].y2-l[cnt].y1)/(l[cnt].x2-l[cnt].x1);
 50         l[cnt].B=(l[cnt].y1)-l[cnt].K*l[cnt].x1;
 51         return ;
 52     }
 53 }
 54 double F(double x){
 55     double ans=0;
 56     for(int i=1;i<=cnt;i++){
 57         if(x>l[i].x2)continue;
 58         if(x<l[i].x1)break;
 59         ans=max(ans,l[i].K*x+l[i].B);
 60     }
 61     for(int i=1;i<=n;i++){
 62         if(x>p[i].x+p[i].r)continue;
 63         if(x<p[i].x-p[i].r)break;
 64         if(!p[i].r)continue;
 65         double now=(x-p[i].x);
 66         ans=max(ans,sqrt(p[i].r*p[i].r-now*now));
 67     }
 68     return ans*2.0;
 69 }
 70 double calc(double l,double r,double fl,double fr,double fm){
 71     return (r-l)/6.0*(fl+fr+4*fm);
 72 }
 73 double simpson(double l,double r,double mid,double fl,double fr,double fm,double s){
 74     double m1=(l+mid)/2.0,m2=(mid+r)/2.0;
 75     double f1=F(m1),f2=F(m2);
 76     double s1=calc(l,mid,fl,fm,f1),s2=calc(mid,r,fm,fr,f2);
 77     if(fabs(s-s1-s2)<eps)return s1+s2;
 78     return simpson(l,mid,m1,fl,fm,f1,s1)+simpson(mid,r,m2,fm,fr,f2,s2);
 79 }
 80 int del[N];
 81 double h[N],r[N],sum[N];
 82 int main(){
 83     scanf("%d%lf",&n,&alp);
 84     alp=tan(alp);
 85     for(int i=1;i<=n+1;i++)scanf("%lf",&h[i]),sum[i]=sum[i-1]+h[i];
 86     for(int i=1;i<=n;i++)scanf("%lf",&r[i]);
 87     for(int i=1;i<=n;i++)p[i].x=sum[i]/alp,p[i].r=r[i];
 88     p[n+1].x=sum[n+1]/alp;p[n+1].r=0;
 89     for(int i=1;i<=n;i++)if(fabs(p[i+1].x-p[i].x)>fabs(p[i+1].r-p[i].r))
 90         solve(p[i],p[i+1]);
 91     sort(l+1,l+cnt+1);
 92     for(int i=1;i<=n+1;i++)
 93         for(int j=1;j<=n+1;j++)if(i!=j&&!del[j])
 94             if(fabs(p[i].x-p[j].x)+p[i].r<=p[j].r){del[i]=1;break;}
 95     for(int i=1;i<=n+1;i++)if(!del[i])d[++tot]=p[i];
 96     n=tot;
 97     for(int i=1;i<=n;i++)p[i]=d[i];
 98     double l=p[1].x-p[1].r,r=p[n].x+p[n].r,mid=(l+r)/2.0;
 99     double fl=F(l),fr=F(r),fm=F(mid),s=calc(l,r,fl,fr,fm);
100     double ans=simpson(l,r,mid,fl,fr,fm,s);
101     printf("%0.2f\n",ans);
102     return 0;
103 }
View Code

我就找到这三道水题，感觉这个也就是用来骗分，掌握了也没有什么坏处。

posted @ 2018-02-24 12:23  Ren_Ivan  阅读(1894)  评论(0)  编辑  收藏  举报