P4207 WA on #9

#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;

const int MAXN=5e2+10;
const double eps=1e-5;

struct Seg{
	double lx,ly,rx,ry;
}s[MAXN];

double f(Seg i,double x){
	double lx=i.lx,rx=i.rx,ly=i.ly,ry=i.ry;
	if(x<lx||x>rx){
		return 0;
	}
	return (x*(ry-ly)+ly*rx-lx*ry)/(rx-lx);
}

struct C{
	double O,r=0;
}c[MAXN];

double f(C i,double x){
	double O=i.O,r=i.r;
	if(x<O-r||x>O+r){
		return 0;
	}
	return sqrt(r*r-(O-x)*(O-x));
}

int num1,num2;

double f(double x){
	double re=0;
	for(int i=1;i<=num1;i++){
		re=max(re,f(s[i],x));
	}
	for(int i=1;i<=num2;i++){
		re=max(re,f(c[i],x));
	}
	return re;
}

double simpson(double l,double r,double fl,double fr){
	double mid=l+(r-l)/2;
	double fm=f(mid);
	return (fl+4*fm+fr)*(r-l)/6;
}

double asr(double l,double r,double ans,double eps,double fl,double fr){
	double mid=l+(r-l)/2;
	double fm=f(mid);
	double l_ans=simpson(l,mid,fl,fm);
	double r_ans=simpson(mid,r,fm,fr);
	if(fabs(ans-l_ans-r_ans)<15*eps){
		return l_ans+r_ans+(l_ans+r_ans-ans)/15;
	}else{
		l_ans=asr(l,mid,l_ans,eps/2,fl,fm);
		r_ans=asr(mid,r,r_ans,eps/2,fm,fr);
		return l_ans+r_ans;
	}
}

double asr(double l,double r,double eps){
	double fl=f(l);
	double fr=f(r);
	return asr(l,r,simpson(l,r,fl,fr),eps,fl,fr);
}

int n;
double d[MAXN],r[MAXN],alp;

int main(){
	cin>>n>>alp;
	alp=1.0/tan(alp);
	for(int i=1;i<=n+1;i++){
		cin>>d[i];
		d[i]+=d[i-1];
	}
	for(int i=1;i<=n;i++){
		cin>>r[i];
		num2++;
		c[num2].r=r[i];
		c[num2].O=d[i]*alp;
	}
	num2++;
	c[num2].O=d[n+1]*alp;
	for(int i=1;i<=n;i++){
		num1++;
		if(fabs(c[i].r-c[i+1].r)<eps){
			s[i].lx=c[i].O,s[i].rx=c[i+1].O;
			s[i].ly=s[i].ry=c[i].r;
			continue;
		}
		double co=-(c[i].r-c[i+1].r)/(c[i].O-c[i+1].O);
		if(fabs(co)>1){
			continue;
		}
		double si=sqrt(1-co*co);
		s[i].lx=c[i].O+c[i].r*co,s[i].ly=si*c[i].r;
		s[i].rx=c[i+1].O+c[i+1].r*co,s[i].ry=si*c[i+1].r;
	}
//	for(int i=1;i<=num1;i++){
//		cout<<c[i].O<<' '<<c[i].r<<' '<<c[i+1].O<<' '<<c[i+1].r<<' '<<s[i].lx<<' '<<s[i].ly<<' '<<s[i].rx<<' '<<s[i].ry<<endl;
//		if(fabs(c[i].O-c[i+1].O)<fabs(c[i].r-c[i+1].r)){
//			cout<<1<<endl;
//		}else{
//			cout<<0<<endl;
//		}
//	}
	cout<<fixed<<setprecision(2)<<asr(c[1].O-c[1].r,max(c[num2].O,c[num2-1].O+c[num2-1].r),eps)*2.0;
	
	return 0;
}