首页 > 代码库 > 【计算几何】【辛普森积分法】UVALive - 7076 - Highway
【计算几何】【辛普森积分法】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; }
【计算几何】【辛普森积分法】UVALive - 7076 - Highway
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。