hdu4498 求曲线长度 数值积分 simpson公式
http://acm.hdu.edu.cn/showproblem.php?pid=4498
simpson公式可以求连续曲线的定积分。
这里注意解ax^2+bx+c=0时a=0,退化为一次的情况。
#include<iostream> #include<cstdio> #include<cstring> #include<cstdlib> #include<algorithm> #include<math.h> #define REP(i,a,b) for(int i=a;i<=b;i++) #define MS0(a) memset(a,0,sizeof(a)) using namespace std; typedef long long ll; const int maxn=1000100; const int INF=1e9+10; const double EPS=0.0000000001; int n; struct Line { double k,a,b; void read() { scanf("%lf%lf%lf",&k,&a,&b); } double y(double x) { return k*(x-a)*(x-a)+b; } };Line L[maxn]; struct Point { double x,y; friend bool operator<(Point A,Point B) { return A.x<B.x; } void out() { printf("x=%.2f y=%.2f\n",x,y); } };Point p[maxn];int pn; bool can(double x,double y) { if(x<=0||x>=100.0) return 0; REP(i,0,n){ if(L[i].y(x)<y-EPS) return 0; } return 1; } void work(double a,double b,double c,Line A) { double dt=b*b-4*a*c; if(dt<-EPS) return; if(fabs(a)<EPS){ if(fabs(b)>EPS){ double x=-c/b; if(can(x,A.y(x))) p[++pn]={x,A.y(x)}; } return; } double x1=(-b-sqrt(dt))/(2*a); double x2=(-b+sqrt(dt))/(2*a); if(fabs(dt)<EPS){ if(!can(x1,A.y(x1))) return; p[++pn]={x1,A.y(x1)}; } else{ if(can(x1,A.y(x1))) p[++pn]={x1,A.y(x1)}; if(can(x2,A.y(x2))) p[++pn]={x2,A.y(x2)}; } } void inter(Line A,Line B) { double a=A.k-B.k; double b=-2*(A.k*A.a-B.k*B.a); double c=A.k*A.a*A.a-B.k*B.a*B.a+A.b-B.b; work(a,b,c,A);///解ax^2+bx+c=0,并取交点 } void get_duan() { double y=100.0; REP(i,1,n) y=min(y,L[i].y(0)); p[++pn]={0,y}; y=100.0; REP(i,1,n) y=min(y,L[i].y(100.0)); p[++pn]={100.0,y}; } void get_inter() { pn=0; REP(i,0,n){ REP(j,i+1,n){ inter(L[i],L[j]);///取两条曲线的交点 } } ///取端点0和100 get_duan(); sort(p+1,p+pn+1); } double F(double x,double k,double a,double b) { double t=2*k*(x-a); return sqrt(1+t*t); } /// simpson求定积分 double simpson(double a,double b,double tk,double ta,double tb) { double c=a+(b-a)/2; return (F(a,tk,ta,tb)+4*F(c,tk,ta,tb)+F(b,tk,ta,tb))*(b-a)/6; } double asr(double a,double b,double EPS,double A,double tk,double ta,double tb) { double c=a+(b-a)/2; double L=simpson(a,c,tk,ta,tb),R=simpson(c,b,tk,ta,tb); if(fabs(L+R-A)<=15*EPS) return L+R+(L+R-A)/15.0; return asr(a,c,EPS/2,L,tk,ta,tb)+asr(c,b,EPS/2,R,tk,ta,tb); } double asr(double a,double b,double EPS,double tk,double ta,double tb) { return asr(a,b,EPS,simpson(a,b,tk,ta,tb),tk,ta,tb); } ///--- bool In(Line L,Point A) { if(fabs(L.y(A.x)-A.y)<EPS) return 1; return 0; } double Len(Point A,Point B) { if(fabs(A.x-B.x)<EPS) return 0; double mx=(A.x+B.x)/2; int cur=-1; REP(i,0,n){ if(!In(L[i],A)||!In(L[i],B)) continue; if(cur==-1||L[i].y(mx)<L[cur].y(mx)) cur=i; } if(cur==-1) return 0; return asr(A.x,B.x,EPS,L[cur].k,L[cur].a,L[cur].b); } double solve() { get_inter();///取交点 /* REP(i,1,pn){ p[i].out(); } */ double res=0; REP(i,1,pn-1){ res+=Len(p[i],p[i+1]);///计算相邻两点的曲线长度 } return res; } int main() { freopen("in.txt","r",stdin); int T;cin>>T; while(T--){ scanf("%d",&n); L[0].k=L[0].a=0;L[0].b=100.0; REP(i,1,n) L[i].read(); printf("%.2f\n",solve()); } return 0; }
没有AC不了的题,只有不努力的ACMER!