hdu4498 数值积分

来源：互联网发布：信用卡账单伪造软件编辑：程序博客网时间：2024/05/21 22:30
写了个牛顿-科茨积分模板，然后求出抛物线再[0,100]之间的所有交点分段积分。
ACcode：
#include<cstdio>#include<cstring>#include<cmath>#include<vector>#include<algorithm>#include<iostream>#include<set>using namespace std;typedef long long LL;typedef double db;const int NS=55;const db eps=1e-8;int n,id;db a[NS],b[NS],k[NS];set<db> ins;set<db>::iterator iter;void cross(){    ins.clear();    for (int i=0;i<n;i++)    for (int j=i+1;j<=n;j++)    {        db A=k[i]-k[j];        db B=-2.0*(k[i]*a[i]-k[j]*a[j]);        db C=b[i]-b[j];        C+=k[i]*a[i]*a[i]-k[j]*a[j]*a[j];        if (fabs(A)<eps)        {            if (fabs(B)>=eps)            {                db x1=-C/B;                ins.insert(x1);            }            continue;        }        db det=B*B-4.0*A*C;        if (det>=eps&&fabs(A)>=eps)        {            db x1=(-B+sqrt(det))/A/2.0;            db x2=(-B-sqrt(det))/A/2.0;            ins.insert(x1),ins.insert(x2);        }    }    iter=ins.begin();    while(iter!=ins.end())    {        if(*iter<eps || *iter+eps>100.0)        {            ins.erase(iter);            iter=ins.begin();            continue;        }        iter++;    }}double func(double x){    db t=2.0*k[id]*(x-a[id]);    return sqrt(1.0+t*t);}double simpson(double a,double b,int n=500){    double ret=func(a)-func(b);    double h=0.5*(b-a)/n,p=a;    for (int i=0;i<n;i++)    {        p+=h; ret+=func(p)*4.0;        p+=h; ret+=func(p)*2.0;    }    return ret*h/3.0;}double cotes(double a,double b,int n=40){    double ret=7.0*(func(a)-func(b));    double h=0.25*(b-a)/n,p=a;    for (int i=0;i<n;i++)    {        p+=h; ret+=32.0*func(p);        p+=h; ret+=12.0*func(p);        p+=h; ret+=32.0*func(p);        p+=h; ret+=14.0*func(p);    }    return ret*h*2.0/45.0;}int main(){    int T;    a[0]=k[0]=0.0,b[0]=100.0;    for (scanf("%d",&T);T--;)    {        scanf("%d",&n);        for (int i=1;i<=n;i++)            scanf("%lf%lf%lf",&k[i],&a[i],&b[i]);        cross();        ins.insert(a[0]),ins.insert(b[0]);        iter=ins.begin();        db l=0,r,ans=0.0;        for (iter++;iter!=ins.end();iter++)        {            r=*iter,id=0;            db yy=100.0,mid=(l+r)/2.0;            for (int i=1;i<=n;i++)            {                db t=k[i]*(mid-a[i])*(mid-a[i])+b[i];                if (t<yy) yy=t,id=i;            }            ans+=cotes(l,r); l=r;        }        printf("%.2lf\n",ans);    }    return 0;}