Simpson公式的应用（HDU 1724/ HDU 1071）

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>

using namespace std;

double coe[4][11];
const double EPS = 1e-8;

int k;
double cal(double x, double *c) {
double ret = c[0];
for (int i = 1; i <= k; i++) ret *= x, ret += c[i];
return ret;
}

inline double cal(double x, double *p, double *q) { return cal(x, p) / cal(x, q);}
inline double cal(double x, double y, double *p, double *q) { return max(cal(x, p, q) - y, 0.0);}
inline double simpson(double y, double l, double r, double *p, double *q) { return (cal(l, y, p, q) + cal(r, y, p, q) + 4 * cal((l + r) / 2, y, p, q)) * (r - l) / 6;}

inline double getpart(double y, double l, double r, double *p, double *q) {
double sum = simpson(y, l, r, p, q);
//cout << l << ' ' << r << ' ' << sum << endl;
if (fabs(sum - simpson(y, l, (l + r) / 2, p, q) - simpson(y, (l + r) / 2, r, p, q)) < EPS) return sum;
return getpart(y, l, (l + r) / 2, p, q) + getpart(y, (l + r) / 2, r, p, q);
}

inline double getarea(double y, double l, double r, double *p, double *q) {
double ret = 0, d = (r - l) / 100;
for (int i = 0; i < 100; i++) {
ret += getpart(y, l + d * i, l + d * (i + 1), p, q);
}
return ret;
}

double dc2(double l, double r, double a, double w) {
double m;
while (r - l > EPS) {
m = (l + r) / 2.0;
//cout << m << ' ' << getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) << endl;
if (getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) > a) l = m;
else r = m;
}
return l;
}

int main() {
//freopen("in", "r", stdin);
//freopen("out", "w", stdout);
double w, d, a;
while (cin >> w >> d >> a >> k) {
for (int i = 0; i < 4; i++) for (int j = 0; j <= k; j++) cin >> coe[i][j];
for (int i = 0; i < 4; i++) reverse(coe[i], coe[i] + k + 1);
//cout << getarea(-5.51389, 0, w, coe[0], coe[1]) - getarea(-5.51389, 0, w, coe[2], coe[3]) << endl;
//cout << cal(3, coe[0], coe[1]) << endl;
printf("%.5f\n", -dc2(-d, 0, a, w));
}
return 0;
}