hdu4498 数值积分

写了个牛顿-科茨积分模板，然后求出抛物线再[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;
}