BZOJ 1502 [NOI2005]月下柠檬树 自适应Simpson积分

题意:链接

方法:自适应Simpson积分

解析:

大半夜刷这道题真是爽歪歪。

求一棵树(圆锥加圆台组成)在平面上的投影的面积。

给定投影角度(0.3 < alpha <= pi/2)。

先来想想圆的投影是什么样子

还是他自己。

再想圆锥投影是什么样子

一个点加一个圆，并且有这个点与该圆的两条切线(该点在圆内部时没有切线)

再想圆台

两个圆，加上两个圆的外公切线组成的一坨图形。

不妨随意画一个。

好难画- -!

大概就转化成这个样子了。

观察这个图形…

轴对称啊- -!

这意味着我们好多东西都不用求了。

比如区间并就转化成求最大值辣。

好吧问题来了，公切线怎么求。

我知道我的做法肯定不太优越…

然而我还是要说2333

首先AC长是圆心距，可求。

AI长是半径差，可求。

所以CI可求。

连接FC，观察△FAC

2*S△FAC=FG*AC=CI*AF

AF为半径，已知。

所以FG可求。

于是AG可求。

A点坐标已知，所以F点坐标已知。

E点，直接相似即可。

代码:

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 1010
#define eps 1e-7
using namespace std;
int n;
double alpha;
struct Point
{
    double x,y;
    Point(){}
    Point(double _x,double _y):x(_x),y(_y){}
    Point operator + (const Point &a)
    {return Point(x+a.x,y+a.y);}
    Point operator - (const Point &a)
    {return Point(x-a.x,y-a.y);}
    Point operator * (double rate)
    {return Point(x*rate,y*rate);}
    double operator * (const Point &a)
    {return x*a.x+y*a.y;}
    double operator ^ (const Point &a)
    {return x*a.y-y*a.x;}
};
struct Line
{
    Point p,v;
    double k,b;
    Line(){}
    Line(Point _a,Point _b):p(_a),v(_b-_a)
    {
        k=((p+v).y-p.y)/((p+v).x-p.x);
        b=p.y-k*p.x;
    }
    double pt(double x)
    {
        return k*x+b;
    }
}line[N<<1];
struct Circle
{
    Point p;
    double R;
    Circle(){}
    Circle(Point _p,double _R):p(_p),R(_R){}
}cir[N];
struct Interval
{
    double l,r;
    Interval(){}
    Interval(double _l,double _r):l(_l),r(_r){}
    friend bool operator < (Interval a,Interval b)
    {
        if(a.l==b.l)return a.r<b.r;
        return a.l<b.l;
    }
}interval[N<<2];
int intervals,lines;
void Get_Lines(Circle a,Circle b)
{
    if(fabs(a.p.x-b.p.x)<fabs(a.R-b.R))return;
    if(a.R==b.R)
    {
        line[++lines]=Line(Point(a.p.x,a.R),Point(b.p.x,b.R));
    }else if(a.R<b.R)
    {
        double height1=sqrt((a.p.x-b.p.x)*(a.p.x-b.p.x)-(b.R-a.R)*(b.R-a.R));
        double height2=b.R*height1/(a.p.x-b.p.x);
        double length1=sqrt(b.R*b.R-height2*height2);
        Point aaaa=b.p+Point(length1,height2);
        Point bbbb=a.p+Point(length1*a.R/b.R,height2*a.R/b.R);
        line[++lines]=Line(aaaa,bbbb);
    }else
    {
        double height1=sqrt((a.p.x-b.p.x)*(a.p.x-b.p.x)-(a.R-b.R)*(a.R-b.R));
        double height2=a.R*height1/(a.p.x-b.p.x);
        double length1=sqrt(a.R*a.R-height2*height2);
        Point aaaaaa=a.p+Point(-length1,height2);
        Point bbbbbb=b.p+Point(-length1*b.R/a.R,height2*b.R/a.R);
        line[++lines]=Line(bbbbbb,aaaaaa);
    }
}
double getf(double x)
{
    double ans=0;
    for(int i=1;i<=n;i++)
    {
        double dis=fabs(x-cir[i].p.x);
        if(dis>cir[i].R)continue;
        double height=sqrt(cir[i].R*cir[i].R-dis*dis);
        ans=max(ans,height);
    }
    for(int i=1;i<=lines;i++)
    {
        if(x>=line[i].p.x&&x<=(line[i].p+line[i].v).x)
        {
            ans=max(ans,line[i].pt(x));
        }
    }
    return ans;
}
double calculate(double fl,double fm,double fr)
{
    return (fl+fr+4*fm)/6;
}
double Simpson(double l,double mid,double r,double fl,double fm,double fr,double s)
{
    double lm=(l+mid)/2,rm=(mid+r)/2;
    double flm=getf(lm),frm=getf(rm);
    double g1=calculate(fl,flm,fm)*(mid-l),g2=calculate(fm,frm,fr)*(r-mid);
    if(fabs(g1+g2-s)<eps)return s;
    else return Simpson(l,lm,mid,fl,flm,fm,g1)+Simpson(mid,rm,r,fm,frm,fr,g2);
}
double h[N],sum[N];
int main()
{
    scanf("%d%lf",&n,&alpha);
    for(int i=n+1;i>=1;i--)scanf("%lf",&h[i]);
    for(int i=n+1;i>=1;i--)sum[i]=sum[i+1]+h[i];
    for(int i=n;i>=1;i--)scanf("%lf",&cir[i].R);
    double l=0x7f7f7f7f,r=0;
    Point thetop(sum[1]*(1/tan(alpha)),0);
    r=max(r,thetop.x);
    for(int i=1;i<=n;i++)
    {
        if(i==1)
        {
            cir[1].p=Point(sum[2]*(1/tan(alpha)),0);
            l=min(cir[1].p.x-cir[1].R,l);
            r=max(cir[1].p.x+cir[1].R,r);
            if(cir[1].p.x+cir[1].R>=thetop.x)continue;
            double length=cir[1].R*cir[1].R/(thetop.x-cir[1].p.x);
            double height=sqrt(cir[1].R*cir[1].R-length*length);
            line[++lines]=Line(Point(cir[1].p.x+length,height),thetop);
        }else
        {
            cir[i].p=Point(sum[i+1]*(1/tan(alpha)),0);
            l=min(cir[i].p.x-cir[i].R,l);
            r=max(cir[i].p.x+cir[i].R,r);
            Get_Lines(cir[i-1],cir[i]);
        }
    }
    double mid=(l+r)/2;
    double fl=getf(l),fr=getf(r),fm=getf(mid);
    double ans=Simpson(l,mid,r,fl,fm,fr,calculate(fl,fm,fr)*(r-l));
    printf("%.2lf\n",2*ans);
}