[BZOJ 1502] [NOI2005] 月下柠檬树 【Simpson积分】

题目链接: 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\n", Ans);
90     return 0;
91 }

 

posted @ 2015-04-02 10:08  JoeFan  阅读(364)  评论(0)  编辑  收藏  举报