自适应辛普森

数学没学好还真吃亏。

首先介绍下辛普森公式（没法帖图片，只能口述），对于一段函数，我们要求面积，除了求反导之外，有一种近似方法，即用辛普森公式。

（来自维基）

a,b分别指函数最左端和最右端，这个公式推导可以去维基上看，他有个直观的应用便是求拟柱体（比如棱柱）体积，f函数便代表截面面积，套一套求体积的公式或许会好理解一点。

仅仅是套公式是不行的，或许这三个点均是特殊点（比如极值），误差那不是一般的大。怎么办呢？我们可是学信息的，一边不行多套几边公式计算机还是办得到的。

下面介绍自适应辛普森

将一段区间二分，再将两端区间求出的面积合并，如此递归下去。为了保证误差足够小，如果两端用上述公式算出面积之和近似等于该区间用上述公式求出面积，便不再递归，返回值。

如此一来，即缩小误差，又节省了在函数图象较平整区域的计算时间。

题目参见：http://hi.baidu.com/cheezer94/blog/item/541e5cd5d761d8c38c1029e6.html

我写的不是很漂亮，卡精度才过的（主要是求f值时，每次我都得扫一遍o（n）），cheezer哲爷写的应该不错

uses math;const eps=0.00001;      eps1=0.01;//type real=extended;var x,y,mid,m,r,l,h:array[1..1000]of real;    n:longint;    ans,w,w0,hh:real;function dist(x1,x2,y1,y2:real):real;inline;begin exit(sqrt(sqr(x1-x2)+sqr(y1-y2)))end;function delt(b,a,c:real):real;inline;begin exit(sqrt(sqr(b)-4*a*c))end;function func(x:real):real;inline;var tmp,cosi,min,a,b,c,d,ne:real;    i:longint;begin func:=maxlongint;min:=maxlongint; for i:=1 to n do  if (-cos(l[i])<=-cos(x))and(-cos(x)<=-cos(h[i])) then begin   tmp:=abs(x-m[i]);   if tmp<eps then min:=mid[i]-r[i]   else begin    tmp:=abs(x-l[i]);    if tmp<eps then min:=sqrt(sqr(mid[i])-sqr(r[i]))    else begin     tmp:=abs(x-h[i]);     if tmp<eps then min:=sqrt(sqr(mid[i])-sqr(r[i]))     else begin      ne:=abs(x-m[i]);      a:=1;b:=-2*mid[i]*cos(ne);c:=sqr(mid[i])-sqr(r[i]);       d:=delt(b,a,c);      cosi:=(-b-d)/(2*a);      if cosi>=0 then min:=cosi     end    end   end;   if min<func then func:=min  end; tmp:=abs(x-pi/2); if (tmp>eps)and(x<pi/2)and(x<>0) then begin  min:=(w-w0)*tan(x);  min:=dist(w,w0,min,0) end else if (x>pi/2)and(tmp>eps) then begin  min:=-w0*tan(x);  min:=dist(w0,0,min,0) end else if x=0 then min:=w-w0; tmp:=abs(x-pi); if min<func then func:=min; if (abs(x-pi/2)>eps)and(x<>0)and(tmp>eps) then begin  min:=hh/tan(x)+w0;  min:=dist(w0,min,hh,0) end else if (x<>0)and(tmp>eps) then begin  min:=hh end; if min<func then func:=min; func:=(1/2)*sqr(func)end;function simpson(l,r,tmp1,tmp2,tmp3:real):real;inline;begin exit((r-l)*(tmp1+tmp3+tmp2*4)/6)end;function count(tmp1,tmp2,tmp3,l,r:real):real;inline;var m,x,tmpl,tmpr,tmp:real;begin m:=(l+r)/2; x:=simpson(l,r,tmp1,tmp2,tmp3); tmpl:=func((l+m)/2);tmpr:=func((m+r)/2); tmp:=simpson(l,m,tmp1,tmpl,tmp2)+simpson(m,r,tmp2,tmpr,tmp3); if abs(x-tmp)<eps1 then exit(tmp) else begin  count:=count(tmp1,tmpl,tmp2,l,m)+count(tmp2,tmpr,tmp3,m,r) endend;procedure init;var i:longint;    mm,ne,tmp1,tmp2,tmp3:real;begin readln(n,w,hh,w0); for i:=1 to n do begin  readln(x[i],y[i],r[i]);  mid[i]:=dist(x[i],w0,y[i],0);  if x[i]>w0 then mm:=arctan(y[i]/(x[i]-w0))  else if x[i]<w0 then mm:=pi-arctan(y[i]/(w0-x[i]))  else mm:=pi/2;  ne:=arcsin(r[i]/mid[i]);  l[i]:=mm-ne;h[i]:=mm+ne;m[i]:=mm end; tmp1:=func(0);tmp2:=func(pi/2);tmp3:=func(pi); ans:=count(tmp1,tmp2,tmp3,0,pi); writeln(ans/(w*hh):0:8)end;beginassign(input,'yc.in');reset(input);assign(output,'yc.out');rewrite(output); init;close(output);close(input)end.