1. 程式人生 > >[BZOJ1502]月下檸檬樹(自適應辛普森積分)

[BZOJ1502]月下檸檬樹(自適應辛普森積分)

成了 typedef 高度 i++ ble 4.0 技術 https SQ

1502: [NOI2005]月下檸檬樹

Time Limit: 5 Sec Memory Limit: 64 MB
Submit: 1387 Solved: 739
[Submit][Status][Discuss]

Description

李哲非常非常喜歡檸檬樹,特別是在靜靜的夜晚,當天空中有一彎明月溫柔地照亮地面上的景物時,他必會悠閑地 坐在他親手植下的那棵檸檬樹旁,獨自思索著人生的哲理。李哲是一個喜愛思考的孩子,當他看到在月光的照射下 檸檬樹投在地面上的影子是如此的清晰,馬上想到了一個問題:樹影的面積是多大呢?李哲知道,直接測量面積是 很難的,他想用幾何的方法算,因為他對這棵檸檬樹的形狀了解得非常清楚,而且想好了簡化的方法。李哲將整棵 檸檬樹分成了n 層,由下向上依次將層編號為1,2,…,n。從第1到n-1 層,每層都是一個圓臺型,第n 層(最上面一 層)是圓錐型。對於圓臺型,其上下底面都是水平的圓。對於相鄰的兩個圓臺,上層的下底面和下層的上底面重合 。第n 層(最上面一層)圓錐的底面就是第n-1 層圓臺的上底面。所有的底面的圓心(包括樹頂)處在同一條與地面垂 直的直線上。李哲知道每一層的高度為h1,h2,…,hn,第1 層圓臺的下底面距地面的高度為h0,以及每層的下底面 的圓的半徑r1,r2,…,rn。李哲用熟知的方法測出了月亮的光線與地面的夾角為alpha。 技術分享圖片
為了便於計算,假設月亮的光線是平行光,且地面是水平的,在計算時忽略樹幹所產生的影子。 李哲當然會算了,但是他希望你也來練練手

Input

第1行包含一個整數n和一個實數alpha,表示檸檬樹的層數和月亮的光線與地面夾角(單位為弧度)。 第2行包含n+1個實數h0,h1,h2,…,hn,表示樹離地的高度和每層的高度。 第3行包含n個實數r1,r2,…,rn,表示檸檬樹每層下底面的圓的半徑。 上述輸入文件中的數據,同一行相鄰的兩個數之間用一個空格分隔。 輸入的所有實數的小數點後可能包含1至10位有效數字。 1≤n≤500,0.3<alpha<π/2,0<hi≤100,0<ri≤100

Output

輸出1個實數,表示樹影的面積。四舍五入保留兩位小數。

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

HINT

Source

[Submit][Status][Discuss]

這個題的正解是各種分類討論求面積,但是可以用辛普森積分騙分。
關於辛普森積分的推導可以看 https://www.zhihu.com/question/47728235
其實這個東西很好卡掉,但是對於很多數據還是可以較為精確地出解的。

首先發現要求的東西是軸對稱的,這意味著我們可以只求x坐標上面的函數區域沒。這個圖像的解析式還是很難求,但是對於每個橫坐標上對應的函數值可以比較輕松地求出,這個時候就可以用Simpson積分了。

普通辛普森積分根本不可能出解,因為僅僅一個二次函數不可能擬合這麽復雜的圖像。但是要分成很多份的話時間復雜度又過高,這個時候就需要引入自適應:如果要對[a,b]積分,可以先對[a,b]擬合,再對[a,(a+b)/2]和[(a+b)/2,b]擬合,如果兩個相差很小則直接退出,否則繼續遞歸。
個人認為這個方法在這裏可以成功很大程度是因為圖像中有很多直線圍城的區域,而自適應Simpson積分在直線積分方面的表現應該是很不錯的。

具體做法看 https://www.cnblogs.com/DaD3zZ-Beyonder/p/5676841.html
上面講的很清楚了,直接作為模板即可。

 1 #include<cmath>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #define rep(i,l,r) for (int i=l; i<=r; i++)
 5 typedef double db;
 6 using namespace std;
 7 
 8 const int N=1010;
 9 double eps=1e-5,inf=1e12;
10 int n,num;
11 db alpha;
12 struct P{ db x,y; P(db X=0,db Y=0){ x=X; y=Y; } };
13 struct Ci{ db r; P c; Ci (P C=(P(0,0)),db R=0):c(C),r(R){}; }C[N];
14 
15 struct Li{
16     P s,t; db k,b;
17     Li (P S=P(0,0),P T=P(0,0)){
18         s=S,t=T; if (s.x>t.x) swap(s,t);
19         k=(s.y-t.y)/(s.x-t.x); b=s.y-k*s.x;
20     }
21     db f(db x){ return k*x+b; }
22 }l[N];
23 
24 int dcmp(db x){ if (fabs(x)<eps) return 0; return (x<0) ? -1 : 1; } 
25 db F(db x){
26     db res=0;
27     rep(i,1,n){
28         db d=fabs(x-C[i].c.x);
29         if (dcmp(d-C[i].r)>0) continue;
30         db len=2*sqrt(C[i].r*C[i].r-d*d);
31         res=max(res,len);
32     }
33     rep(i,1,num) if (x>l[i].s.x && x<=l[i].t.x) res=max(res,2*l[i].f(x));
34     return res;
35 }
36 
37 db calc(db l,db r){ db mid=(l+r)/2; return (F(l)+F(r)+F(mid)*4)*(r-l)/6; }
38 db Simpson(db l,db r,db now){
39     db mid=(l+r)/2,x=calc(l,mid),y=calc(mid,r);
40     if (!dcmp(now-x-y)) return now; else return Simpson(l,mid,x)+Simpson(mid,r,y);
41 }
42 
43 void solve(){
44     db L=inf,R=-inf;
45     rep(i,1,n+1) L=min(L,C[i].c.x-C[i].r),R=max(R,C[i].c.x+C[i].r);
46     rep(i,1,n){
47         db d=C[i+1].c.x-C[i].c.x;
48         if (dcmp(d-fabs(C[i].r-C[i+1].r))<0) continue;
49         db sina=(C[i].r-C[i+1].r)/d,cosa=sqrt(1-sina*sina);
50         l[++num]=Li(P(C[i].c.x+C[i].r*sina,C[i].r*cosa),P(C[i+1].c.x+C[i+1].r*sina,C[i+1].r*cosa));
51     }
52     printf("%.2lf\n",Simpson(L,R,calc(L,R)));
53 }
54 
55 int main(){
56     scanf("%d%lf",&n,&alpha); db h,r;
57     rep(i,1,n+1)
58         scanf("%lf",&h),C[i]=Ci((P((h/tan(alpha))+C[i-1].c.x,0)),0);
59     rep(i,1,n) scanf("%lf",&r),C[i].r=r;
60     solve();
61     return 0;
62 }

[BZOJ1502]月下檸檬樹(自適應辛普森積分)