BZOJ1502: [NOI2005]月下柠檬树【Simpson积分】

dzy posted @ 2013年6月03日 21:07 in BZOJ with tags NOI 计算几何 bzoj , 1273 阅读

开始题目什么意思都没看懂。后来明白了是求一些圆和等腰梯形的面积并。

先根据样例画个图。明确算的是哪一部分。

样例图片

最后要求的就是这些图形的面积并了。怎么求呢?需要用到simpson积分公式。由于我太弱所以就讲的意识流一些。或许不是那么严谨吧。

首先要预处理出这些等腰梯形。他们是由相邻两圆的公切线(如果不存在就不用算了)和两条过切点垂直于x轴的线段组成的。这4个点可以通过相似算出。

重点:simpson(l,r)=(f(l)+4f(mid)+f(r))*(r-l)/6. 其中f(x)当然是一个分段函数.每次O(n)可以算出(枚举所有的等腰梯形和圆)

对于定义在[l,r]上的连续函数,求它和x轴、直线x=l、直线x=r围成的面积S=rsimpson(l,r)。当simpson(l,r)=simpson(l,mid)+simpson(mid,r)时,返回simpson(l,r);否则返回rsimpson(l,mid)+rsimpson(mid,r)。显然,答案就是2*simpson(l,r)。

【拓展】

辛普森的简化公式可以用来算体积V:V=h(a+4b+c)/6。

h是立体(常指拟柱体)的高度,a是下底面积,b是中间截面面积(在一半高度上的截面面积),c是上底面积。

 

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double eps=1e-6;
double x[505],r[505],sx[505],sy[505],tx[505],ty[505],L,R,af,t=0;int n;
double f(double xx){
    double fx=0;
    for(int i=1;i<=n+1;i++) if(fabs(xx-x[i])<=r[i]) fx=max(fx,sqrt(r[i]*r[i]-(xx-x[i])*(xx-x[i])));
    for(int i=1;i<=n;i++)   if(x[i+1]-x[i]-fabs(r[i+1]-r[i])>eps && sx[i]<=xx && xx<=tx[i]){
        if(sy[i]<ty[i]) fx=max(fx,sy[i]+(ty[i]-sy[i])*(xx-sx[i])/(tx[i]-sx[i]));
        else            fx=max(fx,ty[i]+(sy[i]-ty[i])*(tx[i]-xx)/(tx[i]-sx[i]));
	}return fx;
}double simpson(double l,double r){
    return (f(l)+f(r)+4*f((l+r)/2.0))*(r-l)/6.0;
}double rsimpson(double l,double r){
    double mid=(l+r)/2.0;
    if(fabs(simpson(l,r)-simpson(l,mid)-simpson(mid,r))<eps)	return simpson(l,mid)+simpson(mid,r);
    return rsimpson(l,mid)+rsimpson(mid,r);
}int main(){
    scanf("%d%lf",&n,&af);  af=1.0/tan(af);
    for(int i=1;i<=n+1;i++) scanf("%lf",&x[i]),t+=x[i],x[i]=t*af;
    for(int i=1;i<=n;i++)   scanf("%lf",&r[i]);L=x[1],R=x[n+1],r[n+1]=0;
    for(int i=1;i<=n;i++){
        L=min(L,x[i]-r[i]), R=max(R,x[i]+r[i]);
        if(x[i+1]-x[i]-fabs(r[i+1]-r[i])>eps){
            double o=(r[i+1]-r[i])/(x[i+1]-x[i]);
            sx[i]=x[i]-r[i]*o,  sy[i]=sqrt(r[i]*r[i]*(1-o*o));
            tx[i]=x[i+1]-r[i+1]*o,  ty[i]=sqrt(r[i+1]*r[i+1]*(1-o*o));
        }
    }printf("%.2lf\n",2*rsimpson(L,R));	return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Copyright © 2007

Webdesign, tvorba www stránek

Valid XHTML 1.1