UVALive - 6135 Environment Protection 二分+simpson

题目：给你宽为W，深度为D的一块区域，再给你p1,q1,p2,q2；y1=p1/q1,y2=p2/q2，让你确定一个深度d,使得y1y2之间的面积在W*d这个区域内为A.

思路：二分+辛普森

代码：

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<iostream>
#include<algorithm>
#include<ctime>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<list>
#include<numeric>
using namespace std;
#define LL long long
#define ULL unsigned long long
#define INF 0x3f3f3f3f
#define mm(a,b) memset(a,b,sizeof(a))
#define PP puts("*********************");
template<class T> T f_abs(T a){ return a > 0 ? a : -a; }
template<class T> T gcd(T a, T b){ return b ? gcd(b, a%b) : a; }
template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
// 0x3f3f3f3f3f3f3f3f
// 0x3f3f3f3f

const double eps=1e-6;//根据实际情况修改
double arr[4][10];
double W,D,A,mid;
int K;
double F(double x,double *p,double *q){
	//Simpson公式用到的函数
    double fz=0,fm=0;
    for(int i=0;i<=K;i++){
        fz=x*fz+p[i];
        fm=x*fm+q[i];
    }
    return max(mid,fz/fm);
}
double simpson(double a, double b,double *p,double *q){
	double c = a + (b - a) / 2;
	return (F(a,p,q) + 4 * F(c,p,q) + F(b,p,q))*(b - a) / 6;
}
double asr(double a, double b, double eps, double A,double *p,double *q){
	double c = a + (b - a) / 2;
	double L = simpson(a, c,p,q), R = simpson(c, b,p,q);
	if (fabs(L + R - A) <= 15 * eps) return L + R + (L + R - A) / 15.0;
	return asr(a, c, eps / 2, L,p,q) + asr(c, b, eps / 2, R,p,q);
}
double asr(double a, double b, double eps,double *p,double *q){
	return asr(a, b, eps, simpson(a, b,p,q),p,q);
}
int main(){

    while(~scanf("%lf%lf%lf%d",&W,&D,&A,&K)){
        for(int i=0;i<=3;i++)
            for(int j=K;j>=0;j--)
                scanf("%lf",&arr[i][j]);
        double l=-D,r=0,ans;
        while(r-l>eps){
            mid=(l+r)/2;
            double area=asr(0,W,eps,arr[0],arr[1])-asr(0,W,eps,arr[2],arr[3]);
            if(area>=A){
                ans=mid;
                l=mid;
            }
            else r=mid;
        }
        printf("%.5f\n",f_abs(ans));
    }
    return 0;
}