萌新刚学 simpson,求助 mle
查看原帖
萌新刚学 simpson,求助 mle
455558
Imiya楼主2022/4/5 19:36

rt qwq

#include<iostream>
#include<cmath>
using namespace std;
double R;
inline double f(double x){return sqrt(R*R-x*x);}
inline double simpson(double l,double r){
    return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}
double area(double l,double r,double eps,double ans){
    double mid=(l+r)/2;
    double area1=simpson(l,mid);
    double area2=simpson(mid,r);
    if(fabs(area1+area2-ans)<eps)return area1+area2;
    else return area(l,mid,eps,area1)+area(mid,r,eps,area2);
}
struct node{double x,y;};
struct edge{double k,b;};
inline edge New(node x,node y){
    double k=(y.y-x.y)/(y.x-x.x);
    double b=x.y-k*x.x;
    return{k,b};
}
inline node isec(edge l1,edge l2){
    double x=(l2.b-l1.b)/(l1.k-l2.k);
    double y=l1.k*x+l1.b;
    return {x,y};
}
inline node get(double r1,double r2,double dis){
    double C=(r1-r2)/dis;
    double S=sqrt(1-C*C);
    return {C*r1,S*r1};
}
double esp=1e-11;
inline node isec2(double r,double h,double k,double b){
    double derta=4*(k*k*(r*r-h*h)-b*b+r*r-2*k*b*h);
    derta=sqrt(derta+esp);
    double x=max((derta-2*k*b+2*h)/(2*k*k+2),(-derta-2*k*b+2*h)/(2*k*k+2));
    return {x,k*x+b};
}//求圆与直线的交点
const int N=511;
int n;
double alpha,h[N],r[N],ans;
node lp[N],rp[N];
edge lst[N];
void init(){
    cin>>n>>alpha;
    for(int i=1;i<=n+1;i++)cin>>h[i],h[i]/=tan(alpha);
    h[1]=0;
    for(int i=1;i<=n+1;i++)h[i]+=h[i-1];
    for(int i=1;i<=n;i++)cin>>r[i];
}
int main(){
//    freopen("read.in","r",stdin);
    init();
    lp[1]={-r[1],0};
    int p=1;
    for(int i=2;i<=n+1;i++){
        if(h[i]+r[i]<=h[p]+r[p]){
            lp[i]=rp[i]={-1,-1};
            continue;
        }
        lp[i]=get(r[i],r[i-1],h[i-1]-h[i]);
        node t=get(r[i-1],r[i],h[i]-h[i-1]);
        lp[i].x+=h[i];
        t.x+=h[i-1];
        lst[i]=New(lp[i],t);
        rp[p]=isec2(r[p],h[p],lst[i].k,lst[i].b);
        if(rp[p].x<lp[p].x){
            rp[p]=lp[p]=isec(lst[p],lst[i]);
        }
        p=i;
    }
    if(lp[n+1].x==-1)rp[p]={h[p]+r[p],0};
    else lp[n+1]=rp[n+1]={h[n+1],0};
//    for(int i=1;i<=n+1;i++)cout<<lp[i].x<<' '<<lp[i].y<<' '<<rp[i].x<<' '<<rp[i].y<<' '<<h[i]<<endl;
    int m=0;
    for(int i=1;i<=n+1;i++)if(lp[i].x!=-1&&lp[i].y!=-1)lp[++m]=lp[i],rp[m]=rp[i],h[m]=h[i],r[m]=r[i];
    for(int i=1;i<m;i++){
        R=r[i];
        ans+=area(lp[i].x-h[i],rp[i].x-h[i],1e-4,0)*2;
        ans+=(rp[i].y+lp[i+1].y)*(lp[i+1].x-rp[i].x);
    }
    R=r[m];
    ans+=area(lp[m].x-h[m],rp[m].x-h[m],1e-4,0)*2;
    printf("%.2lf",ans);
    return 0;
}
2022/4/5 19:36
加载中...