【计算几何】【辛普森积分法】UVALive - 7076 - Highway

#include<cstdio>
#include<cmath>
using namespace std;
#define EPS 0.00000000001
const double PI=acos(-1.0);
double v0,v1,D,T;
double x0,R,x3,x1,Y1,k,b,A,B,C,x2,y2,k2,b2;
double f(double x)
{
	return k*x+b-(sqrt(R*R-x*x)-D);
}
double area(double l,double r)
{
	return (r-l)/6.0*(f(l)+4.0*f((l+r)/2.0)+f(r));
}
double Sinp(double l,double r)
{
	double m=(l+r)*0.5;
	double t2=area(l,r),t3=area(l,m)+area(m,r);
	if(fabs(t2-t3)<EPS)
	  return t2;
	return Sinp(l,m)+Sinp(m,r);
}

double g(double x)
{
	return -sqrt(R*R-x*x)+D;
}
double are2(double l,double r)
{
	return (r-l)/6.0*(g(l)+4.0*g((l+r)/2.0)+g(r));
}
double Sin2(double l,double r)
{
	double m=(l+r)*0.5;
	double t2=are2(l,r),t3=are2(l,m)+are2(m,r);
	if(fabs(t2-t3)<EPS)
	  return t2;
	return Sin2(l,m)+Sin2(m,r);
}

double w(double x)
{
	return -sqrt(R*R-x*x)-D-(k*x+b);
}
double are3(double l,double r)
{
	return (r-l)/6.0*(w(l)+4.0*w((l+r)/2.0)+w(r));
}
double Sin3(double l,double r)
{
	double m=(l+r)*0.5;
	double t2=are3(l,r),t3=are3(l,m)+are3(m,r);
	if(fabs(t2-t3)<EPS)
	  return t2;
	return Sin3(l,m)+Sin3(m,r);
}

double sqr(double x)
{
	return x*x;
}
int main()
{
	//freopen("a.in","r",stdin);
	int zu=0;
	while(scanf("%lf%lf%lf%lf",&v0,&v1,&D,&T)!=EOF)
	  {
	  	++zu;
	  	x0=-D*v1/v0+v1*T;
	  	R=v0*T;
	  	x3=sqrt(R*R-D*D);
	  	x1=sqr(T*v0-D)/x0;
	  	Y1=sqrt(sqr(T*v0-D)-x1*x1);
	  	k=Y1/(x1-x0);
	  	b=R*sqrt(1.0+k*k)-D;
	  	x0=-b/k;
	  	k2=-1.0/k;
	  	b2=-D;
	  	A=1.0+k2*k2;
	  	B=2.0*(b2+D)*k2;
	  	C=sqr(b2+D)-R*R;
	  	x2=(-B+sqrt(B*B-4.0*A*C))/(2.0*A);
	  	y2=k*x2+b;
	  	if(y2<EPS)
	  	  {
	  	  	printf("Case %d: %.8lf\n",zu,PI*R*R);
	  	  	continue;
	  	  }
	  	double __area=Sinp(x2,x3)+(x0-x3)*(k*x3+b)*0.5;
	  	  	
	  	k=-k;
		b=-b;
	  	A=1.0+k*k;
	  	B=2.0*(b+D)*k;
	  	C=sqr(b+D)-R*R;
	  	x2=(-B+sqrt(B*B-4.0*A*C))/(2.0*A);
	  	y2=k*x2+b;
	  	if(y2+D>(-EPS))
	  	  __area+=(Sin2(x3,x2)+(x0-x2)*(k*x2+b)*(-0.5));
	  	else
	  	  __area+=(Sin2(x3,R)+Sin3(x2,R)+(x0-R)*(k*R+b)*(-0.5));
	  	printf("Case %d: %.8lf\n",zu,__area*2.0+PI*R*R);
	  }
	return 0;
}