数值积分 (辛普森公式 辛普森自适应法则) UVA

double F(double x) {
}
double simpson(double l, double r) {
double m = (l + r) * 0.5;
return (F(l) + 4 * F(m) + F(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double A) {
double m = (l + r) * 0.5;
double L = simpson(l, m), R = simpson(m, r);
if (fabs(L + R - A) <= 15 * eps)
return L + R + (L + R - A) / 15.0;
return asr(l, m, eps * 0.5, L) + asr(m, r, eps * 0.5, R);
}
double asr(double l, double r, double eps) {
return 2.0 * asr(l, r, eps, simpson(l, r));
}
int dblcmp(double x) {
if (x < -eps)
return -1;
return x > eps;
}

#include<cstdio>
#include<cstring>
#include<queue>
#include<vector>
#include<iostream>
#include<string>
#include<sstream>
#include<cctype>
#include<set>
#include<stack>
#include<numeric>
#include<functional>
#include<memory>
#include<deque>
#include<list>
#include<cmath>
#include<fstream>
#include<cstdlib>
#include<climits>
#include<iomanip>
#include<cstring>
#include<memory>
#include<bitset>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<VVI> VVVI;
const double eps = 1e-10;
const int maxn = 1111;
const int mod = 1000000007;
const int N = 1111;
const int M = 111111;
const int INF = 0x7f7f7f7f;
double a, w, h;
double F(double x) {
return sqrt(1 + 4.0 * a * a * x * x);
}
double simpson(double l, double r) {
double m = (l + r) * 0.5;
return (F(l) + 4 * F(m) + F(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double A) {
double m = (l + r) * 0.5;
double L = simpson(l, m), R = simpson(m, r);
if (fabs(L + R - A) <= 15 * eps)
return L + R + (L + R - A) / 15.0;
return asr(l, m, eps * 0.5, L) + asr(m, r, eps * 0.5, R);
}
double asr(double l, double r, double eps) {
return 2.0 * asr(l, r, eps, simpson(l, r));
}
int dblcmp(double x) {
if (x < -eps)
return -1;
return x > eps;
}
int main() {
double D, H, B, L;
int cas;
int n, T;
double L1;
scanf("%d", &T);
for (cas = 1; cas <= T; ++cas) {
scanf("%lf%lf%lf%lf", &D, &H, &B, &L);
if (cas > 1)
puts("");
n = (int) ceil(B * 1.0 / D);
D = B * 1.0 / n;
L1 = L * 1.0 / n;
double len, l, r, m;
l = 0, r = H;
w = D;
while (r - l > eps) {
m = (l + r) * 0.5;
h = m;
a = 4.0 * h / (w * w);
len = asr(0, w * 0.5, 1e-6);
if (dblcmp(len - L1) < 0)
l = m;
else
r = m;
}
printf("Case %d:\n%.2lf\n", cas, 0.0 + H - h);
}
return 0;
}