题目链接: BZOJ - 1502
题目分析
这是我做的第一道 Simpson 积分的题目。Simpson 积分是一种用 (fl + 4*fmid + fr) / 6 * (r - l) 来拟合 fl...fr 的方法。自适应 Simpson 的自适应指的是，如果分成左右两端分别 Simpson 的和与对整段 Simpson 的差值不超过一个 Eps，那么就接受这个值，否则递归下去求两段的 Simpson 值。
代码
1 #include <iostream> 2 #include <cstdlib> 3 #include <cstring> 4 #include <cstdio> 5 #include <cmath> 6 #include <algorithm> 7 8 using namespace std; 9 10 typedef double LF; 11 12 const LF Eps = 1e-8; 13 14 inline LF gmin(LF a, LF b) {return a < b ? a : b;} 15 inline LF gmax(LF a, LF b) {return a > b ? a : b;} 16 inline LF Sqr(LF x) {return x * x;} 17 18 const int MaxN = 500 + 5; 19 20 int n; 21 22 LF Alpha, talpha, Ht, Lp, Rp, Ans; 23 LF A[MaxN], P[MaxN], Rad[MaxN], Ls[MaxN], Rs[MaxN], Lh[MaxN], Rh[MaxN]; 24 25 inline LF f(LF x) 26 { 27 LF ret = 0.0; 28 for (int i = 1; i <= n; ++i) 29 { 30 if (fabs(x - P[i]) < Rad[i]) ret = gmax(ret, sqrt(Sqr(Rad[i]) - Sqr(x - P[i]))); 31 if (x > Ls[i] && x < Rs[i]) ret = gmax(ret, Lh[i] + (Rh[i] - Lh[i]) * (x - Ls[i]) / (Rs[i] - Ls[i])); 32 } 33 return ret; 34 } 35 36 inline LF Simpson(LF l, LF r, LF fl, LF fmid, LF fr) 37 { 38 return (fl + fmid * 4.0 + fr) / 6.0 * (r - l); 39 } 40 41 LF RSimpson(LF l, LF r, LF fl, LF fmid, LF fr) 42 { 43 LF mid, p, q, x, y, z; 44 mid = (l + r) / 2.0; 45 p = f((l + mid) / 2.0); 46 q = f((mid + r) / 2.0); 47 x = Simpson(l, r, fl, fmid, fr); 48 y = Simpson(l, mid, fl, p, fmid); 49 z = Simpson(mid, r, fmid, q, fr); 50 if (fabs(x - y - z) < Eps) return y + z; 51 else return RSimpson(l, mid, fl, p, fmid) + RSimpson(mid, r, fmid, q, fr); 52 } 53 54 int main() 55 { 56 scanf("%d%lf", &n, &Alpha); 57 talpha = tan(Alpha); 58 Ht = 0.0; 59 for (int i = 1; i <= n + 1; ++i) 60 { 61 scanf("%lf", &A[i]); 62 Ht += A[i]; 63 P[i] = Ht / talpha; 64 } 65 Lp = P[1]; Rp = P[n + 1]; 66 for (int i = 1; i <= n; ++i) 67 { 68 scanf("%lf", &Rad[i]); 69 Lp = gmin(Lp, P[i] - Rad[i]); 70 Rp = gmax(Rp, P[i] + Rad[i]); 71 } 72 Rad[n + 1] = 0.0; 73 for (int i = 1; i <= n; ++i) 74 { 75 if (P[i + 1] - P[i] > fabs(Rad[i + 1] - Rad[i])) 76 { 77 Ls[i] = P[i] + Rad[i] * (Rad[i] - Rad[i + 1]) / (P[i + 1] - P[i]); 78 Rs[i] = P[i + 1] + Rad[i + 1] * (Rad[i] - Rad[i + 1]) / (P[i + 1] - P[i]); 79 Lh[i] = sqrt(Sqr(Rad[i]) - Sqr(Ls[i] - P[i])); 80 Rh[i] = sqrt(Sqr(Rad[i + 1]) - Sqr(Rs[i] - P[i + 1])); 81 } 82 else 83 { 84 Ls[i] = -1; 85 Rs[i] = -1; 86 } 87 } 88 Ans = RSimpson(Lp, Rp, f(Lp), f((Lp + Rp) / 2.0), f(Rp)) * 2; 89 printf("%.2lf ", Ans); 90 return 0; 91 }